text:
# 1 Introduction

## 1 Introduction

Despite the progress made in recent years , our understanding of the dynamics of deeply supercooled liquids and the properties of glasses is still far from satisfactory . Although the so-called mode-coupling theory (MCT)  seems to give a very reliable description of this dynamics around the MCT-temperature $`T_c`$, the relaxation dynamics significantly below $`T_c`$, i.e. close to the experimental glass transition temperature $`T_g`$, still waits for a deeper theoretical understanding. The problem is that in this temperature range the mechanism leading to the relaxation of the system seems to be governed by the so-called hopping processes, i.e. the system relaxes in an activated way, and that these processes are not readily describable by MCT in its ideal version. (Although the extended version of MCT , which to some extend takes into account these hopping processes, is able to make precise predictions for the dynamics at intermediate time scales, it looses its predictive power at long times, i.e. the $`\alpha `$relaxation regime.) Although experiments are very useful to give information on particle averaged properties of supercooled liquids and glasses, such as the viscosity or the intermediate scattering function, they are less well adapted to study the system on a very local scale, e.g. on the level of a single particle. This is, however, exactly the information which is needed if one wants to come to an understanding of the relaxation mechanism in these systems if one wants to go beyond the use of simple models, such as a particle moving around in a frozen environment or similar simple pictures. One method which in the past has been proven to be extremely useful to reveal the necessary details of structure and dynamics of complex systems are computer simulations . Since such simulations allow to access at any instant the full microscopic information, they are also ideally suited to learn more about the relaxation mechanisms in supercooled systems and glasses . One problem of such simulations is, however, that the time and length scales accessible are rather limited (typically 100 ns and box sizes of 100Å3) and this has prevented in the past to study, e.g., the equilibrium dynamics of supercooled liquids close to $`T_g`$, where the typical time scales are on the order of 10-1000 s, i.e. 15-17 decades longer than the typical time step of used in a simulation which is on the order of 1/100 of the inverse of the Debye frequency, i.e. 1-10 fs. One possibility to avoid this problem to some extend is to use a sophisticated Monte Carlo algorithm that allows to equilibrate the system even at those low temperatures at which standard molecular dynamics or Monte Carlo algorithms fail to reach equilibration. In recent years various authors have proposed different Monte Carlo algorithms that should be suited to reach the equilibrium state even at low temperatures  and have also been successfully applied to various systems  (see, however, Ref.  for discouraging results). Here we will discuss the application of the so-called parallel tempering (PT) algorithm, proposed by Hukushima et al. , a method that is similar in spirit to the replica exchange algorithm proposed earlier by Swendsen and Wang , to the case of structural glasses and to the Potts-glass. Although for structural glasses with a small number of particles ($`N36`$) the PT method has been successfully tested by Coluzzi and Parisi  we are here interested in using larger systems ($`N=336`$ and $`N=1000`$). In the next sections we will discuss the details of the PT algorithm, then introduce the models we investigated and subsequently discuss the results.

## 2 The parallel tempering algorithm

The PT method proposed by Hukushima et al. can be summarized as follows : 1) We denote the Hamiltonian of interest by $`H=K(𝐩)+𝐄(𝐪)`$, where $`K`$ and $`E`$ are the kinetic and potential energy, respectively, and $`𝐩=(p_1,p_2,\mathrm{},p_N)`$ and $`𝐪=(q_1,q_2,\mathrm{},q_N)`$ are the momenta and coordinates of the particles, respectively. (If this is a spin-like Hamiltonian, we put $`K=0`$.) We now construct a new system consisting of $`M`$ noninteracting subsystems, each composed of $`N`$ particles, with a set of arbitrary particle configurations $`\{𝐪_1,\mathrm{},𝐪_M\}`$ and momenta $`\{𝐩_1,\mathrm{},𝐩_M\}`$. The Hamiltonian of the $`i`$-th subsystem is given by

$$H_i(𝐩_i,𝐪_i)=K(𝐩_i)+\mathrm{\Lambda }_iE(𝐪_i),$$
(1)
where $`\mathrm{\Lambda }_i\{\lambda _1,\mathrm{},\lambda _M\}`$ is a parameter to scale the potential. 2) A molecular dynamics simulation is done for the total system, whose Hamiltonian is given by $`=_{i=1}^MH_i`$, at a constant temperature $`T=\beta _0^1`$. In this way we obtain a canonical distribution $`P(𝐪_1,\mathrm{},𝐪_M;\beta _0)=_{i=1}^MP(𝐪_i;\mathrm{\Lambda }_i\beta _0)\mathrm{exp}[\beta _0_{i=1}^M\mathrm{\Lambda }_iE(𝐪_i)]`$ in configuration space. 3) After each time interval $`\mathrm{\Delta }t_{PT}`$, we attempt to exchange the potential scaling parameter of the $`m`$-th and $`n`$-th subsystem, while $`\{𝐪_1,\mathrm{},𝐪_M\}`$ and $`\{𝐩_1,\mathrm{},𝐩_M\}`$ are unchanged. The acceptance of the exchange is decided in such a way that it takes care of the condition of detailed balance. Here we use the Metropolis scheme, and thus the acceptance ratio is given by

$$w_{m,n}=\{\begin{array}{cc}1,\hfill & \mathrm{\Delta }_{m,n}0\hfill \\ \mathrm{exp}(\mathrm{\Delta }_{m,n}),\hfill & \mathrm{\Delta }_{m,n}>0,\hfill \end{array}$$
(2)
where $`\mathrm{\Delta }_{m,n}=\beta _0(\mathrm{\Lambda }_n\mathrm{\Lambda }_m)(E(𝐪_m)E(𝐪_n))`$. 4) If steps 2) and 3) are repeated for a sufficiently long time this scheme leads to canonical distribution functions $`P(E;\beta _i)`$ at a set of inverse temperatures $`\beta _i=\lambda _i\beta _0`$. To make a measurement at an inverse temperature $`\beta _l`$ one has to average over all those subsystems $`(i1,\mathrm{},M)`$ for which we have (temporarily) $`\beta _l=\lambda _i\beta _0`$.

Although the algorithm presented is correct for general choice of the values of $`m`$ and $`n`$ in Eq. (2), it is advisable to exchange only neighboring subsystems in order to allow for a reasonably high acceptance rate. In addition one has also to choose the values of the coupling constants $`\{\lambda _1,\mathrm{},\lambda _M\}`$ in such a way that neighboring subsystem have an sufficiently large overlap in their distribution of the energy.

We also note that the algorithm we just described shares many properties with the one proposed by Swendsen and Wang in 1986 . These authors already realized that it is very useful to connect dynamically configurations at low temperatures to the ones at high temperatures and also gave a specific method how this can be successfully done in the case of a frustrated spin system.

## 3 Models and Details of the Simulations

In this paper we consider three different types of glassy systems: A realistic model for SiO2, one of the prototypical glass formers. A binary mixture of Lennard-Jones particles, i.e. an example for a simple glass former, and finally a ten state Potts-glass, an example for a spin system with a discontinuous transition from the paramagnetic phase to the spin glass phase.

SiO2 is a prototype for a so-called strong glass former . This means that the temperature dependence of transport quantities like the viscosity or the diffusion constant show an Arrhenius dependence. It is believed that this property is related to the fact that the structure of amorphous silica is a open tetrahedral network that is essentially independent of temperature.

The SiO2 model we use has been proposed by van Beest et al. (BKS)  on the basis of ab initio calculations. In this system two ions of type $`\alpha `$ and $`\beta `$ ($`\alpha ,\beta \{\mathrm{Si},\mathrm{O}\})`$ that are a distance $`r`$ apart interact via the following potential:

$$\varphi (r)=\frac{q_\alpha q_\beta e^2}{r}+A_{\alpha \beta }\mathrm{exp}\left(B_{\alpha \beta }r\right)\frac{C_{\alpha \beta }}{r^6},$$
(3)
where the values of the parameters $`q_\alpha `$, $`A_{\alpha \beta }`$, $`B_{\alpha \beta }`$, and $`C_{\alpha \beta }`$ can be found in Ref. . Previous computer simulations have shown that the BKS model gives a very good description of the static and dynamic properties of real silica . However, so far these tests could be done only at relatively high temperatures, i.e. $`T2750`$K, since below this temperature the relaxation time exceeds the time scale accessible to normal molecular dynamics simulations (in this case 20 ns). This temperature has to be compared with the experimental value of $`T_g`$, 1450K, i.e. so far it has been possible to equilibrate the system only at temperatures two times higher than $`T_g`$.

In the following we will present results in which we used the PT algorithm to equilibrate the system at low temperatures. Between the attempted exchanges between two subsystems we propagated the particles in the isokinetic ensemble at constant volume. The system size was 336 ions and the time step was 1.6 fs. The results presented below are for $`M=32`$ subsystems, the number of time steps for equilibration and production was each $`410^6`$, and the value of $`\mathrm{\Delta }t_{PT}`$ was 1000 time steps. More details on this simulation can be found in Ref. .

The second structural glass is a binary (80:20) mixture of Lennard-Jones particles. If we denote the majority species by A and the minority species by B the interaction between two particles are given by $`\varphi _{\alpha \beta }(r)=4ϵ_{\alpha \beta }[(\sigma _{\alpha \beta }/r)^{12}(\sigma _{\alpha \beta }/r)^6]`$, where $`r`$ is the distance between particles $`i`$ and $`j`$. The interaction parameters are $`\alpha ,\beta \{\mathrm{A},\mathrm{B}\}`$, $`ϵ_{\mathrm{AA}}=1`$, $`ϵ_{\mathrm{AB}}=1.5`$, $`ϵ_{\mathrm{BB}}=0.5`$, $`\sigma _{\mathrm{AA}}=1`$, $`\sigma _{\mathrm{AB}}=0.8`$, and $`\sigma _{\mathrm{BB}}=0.88`$. In the following we will measure length and energy in units of $`\sigma _{\mathrm{AA}}`$ and $`ϵ_{\mathrm{AA}}`$, respectively, (setting $`k_B=1`$) and time in units of $`(m\sigma _{AA}^2/48ϵ_{AA})^{1/2}`$, where $`m`$ is the mass of the particles. In the simulation we used a cubic box, of length 9.4, with periodic boundary conditions, the total number of particles was 1000, and the time step was 0.01732. The number of subsystems was 16, the number of time steps for equilibration and production was each $`510^6`$, and the value of $`\mathrm{\Delta }t_{PT}`$ was 1000 time steps. More details on the simulation can be found in Ref. . In the past many properties of this model have been investigated  and it has been found that its relaxation dynamics becomes very slow at around $`T=0.450.43`$, i.e. its starts to exceed the time scales accessible to normal molecular dynamics simulations, which is on the order of $`10^8`$ time steps. In contrast to the Arrhenius-like increase of the relaxation times, as it is found in SiO2, this system shows at low temperatures an increase which can be fitted well with a power-law, and thus this is considered to be a fragile glass former . The structure of the system resembles the one of randomly closed packed hard spheres and is thus very different from the open network in silica. Last not least it has to be mentioned that the simplicity of the interaction of this model allows to obtain results for this system which have a significantly higher statistical accuracy than the ones for the silica model discussed above, since in the latter one has to calculate numerically expensive long range interactions. This is the reason why in the results discussed below the data for the silica system is quite a bit more noisy than the one for the Lennard-Jones system.

The Potts-glass we consider is an example of a spin glass which shows a discontinuous transition from the paramagnetic phase to a spin glass phase if the number of states is larger than 4 . It has been suggested that such type of models show a qualitatively similar dynamics as structural glasses , and therefore it is of interest to understand their properties in more detail. Here we use the version in which each spin $`\sigma _i`$ can take one of $`q=10`$ different states ($`\sigma _i\{1,2,\mathrm{},10\}`$) and each spin interacts with every other one with an interaction $`J_{ij}`$. Thus the Hamiltonian is:

$$H=\frac{1}{2}\underset{i}{\overset{N}{}}\underset{ji}{\overset{N}{}}J_{ij}.(q\delta _{\sigma _i,\sigma _j}1)$$
(4)
The interactions $`J_{ij}`$ are drawn from a gaussian probability distribution with mean $`(3q)/(N1)`$ and variance $`(N1)^{1/2}`$. It has been show that in the thermodynamic limit this model has a dynamical singularity at $`T_d=1.14`$, slightly above the static singularity at $`T_c=1.13`$ . We have considered system sizes $`N`$ between 32 and 2560  but here we will discuss only results for $`N=320`$. In principle one has of course to average all the results obtained over the quenched disorder, i.e. the interactions $`J_{ij}`$. For the present case we have not done this (i.e. we consider only one realization of the disorder) but we have tested that the results presented below are independent of the choice of $`J_{ij}`$. The number of subsystems was $`M=16`$ and the exchange time $`\mathrm{\Delta }t_{PT}`$ was 10 Monte Carlo steps per spin. Below we will compare the results of the relaxation dynamics with the PT algorithm with the one of standard single spin flip Monte Carlo scheme. In the latter case we used the Metropolis criterion to accept or reject a move. More details on this investigation can be found in Ref.

## 4 Results

In this section we discuss the results, i.e. we compare the efficiency of the PT algorithm with the more conventional methods (standard molecular dynamics and Metropolis Monte Carlo) to propagate the system through configuration space. First we present our findings for the SiO2 system, then for the Lennard-Jones system, and finally for the Potts-glass. In all cases we first used the PT algorithm for a sufficiently long time to allow all subsystems to equilibrate. This equilibration was tested by comparing quantities like the energy or the specific heat obtained from such runs with the results from similar runs done with a standard method (molecular dynamics, single flip Monte Carlo).

From the setup of the PT algorithm it follows that each subsystem makes a random walk in temperature space. A rough estimate for the time needed until a new low-temperature configuration is produced can be obtained by looking at this random walk. For the case of SiO2 a typical example is shown in Fig. 1, where we plot the temperature of one subsystem, or rather the value of the coupling constant $`\mathrm{\Lambda }_i`$

(see (1)), as a function of time. The lowest temperature ($`p=1`$) is 2750K and the highest one ($`p=32`$) is 3922K. From this figure we recognize that it takes the system a bit more than 1000ps to generate a new independent configuration, i.e. to go from $`p=1`$ to $`p=M=32`$ and back. As we will see below this is significantly shorter than the time it takes in a standard molecular dynamics simulation. We also mention that the computational cost is of course also $`M`$ times (=number of subsystems) higher than in a conventional scheme. However, it has to be remembered that we also obtain $`M`$ independent configurations and not just one (this holds for each temperature!). Therefore if one needs many independent configurations, as it is usually the case in order to obtain reliable averages, the PT algorithm will pay off.

From this run we can calculate the distribution of the potential energy at the various temperatures. These distributions are shown in

Fig. 2. (In order to avoid overcrowding of the figure we show only every second temperature and the highest one.) We see that in order to obtain a sufficiently high acceptance ratio it is necessary to have a good overlap of neighboring distributions. A smaller overlap will lead to a smaller acceptance probability and hence the random walk will take longer. On the other hand a smaller value of $`M`$ will allow the random walk to go faster from low temperatures to high temperatures. It is presently not clear what the optimal choice is since this will depend on the details of the system. More results on this can be found in Ref. .

At the beginning of this section we have mentioned that one possibility to test whether or not the system has reached equilibrium is to compare various observables with the ones obtained with a conventional simulation method. This approach is, however, only possible for those temperatures at which it is feasible to equilibrate the system also with one of the latter methods. A different possibility is to use the data for the energy distribution for the different subsystems and to test whether

it is possible to reweight neighboring distributions to one common temperature . If we denote the distribution of the energy at the different coupling constants $`T_i=1/(\beta _0\lambda _i)`$ by

$$P_i(E)=P(E;T_i)$$
(5)
we should have the identity

$$P_i(E;T_j)=\frac{P_i(E)\mathrm{exp}[(\lambda _i\lambda _j)\beta _0E]}{𝑑E^{}P_i(E^{})\mathrm{exp}[(\lambda _i\lambda _j)\beta _0E^{}]}$$
(6)
for all pairs $`i`$ and $`j`$. Note that in general this identity holds only if the different subsystems are in equilibrium and hence it can be used to check whether the total system is in equilibrium or not.

That for the PT run such a reweighting does indeed lead to a nice collapse of the different $`P_i(E;T_j)`$ onto one master curve is shown in Fig. 3, where we plot these functions for the case $`j=16`$, which corresponds to $`T=3273`$K, for $`i=4,\mathrm{},15`$ and $`i=17,\mathrm{},27`$. Thus from this plot we have evidence that the system has indeed equilibrated within the time span of the simulation.

Having checked that the algorithm does indeed allow to equilibrate the system even at low temperatures it is of course important to see how efficiently this is done. One possibility to measure this is to calculate the mean squared displacement (MSD) of a tagged particle,

$$\mathrm{\Delta }r^2(t)=|𝐫_j(t)𝐫_j(0)|^2.$$
(7)
The time dependence of $`\mathrm{\Delta }r^2(t)`$ is shown in Fig. 4 for the case of silicon for various temperatures (dashed lines) 111Note that the MSD from the PT simulation was calculated by starting at configurations that were not yet in equilibrium. Thus the equilibrium MSD will be slightly different for times smaller than the equilibration time. We have checked, however, that at long times the shown curves are identical to the one in equilibrium.. Also included in the graph are the MSD obtained from a standard microcanonical run of the same system at similar temperatures (solid lines) . From this figure we recognize that at low temperatures and long times the MSD from the PT is larger by about a factor to 100, thus demonstrating that this type of dynamics is significantly faster than the conventional one.

The silica system we just considered was a glass former whose structure is given by an open tetrahedral network and whose dynamics at low temperatures shows an Arrhenius dependence . The Lennard-Jones system we consider next has very different properties in that its structure is rather similar to a dense packing of hard spheres and its dynamics shows a temperature dependence which is stronger than a Arrhenius law . Therefore it is of interest to see whether the PT method also works for this kind of system, which is done now.

In Fig. 5 we show the distribution of the energy for all 16 subsystems. From the curves we see that in this case it is possible to obtain much better data than for the silica system, since the numerical demand is quite a bit smaller in the former type of system (due to the short range nature of the interactions). That also in this case the PT dynamics is indeed able to equilibrate the system is demonstrated in Fig. 6a where we show the same distribution functions

reweighted to the temperature $`T=0.506`$. From this figure we see that all neighboring distributions collapse nicely onto the one for $`T=0.506`$, thus giving evidence that the system is indeed in equilibrium. In order to check whether this type of test is indeed sufficiently sensitive to detect whether or not the system is in equilibrium we have also made a run with standard molecular dynamics (at constant “temperature” 0.45) using as starting configuration the same as we used for the PT at that temperature. The length of this simulation was the same as the one for the PT and from earlier simulations it is known  that this time is not sufficient to equilibrate the system via conventional molecular dynamics. From Fig. 6b one sees that the curves stemming from the low temperatures do not fall on the master curve obtained from the higher temperatures. Thus we conclude that this way to analyze the data is indeed able to detect whether or not the system is in equilibrium and hence we have good evidence that Fig. 6a shows that the PT method has equilibrated the system.

Also in this case we judge the efficiency of the PT method by calculating the time dependence of the mean squared displacement of the A particles and compare it with the one obtained from standard molecular dynamics simulations. In Fig. 7 we show the two curves for the lowest temperature, $`T=0.45`$. It is recognized immediately that also for the Lennard-Jones system the PT algorithm leads to a much faster propagation of the system through configuration space in that at long times the MSD from the PT is around 100 times larger than the one from the standard molecular dynamics.

Finally we investigate the efficiency of the PT algorithm for the third type of system, the fully connected Potts-glass. We have done two PT simulations with $`T=0.9`$ and $`T=0.7`$ as the lowest temperatures. Also for this system we have checked that the distribution functions for the energy, $`P_i(E)`$, obtained in the different subsystems can be reweighted onto one master curve, which is again evidence that the PT algorithm does indeed equilibrate the system. In order to see how fast this is done we have calculated the normalized spin-autocorrelation function $`C(t)`$:

$$C(t)=\frac{1}{N(11/q)}\underset{i}{\overset{N}{}}(\delta _{\sigma _i(t)\sigma _i(0)}1/q).$$
(8)
The time dependence of this quantity is shown in Fig. 8 at the two mentioned temperatures. Using configurations which have been equilibrated with the PT algorithm as starting configurations, we also did conventional single spin flip Monte Carlo simulations at the same temperature. The resulting autocorrelation functions are included in the figure also. From the different curves we recognize that at the higher temperature the PT method leads to a relaxation which is more than one decade faster than the one of the standard Monte Carlo procedure. This factor has increased to more than 100 at the lower temperature and we see that at this temperature the equilibration of the system with the standard Monte Carlo method becomes hardly feasible. We have also checked that this type of speedup is typical in that it does not depend on the realization of the disorder, i.e. the bonds $`J_{ij}`$ in the Hamiltonian (4. Thus we conclude that also for this glass model the PT algorithm leads to a much quicker equilibration of the system than a standard simulation method.

## 5 Summary

In this article we have tested to what extent the parallel tempering algorithm is useful to equilibrate glassy systems. Whereas in the past it has been shown that this method works well for lattice systems, such as Ising spin systems , we now focused on off-lattice models. In particular we investigated a model for SiO2, a network glass-former whose transport properties at low temperatures show an Arrhenius dependence, and a binary Lennard-Jones system, a glass former which is structurally similar to a hard sphere system and whose temperature dependence of the relaxation times shows at low temperatures a strong deviation from the Arrhenius law. We have found that for both types of glass-formers the PT method leads to a significant ($`O(10^2)`$) acceleration of the relaxation dynamics and that this factor increases even more with increasing temperature. These results are confirmed by our simulation of a fully connected Potts-glass, a frustrated system which is believed to share many similarities with structural glasses.

It must be expected that there exists a temperature $`T_{g,PT}`$ below which also the PT method is not able to equilibrate the system, thus hindering one to investigate the equilibrium dynamics below this effective glass transition temperature. What the value of $`T_{g,PT}`$ is, how it depends on the type of system, the system size, or the PT exchange time $`t_{PT}`$, is presently not clear and thus has to be investigated in more detail. However, already now we see that this “new” algorithm (its main features dates back to 1986! ) will allow that computer simulations probe the static and dynamics properties of glassy systems at temperatures which are well below the ones that have been accessible so far.

We thank R. H. Swendsen for helpful discussions and gratefully acknowledge the financial support by DFG under SFB 262, the Grants in Aid for Scientific Research from the Ministry of Education, Science, Sports and Culture of Japan, the NIC in Jülich, and the Human Genome Center, Institute of Medical Science, University of Tokyo.


format: text

text:
# 1 Introduction

## 1 Introduction

The recent development of geometry has shown the importance of non-associative algebraic structures such as quasigroups, loops and odules. For instance, it is possible to say that the non associativity is the algebraic equivalent of the differential geometric concept of curvature. The corresponding construction may be described as follows. In a neighborhood of an arbitrary point on a manifold with an affine connection one can introduce the geodesic local loop which is uniquely defined by means of the parallel translation of geodesics along geodesics \[K, S1, S2\]. The family of local loops constructed in this way uniquely defines a space with affine connection, but not every family of geodesic loops on a manifold defines an affine connection. It is necessary to add some algebraic identities connecting loops in different points. Later, the additional algebraic structures (so called geoodular structures) were introduced and the equivalence of the categories of geoodular structures and of affine connections was shown by Sabinin \[S3, S4\]. The main algebraic structures arising in this approach are related to non associative algebra and theory of quasigroups and loops.

In our paper the new mathematical theory, non-associative geometry, which may help us to understand the discrete structure of spacetime, is introduced. The key point is to give an algebraic (non-local) description of manifold, which may be used in continuous and discrete cases. Non-associative geometry provides the unified algebraic description of continuous and discrete spacetime and admits all basic attributes of spacetime geometry including generalized Einstein’s equations.

## 2 What is the non-associative geometry?

Here we survey, in brief, algebraic foundations of non-associative geometry due to L.V. Sabinin (see on the matter \[S5, S6, S7, S8\] ).

###### Definition 2.1.

Let $`Q,\mathbf{}`$ be a groupoid with a binary operation $`(a,b)a\mathbf{}b`$ and $`Q`$ be a smooth manifold. Then $`Q,\mathbf{}`$ is called a quasigroup if the equations $`a\mathbf{}x=b,y\mathbf{}a=b`$ have a unique solutions: $`x=a\backslash b`$, $`y=b/a`$. A loop is a quasigroup with a two-sided identity, $`a\mathbf{}e=e\mathbf{}a=a,aQ`$. A loop $`Q,\mathbf{},e`$ with a smooth functions $`\varphi (a,b):=a\mathbf{}b`$ is called a smooth loop.

Let $`Q,\mathbf{},e`$ be a smooth local loop with a neutral element $`e`$. We define

$$L_ab=R_ba=a\mathbf{}b,l_{(a,b)}=L_{a\mathbf{}b}^1L_aL_b,$$
(1)
where $`L_a`$ is a left translation, $`R_b`$ is a right translation, $`l_{(a,b)}`$ is an associator.


text:
# The (in)stability of global monopoles revisited

## I Introduction

Global monopoles have been investigated for years as possible seeds for structure formation in the Universe . Although they appear to be ruled out by the latest cosmological data , their appearance in condensed matter –and other– systems and their peculiar properties make them worthy of investigation. These objects have divergent energy, due to the slow fall-off of angular gradients in the fields, which has to be cut-off at a certain distance R (in practice, the distance to the nearest monopole or antimonopole) and has two important consequences, in particular for cosmology. First, the evolution of a network of global monopoles is very different from that of gauged monopoles, as long-range interactions enhance annihilation to the extent of eliminating the overabundance problem altogether . Second, their gravitational properties include a deficit solid angle , which makes them rather exotic.

The stability of global $`O(3)`$ monopoles has been the subject of some debate in the literature . In this paper we try to settle the issue by: a) analysing the axial perturbation equations in the limit where the cut-off is taken to infinity, and b) proving that the energy barrier between the monopole and the vacuum (meaning the extra energy required by the monopole to reach an unstable configuration that decays to the vacuum) is finite. It is somewhat surprising for different topological sectors to be separated by finite energy barriers, but in this case it is a consequence of the scale invariance of gradient energy on two dimensional surfaces ($`r=`$ constant), and therefore independent of the details of the scalar potential.

## II The Model

We consider the simplest model that gives rise to global monopoles, the $`O(3)`$ model with lagrangian:

$$L=\frac{1}{2}_\mu \mathrm{\Phi }^a^\mu \mathrm{\Phi }^a\frac{1}{4}\lambda (|\mathrm{\Phi }|^2\eta ^2)^2a=1,2,3.$$
(1)
$`\mathrm{\Phi }^a`$ is a scalar triplet, $`|\mathrm{\Phi }|\sqrt{\mathrm{\Phi }^a\mathrm{\Phi }^a}`$ and $`\mu =0,1,2,3`$. The $`O(3)`$ symmetry is spontaneously broken to $`O(2)`$, leading to two Goldstone bosons and one scalar excitation with mass $`m_s=\sqrt{2\lambda }\eta `$. The set of ground states is the two-sphere $`|\mathrm{\Phi }|=\eta `$ and, since $`\pi _2(S^2)=𝐙`$, there are field configurations with non-trivial topological charge. One such configuration with unit winding is the spherically symmetric monopole,

$$\mathrm{\Phi }^a=\eta f(r)\frac{x^a}{r},$$
(2)
where $`f(0)=0`$ and $`f(r\mathrm{})=1`$. Its asymptotic behaviour is $`f(r0)\alpha r`$, $`\alpha 0.5`$ and $`f(r\mathrm{})11/r^23/2r^4`$, as can be seen from the e.o.m. of $`f(r)`$

$$f^{\prime \prime }+\frac{2}{r}f^{}\frac{2}{r^2}ff\left(f^21\right)=0.$$
(3)
The two parameters ($`\eta ,\lambda `$) appearing in the lagrangian can be absorbed by the rescaling $`\mathrm{\Phi }^a\stackrel{~}{\mathrm{\Phi }}^a=\mathrm{\Phi }^a/\eta ,x^\mu \stackrel{~}{x}^\mu =\sqrt{\lambda \eta ^2}x^\mu ,`$ which amounts to choosing $`\eta `$ as the unit of energy and the inverse scalar mass as the unit of length (up to a numerical factor). Note, however, that the energy of a configuration with non-trivial winding such as (2) is (linearly) divergent with radius, due to the slow fall-off of angular gradients and has to be cut off at $`r=R`$, say. Unlike $`\eta `$ and $`\lambda `$, the (rescaled) cut-off is an important parameter which could affect the dynamics of solutions with non-trivial topology. Dropping tildes:

$$E=_0^R\frac{1}{2}_\mu \mathrm{\Phi }^a^\mu \mathrm{\Phi }^a+\frac{1}{4}(|\mathrm{\Phi }|^21)^2.$$
(4)
Since the energy diverges, Derrick’s theorem does not apply in this case, and in  it was shown that the global monopole is stable towards radial rescalings. On the other hand the question of stability with respect to angular perturbations has led to some discussion in the literature after Goldhaber  pointed out that the ansatz

$`\mathrm{\Phi }^1`$ $`=`$ $`F(r,\theta )\mathrm{sin}\overline{\theta }(r,\theta )\mathrm{cos}\phi `$ (5)
$`\mathrm{\Phi }^2`$ $`=`$ $`F(r,\theta )\mathrm{sin}\overline{\theta }(r,\theta )\mathrm{sin}\phi `$ (6)
$`\mathrm{\Phi }^3`$ $`=`$ $`F(r,\theta )\mathrm{cos}\overline{\theta }(r,\theta ),`$ (7)
which describes axially symmetric deformations of the spherical monopole (2), leads to the following expression for the energy after a change of variables $`y=\mathrm{ln}\mathrm{tan}(\theta /2`$):

$`E={\displaystyle _0^{2\pi }}d\phi {\displaystyle _{\mathrm{}}^{\mathrm{}}}dy{\displaystyle _0^R}dr{\displaystyle \frac{1}{2}}\left[\rho _1+r^2\mathrm{sech}^2(y)\rho _2\right]`$ (8)
$`\rho _1=F_y^2+F^2[\mathrm{sin}^2(\overline{\theta })+(\overline{\theta }_y)^2]`$ (9)
$`\rho _2=F_r^2+F^2(\overline{\theta }_r)^2+{\displaystyle \frac{1}{2}}(F^21)^2`$ (10)
with $`F_r_rF`$, etc. The term in brackets in $`\rho _1`$ is identical to the energy of the sine-Gordon soliton, so translational invariance in $`y`$ implies that configurations with

$$F(r,y)=f(r),\mathrm{tan}\left(\overline{\theta }(r,y)/2\right)=e^{y+\xi },\xi =\mathrm{const}$$
(11)
have the same energy as (2) (which corresponds to $`\xi =0`$). On a given $`r=`$const. shell, the effect of taking $`\xi \mathrm{}`$ is to concentrate the angular gradients in an arbitrarily small region around the north pole. When the gradient energy is inside a region of size comparable to the inverse scalar mass, it is energetically favourable to “undo the knot” by reducing the modulus of the scalar field to zero and climbing over the top of the mexican hat potential. Unwinding is estimated to occur at a critical value of $`\xi `$ (say $`\xi _0`$) whose dependence with $`r`$ far from the core is logarithmic, $`\xi _0\mathrm{ln}r+\mathrm{const}.`$ Numerical simulation on individual shells closer to the core gives $`\xi _0a_0+b_0\mathrm{ln}r`$ with slowly varying $`b_01`$ and $`a_01.3`$. Unwinding is expected if $`\xi (r)>\xi _0(r)`$.

## III Estimation of the energy barrier

As explained in , the shift (11) creates a tension pulling the monopole core, and the apparent unwinding (which starts in the inner shells) is only a manifestation of the core’s translation. In order to stop the motion of the core, we consider a hybrid configuration such that the monopole core remain unperturbed and the unwinding occur in the outer shells. This is achieved by taking $`\xi =0`$ for, say, $`r<r_1`$, a string-like configuration for $`r>r_2`$, and some continuous interpolation in between. One such configuration (see figure 1a) would be:

$$\widehat{F}(r,y)=f(r),\widehat{\xi }(r)=\{\begin{array}{cc}0\hfill & r<r_1\hfill \\ c\left(1\frac{r_1}{r}\right)\hfill & r_1<r<r_2\hfill \\ a+b\mathrm{ln}(r)\hfill & r_2<r\hfill \end{array}$$
(12)
with $`a>a_0`$, $`b>b_0`$, and $`c`$ given by continuity of $`\widehat{\xi }(r)`$.

If the energy (i.e. mass) inside $`r_2`$ is large enough there will be no appreciable motion of the monopole core (because the tension of the string is constant $`4\pi `$). For instance, simulations with $`r_1=3`$, $`r_2=6`$, $`b=1`$ show different behaviour depending on $`c`$: for $`c<\mathrm{\hspace{0.33em}2.8}`$ the core translates, but for $`c>\mathrm{\hspace{0.33em}2.8}`$, the unwinding happens far from the core, which remains fixed. In the latter case, the decay of the string can also be understood as a monopole-antimonopole pair creation with the new monopole appearing at $`r=R`$ and the antimonopole appearing near $`r=r_2`$. The equations solved in the simulations were

$$\mathrm{}\mathrm{\Phi }^a+\mathrm{\Phi }^a(|\mathrm{\Phi }|^21)+\gamma \dot{\mathrm{\Phi }}^a=0,$$
(13)
with a dissipative term added to make the integration faster. Simulations using different $`\gamma `$ show no appreciable difference. The equations were integrated using cylindrical coordinates ($`\rho ,\phi ,z`$) in a $`200^2`$ grid using explicit Runge-Kutta method with step size control . The output can be seen in fig. 1b, where the potential energy is plotted at various times, confirming that the configuration (12) is unstable and decays to the vacuum.

We will now show that the extra energy required to reach this unstable configuration (12) from the spherical monopole configuration (2) is finite.

Consider the set of configurations with $`F=f(r)`$ and $`\xi =\xi (r)`$. From (10) the difference between the energy of any such configuration ($`E[\xi ]`$) and the energy of the spherically symmetric unperturbed monopole ($`E[0]`$) is:

$`E[\xi ]E[0]`$ $`=`$ $`\left[{\displaystyle _0^r_{}}+{\displaystyle _r_{}^R}\right]\mathrm{d}r\left[{\displaystyle \frac{r^2}{2}}f^2(r)\xi _r^2(r)I[\xi (r)]\right]`$ (14)
$`I[\xi ]`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}dy\mathrm{sech}^2(y)\mathrm{sech}^2(y+\xi (r)).`$ (15)
$`I[\xi ]`$ is bell-shaped: from a maximum value $`I[0]=4/3`$, it rapidly falls to zero for $`|\xi |>\xi _{}5`$ as $`16(|\xi |1)exp(2|\xi |)`$. $`r_{}`$ is the radius at which $`\xi (r_{})=\xi _{}`$. The first integral is clearly finite. The second can be estimated using the asymptotic form of $`I[\xi ]`$ and is negligible for $`\xi =a+\mathrm{ln}r`$ even in the limit $`R\mathrm{}`$, where one gets $`16\pi e^a\xi _{}e^\xi _{}`$.

Moreover, there is a continuous path in configuration space connecting the configurations (2) and (12) such that $`E[\xi ]E[0]`$ remains finite along the whole path: first increase $`\xi `$ in the outer shells $`r>r_{}`$ using Goldhaber’s deformation ($`\xi _r=0`$) until $`\xi `$ reaches $`\xi _{}`$; then adjust the radial dependence to match (12).

Thus, we have shown that the extra energy required by the monopole (2) to go over the energy barrier and decay to the vacuum is finite even as $`R\mathrm{}`$; moreover, since the monopole energy grows with R, the ratio $`(E[\widehat{\xi }]E[0])/E[0]0`$ as $`R\mathrm{}`$.

## IV Stability to small perturbations with axial symmetry

We now turn to the classical stability of (2) by considering small perturbations parametrised by

$`F=f(r)+\delta (t,r,y)`$ (16)
$`\mathrm{tan}(\overline{\theta }/2)=(1+\xi (t,r,y))e^y`$ $`.`$ (17)
Neglecting quadratic terms gives $`\mathrm{sin}\overline{\theta }\mathrm{sin}\theta +\xi \mathrm{sin}\theta \mathrm{cos}\theta `$, $`\mathrm{cos}\overline{\theta }\mathrm{cos}\theta \xi \mathrm{sin}^2\theta `$. Introducing $`Xf\xi `$,

$`\mathrm{\Phi }^1`$ $``$ $`(f\mathrm{sin}\theta +\delta \mathrm{sin}\theta +X\mathrm{sin}\theta \mathrm{cos}\theta )\mathrm{cos}\phi `$ (18)
$`\mathrm{\Phi }^2`$ $``$ $`(f\mathrm{sin}\theta +\delta \mathrm{sin}\theta +X\mathrm{sin}\theta \mathrm{cos}\theta )\mathrm{sin}\phi `$ (19)
$`\mathrm{\Phi }^3`$ $``$ $`(f\mathrm{cos}\theta +\delta \mathrm{cos}\theta X\mathrm{sin}^2\theta ),`$ (20)
which shows that the correct boundary conditions are that $`\mathrm{sin}\theta (\delta +X\mathrm{cos}\theta )`$ should vanish on the $`z`$ axis. Note that $`\delta (0)`$ and $`X(0)`$ need not vanish. An infinitesimal translation of the monopole in the $`z`$ direction corresponds to $`\delta =f_r(r)\mathrm{cos}\theta ,X=f(r)/r`$, and both $`f_r(r)`$ and $`f(r)/r`$ tend to $`0.5`$ as $`r0`$ (see fig 2a ). There is no zero mode associated with global rotations since these have been factored out in the ansatz (7).

As usual, the perturbation equations

$`0=r^2\mathrm{}\delta +2\delta +r^2[3f^21]\delta +2X_y4X\mathrm{tanh}y`$ (21)
$`0=\mathrm{sech}^2yr^2\mathrm{}X+\mathrm{sech}^2y[r^2(f^21)+2]X+2\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}yX_y2\delta _y`$ (22)
reduce to an eigenvalue problem in $`\omega ^2`$ for perturbations of the form $`\delta =e^{i\omega t}\widehat{\delta }(r,y)`$, $`X=e^{i\omega t}\widehat{X}(r,y)`$. Eigenfunctions $`\widehat{\delta },\widehat{X}`$ with negative eigenvalues $`\omega ^2<0`$ correspond to instabilities. Dropping hats and defining $`u=\mathrm{tanh}y`$

$`R_1\delta +_u[(u^21)_u\delta ]2_u[(u^21)X]=\omega ^2r^2\delta `$ (23)
$`R_2X+_u^2[(u^21)X]2_u\delta =\omega ^2r^2X`$ , (24)
where $`R_1`$, $`R_2`$ are radial operators (see fig 2b):

$`R_i=_r(r^2_r)+V_i(r),`$ $`i=1,2`$ (25)
$`V_1(r)=r^2(3f^2(r)1)+2`$ (26)
$`V_2(r)=r^2(f^2(r)1)`$ $`.`$ (27)
Using Legendre polynomials and changing variables: $`\chi :=_u[(1u^2)X]=\chi _l(r)P_l(u),\delta =\delta _l(r)P_l(u)`$ the equations for different values of $`l`$ decouple, giving

$`R_1\mathrm{\Delta }_l+x^2\mathrm{\Delta }_l+2x\chi _l`$ $`=\omega ^2r^2\mathrm{\Delta }_l`$ (28)
$`R_2\chi _l+x^2\chi _l+2x\mathrm{\Delta }_l`$ $`=\omega ^2r^2\chi _l`$ (29)
where we introduced $`x=\sqrt{l(l+1)}`$ and $`\mathrm{\Delta }_l=x\delta _l`$. In order to get (29) we multiplied the $`X`$ eqn. in (24) by $`(1u^2)`$ and differentiated w.r.t. $`u`$, so there may be spureous solutions; in particular, $`l=0`$ corresponds to angular perturbations that are singular on the $`z`$ axis. These are not physical, and will be discarded. But if there is no solution of (28,29) with negative $`\omega ^2`$ there will be no instability in the original problem (22).

Our task is to find the solution to (28,29) with the minimum value of $`\omega ^2`$ over all admissible perturbations and all $`l1`$. We know one $`l=1`$ solution, the translational zero mode ($`\stackrel{~}{\mathrm{\Delta }}_1=\sqrt{2}f_r(r)`$, $`\stackrel{~}{\chi }_1=2f(r)/r`$). Goldhaber’s deformation ($`\chi _1=f(r)`$, $`\mathrm{\Delta }_1=0`$) is also $`l=1`$ but is not a solution of (28, 29), and it can be shown that there are no instabilities with $`\chi _1(r\mathrm{})\mathrm{const}`$.

Let us first consider normalizable perturbations. Note that, for each $`l`$, the equations (28,29) can be obtained by functional variation from

$$E_ldr\left[r^2\left(\mathrm{\Delta }_r^2+\chi _r^2\right)+(V_1+x^2)\mathrm{\Delta }^2+(V_2+x^2)\chi ^2+4x\mathrm{\Delta }\chi \right]=\omega ^2drr^2[\mathrm{\Delta }^2+\chi ^2].$$
(30)
The lowest value of $`\omega ^2`$ can be found minimising $`E_l`$ over normalized functions, $`drr^2[\mathrm{\Delta }^2+\chi ^2]=1`$ and over all $`l1`$. However the minimum must be in the $`l=1`$ sector since, for all $`l>1`$ and for given $`\mathrm{\Delta }`$,$`\chi `$, $`(E_lE_1)`$ is a sum of squares with positive coefficients:

$$E_lE_1=dr[A_{l,+}(\mathrm{\Delta }+\chi )^2+A_{l,}(\mathrm{\Delta }\chi )^2],$$
(31)
where $`A_{l,\pm }[x^22\pm 2(x\sqrt{2})]/2>0`$ for $`l>1`$.

In order to investigate the $`l=1`$ sector, we rewrite $`E_1`$ using arbitrary functions $`G(r)`$, $`H(r)`$ and $`K(r)`$:

$`E_1={\displaystyle _0^{\mathrm{}}}\mathrm{d}r[(r\mathrm{\Delta }_r+{\displaystyle \frac{G}{r}}\mathrm{\Delta })^2+(r\chi _r+{\displaystyle \frac{H}{r}}\chi )^2+2\sqrt{2}({\displaystyle \frac{\mathrm{\Delta }}{K}}+K\chi )^2+(V_1+2{\displaystyle \frac{2\sqrt{2}}{K^2}}+G_r{\displaystyle \frac{G^2}{r^2}})\chi ^2`$ (32)
$`+(V_2+22\sqrt{2}K^2+H_r{\displaystyle \frac{H^2}{r^2}})\mathrm{\Delta }^2]G(r)\mathrm{\Delta }^2(r)|_0^{\mathrm{}}+H(r)\chi ^2(r)|_0^{\mathrm{}}`$ $`.`$ (33)
Choosing $`G(r)=r^2_r\mathrm{ln}\stackrel{~}{\mathrm{\Delta }}_1`$, $`H(r)=r^2_r\mathrm{ln}\stackrel{~}{\chi }_1`$, $`K(r)^2=\stackrel{~}{\mathrm{\Delta }}_1/\stackrel{~}{\chi }_1`$, the coefficients of $`\mathrm{\Delta }^2`$ and $`\chi ^2`$ in (33) vanish identically by virtue of (3). For large r, $`G`$ and $`Hr`$, thus, $`E_10`$ for all functions ($`\mathrm{\Delta }_1`$, $`\chi _1`$) that decay faster than $`1/\sqrt{r}`$ as $`r\mathrm{}`$. This proves that all normalizable perturbations have $`\omega ^20`$. The above argument and a host of numerical simulations strongly suggest that non-normalizable perturbations are at best compatible with $`\omega ^2=0`$, as can be verified directly from the equations for perturbations that fall to zero like a power of r, but in this case we have no analytic proof.

## V Discussion

In this paper we have derived and analysed the axial perturbation equations of $`O(3)`$ monopoles and proved that, contrary to statements in the literature, $`O(3)`$ monopoles are perturbatively stable (or neutrally stable) to infinitesimal, axially symmetric, normalizable (or power-law decay) perturbations. We have also proved that the energy barrier between topological sectors is finite, irrespective of the details of the scalar potential. This feature is specific of global monopoles; global vortices in two dimensions, whose energy grows logarithmically with radius, and gauge monopoles, whose energy is finite, are separated from the vacuum by an energy barrier growing linearly with the cut-off. But the energy of global monopoles is dominated by two-dimensional (angular) gradients far from the core, and these can be deformed with no energy cost.

One would naively expect thermal fluctuations with $`KTE[\xi ]E[0]`$ to cause the monopole to decay. As far as we know, this effect is not seen in “cosmological” numerical simulations of global monopole networks, but perhaps this is not so surprising. First, we worked in flat space. Second, the range of scales introduced by the expansion of the Universe forces drastic approximations on cosmological simulations; in particular, the sigma model approximation, which is widely used, sets the field on the vacuum manifold everywhere, so unwinding and pair creation events are not resolved by the grid. Finally, cosmological simulations do not include thermal effects, and it has been shown that full thermal simulation across the phase transition  gives qualitatively different results. Our results provide further evidence that a more careful analysis of global monopole networks may be required.

## VI Acknowledgments

We thank M. J. Esteban for pointing out an error in an earlier version of this paper, and M. Escobedo, R. Durrer, R. Gregory, A. Vilenkin, K. Kuijken, I.L. Egusquiza, J.M. Aguirregabiria, L. Perivolaropoulos and A. Lande for conversations. We acknowledge support from grants CICYT AEN99-0315 and UPV 063.310-EB225/95.


text:
# 1 Introduction

## 1 Introduction

When a sunspot group appears on the surface of the Sun, it reduce the brightness locally. Does this lead to an actual reduction of the Sun’s brightness as measured on earth (the so-called irradiance)? Or is the ‘missing flux’ re-emitted elsewhere on the solar surface? Faculae are bright areas, visible especially near the limb of the Sun. Do they contribute to the irradiance? Does the irradiance vary only because of spots and faculae, or are there other effects as well, and can these be large enough to influence climate on earth?

Some of these questions have been answered by actual measurements with space-base radiometers over the past 20 years (figure 1). The data show a clear variation, obviously correlated with the solar cycle: the solar magnetic field does indeed cause the Sun’s brightness to vary a bit, by about 0.1% between minimum and maximum of the cycle. The record shows prominent sharp dips; these are associated with the passage of individual sunspot groups across the disk. When a large group appears, the Sun actually does become dimmer by a tiny bit. The smoothed data (thick line in the figure) shows, however, that the average trend is the opposite: at maximum activity, the Sun is brighter.

The level of variation seen in figure 1 is not large, and by general opinion in the climate modeling community, insufficient to cause variations in the earth’s climate, by a factor of 10 or so (see the contributions in Nesme-Ribes et al., 1994, Friis-Christensen et al. 2000). Even if the Sun were to settle into a state of activity as seen during the minima in figure 1, the average reduction of irradiance by 0.05% would have negligible effects, and probably would not explain events like the so-called ‘little ace age’111A period during the 17th century when the climate in Europe was colder by a degree C or so. This cold spell appears to have been limited mostly to Western Europe, however..

For climate modeling, the record of figure 1 is too short. One would like to know what the irradiance has been in the past, during the last 10000 years, for example, so comparisons can be made with climate records such as have been obtained from ice cores in Greenland and the Antarctic (see the contributions Friis-Christensen et al. 2000). Is it conceivable, for example, that figure 1 shows only a fraction of what the solar magnetic field is capable of doing, and that variations of, say, 10 time larger amplitude could have been caused by the solar magnetic field, or something else that has not yet been discovered on the Sun? These are questions that can excite the imagination, including that of funding agencies.

In the absence of reliable independent indicators of the past brightness of the Sun, it is of some importance to understand how the variation seen in figure 1 come about. If theoretical understanding is sufficient, and accepted by the community, it may be possible to assess with more certainty whether irradiance variations larger than those seen so far are likely to have occurred in the past, and perhaps influence the climate.

Until the advent of helioseismology (see Dziembowski, elsewhere in this volume), our knowledge of the internal structure of the Sun was based on the theory of stellar structure and evolution. Since this theory is reasonably simple, confidence in the models was high even before helioseismological tests became possible. The main uncertainty had to do with the exact values of the opacity of stellar plasma. Increasingly accurate values for this were computed from first principles (atomic physics) by extensive numerical work. The basic correctness of the theory was then confirmed by helioseismology, to such an extent even that one is now confident that the ‘solar neutrino problem’ has its cause in neutrino physics, not inaccuracies in the solar models (e.g. Bahcall 1996, Richard et al. 1996, Brun et al. 1998).

Though these solar models are spherical and only slowly evolving in time, the theory behind them also apply to time dependent and non-spherical disturbances. In this form, the theory makes definite predictions about the thermal behavior of the Sun, that is, the evolution in time and position in the star of perturbations that might be related to the observed irradiance variations. In the following I present my view how this theory can be used to locate the most plausible source of irradiance variations. After this, I discuss how the new heliosmological measurements currently being made can be used to test the predictions of proposed mechanisms of irradiance variation.

Conceptually, one can separate the thermal effects of magnetic fields into three types:

1. ‘Sources and sinks’. The generation of a magnetic field involves the conversion of energy of motion into magnetic energy. Since the motions in the solar envelope are thermally driven, this ultimately means conversion of thermal into magnetic energy: building up a magnetic field produces a thermal sink somewhere. The opposite happens when the field decays: magnetic energy is converted into heat. These thermal effects exist only during changes in the magnetic energy content of the envelope.

2. ‘Shadows’, or changes in the heat transport coefficient. Magnetic fields interfere with convection, causing a reduction in the efficiency of heat transport in the envelope. In contrast to (1), these changes last as long as the magnetic field itself is present.

Both these sources and sinks, and changes in the transport coefficient cause thermal perturbations, varying with magnetic activity, which propagate through the envelope and cause variations in surface energy flux. Related to the second class of perturbations are:

3. The effects of magnetic fields at the surface of the star. Sunspots, being dark, radiate less than the surrounding photosphere, while the small elements that make up plages and the network have an excess emission. In addition, it is conceivable that the magnetic elements have an indirect effect surface by modifying the convective flow in their surroundings slightly (for which there is some observational evidence, see section 3.6).

## 2 Time scales

The response of the Solar envelope to thermal perturbations is not governed by a single time scale, but by a wide range of time scales. The longest of these is the thermal time scale of the Sun as a whole, called the Kelvin-Helmholtz time scale. It is about $`10^7`$ years, namely the thermal energy content of the Sun, of the order $`10^{41}`$ Joule, divided by its luminosity of $`\mathrm{4\hspace{0.17em}10}^{26}`$ W. Hence, on time scales of human interest, the Sun has an extremely large thermal inertia. If the central heat source of the Sun were switched off, the internal structure and the luminosity would start to change only on this long time scale.

The Kelvin-Helmhotz time scale is a global time scale for the Sun as a whole; more generally, we can define the thermal time scale $`\tau _\mathrm{t}`$ as a function of depth:

$$\tau _\mathrm{t}(z)U(z)/L(z)\frac{1}{L}_{Rz}^R4\pi r^2udr,$$
(1)
where $`L`$ is the luminosity at depth $`z`$, $`U`$ the thermal energy of the envelope down to a depth $`z`$, and $`u`$ the thermal energy per unit volume, approximately (for an ideal gas of constant $`\gamma `$) given by $`u=P/(\gamma 1)`$. This is the time scale on which the structure of the envelope, and the observed luminosity, would start changing when the heat flux in the star were interrupted, by some magical means, at depth $`z`$. Some rough values for this quantity are $`\tau _\mathrm{t}10^5`$yr at $`z=\mathrm{2\hspace{0.17em}10}^5`$km (depth of the convection zone), 10 yr at $`\mathrm{20\hspace{0.17em}000}`$km (the size of a supergranule), 10 hrs at $`2000`$km (size of a granule). This shows that the thermal time scale depends rather strongly on depth in the Sun. As a result, the thermal response of the Sun also depends critically on the location of the disturbance. Near-surface disturbances have much larger effects than ones orginating near the base of the convection zone.

A second kind of time scale involved in thermal readjustments is the diffusive time scale. This is the time scale on which differences in entropy between different parts of the convection zone are ironed out by the convective flows. This is easily computed in a ‘mean field’ or diffusion model for convection, in which the local convective heat flux is assumed to be a function of only the local conditions (thermodynamic state of the gas and the temperature gradient). The validity of this model is questionable, since it is now known that stellar envelope convection is in fact extremely non-local, with almost all flows generated by cooling at the stellar surface (e.g. Nordlund 1986, Stein and Nordlund 1989, for reviews see Nordlund and Stein 1996, Spruit 1997). The diffusion model has the advantage, however, that it is simple to apply, and is (unfortunately) still almost universally accepted. For these historical and practical reasons, I still use it for the estimates in this section. Estimates in the context of the new nonlocal view of envelope convection have not yet been developed to the same level of detail. I return to this in section 4, where the examples of heat flux blocking by sunspots and ‘thermal shadows’ are considered in the nonlocal convection picture. I argue there that the new nonlocal view only strengthens the conclusions based on the diffusion model.

For small perturbations, the heat flux in the local model is proportional to the entropy gradient, and time dependent perturbations are governed by a diffusion equation. The ‘turbulent’ diffusion coefficient in this equation is $`\kappa _\mathrm{t}\frac{1}{3}l_\mathrm{c}v_\mathrm{c}`$ where $`l_\mathrm{c}`$ and $`v_\mathrm{c}`$ are the convective length scale and velocity. In mixing length models of the solar convective envelope, this quantity varies only weakly, at a value of the order $`10^{13}`$ cm2/s. In this model thermal inhomogeneities (more precisely: entropy inhomogeneities) of length scale $`d`$ are smoothed by turbulent diffusion on a time scale

$$\tau _\mathrm{d}=d^2/\kappa _\mathrm{t}.$$
(2)
For $`d=\mathrm{2\hspace{0.17em}10}^5`$km this is about 1 yr, for $`d=2000`$km about 1 hr. Comparing $`\tau _\mathrm{t}`$ and $`\tau _\mathrm{d}`$, we see that they are of similar magnitude close to the surface (to be precise: in the surface boundary layer where convection is not efficient enough to keep the stratification close to adiabatic). In deeper layers, the thermal time scale is much longer than the diffusive time scale, by a factor of up to $`10^5`$.

The two time scales measure different types of thermal adjustment process. These same processes appear in the thermal behavior of, say, a chunk of metal, aluminum for example, heated from the inside and suspended in space. The thermal time scale is the time scale on which its temperature adjusts to a change in the heat input, such that the heat radiated from the surface into space balances the heat input again. It is determined by the heat capacity \[$`U`$ in eq. (1)\] and the power level ($`L`$). The time scale on which different parts of the chunk equilibrate to the same temperature is governed by a different process, namely thermal conduction (the equivalent of the turbulent diffusion in the Sun). The diffusion (conduction) time scale is much shorter than the thermal time scale, because of the large heat conductivity of Al. In the Sun, it is the very large turbulent diffusivity in the bulk of the convection zone that causes the very short diffusive time scale compared with the thermal time scale.

How do these different time scales come into play when the convection zone is thermally perturbed by, say, the storage of energy in a growing magnetic field? Such perturbations can be computed in detail, either by numerical methods (Endal et al. 1985) or more analytically. We can, for example, consider the initial value problem in which a perturbation is allowed to evolve in time by heat transport in the convection zone. In general this evolution has components on all the time scales of the problem, including the very long thermal time scale. Detailed calculations of this problem are discussed briefly in the next subsection (for a more complete analysis, see Spruit, 1982ab, 1991, Arendt, 1992).

### 2.1 Modes of thermal relaxation

In the mixing length approximation for convective energy transport the heat flux is given by:

$$𝐅=\kappa _\mathrm{t}\rho TS,$$
(3)
where $`𝐅`$ is the heat flux, $`\kappa _\mathrm{t}`$ the turbulent diffusivity, $`\rho `$ the density, $`T`$ the temperature and $`S`$ the entropy. The energy equation, in the absence of sources, is given in terms of the entropy by

$$\rho T\frac{\mathrm{d}S}{\mathrm{d}t}=\mathrm{div}𝐅.$$
(4)
As an example, consider the one dimensional problem in which all quantities depend on depth $`z`$ only. In this case, combining (3) and (4), and neglecting flows we get

$$\frac{S}{t}=\frac{F}{\rho TH}+\kappa \frac{S}{z^2},$$
(5)
where $`H=\mathrm{ln}\rho T/z^1`$ is the pressure scale height. For perturbations in which the left hand side is balanced by the first term on the right, the time scale is of the order $`F/(c_p\rho TH)`$, which is just a thermal time scale (cf. 1). For modes in which the left hand side is balanced by the second term, we evidently have an ordinary diffusion problem, hence these modes evolve on diffusive time scales. The first term has appeared as a consequence of the stratification of pressure, ie. the inhomogeneity of the convection zone. In a homogeneous medium (as in a chunk of Al), the thermal time scale comes in only through the radiating surface boundary condition (which one could also consider as a form of inhomogeneity). If $`T_s`$ is the surface temperature, this condition is

$$\sigma T_s^4=F.$$
(6)
Since the perturbations in entropy outside spots and the magnetic elements making up the small scale field are small, the reponse of the convection zone can be computed by linearizing eqs. (3),(4),(6). An arbitrary initial perturbation may then be decomposed in the standard way in terms of the eigenmodes of the problem. These are found by setting $`\delta S=e^{\eta t}f(𝐫)`$, and solving for $`f(𝐫)`$ with the decay rate $`\eta `$ as eigenvalue (the problem thus defined has only decaying solutions). Specializing to the one dimensional case ($`f(z)`$) one finds a series of modes ordered by the number of nodes $`n`$ of $`f`$. The fundamental ($`n=0`$) has a decay rate of the order $`\eta _0\tau _t^1`$ that is, this mode decays on the thermal time scale. The higher modes decay with rates of the order $`\eta _nn^2\tau _d^1`$, they decay on diffusive time scales.

The extreme separation of time scales, $`10^5`$yr for the thermal time scale and 1yr for the longest diffusive time scale (with $`n=1`$) leads to a simple picture. The high heat conductivity tends to keep the convection zone entropy uniform, so that sources and sinks of energy are shared by a large part of the convection zone. This large part of the convection zone has a large heat capacity, however, so that sources and sinks tend to spend their efforts in slow heatings and coolings of large amounts of mass, rather than propagating to the surface.

How do these different time scales come into play when the convection zone is thermally perturbed by, say, the storage of energy in a growing magnetic field? The evolution of such perturbations can be computed in detail, either by numerical methods (Endal et al. 1985, Gilliland, 1988) or more analytically. We can, for example, consider the initial value problem in which a perturbation is allowed to evolve in time by thermal transport in the convection zone. In general this evolution has components on all the time scales of the problem, including the very long thermal time scale. Formal aspects of this problem have been discussed elsewhere (Spruit 1982ab, 1991, Arendt 1992). In the following, the basic conclusions of these analyses are summarized.

## 3 Expected level of luminosity variations

For quantitative estimates, the strength of the field and its filling factor in the convection zone have to be specified. Assume that we have a layer of field with strength of the order of $`\mathrm{10\hspace{0.17em}000}`$G (equipartition with the convective flows as estimated by a mixing length model), one scale height deep, near the base of the convection zone (where most of the magnetic flux is probably located). I summarize here some results, discussed in greater detail elsewhere (Spruit 1991).

### 3.1 sources-and-sinks

If the energy needed to build up the solar magnetic field during one half of the solar cycle is taken out of the thermal energy near the base of the convection zone, the calculations show that a surface luminosity variation of only $`\delta L/L10^7`$ results. This is due to the very large heat capacity of the lower convection zone. The effect is stronger if the source of the magnetic field is assumed to be closer to the surface, but is still much smaller than the observed effect. Recent models for the emergence of magnetic flux from the base of the convection zone (D’Silva and Choudhuri 1993, D’Silva and Howard 1993, Caligari et al. 1998) indicate that the actual field strength at the base of the convection zone is probably about 10 times higher than the equipartition estimate, of the order $`10^5`$ G. The magnetic energy per unit of magnetic flux is then also 10 times higher, but this still does not lead to a significant luminosity effect.

### 3.2 shadows

If magnetic fields interfere with convection, a thermal perturbation develops as well. If at some depth below the surface there is an area where magnetic fields reduce the convective flows, the surroundings will be affected thermally, and one might hope to see a ‘thermal shadow’ at the surface. This can again be computed in the diffusion assumption for convection (Spruit, 1977). The amplitude of the effect depends crucially on the ‘covering factor’. If a reduction of convective efficiency is assumed that uniformly covers a horizontal surface at some depth $`z`$, magnetic fields comparable to equipartition with convective flows can have stronger effects than the source-and-sink perturbations (Gilliland, 1988). A measurable effect, however, is predicted only if the field is located close below the surface. At face value, this would make changes in convective efficiency a reasonable candidate for irradiance variations.

In practice, the effect is very strongly reduced, however, if the change in convective efficiency does not extend over an entire horizontal surface. If there are ‘holes’ in this cover, the effect is much smaller. This is because the turbulent heat conductivity in the convection zone is so high that the heat flux is easily ‘shunted’ past blocking objects below the surface (Spruit, 1977). Compare this with efforts to thermally insulate homes. To be effective, it is not enough that most of the building’s outer surface is very well insulated. A more modest insulation but without any leaks is more effective. For this reason, thermal shadows of subsurface fields are probably not very important (see also section 4.3).

### 3.3 surface effects

By far the most effective way in which a magnetic field influences the irradiance is by its effect on the net surface emissivity. The reduced emission from a sunspot area shows up directly in the irradiance records as a dip tracking the passage of the area across the disk. The only complication is that one might expect, depending on one’s prejudices, that part of this reduction could be compensated by a brightening elsewhere. For example in the form of a ‘bright ring’ surrounding the spot. Evidence of such bright rings is absent for most spots that have been studied for this effect, and where observation of a ring is reported, it compensates for only a small fraction of the spot blocking (e.g. Rast et al., 1999). Brightening of the photosphere around a spot is hard to measure because it is usually masked by facular emission, which is a different effect. Such faculae are part of all active regions. Though their excess emission is rougly of the same order as the spot deficit, this is a coincidence that does not hold any more on closer inspection. For example, they are present in the same amount in active regions with or without sunspots, and their contribution to the irradiance variation does not correlate very well with spot areas (Foukal and Lean 1986, Lean et al. 1998).

The absence of bright rings is understood in terms of a turbulent diffusion model for the heat flux in the convection zone (Spruit 1977, 1982b, Foukal et al. 1983, Chiang et al. 1984). The ‘blocked heat flux’, for the most part, does not reappear elsewhere on the surface, but stays inside the convection zone, being stored/released on the very long thermal time scale of the convective envelope. This conclusion holds, in the diffusion model, as long as the blocking effect of the spot extends to a depth of at least 1000 km, a mild requirement given that the observed Wilson depression of the umbra of a spot is already of the order of 500km.

### 3.4 spot blocking: stove top analogy

The absence of bright rings in the diffusion model can be understood by a kitchen analogy (for an improved, but more qualitative, model based on the new nonlocal view of convection, see section 4.1). Imagine a stove top with an electric heater plate. Assume this plate to be made of a massive block of a metal with a high thermal conductivity, Aluminum or Copper, say. Such a block, like the solar convection zone, has two time scales, the thermal time scale, given by the heat content divided by the heating rate, and the thermal conduction time scale given by the thermal diffusion time across the block. These are different physical quantities and their values can be quite different. If the thermal conductivity is large, the diffusion time scale can be much shorter than the thermal time scale. In the solar convection zone this is the same except that the ratio of time scales ($`10^5`$) is much more extreme.

Suppose the heater plate has been on and has reached a steady temperature, so that heat input balances the conduction/convection losses at the top. Suppose also that the ratio of time scales is sufficiently large. Then put a piece of thermally insulating material on top, covering part of the surface of the plate (as a model for the effect of a spot). Where does the blocked heat flux go? The answer is of course that it does not go anywhere, the blocked heat stays in the block, heating it up slowly, instead of being radiated at some other part of the surface. On account of the large thermal conductivity, the temperature in the block stays uniform, and almost no ‘bright ring’ appears around the spot. When the insulating spot is removed after a while, the heat flux returns to its old value instantaneously. The temperature has increased only slightly during the blocking phase, and this slight increase disappears again on a thermal time scale.

### 3.5 excess emission by small scale magnetic fields

Apart from the highly visible sunspots, the Sun has a small scale magnetic field consisting of a very large number of tiny magnetic elements. They are called faculae if their density is high, and ‘network’ at lower density, and collectively contain much more magnetic flux than sunspots. They are visible especially near the solar limb, showing that they radiate more isotropically than the limb-darkened normal solar atmosphere.

This excess emission has two components, one originating in somewhat higher atmospheric layers, in the form of line- and UV-continuum emission. This component probably is due to dissipation of mechanical (wave) energy carried by the magnetic field. In addition, there is a component seen in white light, originating in the surface layers of the magnetic structures. This second component is understood theoretically as a radiative transfer effect (Spruit 1977): the radiating surface in a magnetic element is a little ‘dimple’ in the surface. The large number of small magnetic elements effectively corrugates the surface, and makes it radiate both more effectively and more isotropically, just like a rough surface has a higher emissivity than a polished surface. This effect increases the radiative flux from any part of the surface where small scale magnetic fields are present. It is a net effect (predicted before observation, Spruit 1977) that is not balanced by a deficit somewhere else on the surface. Since the magnetic elements are so small, they are hardly resolved in the observations. Recent high-resolution observations (Topka et al. 1997), however, verify the picture predicted by flux tube models (Spruit 1976 figure 11, Steiner et al. 1997) for the appearance of small magnetic elements.

Since the effect of small scale fields is an increased emissivity, it increases the cooling rate at the surface.

### 3.6 changes in granulation in active regions

The changes of surface emissivity by magnetic fields discussed above are due to the different radiative properties of small magnetic elements (‘flux tubes’). More indirect effects may play a role as well. The shape of granules appears to be different in magnetic regions (Macris and Roesch 1983, Muller and Roudier 1984, Muller 1986, Title et al. 1992). They are smaller, more irregular, and the measured flow speeds are lower. These effects are present in addition to the magnetic elements themselves, which also change the appearance of granulation by filling in the intergranular lanes, especially when seen at low spatial resolution. The changes may be the result of the geometric constraints the magnetic tubes put on the convective flow outside them.

Since the flow appears to be different, it would seem possible that the heat flux it carries is also different. This might contribute to the observed solar cycle variations of irradiance (Muller 1986, Kuhn et al. 1988). Limits on this contribution can be put by observations of the colors of the solar cycle variation signal. Solanki and Unruh (1998) find that the wavelength dependence of the irradiance signal argues against a dominant contribution from small changes of surface temperature over a larger areas, but do not quote quantitative limits. In order to contribute in the right sense to the solar cycle variation, the constraints imposed on the flow by the magnetic flux elements would have to lead to an increase in the heat flux carried by granulation. Direct (spatially resolved) measurements of a heat flux change in granulation in magnetic regions are probably difficult, since the effect would be small and hard to separate from the enhanced emission from the magnetic elements.

Changes in surface heat flux by modified granulation, if they exist, will affect the convective envelope in the same way as the excess emission from the small scale magnetic field. They can also be represented by an effective change of the emissivity of the solar surface. In particular, one does not expect these changes to be ‘compensated’ by opposite changes elsewhere on the surface (except, as before, on the $`10^5`$yr thermal time scale of the envelope).

## 4 Improving on the turbulent diffusion picture

In the results quoted a diffusion model for convection was used. It assumes that convection can be modeled by a turbulent viscosity for momentum and a turbulent diffusivity for the transport of heat. While this was a simple and somewhat justifiable model of convection in the absence of detailed knowledge of the convective flows in a stellar envelope, we now know that it does not represent stellar convective envelopes well. Numerical simulations (Nordlund 1982, 1985ab, 1986, 1990, Nordlund and Dravins 1990, Nordlund and Stein 1990, 1991, 1996, Stein and Nordlund 1989, 1991, 1998, Steffen et al. 1989, Steffen 1993, Ludwig et al. 1999) show an extremely nonlocal picture. The convective flows are driven almost exclusively by cooling at the surface, with narrow fast moving downdrafts between slow almost isentropic upflows (for a discussion see Spruit, 1997). In the present context, the most important property is that the convective flow at all depths is driven by cooling at the surface rather than by a local overturning process. The material in the cool downdrafts survives to large depths below the surface with little mixing into the upflows.

Given this extremely nonlocal picture, it is appropriate to ask how the thermal behavior of perturbations in the convection zone can be estimated without appealing to a turbulent diffusion model. I do this below for three of the effects discussed, the spot blocking, the facular emission, and the ‘thermal shadow’ problems. In all three cases, the conclusion is that thermal perturbations of the normal solar surface are even smaller than in the diffusion picture. This has to with the fact that the flows in the new non-local picture are caused by cooling at the surface, and that most of the changes are one-way (downward, away from the visible surface). The diffusion picture on the other hand communicates perturbations equally well in all directions.

### 4.1 spot blocking

Below the spot (modeled as a region of reduced heat transport efficiency extending to some depth below the surface) the upflows have exactly the same temperature as upflows in the unspotted surroundings at the same level, namely that given by the entropy at the base of the upflows. In this sense, there is no ‘pile up of heat below the spot’. Because of the reduced heat loss at the surface, however, the downflows below the spot will be less vigorous. The unspotted surface notices nothing of the spot’s presence (except for an extremely narrow ring where lateral radiative exchange takes place, and except for the presence of a moat flow, see below). It continues to cool upwellings into downdrafts as before, since the entropy in the upflows has not changed. Thus, we expect again that bright rings will be absent around spots, but the reason is even simpler than in the diffusion model. The spot is a region at the surface where less heat is radiated away, and this is now independent of the depth of the spot below the surface (in contrast to the diffusion model, where the spot has to extend to a minimum depth of 1000 km for the explanation to work).

Still, one may wonder what happens to the amount of heat generated in the solar interior that now fails to be emitted at the surface. This part of the problem is the same as in the diffusion model: the imbalance causes a secular increase of the entropy in the entire convection zone until a new thermal equilibrium is reached. Because of the very long thermal time scale of the convection zone ($`10^5`$ yr), the effect is negligible on observable time scales. In a steady state, when the average number of spots does not change, the convection zone does not heat up, because its mean temperature is higher than it would be without spots (Spruit and Weiss, 1986). Episodes of larger than average spot coverage cause heating, those of less than average spot coverage cause cooling on this time scale.

### 4.2 facular emission

The lower internal pressure in a small magnetic element (‘tube’) in the photosphere causes it to be more transparent: a surface of unit optical depth as seen from earth has a little dimple at the loaction of the small magnetic element. Its side walls radiate an additional heat flux (Spruit 1976, 1977). This heat flux is supplied to the element by convection, i.e. by the granulation surrounding it, and conversely the side wall, by the extra radiation it emits, cools the surroundings. In a simple minded diffusion model for convection, this causes a dark ring around the element, which compensates for a part of the extra emission from the element. The cooling, however, also makes the photosphere next to the tube heavier, so that a circulation is set up: down along the tube walls, and towards the element at the surface. This effect has been studied with numerical simulations by Knölker et al. (1991, see also Steiner et al. 1998), who find a very fast downflow along the tube wall (km/s) due to the excess cooling. In a diffusion model, the information of the low temperature in the downdrafts diffuses back to the surface. In reality, the cool downdraft drops below the visible surface because it is heavy, and only a negligible fraction of it finds its way back to the surface. The excess emission is therefore expected to be stronger than in the diffusion model, and not compensated by a dark ring.

### 4.3 thermal shadows

The thermal shadow expected from a blocking object below the surface (Spruit, 1977), which is a small effect already in the turbulent diffusion model, is likely to be even smaller in the new picture of stellar envelope convection. Since all upward flowing gas has nearly the same entropy, a blocking object below the surface has no effect on what happens above it, as long as there is enough room above it for upflows to comfortably reach the solar surface. It may interfere somewaht with downflowing gas, but this has no consequence for the surface temperature above the blocking object, since only a small fraction of the downward moving gas is carried back up to the surface. As before, this is a consequence of the non-local nature of the convective flow.

## 5 Helioseismological tests

In the above I have concentrated on surface effects, and argued that temperature anhancements with origins deep inside the convection zone are less likely to contribute. Such deep origins have been advocated by Kuhn et al. (1988), and Kuhn and Stein (1996) who propose that a temperature enhancement with origin at the base of the convection zone causes a temperature enhancement of 0.5K in the active latitudes.

Potentially exhausing theoretical debates for or against these models can probably be circumvented since the models make different predictions for the associated perturbations below the surface that can be tested by observation. The excess radiation from small scale magnetic fields comes about through an increased surface emissivity, associated with the increased ‘surface roughness’. This increases the average cooling rate at the surface and creates somewhat stronger and/or cooler downdrafts. It therefore causes slightly lower horizontally averaged temperatures below the surface. The two models thus predict opposite thermal effects below the surface.

The prospects for a test of these predictions are good, since helioseismological data are now quite sensitive, with the data obtained by SOHO/MDI and the GONG network (Kosovichev 1996, Kosovichev and Schou 1997, Duvall et al. 1998, Duvall and Kosovichev 1999). Changes in the p-mode frequencies related to the solar cycle were reported by e.g. by Libbrecht and Woodard (1990), Dziembowski et al. (1998), and Dziembowski (this volume). The dependence of the signal on wavenumbers $`(l,m)`$ shows that most of the effect is due to a change in wave propagation in a layer very close to the solar surface, at the latitudes of magnetic activity. The sense of the effect, an increase of the frequencies, is the opposite (Goldreich et al. 1991) of what would be expected if the temperature in these latitudes were higher. Though a higher temperature increases the propagation speed (increasing the frequencies), it also expands the solar envelope at these latitudes, by vertical hydrostatic balance. This lengthens the path traveled and decreases the mode frequencies. The net effect turns out to be dominated by the the second effect, the path length increase, because it is linear in $`T`$, while the propagation speed increases only as $`T^{1/2}`$.

In addition to the thermal effect, the presence of the magnetic field also has an effect, since it increases the ‘stiffness’ of the gas (in terms of propagation, the relevant speed is the fast mode speed, which is larger than the sound speed). Since the contribution of this effect scales with the ratio of magnetic to gas pressure, it is important only in a thin layer near the surface. The observed changes in mode frequency can thus be due either to a somewhat lower temperature, or the effect of the magnetic field in the atmosphere. The observation that most of the effect is localized at the surface argues for a major contribution from the magnetic stiffness effect (Goldreich et al. 1991, Woodard et al. 1991). In any case, the observed sign of the effect is incompatible with enhanced temperatures being the main change in envelope structure in the active latitudes.

A second way of measuring variations in propagation conditions below the solar surface is through time-distance seismology (Duvall et al. 1996, Kosovichev 1996, Braun et al. 1997), in which one does not measure the mode frequencies but the wave travel times between points on the surface (much like seismology is used to measure irregularities in the earth’s interior). These measurements show shorter travel times in the active regions (Duvall et al. 1998, Duvall and Kosovichev 1998, Giles et al. 1997, Bogdan et al. 1998, Goode and Dziembowski 1999), and downflows below the surface. The shorter travel times are compatible with a contribution of the magnetic field to the propagation speed.

The lower temperatures expected by enhanced surface cooling in the small scale magnetic field has further observational consequences. In the absence of rotation, the increased downward buoyancy would just cause a downdraft below active regions. Associated with it a horizontal flow converging on the active region is expected. Effects like these are beginning to be explored by time-distance seismology.