text:
# Constraints on a planetary origin for the gap in GM Aurigae’s protoplanetary disc

## 1 Introduction

The detection of extrasolar planets and planetary systems (e.g., Mayor & Queloz 1995) has invigorated studies of planetary formation (Lissauer 1993; Boss 1998, 2000), migration (e.g., Lin & Papaloizou 1986; Lin et al. 1996; Trilling et al. 1998; R. Nelson et al. 2000), and interaction with a protoplanetary disc (Lin & Papalizou 1979a, b; Artymowicz & Lubow 1994). Recent theoretical work has investigated mechanisms for halting planetary migration and explaining the statistical distribution of the orbital radii of extrasolar planets (Armitage et al. 2002; Kuchner & Lecar 2002; Trilling, Lunine & Benz 2002). Theoretical models for the interactions of (proto)planets with a protoplanetary disc show that large low density regions (gaps) can be created within the disc (Lin & Papaloizou 1979a, b; Artymowicz & Lubow 1994). The detection of these gaps would provide strong support for planetary formation theories, but current observations do not have the spatial resolution to directly image the au sized gaps that are predicted. Future instrumentation, such as ALMA, may allow direct detection of gaps (Wolf, Henning & Kley 2002), but for now we are restricted to spectral techniques for inferring the presence of disc gaps.

Spectral energy distributions (SEDs) from protoplanetary discs arise from the thermal reprocessing of starlight and dissipation of viscous accretion luminosity (Lynden-Bell & Pringle 1974; Adams, Lada & Shu 1987; Kenyon & Hartmann 1987). The resulting SEDs display characteristic infrared excesses that are sensitive to the size, shape, and mass of the disc. Disc gaps may manifest themselves in SEDs by removing material that would normally have contributed to the SED at a given wavelength: gaps close to the star result in smaller than usual near-IR excesses, while gaps at large disc radii will cause distortions to the SED at longer wavelengths (Beckwith 1999). The classical T Tauri stars GM Aur and TW Hya have SEDs that indicate the inner regions of their discs have been cleared of material out to around 4 au from the central stars (Koerner, Sargent, & Beckwith 1993; Chiang & Goldreich 1999; Schneider et al. 2002; Calvet et al. 2002).

There are a number of possible mechanisms for clearling such gaps. Photoevaporation by ultraviolet radiation, either from the central star or from surrounding radiation, can produce a disc wind at radii where the sound speed exceeds the escape velocity (Clarke et al. 2001). The inner regions of the disc will then be cleared on a viscous timescale once the wind mass loss rate is comparable to the accretion rate. This, however, generally requires an accretion rate two orders of magnitude lower than the $`10^8M_{}\mathrm{yr}^1`$ accretion rate expected for GM Aur (Gullbring et al. 1998). Clarke et al. (2001) also assumed a constant ionizing flux. Matsuyama et al. (2002) show that the clearing of the inner disc is even less likely if the decrease in ionizing flux with decreasing accretion rate is included in the calculation. A. Nelson et al. (2000) also suggest that viscous dissipation and the dissipation of gravitational instabilities may increase the midplane temperature of the inner regions of the disc to above the grain destruction temperature. This would then manifest itself as a reduction in the flux of radiation with wavelengths between a few and ten microns. Their simulation, however, assumed a disc luminosity ($`0.5L_{}`$) an order of magnitude greater than that calculated by Gullbring et al. (1998) for GM Aur ($`0.071L_{}`$).

It has been suggested (Calvet et al. 2002) that the inner regions of TW Hydra’s disc could be cleared by tidal interactions with an embedded planet (Lin & Papalizou 1979a, b). In this paper we consider if the gap in GM Aur’s disc could be formed in the same way. We use a dynamical model of an approximately Jupiter mass planet in a protoplanetary accretion disc, with parameters appropriate for GM Aur (e.g., Koerner, Sargent & Beckwith 1993; Simon, Dutrey & Guilloteau 2000; Schneider et al. 2002), to investigate if the resulting disc structure can reproduce the observed SED. We also consider other planet masses to test if the SED can be used to constrain the planet mass. Estimates for GM Aur’s age vary from $`1.5`$ Myr (Beckwith et al. 1990; Siess, Forestini & Bertout 1999) to $`10`$ Myr (Simon & Prato 1995; Hartmann 2002). There has therefore been sufficient time for an embedded planet to have formed via core accretion (Lissauer 1993) or via gravitational collapse (Boss 1998, 2000). The observation of a much younger system (e.g., $`<1`$ Myr) with an SED signature indicating the presence of a disc gap would, if the gap could be shown to be due to the presence of an embedded planet, suggest that either core accretion occurs on a much shorter timescale than expected or that planet formation can occur via an alternative mechanism such as gravitational collapse. In section 2 we describe and present the results of our dynamical simulations, in § 3 we compare model SEDs with the observed SED of GM Aur, and we summarize our results in § 4.

## 2 Dynamical Simulations

### 2.1 Dynamical Model and Initial conditions

The simulations presented here were performed using smoothed particle hydrodynamics (SPH), a Lagrangian hydrodynamics code (e.g., Benz 1990; Monaghan 1992). We initially assume a planet with mass $`M_p=1.7M_J`$ orbits a star, with mass $`M_{}=0.85M_{}`$, at a radius of $`r_p=2.5`$ au. Both the central star and planet are modelled as point masses onto which gas particles can accrete if they approach to within the sink radius (Bate et al. 1995). The planet orbits within a circumstellar disc of mass $`M_{\mathrm{disc}}=0.047M_{}`$, modelled using 300000 SPH gas particles distributed such as to give a surface density profile of $`\mathrm{\Sigma }r^1`$ (D’Alessio et al. 1998). The disc extends from $`r_{\mathrm{in}}=0.25`$ au to $`r_{\mathrm{out}}=300`$ au. The disc temperature varies as $`Tr^{1/2}`$, as expected for discs heated by starlight (D’Alessio et al. 1998), and we assume hydrostatic equilibrium in the vertical direction, giving a central density profile of $`\rho r^{2.25}`$. For the SPH simulation, the disc is assumed to be locally isothermal and the above temperature profile is maintained throughout the simulation. The parameters chosen are appropriate for GM Aur based on calculations of the stellar mass (Simon, Dutrey & Guilloteau 2000) and the disc mass and size (Schneider et al. 2002).

The vertical scale height, $`H`$, of the disc is given by $`H^2=c_s^2r^3/GM_{}`$ where $`c_s`$ is the local sound speed, $`r`$ is the radius, $`G`$ is the gravitational constant, and $`M_{}`$ is the stellar mass. The temperature is normalised such as to give a scale height of 10.5 au at 100 au. The scale height at 2.5 au is then 0.1 au, giving an $`H/r`$ value of 0.04. This choice of scale height was based on an SED model using an analytic disc geometry and is similar to that obtained by Schneider et al.(2002). In using the above form for the scale height, we assume that the dust has the same scale height as the gas. Although this may not necessarily be the case, the need to use a dust scale height that has the same form as a standard equilibrium disc (e.g., Pringle 1981) would seem to make this a reasonable assumption.

Using 300000 SPH particles the vertical profile of the disk can be modelled accurately within $`3`$ scale heights of the disc midplane. For higher altitudes, the density is calculated analytically using (e.g., Pringle 1981)

$$\rho (z,r)=\rho _c(r)\mathrm{exp}(z^2/2H^2)$$
(1)
where $`z`$ is the height above the disc midplane, and $`\rho _c(r)`$ is the central density a distance $`r`$ from the central star, determined using the dynamical model. Very little mass is contained in these regions of the disc and is hence unimportant for the dynamical calculation. It is, however, extremely important for calculating the radiative transfer through the disc (see D’Alessio et al. 1998).

Self-gravity is included in the calculation and both the point masses and the gas use a tree to calculate gravitational forces (Benz et al. 1990). The inclusion of self-gravity means that the central star and planet are both free to move under the gravitational influence of the disc. A great saving in computational time is obtained by using individual, particle time-steps (Bate et al. 1995; Navarro & White 1993) which are limited by the Courant condition and a force condition (Benz et al. 1990).

### 2.2 Planet-disc interactions

A planet orbiting within a circumstellar disc will tidally interact with the disc. Disc material exterior to the planet will gain angular momentum from the planet and will move to larger radii, while material interior to the planet will lose angular momentum and move to smaller radii (Lin & Papaloizou 1979a, b; Artymowicz & Lubow 1994). Viscous processes within the disc will tend to act against the tidal forces and the formation of a gap will depend on the tidal forces overcoming these viscous processes (Lin & Papaloizou 1979a). This depends on the relative magnitudes of the mass ratio of the planet and the star, $`q`$, and the Reynolds number within the disc, $``$. If $`q>^1`$ a gap will form, while if $`q^2>^1`$ the disc will be truncated at the outer Lindblad resonance. The Reynolds number is given by $`\mathrm{\Omega }r^2/\nu `$ where $`\mathrm{\Omega }`$ is the angular frequency, $`r`$ is the radius, and $`\nu `$ is the turbulent viscosity in the disc generally described as an alpha viscosity (Shakura & Sunyaev 1973).

For the initial simulation ($`M_p=1.7M_J`$, $`M_{\mathrm{disc}}=0.047M_{}`$) considered here $`q=0.002`$, and the Reynolds number, determined from the artificial viscosity included in the code (Murray 1996), is $`20000`$. Therefore $`q>^1`$ but $`q^2<^1`$. A gap should form with the gap edge at the outer Lindblad resonance located, for a planet orbiting at $`2.5`$ au, at a radius $`r=4`$ au. However, for the mass ratio and Reynolds number in this initial simulation, the disc should not be strongly truncated at the outer Lindblad resonance. The Lindblad resonance can be determined from $`\mathrm{\Omega }(r)=\omega m/(1+m)`$ where $`\omega `$ is the angluar frequency of the orbiting planet and, for the outer Lindblad resonance, $`m=1`$ (Binney & Tremaine 1987; Lin & Paploizou 1979a). The angular frequency at the outer Lindblad resonance is therefore $`\mathrm{\Omega }(r)=\omega /2`$. The assumption of a planet orbiting at $`2.5`$ au was based on SED calculations using an analytic disc geometry and the results of Schneider et al.(2002) which suggest that the best fit to the SED occurs for gap extending to a radius of $`4`$ au (i.e., the radius of the outer Lindblad resonance for a planet orbiting at $`2.5`$ au).

Figure 1 shows the azimuthally averaged midplane density at the beginning (dashed line) and end (solid line) of the simulation. After $`1820`$ years the planet orbiting at $`2.5`$ au has produced a gap in the disc out to $`r4`$ au. The gap edge is relatively soft, consistent with the relative magnitudes of the mass ratio ($`q`$) and Reynolds number ($``$). The softness of the gap edge is likely to have an impact on the radiation transfer through the disc. As will be considered in a later section, to produce a gap with a significantly harder edge would require a more massive companion or a significantly less viscous disc. The mass within $`5`$ au has reduced from $`4\times 10^4M_{}`$ at $`t=0`$ to $`4\times 10^5M_{}`$ at $`t=1820`$ years. The reduction in the number of SPH particles respresenting the gas within $`5`$ au requires the use of equation 1 to determine the vertical density at distances greater than $`1`$ scale height above the disc midplane.

If the planet clears sufficient material from within its orbit, torques from material outside its orbit should cause it to eventually lose angular momentum and spiral inwards. However, the planet in this simulation is approximately four times as massive as the disc material contained within its orbit and hence we expect it to suffer minimal orbital migration (Syer & Clarke 1995; R. Nelson et al. 2000).

Accretion onto the central star and the orbiting planet continues throughout the simulation despite the formation of the gap, consistent with simulations performed by Artymowicz & Lubow (1996). The central star has increased from it’s initial mass of $`0.85M_{}`$ to a mass of $`0.85019M_{}`$ after $`1820`$ years, while the orbiting planet has increased from $`1.7\times 10^3M_{}`$ to $`1.88\times 10^3M_{}`$ in the same time period. The accretion rate at the end of the simulation (averaged over the last 100 years) was $`5\times 10^8M_{}\mathrm{yr}^1`$ and $`1\times 10^8M_{}\mathrm{yr}^1`$ for the central star and orbiting planet respectively. An accurate calculation of the accretion rate in our simulation is limited by the mass resolution, since each SPH particle has a mass of $`1.6\times 10^7M_{}`$, and by our simplistic accretion process. We do not include the rotating magnetosphere of the T Tauri star which is likely to truncate the disc at the corotation radius and force the accretion flow to follow the stellar magnetic field lines (Ghosh & Lamb 1978). The accretion rate from the model is however similar to that expected for GM Aur (Gullbring et al. 1998). That mass is flowing through the gap and accreting onto the central star and planet is also consistent with other simulations (Artymowicz & Lubow 1996; Lubow et al. 1999).

The disc velocities also remain very close to Keplerian despite the presence of an orbiting planet that perturbs the potential. This is consistent with CO observations (Koerner, Sargent & Beckwith 1993; Dutrey et al. 1998) which suggest that GM Aur is surrounded by a Keplerian circumstellar disc inclined at $`50^{}`$.

## 3 Radiation Transfer Models

We calculate the SED for our simulations using the Monte Carlo radiative equilibrium technique of Bjorkman & Wood (2001). The output of our radiative transfer simulation is the disc temperature structure and the emergent SED. The disc heating is assumed to be from starlight and the Monte Carlo “photons” are tracked within a two dimensional grid in $`r`$ and $`\theta `$. Our radiation transfer code calculates disc temperature and does not assume a vertically isothermal structure as in the SPH simulation. The input stellar spectrum is a 4000K, $`\mathrm{log}g=3.5`$ Kurucz model atmosphere (Kurucz 1994) and the stellar radius ($`R_{}`$) is taken to be 1.75 $`R_{}`$. The circumstellar dust opacity is taken to be that described in Wood et al. (2002a) which fits the HH 30 IRS SED and was also used to fit the GM Aur SED (Schneider et al. 2002). The dust model has a shallower wavelength dependent opacity than “standard” ISM dust (e.g., Mathis, Rumpl, & Nordsieck 1977) and a larger grain size as inferred for many other classical T Tauri stars (e.g., Beckwith et al. 1990; Beckwith & Sargent 1991).

We find that the dearth of near-IR excess in GM Aur’s SED cannot be explained by a dust composition different to that used here. Using an ISM mixture and a disc without a substantial gap still results in a near-IR excess not observed in GM Aur. Our results are also in agreement with the Calvet et al. (2002) modelling of TW Hya. In their model they do use different dust types in the inner and outer disc, but still require a very low density inner disc to reproduce the near-IR emission observed in TW Hya. Further, ISM grains yield too steep a slope to the sub-mm SED (see also Beckwith & Sargent 1991). We therefore conclude that the dearth of near-IR excess emission from GM Aur can only be explained by an evacuated (low density) inner region.

Circumstellar discs absorb and reprocess stellar photons in a layer which, depending on the disc optical depth, can lie four to five scaleheights above the disc midplane (D’Alessio et al. 1999). Our SPH simulations cannot resolve this low density material, so to provide a density grid for our radiation transfer simulations we have extrapolated the density, using Equation (1), to large heights above the disk midplane. The density calculated from the SPH simulations was azimuthally averaged and then interpolated onto a two dimensional grid in $`r`$ and $`\theta `$. The grid has 100 cells for $`0.25\mathrm{au}r300\mathrm{a}\mathrm{u}`$ spaced proportional to $`r^3`$ and 400 cells linearly spaced for $`0\theta \pi `$. This gridding allows us to adequately resolve the SPH density for our SED calculations.

### 3.1 GM Aur SED

Figure 2 shows the GM Aur SED data. The dotted line shows the input model atmosphere, and the open squares are data points obtained using various instruments (e.g., Cohen & Kuhi 1979; Kenyon & Hartmann 1995; Weaver & Jones 1992; Beckwith & Sargent 1991). The solid line shows the SED calculated using the azimuthally averaged SPH density as input to the Monte Carlo radiation transfer calculation. The SPH simulation was run for $`450`$ orbits of the $`1.7`$ Jupiter mass planet at $`2.5`$ au (1820 years), giving sufficient time for the density profile to reach a steady state. The disc structure is also approximately azimuthally symmetric since the viscous timescale is significantly greater than the planet’s orbital period. The inclination of the disc was taken to be $`50^{}`$, consistent with the value of $`54^{}\pm 5^{}`$ determined by Dutrey et al. (1998), and with the recent scattered light models of Schneider et al. (2002). The model SED provides a good fit to the observations, particularly the lack of near-IR excess and the level and slope of the long wavelength continuum. The dashed line in Figure 2 is the SED for an equivalent disc in which a gap is not present and in which the inner edge of the disc is taken to be at $`7`$ stellar radii. Comparison between the solid line and dashed line clearly shows that the effect of the disc clearing is to reduce the excess emission at short wavelengths ($`210`$ microns) and slightly increase the excess at longer wavelengths ($`1030`$ microns).

This is consequence of the stellar flux being reprocessed at the disc edge. In the simulation with $`R_{\mathrm{min}}=7R_{}`$ (dashed line in Figure 2) stellar photons are reprocessed in high temperature ($`>1000`$ K) regions of the disc giving a near-IR excess. When a gap is present the photons are reprocessed further out in the disc where the temperature is lower, decreasing the near infrared flux and increasing the flux at longer wavelengths. A fit to the data is only possible if the gap extends all the way to the inner disc boundary. If the planet has only had sufficient time to clear a ring of material, leaving the inner disc almost unchanged, this inner material should remain hot, resulting in $`2`$ and $`3`$ micron fluxes similar to that for the dashed line in Figure 2, and higher than that observed for GM Aur (open squares).

### 3.2 Companion (planet) mass and viscosity

Although we are able to reproduce the observed SED with a $`1.7`$ Jupiter mass planet orbiting at $`2.5`$ au, the exact form of the density profile depends on the relationship between the mass ratio, $`q`$, and the disc Reynolds number, $``$. A weakness of SPH is the difficulty in determining exactly the viscosity and hence Reynolds number. In light of this, we have performed a number of additional simulations considering substellar companions (planets) of various masses, each of which was run for at least $`250`$ orbits of the companion around the central star ($`>1000`$ years). In these simulations the disc was modelled between $`0.25`$ and $`100`$ au using $`150000`$ SPH particles with a surface density profile such that the total disc mass would be $`0.047M_{}`$ were the disc extended to $`300`$ au. Using $`150000`$ SPH particles to simulate the inner $`100`$ au of the disc should give approximately the same resolution as the $`300000`$ particles used to simulate the disc out to $`300`$ au. The Reynolds number in the simulations using $`150000`$ SPH particles should therefore be very similar to the Reynolds number in a $`300000`$ particle simulation (Murray 1996). The companion masses considered for these additional simulations were $`0.085`$, $`21`$, and $`43M_J`$. In simulating companions with relatively large masses ($`21`$ and $`43M_J`$) we are not necessarily suggesting that one would expect to observe substellar companions with such masses, but simply illustrating how the SED can vary with mass. Strictly speaking, such massive objects would generally be regarded as brown dwarfs, and we therefore refer to the $`21`$ and $`43M_J`$ companions as substellar objects rather than planets.

The process of star formation is known to preferentially produce binary stellar systems rather than single stars (Duquennoy & Mayor 1991). On the main sequence, however, the mass ratio of close binaries is rarely so extreme as to pair a Solar mass star with a brown dwarf (Marcy & Butler 2000; Halbwachs et al. 2000). If this trend also exists during the pre-main-sequence phase, one would expect disc-embedded planetary mass objects to be more common than embedded close brown dwarfs. The formation mechanism of close binaries is, however, unclear with some simulations (Bate, Bonnell & Bromm 2002) suggesting substantial evolution of binary orbital parameters while a disc is present. The possibility of observing substellar mass companions at this early stage could therefore be an important test of different binary formation models.

The Reynolds number, $``$, in the disc is calculated to be $`20000`$. The $`0.085`$ and $`1.7`$ Jupiter mass planets therefore have $`q>^1`$ but $`q^2<^1`$ while the $`21`$ and $`43`$ Jupiter mass substellar objects have $`q>^1`$ and $`q^2>^1`$. The two planets ($`0.085`$ and $`1.7M_J`$) should therefore produce gaps that are not strongly truncated at the outer Lindblad resonance while the two substellar objects ($`21`$ and $`43M_J`$) should truncate the disc strongly at the outer Lindblad resonance (Lin & Papaloizou 1979a; Syer & Clarke 1995). For a Reynolds number larger than used here, a planetary mass companion ($`<10M_J`$) would satisfy easily the latter condition ($`q>^1`$ and $`q^2>^1`$). The results we obtain for the companions referred to as substellar objects would therefore also be applicable to planetary mass objects in a disc with a Reynolds number larger than in the simulation presented here.

Figure 3 shows the azimuthally averaged midplane densities for the three additional masses together with the azimuthally averaged midplane density profile for the $`1.7M_J`$ planet, between radii of $`0.25`$ and $`25`$ au. As expected the $`0.085`$ and $`1.7`$ Jupiter mass planets produce gaps with soft edges while the $`21`$ and $`43`$ Jupiter mass objects truncate the disc strongly at the outer Lindblad resonance. The position of the gap edge varies slightly for the different masses, but in all cases is close to the outer Lindblad resonance at $`4`$ au.

Figure 4 shows the SEDs calculated using the azimuthally averaged density profiles for the four masses as input to the radiation transfer calculation. Since the simulations considering masses of $`0.085`$, $`21`$, and $`43M_J`$ were only performed out to $`100`$ au, the density profiles were extrapolated to $`300`$ au. The solid and dotted lines are for the $`1.7M_J`$ and $`0.085M_J`$ planets, while the dashed and dash-dot lines are for the $`43M_J`$ and $`21M_J`$ substellar objects. The SEDs for the $`1.7`$ and $`0.085M_J`$ planets are very similar, as are the SEDs for the two substellar objects ($`21`$ and $`43M_J`$), which are, in fact, very difficult to separate. Figure 5 shows a detailed view of the computed SEDs between $`5`$ and $`200`$ microns with the line styles the same as in Figure 4. Figure 5 also shows $`\pm 25`$ % error bars on the data points (open squares). This choice of error was based on worst case 1 $`\sigma `$ error estimates (Cohen & Kuhi 1979; Weaver & Jones 1992; Kenyon & Hartmann 1995). Apart from the SED flux at $`12`$ microns due to the $`0.085M_J`$ planet (dotted line) being slightly higher than that observed, the SEDs computed from the density profiles fit the $`2`$, $`3`$, and $`12`$ micron data points well. The flux between $`3`$ and $`12`$ microns varies quite considerably for the different masses, suggesting that more observations in this wavelength range could allow for the determination of the companion mass or, since the Reynolds number in the disc is unlikely to be accurately known, the relationship betwen the mass ratio and the Reynolds number. The SEDs due to the $`21`$ and $`43M_J`$ substellar objects also have a flux at $`25`$ microns that is considerably higher than that observed. Although not definitive, this may suggest that if the GM Aur SED is due to an embedded companion (planet) orbiting at $`2.5`$ au, the mass must be such as to produce a gap that is not strongly truncated at the outer Lindblad resonance, i.e. $`q>^1`$ but $`q^2<^1`$. This also appears to be the case for TW Hydra as Calvet et al. (2002) require low density material within 4 au to fit the $`210`$ micron SED flux.

At the distance of GM Aur ($`140`$ pc), a few Jupiter mass planet orbiting at $`2.5`$ au will produce a ‘wobble’ in the central star with an angular displacement of $`0.1`$ milliarcsec, and a period of $`4`$ years. Future high precision astrometric instruments should be sufficiently sensitive to make such an observation. Observations of this kind have been suggested for detecting both the presence of embedded planets (Boss 1998) and for detecting instabilities in self-gravitating protoplanetary discs (Rice et al. 2002).

Although it is difficult in SPH both to determine and set the viscosity, by performing a number of simulations with various companion masses we are able to determine how the density profile, and resulting SED, varies for various values of the mass ratio to Reynolds numbers ratio. To further test the relationship between planet mass and viscosity we have also performed a high resolution SPH simulation using 150000 SPH particles to simulate the inner $`25`$ au of the accretion disc. This increases the Reynolds number by a factor of 2 and allows the inner few au of the disc to be represented by more SPH particles once the gap has formed. Figure 6 shows the midplane disc density for the simulation using 300000 particles to represent the entire disc (solid line) and the simulation using 150000 particles to respresent the inner 25 au (dashed line). In both cases a planet mass of $`1.7M_J`$ was used. Although the density profiles differ, this is entirely consistent with the increased Reynolds number of the high resolution simulation. Extrapolating the high resolution simulation to 300 au and using the resulting density profile as input to the radiation transfer produces an SED that, like the previous $`0.085`$ and $`1.7M_J`$ simulations, fits the data remarkably well.

## 4 Summary

We have presented an SPH simulation of a planet opening up a gap in GM Aur’s protoplanetary disk. For a $`1.7M_J`$ planet orbiting at 2.5 au, the disc density within 4 au is decreased by three orders of magnitude. The almost Keplerian velocity structure matches that derived from CO observations. Further, our radiation transfer simulations of a passively heated disk show that the disc structure from the SPH simulation reproduces the observed SED of GM Aur.

Additional simulations of companions (planets) with masses of $`0.085`$, $`21`$, and $`43M_J`$ were also performed. The structure of the resulting gaps depend on the relationship between the mass ratio and the Reynolds number in the disc. For the Reynolds number in our simulations, the two planets ($`0.085`$ and $`1.7M_J`$) produce gaps with relatively soft edges, while the two substellar objects ($`21`$ and $`43M_J`$) truncate the disc strongly at the outer Lindblad resonance ($`4`$ au). As with the $`1.7M_J`$ planet, the SED computed using the density profile due to the $`0.085M_J`$ planet fits the data points extremely well. The SEDs computed using the density profiles due to the two substellar objects ($`21`$ and $`43M_J`$) fit most data points well but have quite a significant excess at $`25`$ microns. Although not definitive, this may suggest that the mass must be such as to produce a gap that is not strongly truncated at the outer Lindblad resonance. Further observations between $`3`$ and $`20`$ microns, where there is a considerable difference in the computed SEDs, but currently no data points, will be needed to determine if this is indeed the case. Future SIRTF observations will provide the required measurements to test our predictions.

Since the density profile depends on the relationship between the mass ratio and the Reynolds number, additional information would be required to determine the actual planet mass. A few Jupiter mass planet orbiting at 2.5 au will produce a ‘wobble’ in the central star, at the distance of GM Aur ($`140`$ pc), with an angular displacement of $`0.1`$ milliarcsec, and a period of $`4`$ years. Future high precision astrometric instruments should be sufficiently sensitive to make such observations. Not only would this establish if a planet is indeed orbiting at $`2.5`$ au but should also allow for the determination of the planet mass from which the Reynolds number, and hence the turbulent viscosity, could be estimated.

The transfer of angular momentum between the embedded planet and the disc material is ultimately expected to cause the planet to migrate inwards (Lin & Paploizou 1986). Simulations (R. Nelson et al 2000) suggest that the migration timescale is of order $`10^5`$ yrs for a Jupiter mass planet starting at 5 au. Since as many as 15 percent of Solar type stars are expected to have planets within $`5`$ au, that there are only currently two T Tauri stars with SED signatures of close-in massive planets is somewhat surprising. This could reflect the relatively sparse observations in the relevant wavelength range ($`312`$ microns). Alternatively, fully evacuated inner holes may form in only a fraction of classical T Tauri stars with planets. An annular gap with a sizeable inner disc would lead to less obvious signatures in the SED.

## acknowledgements

We acknowledge financial support from a PPARC standard grant (WKMR) and Advanced Fellowship (KW). BW would like to acknowledge support from NASA’s Long Term Space Astrophysics (LTSA) Research Program, NAG 5-8412.


format: text

text:
# Exploring Valleys of Aging Systems: The Spin Glass Case

## I Introduction

Physical properties of glassy systems quenched from a high temperature slowly change with the time or age $`t_w`$ elapsed since the quench. For any $`t_w`$ well below the true equilibration time, the dynamics deceptively appears to be stationary when observed on time scales shorter than $`t_w`$. Experimental and numerical evidence for the presence of a quasi-stationary fluctuations regime for $`tt_w`$ followed by non-equilibrium drift for $`tt_w`$, stems e.g. from measurements of conjugate linear response and autocorrelation functions Reim86 ; Refrigier87c ; Andersson92  which obey, respectively violate, the fluctuation-dissipation theorem in the two regimes. In general, the (apparent) age of the system can be deduced from the decrease in the rate of change of macroscopic averages. This apparent age can however be ‘reset’ to an earlier value by applying a perturbation of short duration, as e.g. a temperature pulse, which thus rejuvenates the system.

Non-equilibrium memory effects such as aging and rejuvenation were first noticed and studied in a spin glass context Granberg88 ; Granberg90 ; Schultze91 ; Hoffmann97 , but are now observed in a variety of glassy systems Kityk02 ; Nicodemi01 ; Sibani01 ; Bureau02 ; Hannemann02 . Expanding a previous brief exposition sibani02a , we describe and test a statistical approach to landscape explorations designed to find the generic properties of the energy landscape which produce these effects.

Our main point is that the dynamical events marking the transition between the quasi-equilibrium and the off-equilibrium dynamical regimes are triggered by the attainment of energy values of record magnitude. This immediately allows a description of the non-equilibrium dynamics in terms of a log-Poisson process Sibani93a ; sibanidall03a , i.e. a stochastic process which is homogeneous when viewed on a logarithmic time scale.

The sequel is organized as follows: the next section introduces the landscape exploration method. Section III briefly explains the relevant properties of the log-Poisson statistics used throughout the paper as a semi-analytical description of non-equilibrium dynamics. Section V presents results of an extensive application of the method to spin glass systems, demonstrating the viability of the method and the usefulness of the log-Poisson description. Section VI puts the results in a broader perspective, with special reference to coarse-grained mesoscopic models of configuration space. In particular, we discuss the connection between aging in thermalizing systems and in dissipative driven system and biological evolution. Section VII is a summary and an outlook.

As the nature of true equilibrium and the final stages of thermal relaxation are only weakly related to the presented data, no theoretical consideration is given to these aspects in the paper.

## II The Method

In models of driven dissipative systems with multiple attractors Sibani01 ; Sibani93a ; Coppersmith87 , marginally metastable attractors with an *a priori* negligible statistical weight are nevertheless those typically selected by the dynamics. In such systems, this mechanism underlies memory and rejuvenation effects analogous to those observed in the thermalization of e.g. spin glasses Jonason98  after a quench. One can speculate that similar mechanisms could generally be present in glassy systems with an extensive number of metastable attractors. However, the issue of attractor selection is not explicitly considered in widely used landscape exploration methods such as the Stillinger-Weber approach Stillinger83 ; Nemoto88 ; Becker97 ; Crisanti02 ; Mossa02 , which study a set of local energy minima (inherent states) generated by quenches. The same applies to exhaustive landscape exploration techniques Sibani93 ; Sibani98 ; Schon96  and studies of the real space morphology of low-energy excitations by techniques requiring quenches Andersson96 , genetic algorithms Palassini99  and energy minimization of excitations of fixed volume Houdayer00 ; Lamarcq02 .

The ability of small external perturbations to induce strong rejuvenation and memory effects in complex dynamics strongly suggests that any probe introducing extraneous elements in the dynamical evolution might at the same time yield a biased picture of the energy landscape. In other words: are the attractors identified also those which would be selected by e.g. the unperturbed thermalization process after a deep quench or any other dynamical evolution of interest? Analyzing the regions of state space surrounding intrinsic states provides valuable information about the quasi-equilibrium (fluctuation) dynamics in the energy landscape, but it does not tell the whole story. The question of how to properly describe the non-equilibrium process of ‘selecting’ the metastable states is still open. This question motivates the present approach which is based solely on statistical information collected during an *undisturbed aging process*.

Conventionally, a valley is a connected neighborhood of configuration space which supports a state of approximate thermal equilibrium centered on a local energy minimum. During the time a trajectory ‘resides’ in a valley, the energy and other physical quantities fluctuate around a fixed average, and the state of lowest energy is often revisited; the dynamics has a recurrent character. Non-equilibrium events, henceforth ‘quakes’, move the system irreversibly from the neighborhood of the initial local energy minima of high value and into progressively deeper valleys. A sequence of such events seldom or never revisits the same configurations and has a transient character.

As widely recognized, the lack of time translational invariance in aging systems stems from the dynamical in-equivalence of the valleys visited. Consider therefore a landscape with multiple valleys of varying degrees of metastability, or depth. On time scales larger than the residence time of the deepest (i.e. most stable) valley seen up to the age $`t`$, all valleys shallower than this valley are, by definition, unstable and hence irrelevant for the non-equilibrium dynamics. The interesting dynamical objects are thus the valleys deeper than the deepest valley seen up to time $`t`$. This points to energy records as prospective markers of the non-equilibrium events.

We shall use the term ‘energy barrier’ to denote the difference between the energy of the current state and the lowest seen or best-so-far energy minimum. The lowest energy minimum value and the highest energy barrier observed up to time $`t`$ will identify the valleys as they successively appear in the landscape.

Our classification procedure of the undisturbed dynamics is as follows: We save the minima and barriers of record values encountered and the times at which they occur. We stipulate that a valley is entered at time $`t`$ if the barrier record achieved at time $`t`$ happens to be the last one prior to the attainment of a energy minimum record. To leave a valley, a barrier record must again be followed by a record in the lowest energy. Whenever several minima records are found between two barrier records, we only keep the latest, and therefore deepest, minimum.

In short, we operationally identify the valleys encountered by a series of minima records with at least one barrier event between them. Note that this selection procedure must be performed retrospectively, since it is impossible to know ‘on the fly’ whether a new valley has been entered or not. We also stress that the barrier records will not necessarily coincide with the lowest barrier separating two consecutive energy minima records, as the dynamics, due to entropic effects, is not likely to follow the path of lowest energy, an effect noticed by Wevers et al. Wevers99  in the landscape of metastable ionic compounds.

The discovery of a new record in low energy is a non-equilibrium event. However, by no means does it imply that the internal dynamics in a valley is entirely equilibrium-like. Several sub-valleys are typically explored before the energy minimum is encountered which eventually remains as the lowest state within the valley. Only then does the dynamics acquire the recurrent, fluctuation-like nature which is characteristic of quasi-equilibrium.

Resetting the highest barrier to zero at an arbitrary point in time produces numerous barrier records but no new valleys before a lower energy value is again recorded. However, resetting both the record energy and barrier values to zero may result in a series of new records being registered, which describe sub-valleys within the valley originally explored. By repeating the simulation with the exact same random numbers, this procedure allows one to take a closer look at the internal structure of a valley if the resetting is done at the time of entry.

The method presented is generally applicable, easily implemented, and does not add much to the total runtime of a simulation. If the energy landscape explored is simple, e.g. if it contains one large, structureless valley or if it is perfectly periodic, our scheme only detects a single valley, since degenerate and hence dynamically equivalent minima are appropriately lumped together. In other words, our method produces simple results when used on simple systems. In the following sections we show that non-trivial results are indeed the outcome when complex landscapes are explored.

## III Log-Poisson statistics

In this section we introduce and motivate the log-Poisson statistics used through the paper as an idealized analytical description of non-equilibrium events—such as the quakes which lead to new valleys. We mainly focus on the consequences and predictions for various quantities of physical interest, leaving the empirical justification of the formalism to the next section, where the data are also presented.

A log-Poisson description of complex dynamics was first introduced in connection with a model of Charge Density Waves Sibani93a , and later used to explain macro-evolutionary patterns from the fossil record Sibani95 ; Sibani98a . It describes Sibani99a  the coarse grained dynamics of a population evolving in a NK landscape Kauffman87 , and there are indications that it might also apply to far more realistic models of evolution Hall02 .

The familiar Poisson process with the time argument replaced by its logarithm is in short denoted log-Poisson:

$$P_k(t)=\frac{(\alpha \mathrm{ln}t)^k}{k!}t^\alpha ,t1.$$
(1)
As shown in Sibani93a ; Sibani98a , the probability that $`k`$ records occur in a sequence of $`t`$ random numbers is given by Eq. (1), independently of the underlying process generating the numbers. Log-Poisson statistics implies that the tempo at which the events occur falls off as $`1/t`$. Switching to $`\mathrm{log}t`$ as the independent variable gives a constant (logarithmic) rate of events, and restores time homogeneity. Several other interesting mathematical properties of log-Poisson processes are listed below together with their physical implications for the spin glass problem at hand.

The probability that $`k`$ ‘events’ occur between $`t_w`$ and $`t_w+t`$ is sibanidall03a

$$P_k(t_w,t_w+t)=\frac{1}{k!}\left[\alpha \mathrm{ln}\left(\frac{t_w+t}{t_w}\right)\right]^k\left[\frac{t_w+t}{t_w}\right]^\alpha .$$
(2)
Consider a function $`c(m,m+k)`$ describing the effect of $`k`$ events, for example the overlap between the configurations of lowest energy in valleys $`m`$ and $`m+k`$. To connect with the thermal correlation $`C(t_w,t_w+t)`$, we can utilize a cartoon rendering of the coarse-grained non-equilibrium dynamics as a one dimensional walk on a set of states corresponding to the minima of the valleys, which is similar to e.g. one of the basic assumptions of Ref. rinn00 . We assume, in our terminology, that the dynamical trajectories of an aging system dwell at the bottom of the ‘current’ valley for a certain residence time, until a quake instantaneously moves them to the bottom of the next valley. This coarse grained picture neglects the internal structure of the valleys, which is questionable since the time spent searching for the bottom state within a valley is of the same order of magnitude as the full residence time, see e.g. Fig. 7. Nevertheless, the cartoon has the virtue of simplicity, and allows us to write

$$C(t_w,t_w+t)=\underset{m,k}{}P_m(t_w)P_k(t_w,t_w+t)c(m,m+k).$$
(3)
which is a function of $`t/t_w`$ if and only if $`c(m,m+k)`$ is independent of $`m`$. This *pure* or *full* aging form has very recently been shown to describe aging in real spin glasses, if cooling rate effects through the critical temperature are accounted for rodriguez02 . The assumption that $`c(m,m+k)`$ only depends on a single argument has been checked separately, and an additional dependence on the valley index $`m`$ has been found, which entails a deviation from $`t/t_w`$ scaling dall03b  for the numerical models.

For later convenience we finally note that, if $`t_k`$ marks the time of the $`k`$’th event in a log-Poisson process occurring after $`t`$, the ‘log-waiting time’ $`\mathrm{log}(t_k)\mathrm{log}(t)`$ is exponentially distributed, in perfect analogy to the usual Poisson process. The same exponential distribution also describes the series of *independent* variables $`\mathrm{log}(t_k)\mathrm{log}(t_{k1})`$. This property of the log-waiting time distribution is very well fulfilled by our data.

## IV Models and Dynamics

As glassy dynamics is insensitive to many details of the interactions, computational convenience is a prime criterion for choosing a test model. Spin glass models are comparatively easy to simulate and have been investigated experimentally and numerically for more than twenty years. The lack of a comprehensive and coherent picture of their dynamics and statics furthermore endows them with considerable intrinsic interest.

This section demonstrates that a simple and consistent geometrical picture of the spin glass energy landscape is obtained with our ‘non-invasive’ exploration method. Current spin glass issues are mentioned as needed, while a more complete discussion can be found in Section VI.

We consider $`N`$ Ising spins, where the energy of a spin configuration $`𝐬=\{s_1,\mathrm{},s_N\}`$ is given by

$$E(𝐬)=\frac{1}{2}\underset{i,j}{}J_{ij}s_is_j,$$
(4)
where the couplings $`J_{ij}`$ are symmetric, independent Gaussian variables of unit variance. To test the importance of the topology of the system, we have simulated lattices with periodic boundary conditions in both $`2d`$, $`3d`$ and $`4d`$ ($`N=L^d`$), where $`J_{ij}`$ are non-zero only if $`i`$ and $`j`$ are lattice neighbors. Additionally, we have placed the spins in a $`k`$-regular random graph, i.e. each spin interacts with exactly $`k`$ spins chosen at random.

We use single spin flip dynamics coupled with a rejectionless algorithm, the Waiting Time Method Dall01 . The latter generates a sequence of moves equal in probability to the sequence of *accepted* moves in the Metropolis algorithm. Thus, the results of this paper can also be obtained with the standard Metropolis algorithm, albeit at the price of considerably longer run times. The ‘intrinsic’ and size independent time variable $`t`$ used throughout corresponds to the number of Monte Carlo (lattice) sweeps in the Metropolis as well as to the physical time of a real experiment. We refer to Dall01  for a detailed account of the Waiting Time Method. In all our runs, we conventionally skip the data within the first $`10`$ time units in order to let the system settle down from the random $`T=\mathrm{}`$ initial configuration.

## V Results

Figure 1 shows the average number of valleys $`\overline{n_V}(t)`$ (which is equal to the number of quakes leading from a valley to another) observed in $`3d`$ at various low temperatures. Neglecting the earliest part of the simulations where the memory of the initial random spin configuration is still present, $`\overline{n_V}(t)`$ has a logarithmic shape. As the ratio of the variance $`\sigma _{n_V}^2(t)`$ to $`\overline{n_V}(t)`$ is constant and close to one (see Fig. 2 of ref. sibanidall03a ), the statistics of valleys is essentially a $`\mathrm{log}`$-Poisson process, as expected in a record induced dynamics Sibani93a . As this holds true for barrier records as well (not shown), the ratio of the number of barrier records to minima records is therefore constant on average. This non-trivial geometrical feature of the spin glass landscape has not previously been noted, and it provides a link to an analytical description of the non-equilibrium dynamics as a log-Poisson process Sibani93a ; sibanidall03a .

The barely perceptible curvature seen in Fig. 1 for large $`t`$ decreases systematically as the system size increases at a fixed temperature, as shown in the insert. The sub-logarithmic form of $`\overline{n_V}(t)`$ means that the likelihood that a record in low energy follows a barrier record decreases very slowly, but systematically as the system ages. Taking into consideration that the curvature seems to vanish in the limit of a very large system (see Fig. 2 of ref. sibani02a , which shows that the slope of the straight part of $`\overline{n_V}(t)`$ seems to be proportional to $`\mathrm{log}(N)`$), we surmise that its presence reflects the increasing difficulty in finding new low energies as the ground state is approached.

In summary, the series of quakes moving the system from one valley to the next can be meaningfully idealized into a log-Poisson process discussed in the previous section, if one disregards the curvature of the data. How record events are distributed in time is insensitive to the properties of the stochastic process from which the records are drawn Sibani93a . In our case Figs. 1 and 5 have a weak temperature dependence, as opposed to the strong $`T`$ dependence of the underlying fluctuations. Finally, we note that identical results with respect to $`\overline{n_V}(t)`$ are also found in $`2d`$, $`4d`$ and random graphs.

The rest of this section deals with two scaling plots similar to those presented in ref. sibani02a  (Figs. 2 and 4), as well as new results concerning the scaling of Hamming distances (Fig. 3) between the local minima of contiguous valleys. Finally, the statistics of the residence time and the fraction thereof spent before hitting the lowest energy state in the valley is analyzed, along with the correlations between the minima of the valleys.

### V.1 Scaling of barriers

The data shown in this section are averages over many thousands of realizations of the couplings $`J_{ij}`$ and cover a wide range of low temperatures and system sizes. The raw data for each $`T`$ and $`N`$ have very little scatter, and the main sources of error are systematic. For example, our finite runtimes of $`t=10^6`$, combined with the very broad distribution of residence times in the valleys, bias the average energy of valleys discovered late in the process, since only the ’faster’ trajectories are able to explore these. The data are parametrized by the valley index $`i`$.

In our plots the scaling of the ordinate is given in the corresponding labels. The value of the valley index $`i`$ is shifted from one data set to another by up to one unit. This corresponds to a multiplicative shift of the age, and compensates for the arbitrariness of skipping data within the first 10 time units after the quench, irrespective of temperature and system size.

Figure 2 is a scaling plot of the energy barrier $`B_i`$ which on average must be surmounted in order to leave the $`i`$’th valley. The system is a $`k`$-regular random graph with $`k=6`$, i.e. the same number of links per spin as in $`3d`$. $`B_i`$ is seen to scale in a simple way with the temperature and the size of the system: $`B_i(N,T)=B_i(N^aT^b)`$, where $`a=0.38(2)`$ and $`b=1.8(1)`$. Interestingly, the very same values of $`a`$ and $`b`$ are found when scaling the barrier heights in $`2d`$, $`3d`$ and $`4d`$ lattices as well (see e.g. Fig. 3 in ref. sibani02a ). So, in addition to the non trivial fact that a straightforward two parameter fit works very well, the values of the scaling exponents seem universal. As the next section will reveal, this universality is also found when looking at energy minima and the Hamming distance between them.

The close to linear shape of the scaling plot shows that the barrier heights between valleys gradually increase. Again, we stress that the barriers discussed are those found along an unperturbed trajectory, end hence unlikely to be the lowest barriers between energy minima in contiguous valleys.

Since the exponent $`a`$ is positive, the quakes do not remain localized to a finite set of spins in the macroscopic limit. Some information on the size and shape of the quakes can be inferred from the scaling law for the energy barrier, if we assume that quakes correspond to the motion through the system of a generalized domain wall. Letting $`m`$ be the number of spins typically involved in a quake, we write $`Bm^x`$ for some exponent $`x>0`$. Were the barrier energy the outcome of a fluctuation, i.e. the sum of $`m`$ contributions of random sign, the behavior would be diffusion-like, with $`x=1/2`$. However, since the barrier crossing process favors low-energy states, we expect a lower barrier energy, i.e. $`x<1/2`$. As $`mN^{0.38/x}`$ and $`mN`$, we conclude that $`x[0.38,0.5]`$. In summary, quakes have a fractal shape with exponent $`0.38/x`$, and are more space filling than a conventional domain wall in $`d4`$.

An explanation of the strongly non-Arrhenius $`T`$ dependence of the exponent $`b`$ involves entropic effects and the connection between the quasi-equilibrium dynamics in configuration space and in real space. The qualitative argument given below predicts that the barriers grow linearly with the valley index, and that $`b=2`$ which is close to the actual value of $`b=1.8`$. The discrepancy might arise because we neglect that domain sizes are distributed quantities.

In real space, the relevant objects are connected domains Andersson96 ; kisker96  of thermally equilibrated spins which slowly grow in a sea of frozen spins reaching, on average, a linear length scale $`\mathrm{}(t)`$ on a time scale $`t`$. In the progressively longer quiescent periods between the quakes, these domains do not interact, and their quasi-equilibrium properties are therefore determined by the local density of states $`𝒟(E)`$ pertaining to the configurations accessible to the spins constituting each domain. $`𝒟(E)`$ has been investigated by means of the lid-algorithm Sibani93  for a number of different glassy systems Sibani98 ; Sibani94 ; Sibani99 ; Schon98 ; Schon00 , and has consistently been found to have a close to exponential shape: $`𝒟(E)\mathrm{exp}(E/T_0(\mathrm{}))`$. The energy scale $`T_0(\mathrm{})`$ is an upper bound for the temperature at which metastability can hold and decreases monotonically with the linear size of the system considered Sibani98 . While the results quoted pertain to systems of fixed size and shape, preliminary investigations confirm that the same behavior applies to spin domains of size $`\mathrm{}(t)`$ which grow within a large system. In this case $`T_0(\mathrm{})`$ becomes a slowly decreasing function of time through the time dependence of $`\mathrm{}(t)`$, and asymptotically approaches $`T`$ from above. We anticipate that the marginal stability of the valley visited implies that a domain is typically close to its maximal size, i.e. $`T_0(\mathrm{}(t))T`$.

Returning to $`B_i`$, we recall that energy barriers delimiting valleys are extremal values in a series of $`𝒪(t)`$ independent outcomes in a system of age $`t`$, and that each attempt sees energy $`E`$ with probability

$$P_{eq}(E)\mathrm{exp}\left(E\frac{T_0T}{T_0T}\right).$$
(5)
It follows Leadbetter83  that the typical energy barrier scales as

$$B\frac{T_0T}{T_0T}\mathrm{log}t.$$
(6)
Since the valley index grows linearly with $`\mathrm{log}t`$, the linear dependence on the valley index seen in Fig. 2 is recovered. Finally, since $`T_0T+ϵ`$ for some small $`ϵ`$, one obtains $`b2`$ as anticipated.

### V.2 Scaling of Hamming distances and energies

The fact that contiguous valleys contain rather different low energy configurations is shown in Fig. 3, which depicts the scaling of the Hamming distance $`H`$ (number of spins which differ in their orientation) separating the configurations of lowest energy at the bottom of the two valleys in $`3d`$ lattices. Again, a simple scaling relation of the kind $`H_i(N,T)=H_i(N^aT^b)`$ is found in the $`d`$-dimensional lattices simulated as well as in regular graphs. The value of $`a=0.95(2)`$ holds universally, and it indicates that the number of spins involved is more or less an extensive quantity. The temperature exponent $`b`$ seems to depend slightly on the topology; it has a slightly different value for very low temperatures such as $`T=0.3`$. In short, near-perfect data collapse does not seem possible for the whole range of low temperatures simulated.

We must emphasize that data for large values of the valley index must be interpreted cautiously when scaling the quantities in Figs. 24: Since some runs only reach a few valleys, averages for larger valley indices become biased. This seems particularly important when looking at Hamming distances. Thus, we only plot data for valleys reached by at least $`95\%`$ of all runs. By lowering this threshold to, say, $`80\%`$, one will clearly see that the extra data points lie below the master curve.

The strong temperature dependence of $`H`$ in Fig. 3 is remarkable considering that the energies of the states involved are virtually independent of temperature, as implied by Fig. 4 and as one would expect for actual minima. It follows that the observed minima are nearly degenerate, which is also expected in spin glasses. The $`T`$ dependence of $`H`$ is likely due to the fact that at higher $`T`$ the barriers overcome are higher, and hence the number $`m`$ of spins involved in a quake is larger. In the simplest scenario where the ‘downhill’ part of the quake does not substantially change $`m`$, one can assume $`mHN^{0.95}`$. The energy of the barrier state, which scales as $`BN^{0.38}`$, must be carried by these $`m`$ spins. With $`Bm^xN^{0.95x}`$, we find $`x0.38/0.95=0.4`$, in agreement with the considerations in the previous section.

The energy difference $`\mathrm{\Delta }_i`$ between the state of lowest energy in the first and the $`i`$’th valley is plotted in Fig. 4 for $`2d`$ and $`3d`$ lattices. $`\mathrm{\Delta }_i(N,T)=\mathrm{\Delta }_i(N^aT^b)`$ produces a nice data collapse, with $`a=0.87(3)`$ and $`b=0.2(1)`$. Again, these exponents are found in all topologies simulated, as emphasized by the $`2d`$ data in the insert. The $`2d`$ data bends slightly more than the $`3d`$ data in the main panel, preventing the ultimate data collapse. This is true in general of Figs. 24: although the same scaling form holds across $`2d`$, $`3d`$, $`4d`$, and $`k`$-regular random graphs, they have slightly different curvatures.

The linear trend in Fig. 4 means that the lowest energy decreases (almost) logarithmically with the age, a feature already implicit in the early investigations by Grest et al. Grest86 . It also implies that the energy difference between neighboring valleys remains approximately constant along a trajectory. This means that the energy ‘gain’ induced by a quake per participating spin is of order $`N^{0.87}/N^{0.95}𝒪(1)`$ for the range of system sizes investigated. In other words, the height of the barrier which must be overcome to enter a new valley as well as the overlap between them depend strongly on the temperature, while the amount of energy gained is nearly independent of $`T`$. Hence, lowering the temperature only slows down the process of finding *similar* valleys. Finally, we mention that the starting point $`E_1`$ when measuring $`\mathrm{\Delta }_i`$ is only slightly $`T`$-dependent. If this was not the case, we could not claim the similarity of valleys as stated above.

To sum up, the plots of Figs. 24 tell us that barriers, Hamming distances, and energies of the minima of the valleys can all be scaled with respect to size and temperature in a simple way. Furthermore, the scaling exponents are universal; the same values of the latter are found in $`d`$-dimensional euclidean lattices as well as in regular random graphs.

### V.3 Residence time distribution and superaging

The distribution of time spent in ‘traps’ or valleys of the energy landscape has, to the best of our knowledge, never been measured by others in simulations of spin glass models. The assumed form of this residence time distribution enters heuristic scenarios of spin glass relaxation Bouchaud92 ; Vincent95 , as well as the log-Poisson description of non-equilibrium relaxation sibanidall03a ; dall03b .

Thus, we consider the probability $`R(tt_w)`$ that the residence time $`t_r`$ in a valley entered at age $`t_w`$ be less than $`t`$. As with the scaling plots in the previous sections, a clear picture emerges which lends support to the usefulness of the method for identifying valleys described in Section II.

Before we discuss the numerical results in any detail, we present some theoretical remarks on its expected form for a log-Poisson process. Assuming a pure log-Poisson process for $`t_w>1`$, the ‘log waiting time’ $`\mathrm{ln}(t_r+t_w)\mathrm{ln}(t_w)`$ is exponentially distributed with average $`1/\alpha `$, where $`\alpha =d\overline{n_V}/d\mathrm{ln}t`$ is the logarithmic slope of the curves in Fig. 1. For any $`x>0`$, the probability of $`t_r>t_w(e^x1)`$ is therefore $`e^{\alpha x}`$. Taking $`t=e^x1`$, the age-scaled residence time probability distribution is

$$R(tt_w)=1\left(\frac{t_w+t}{t_w}\right)^\alpha .$$
(7)
This results also follows directly from Eq. (2) by noting that $`R(tt_w)=1P_0(t_w,t_w+t)`$. Since $`\alpha >1`$ sibanidall03a , the distribution has the finite average $`t_r=t_w/(\alpha 1)`$. Equation (7) has been shown to fit the empirical residence remarkably well sibanidall03a .

The $`t/t_w`$ scaling form is often called pure aging, as opposed to super- or sub-aging where the scaled variable is $`t/t_w^\mu `$, with $`\mu >1`$ or $`\mu <1`$, respectively. Beside Fig. 5, superaging of the thermal correlation $`C(t_w,t_w+t)`$ has been observed in numerical investigations of spin glasses berthier02 , while most experimental data show sub-aging vincent97 ; bouchaud99 . Very recent experimental work rodriguez02  shows that subaging is a consequence of the finiteness of the cooling rate. In the limit of a large cooling rate $`\mu `$ approaches one and pure aging is obtained.

The empirical data of Fig. 5 are for random graphs of $`N=16000`$ spins at $`T=0.5`$. As with all the scaling plots presented so far, the quality of the fit is equally convincing in $`d`$-dimensional lattices. To study the age dependence of $`R(tt_w)`$, the valleys entered are analyzed separately for a broad range of waiting time $`t_w=10^2,\mathrm{},10^5`$. As the insert shows, a $`t/t_w`$ scaling is a fair approximation, neglecting the small but systematic drift to the right as $`t_w`$ grows.

The main panel reveals that an excellent data collapse can be obtained with

$$R(tt_w)=R(t/t_w^\mu ),$$
(8)
where $`\mu =1.055`$. The growth of the residence time with $`t_w^\mu `$ shows that the valleys explored become more stable as the system grows older, which concurs with the growth law for the barriers given in Fig. 2. Physically, this means that, compared to the idealized log-Poisson case, the residence time in a valley grows even faster than the age. This deviation from log-time homogeneity is likely due to the already mentioned fact that it becomes relatively harder to find a new valley as the ground state is approached dall03b .

Allowing $`\mu `$ to depend on $`T`$, the scaling works well for all temperatures up to $`T=0.8`$. $`d\mu /dT>0`$ is observed dall03b , i.e. the superaging effect becomes more pronounced as $`T`$ increases. Noting that simulations slightly below the critical temperature, i.e. at the high end of the low temperature phase, are able to explore lower energy regions in the same span of time, this strengthens our hypothesis that deviations from $`t/t_w`$ scaling are a reverberation of the finiteness of the ground state energy. As such, these deviations can be expected to become less important the larger the system is, as they indeed do dall03b .

### V.4 Correlations

We have investigated the overlap between configurations $`m`$ and $`m+k`$ for $`m5`$ and $`k12`$ and found that the form

$$c(m,m+k)=c_{\mathrm{}}(m)+(1c_{\mathrm{}}(m))e^{\gamma (m)k}$$
(9)
fits the data for our range of low temperatures in euclidean lattices as well as in random regular graphs. Only in $`2d`$ do we find significant deviations, as expected considering that the aging is interrupted in this case kisker96 . As an example, Fig. 6 shows $`c(5,5+k)`$ in $`4d`$ models. A similar plot for $`c(1,1+k)`$ in $`3d`$ systems can be found in Fig. 4 of Ref. sibanidall03a , where the connection to the experimentally available non-equilibrium exponent $`\lambda `$ is verified.

If the limiting value for $`k\mathrm{}`$ of $`c(m,m+k)`$, $`c_{\mathrm{}}(m)`$, were independent of $`m`$, it would coincide with the Edwards-Anderson order parameter. Empirically, we find a small $`m`$ dependence, for which we have no physical interpretation. The exponent $`\gamma `$ also has a small and almost linear $`m`$ dependence. As Eq. (9) remains a heuristic approximation of limited value, we have not pursued the size and temperature dependence of the parameters involved in Eq. (9). Still, we find it interesting that a simple exponential parametrization in the number of quakes $`k`$ accurately describes the data.

### V.5 Hitting time for the minima

Having argued that our method of partitioning the sampled states into valleys leads to consistent results and provides useful insight into the coarse structure of complex energy landscapes, it is natural to take a first look at the rich *internal structure* of the valleys. The existence of such a structure is implied by the wide distribution of residence times, which requires a matching distribution of internal energy barriers. Direct evidence is presented in this section, where we consider the ’hitting time’ $`t_{hit}`$ elapsing from the time of entry to the time where the state of lowest is encountered.

Consider $`\tau =t_{hit}/t_r`$, the fraction of the residence time spent ’searching’ for the global minimum. We expect $`\tau `$ to be close to zero in a structureless valley, where the global minimum is reached soon after entry, and close to one in the opposite limit of a rugged valley with many internal minima.

In measuring the distribution of $`\tau `$, we keep track of the age $`t_w`$ at the time of entry, by selecting entry times for a range of different $`t_w`$. From the outset, the empirical distribution $`G(\tau t_w)`$ could depend on $`t_w`$, but, as shown by Fig. 7, there is practically no dependence, except for $`T`$ close to the critical temperature.

The right-skewed form of $`G`$ for lower $`T`$ values (the picture is the same for any $`T0.6`$) implies that by far the greatest part of $`t_r`$ is spent before hitting the lowest minimum: the system explores several sub-valleys, each identified by its own local minimum, while it slides towards states of lower energy. Equilibrium-like fluctuations can only occur after the state of lowest energy which remains at the bottom of the valley is hit. This time interval is small compared to the total residence time, but nevertheless grows steadily with the valley index.

Since the measured shape of $`G`$ at low $`T`$ is almost $`t_w`$ invariant, the distributions of both $`t_r`$ and $`t_{hit}`$ must both scale as $`t_w^\mu `$. We conclude that the process of hitting the lowest minimum in a valley at low $`T`$ is limited by internal barrier crossing events, and that these barriers themselves grow as the age increases and the valleys become larger. The landscape geometry is thus invariant under a time dilation, and hence self-similar, in agreement with the general properties of a log-Poisson process. The temperature independence of the distribution at low $`T`$ leads to the same conclusion: Lowering the temperature means exploring smaller valleys, which nevertheless retain the same internal structure within a range of low temperatures.

As the temperature gets closer to the critical temperature, a waiting time dependence appears: $`G(\tau )`$ remains peaked around $`\tau 0.9`$, but smaller values of $`\tau `$ become gradually more probable. The corresponding flattening of the distribution of $`G(\tau )`$ for $`T=0.8`$, a high temperature, is shown in the insert of Fig. 7. Since the number of potentially new valleys diminishes in time because of the lack of new low energy records, $`G`$ must eventually become highly *left*-skewed. While only a very slow trend towards this situation can be seen in the main panel, the data for the higher temperature in the insert seems to be heading in that direction, in accordance with equilibration happening much faster at high $`T`$ than low $`T`$.

We end this section by noting that, like most other quantities presented in this paper, similar results for $`G(\tau t_w)`$ are found across all systems simulated.

## VI Coarse grained landscape description

Having argued that log-Poisson statistics is a (slightly) idealized description which leads to pure aging and other dynamical features found in experiments, we turn to the physical mechanism behind the selection of the attractors and, more broadly, consider the implications of our results for pertinent mesoscopic models of complex landscapes. We do not attempt to deal with domain growth and other real space issues in any detail. These aspects, though very important for a complete understanding of complex relaxation, are too strongly connected to quasi-equilibrium issues which lie beyond the scope of this paper.

The applicability of log-Poisson statistics to *barrier* records only implies that the trajectory fully decorrelates between successive events, as can be expected for an activated process in a landscape with many local minima. The new and important information about the landscape geometry lies in the fact that *barrier records are required in order to find new minima records*. This explains why the log-Poisson statistics is equally relevant for barriers as for minima. Secondly, and more importantly, it implies that the least stable, or *marginally* stable among the available attractors are those selected. Indeed, if at some stage, very deep minima were to follow a shallow barrier, subsequent barrier records would not produce new states of record-low energy, and the log-Poisson description would fail.

The close match between the depth of a valley and the magnitude of the barrier record giving access to it implies a quasi-continuum of available attractors and was first noticed in connection with the noisy dynamics of a driven dissipative system Sibani93a , where the phenomenon was dubbed noise adaptation. That complex memory behavior is linked to marginal stability of metastable attractors has long been known for noiseless models of Charge Density Waves Coppersmith87 ; Tang87 . The extremely simple automaton model of Tang et al. Tang87  describes a sheet of elastically coupled ‘balls’ driven along a sinusoidal potential by a pulsed external force. A recent study of a noisy version of this model Sibani01  shows that its non-equilibrium dynamics is described by log-Poisson statistics and that the age of the system can be reset by a change of the elastic coupling constant.

The physical origin of the bias towards shallow attractors in thermal glassy dynamics is likely to be entropic, i.e. simply the fact that shallow attractors vastly outnumber deeper ones, in line with the general observation that quenches usually produce poor minima. To support such bias for a range of low temperatures, the density of energy minima must dwarf the Boltzmann factor and hence increase at least exponentially with the energy. This concurs with the outcome of numerical exhaustive investigations of the local configuration space structure of different glassy systems Sibani93 ; Sibani98 ; Sibani94 ; Schon98 ; Schon00 ; Klotz98 ; Schon02  which were performed with the lid method Sibani93 ; Sibani99 . In all cases the local density of states and the local density of minima are nearly exponential functions of the energy.

An exponential density of states in connection with activated dynamics implies a dynamical glass transition. This exact feature was built into the tree model of complex relaxation proposed and analyzed by Grossmann et al. Grossmann85  and later studied in more detail in Hoffmann85 ; Sibani87 . It is also incorporated in the even simpler trap model of Bouchaud Bouchaud92  which, nonetheless, is very different from tree models in one important respect: trees have the lowest connectivity possible for a connected set, while each trap of the trap model is connected to all others.

Log-Poisson statistics only applies as long as new and gradually more stable valleys remain available to the dynamics. The simplest way of modeling inequivalent valleys is through a hierarchy of energy barriers separating degenerate states, which are organized in either a linear array Joh96  or in a tree graph. Beside the models already mentioned, the latter approach is followed in Sibani87 ; Joh96 ; Schreckenberg85 ; Sibani86 ; Sibani89 ; Hoffmann90 . Many important features of complex relaxation can be reproduced in this model, but not the fact that in many systems, including spin glasses, the energy decreases logarithmically with the age. This can be achieved by introducing non-degenerate local minima, as done in the so-called LS tree Hoffmann97 ; Sibani91 . The minima have energies which on average decrease linearly with the size of the barrier overcome, i.e. in the same overall fashion as Fig. 4.

In Bouchauds model, trap energies are exponentially distributed and hence non-degenerate. That deeper minima are gradually explored is a statistical consequence of the (assumed) infinite average of the residence time in a trap. If traps and valleys can be identified, the results of Section V.3 are at variance with this interpretation. The average residence time, which equals $`t_w`$ for $`\alpha =2`$, is, in practice, slightly lower than $`t_w`$. This is again reminiscent of the situation encountered in tree models Hoffmann97 ; Sibani87 ; Schreckenberg85 ; Sibani86 ; Sibani89 ; Hoffmann90 ; Sibani91 . The distribution of barriers in tree models has a lower cut-off, unlike the fractal description of configuration space of Dotsenko Dotsenko85 ; Vincent91 , which better captures the growing importance of gradually smaller barriers as the temperature decreases.

A last important issue is the connection between the energy barrier separating two configurations and the distance between them. For spin glasses the relevant metric is the Hamming distance, which, according to Figs. 3 and 4, on average bears a linear relationship to the energy barrier. A similar result was found in numerical work on the SK model Sherrington75  by Nemoto88 ; Vertechi89 , and by the lid-method (i.e. exact exhaustive enumeration) in Sibani94  for short range spin glasses. Since the latter investigations deal with the small scale structures inside a ‘pocket’, while the present ones are concerned with the large scale structures explored by the non-equilibrium dynamics, the agreement in their outcome is further evidence of a self-similar landscape structure. Interestingly, the *largest* distance which can be achieved for a fixed energy barrier grows exponentially with the barrier Sibani98 . A linear relationship between energy and Hamming distance is assumed in tree models of aging dynamics, see e.g. Sibani89 , and also in the barrier model Joh96 . This is, however, not a crucial assumption for aging, and other types of functional relationships can also be utilized Hoffmann97  successfully.

## VII Summary and outlook

In this paper we have presented a general ‘non-invasive’ statistical method for complex energy landscape exploration, especially designed to provide information on the metastable states actually explored by an unperturbed aging process following a quench. The method has been thoroughly tested on Ising spin glasses, and the results obtained both match and extend the established knowledge about spin glasses. In particular, most quantities investigated obey simple scaling laws with universal scaling exponents. The view which emerges is that non-equilibrium aging dynamics is steered by energy barrier *records*, which are the only events capable of opening the route to new valleys. Since this description has previously been shown to apply to driven dissipative models and to evolution modeling, a possible unified theory for non-equilibrium glassy dynamics seems within reach.

A more complete picture can be obtained by looking more closely into the fluctuation dynamics which may differ across different systems. For spin glasses we have argued that real space domains of (pseudo) thermalized spins relax independently as long as the system as a whole remains in the same valley. The present method opens the possibility of identifying the quasi-equilibrium clusters as defined by the dynamics itself: These clusters are separated by a backbone of spins whose orientation remains fixed within each valley and changes slowly from one valley to the next, as seen in Fig. 3. After the lowest energy state in the valley has been hit and before the next valley is entered, there is no drift towards the global minimum. Hence, the spins fall into two categories only: those which are frozen, and those which fluctuate in a quasi-equilibrium fashion. The quasi-equilibrium clusters can thus be extracted and their statistical properties, such as e.g. the density of states, can be studied for each cluster separately.

Acknowledgments: This project has been supported by Statens Naturvidenskabelige Forskningsråd through a block grant and by the Danish Center for Super Computing with computer time on the Horseshoe Linux Cluster. We are grateful to J. Christian Schön for many discussions, and, in particular, for pointing out a flawed mathematical argument in an early version of this work.


text:
# Triangular and 𝑌-shaped hadrons with static sources

(Institute of Theoretical and Experimental Physics,
117218, B.Cheremushkinskaya 25, Moscow, Russia)

## Abstract

The structure of hadrons consisting of three static color sources in fundamental (baryons) or adjoint (three-gluon glueballs) representations is studied. The static potentials of glueballs as well as gluon field distributions in glueballs and baryons are calculated in the framework of field correlator method.

1. The study of the gluon field structure in hadrons is important for the deep understanding of the strong interaction physics. Interesting results on the baryon flux structure were obtained recently in the maximally abelian gauge on the lattice . We demonstrate below how to perform calculations in gauge invariant way using the field correlator method relyed directly on QCD.

The static potential in baryon was considered in detail in  at this conference, and we will extensively refer to this talk in what follows. We start with the calculations of static potentials in three-gluon glueballs. This part of the talk presents the results from . In the rest of the talk the gluon field distributions in hadrons are discussed on the base of papers , , .

In contrast to baryons, the gauge-invariant extended wave function of glueballs may have both $`Y`$-type and triangular structure. In the former case it may be written as

$$G_Y^{(f)}(x,y,z,Y)=f^{abc}g_a(x,Y)g_b(y,Y)g_c(z,Y)\text{ or}$$
(1)
$$G_Y^{(g)}(x,y,z,Y)=g^{abc}g_a(x,Y)g_b(y,Y)g_c(z,Y).$$
(2)
Here $`g_a(x,Y)=g_a(x)\mathrm{\Phi }^{ab}(x,Y)`$ denotes the extended gluon operator, $`g_a(x)`$ is the valence gluon operator and $`\mathrm{\Phi }^{ab}(x,Y)=(\mathrm{P}\mathrm{exp}ig_Y^xA_\mu 𝑑z_\mu )^{ab}`$ is the parallel transporter or Schwinger line in the adjoint representation. The coordinates $`x,y,z`$ and $`Y`$ in (1) apply to the valence gluons and string junction positions respectively, $`f^{abc}`$ and $`g^{abc}`$ denote adjoint antisymmetric and symmetric symbols.

The wave function of the triangular glueball has the form

$$G_\mathrm{\Delta }(x,y,z)=G_\alpha ^\beta (x)\mathrm{\Phi }_\beta ^\gamma (x,y)G_\gamma ^\delta (y)\mathrm{\Phi }_\delta ^\epsilon (y,z)G_\epsilon ^\rho (z)\mathrm{\Phi }_\rho ^\alpha (z,x),$$
(3)
where $`G_\alpha ^\beta (x)=g^a(x)(t^a)_\alpha ^\beta `$, $`t`$ is the generator of $`SU(3)`$, and $`\mathrm{\Phi }_\alpha ^\beta `$ is the parallel transporter in fundamental representation.

According to (1), (2), a Wilson loop of $`Y`$-type glueball has the same structure as in the case of baryon, see eq. (2) and Fig. 1 of . The difference is that the adjoint transporters should be substituted for the fundamental ones and symbols $`f^{abc}`$ or $`g^{abc}`$ for $`ϵ_{\alpha \beta \gamma }`$.

A Wilson loop of triangular glueball, induced by (3) in the limit of large $`N_c`$, is shown in Fig. 1. The trajectories of valence gluons are replaced here with two fundamental (quark) trajectories and the Wilson loop is reduced to three disconnected rectangular meson loops.

The static potential in baryon calculated in the field correlator method (MFC) is proportional to two-point gluon field strength correlators in fundamental representation, see e.g. (5), (6) from . Therefore it is proportional to the fundamental Casimir operator, $`C_F=4/3`$, and the ratio of $`Y`$-type glueball and baryon potentials equals to the ratio of corresponding Casimir operators,

$$\frac{V_Y^{(G)}}{V_Y^{(B)}}=\frac{C_A}{C_F}=\frac{9}{4},$$
(4)
where $`C_A=3`$ is the adjoint Casimir operator and nondiagonal (interference) terms are neglected. The diagonal part of the baryon potential in equilateral triangle has the form

$$V_Y^{(B)}(R)=\frac{6\sigma }{\pi }\left\{R_0^{R/T_g}𝑑xxK_1(x)T_g\left(2\frac{R^2}{T_g^2}K_2\left(\frac{R}{T_g}\right)\right)\right\}3V^{(M)}(R),$$
(5)
where $`R`$ is the distance from quarks to string junction, $`K_1,K_2`$ are McDonald functions, $`\sigma `$=0.18 GeV2 and $`T_g=0.12÷0.2`$ fm are the string tension and the gluon field correlation length, and $`V^{(M)}(R)`$ is the static potential in meson.

Now we turn to the triangular glueball, whose Wilson loop according to Fig. 1 is the product of three meson loops. Therefore, the potential $`V_\mathrm{\Delta }^{(G)}`$, which is the logarithm of the Wilson loop, is the sum of three meson potentials. In the case of equilateral triangle with the side $`r=\sqrt{3}R`$

$$V_\mathrm{\Delta }^{(G)}(r)=3V^{(M)}(r).$$
(6)
The perturbative potential for three-gluon glueballs reads as

$$V_{(\mathrm{adj})}^{\mathrm{pert}}(r)=\frac{3}{2}\frac{C_2(\mathrm{adj})\alpha _s}{r}.$$
(7)
The behavior of total potentials in $`Y`$-type and triangular glueballs are shown in Fig. 2. One can see that they are very close up to distances $`r0.6`$ fm. Note that the interference terms are to change this picture. They will be considered in subsequent publications.

2. We proceed now to the study of the gluon field structure using the connected probe , . The connected probe consists of the probe plaquette joined to the Wilson loop by parallel transporters and forms the frame with the current in four-dimentional euclidian space, see Fig. 3. When the plaquette size is small enough, the connected probe affords to calculate the color-integrated gluon field in hadron, using the Wilson loop $`𝒲^{}`$ of the latter,

$$F_{\mu \nu }^H(x)=ig\frac{𝒲_{}^{}{}_{\beta }{}^{\alpha }(x_0)\mathrm{\Phi }_\gamma ^\beta (x_0,x)(F_{\mu \nu }^a(x)t^a)_\delta ^\gamma (\mathrm{\Phi }^+)_\alpha ^\delta (x,x_0)}{𝒲^{}}.$$
(8)
The following expression is valid in the bilocal approximation of MFC in the case of mesons ,

$$F_{\mu \nu }^M(x)=_S𝑑\sigma _{\rho \sigma }(x^{})𝒟_{\rho \sigma ,\mu \nu }(x^{}x),$$
(9)
where the integration is taken over the surface of Wilson loop $`S`$ and $`𝒟`$ denotes the gauge-invariant bilocal gluon field strenth correlator, see e.g. (5) of  for details. Using the MFC parametrization of correlators , we calculate that the electric component of the $`F^M`$ has the form

$$𝐄^M(𝐱,𝐑)=𝐧\frac{2\sigma }{\pi }_0^{R/T_g}𝑑l\left|l𝐧\frac{𝐱}{T_g}\right|K_1\left(\left|l𝐧\frac{𝐱}{T_g}\right|\right),$$
(10)
while the magnetic one is absent.

The field $`𝐄^M(𝐱,𝐑)`$ (10), which we call the nonperturbative background gluon field in meson, is the force acting on the probe located at the point $`𝐱`$ while the quark situated at zero and antiquark at the point $`𝐑`$. It is related to the field acting on the quark as follows.

The force $`𝐅^M`$ acting on the quark in meson is defined by the static potential $`V^{(M)}`$ (5),

$$𝐅^M(𝐫)=\mathbf{}V^{(M)}(𝐫).$$
(11)
One can verify using (5), (10), (11) that the following relations are valid ,

$$𝐄^M(0,𝐑)=𝐄^M(𝐑,𝐑)=𝐅^M(𝐑),$$
(12)
$$𝐄^M(𝐑,2𝐑)=𝐄^M(0,𝐑)+𝐄^M(𝐑,𝐑)=2𝐅^M(𝐑).$$
(13)
The relation (12) mean that when the locations of probe and quark (antiquark) coinside, the probe recombines with the latter, and that is why it is affected by the same force $`𝐅^M`$. One can conclude from the relation (13) that when the probe is situated in the middle of two sources, it interacts with both the quark and antiquark, and the total field becomes twice as big. If the point $`𝐱`$ is located on the line connecting the quark and antiquark, it is easy to check using (5), (10) that the generalization of (12), (13) has the form

$$𝐄^M(𝐱,𝐑)=𝐅^M(𝐱)𝐅^M(𝐱𝐑).$$
(14)
According to (11), the force $`𝐅^M(𝐑)`$ acting on the quark increases linearly with the slope $`2\sigma /(\pi T_g)`$ at small distances $`R\begin{array}{c}<\hfill \\ \hfill \end{array}T_g`$ and saturates with the value $`\sigma `$ at $`R\begin{array}{c}>\hfill \\ \hfill \end{array}0.7`$ fm.

The field $`𝐄^M`$ in the saturated regime acquires the universal profile $`E^{\mathrm{string}}`$, which does not depend on the quark-antiquark separation,

$$E^{\mathrm{string}}(\rho )=2\sigma \left(1+\frac{\rho }{T_g}\right)\mathrm{exp}\left(\frac{\rho }{T_g}\right),$$
(15)
where $`\rho `$ is the distance from the quark-antiquark axis. We apply now the Gauss law in the form

$$𝐄𝑑𝐒=4\pi C_F\alpha _s$$
(16)
to the saturated field (15) and get the parameter relation

$$3\sigma T_g=C_F\alpha _s.$$
(17)
Taking the freezing value of the strong coupling $`\alpha _s0.4`$  and the phenomenological value of the string tension $`\sigma =0.18`$ GeV2, we determine from (17) the reasonable gluon correlation length value $`T_g0.2`$ fm (see ).

Using the Maxwell equation with magnetic currents,

$$\mathrm{rot}𝐄^M=𝐣_{\mathrm{magn}},$$
(18)
we obtain that currents corresponding to the saturated field (15) form closed circles around the quark-antiquark axis, with the density

$$j_{\mathrm{magn}}(\rho )=\frac{E^{\mathrm{string}}(\rho )}{\rho }=\frac{2\sigma \rho }{T_g^2}\mathrm{exp}\left(\frac{\rho }{T_g}\right).$$
(19)
One can verify using (9), (18) and nonabelian Bianchi identity that magnetic currents arise due to the three-point field strength correlator, which describes the emitting of the color-magnetic gluon field by the color-electric one.

The detailed study of these and other properties of the field $`𝐄^M`$ will be given elsewhere.

3. When the field distribution of the quark-antiquark pair is known, it is straightforward to calculate the field of triangular glueball (see Fig. 1). We join the connected probe to each of quark-antiquark loops and arrive at the expression

$$𝐄_\mathrm{\Delta }^{(G)}(𝐱,𝐫^{(1)},𝐫^{(2)},𝐫^{(3)})=\underset{i=1}{\overset{3}{}}𝐄^M(𝐱𝐫^{(i)},𝐑^{(i)}),$$
(20)
where $`𝐫^{(i)}`$ is the position of $`i`$-th valence gluon and $`𝐑^{(i)}=𝐫^{(i+1)}𝐫^{(i)}`$. The surface defined by the condition $`|𝐄_\mathrm{\Delta }^{(G)}(𝐱)|=\sigma `$ at gluon separations 1 fm is shown in Fig. 4. Values of parameters $`\sigma =0.18`$ GeV2, $`T_g=0.12`$ fm are used. The surface shown in the figure goes through the valence gluon locations. Indeed, one can verify using (20) that two forces of the value $`\sigma `$ having the angle $`2\pi /3`$ between them add to give the resulting force of the value $`\sigma `$.

The same procedure can be applied to the system of arbitrary number of mesons to obtain the field distributions in the leading order of $`1/N_c`$. The next order corrections can also be considered within the approach.

4. The problem of baryon (or $`Y`$-type glueball) field calculations is more complicated. When we join the probe plaquette to the trajectory of the first quark and calculate the field distribution (8) using the bilocal approximation of MFC, we arrive at the following relation for the electric field,

$$𝐄_{(1)}^B(𝐱,𝐑^{(1)},𝐑^{(2)},𝐑^{(3)})=𝐄^M(𝐱,𝐑^{(1)})\frac{1}{2}𝐄^M(𝐱,𝐑^{(2)})\frac{1}{2}𝐄^M(𝐱,𝐑^{(3)})$$
(21)
(the magnetic field is absent just as in the meson case). The vector $`𝐑^{(i)}`$ here is directed from the string junction to the $`i`$-th quark, $`𝐄^M`$ is defined in (10). When we join the probe to the second and third quark trajectories, we will get analogous formulas with the corresponding index transpositions.

To obtain symmetric field distributions, we are to sum squares $`(𝐄_{(i)}^B)^2`$,

$$(𝐄^{(B)})^2=\frac{2}{3}\left((𝐄_{(1)}^B)^2+(𝐄_{(2)}^B)^2+(𝐄_{(3)}^B)^2\right),$$
(22)
The surface given by the condition $`|𝐄^{(B)}(𝐱)|=\sigma `$ at quark separations 1 fm is shown in Fig. 4. It goes through the quark positions and has small convexity near the string junction.

The ratio of the energy density at the string junction position, $`w_{\mathrm{sj}}`$, to the density in the middle of the string with the saturated profile, $`w_{\mathrm{string}}`$, has to be equal to the corresponding field squares and according to (22), (15)

$$\frac{w_{\mathrm{sj}}}{w_{\mathrm{string}}}=\frac{(𝐄^{(B)}(0))^2}{(𝐄^{\mathrm{string}}(0))^2}=\frac{9}{8}.$$
(23)
To conclude, hadron static potentials and the gluon field distributions, which ensure the confinement of color sources, were calculated using the field correlator method. It is stressed that the mechanism of confinement is the emitting process of the color-magnetic gluon field by the color-electric one in the nonperturbative vacuum. The static potentials of three-gluon glueballs were reduced in the first approximation to the appropriate sum of quark-antiquark potentials. The nonperturbative gluon field induced by the static quark-antiquark pair was calculated using the connected probe, and its classical abelian properties were considered. Triangular gluon field distributions were calculated for glueballs and $`Y`$-type ones for baryons.

An extension of the method to the arbitrary number of colors is straightforward. It is also applicable to the nuclear structure study.

This work has been supported by RFBR grants 00-02-17836, 00-15-96786, and INTAS 00-00110, 00-00366.