text:
# History effect in inhomogeneous superconductors

## I Introduction

Recently the so-called history effect (HE) of vortex matter has been frequently observed by a variety of methods such as electric transport measurement, dc magnetization hysteresis, ac susceptibility and torque techniques, and drawn much more attention . HE refers to that critical current density jc of a superconductor has different values at a same field (or temperature) for either a thermal, magnetic or current cycles. The earlier observed hysteresis of R-T curves in the thermal cycle in fact is a kind of HE caused by first order phase transition of vortex matter (melting or freezing) . Different mechanisms have been proposed to account for HE . And a typical explanation contains two assumptions: relatively weak flux pinning and uniform small quenched disorders. Then HE is explained with a theory of disorder-order transition of vortex matter. The disordered vortex matter pinned more strongly (higher jc) is supercooled to low temperature where it is in metastable state and then a process such as a bias current anneals the metastable disordered state into an ordered and stable one with lower jc. An important point is that flux creep in the models is omitted for the low temperatures. In this case, the property of vortex matter is governed only by the competition between interactions of vortices-vortices and vortices-quenched disorders.

However, a different kind of HE has been observed . The hysteresis loop of V-I curve of polycrystalline Bi2-xPbxSr2Ca2Cu3Oy (Bi2223) at liquid nitrogen temperature was measured there. See also Fig.6 in this paper. It is clear that such a V-I loop is a new kind of HE and very different from the V-I loop in low temperature superconductors reported before . The difference between these two kinds of HE is that the direction of the V-I loop is anti-clockwise in reference 1 while clockwise in reference 12 and this paper. Despite the anti-clockwise loop has been explained very well, the clockwise one has not yet been understood up to our knowledge and interested us very much.

Note that this new kind of HE takes place in a system different from those with uniform and weak quenched disorders at low temperatures. First, at higher temperatures, such as liquid nitrogen temperatures, thermal fluctuation should be important and flux creep could influence the property of vortex matter apparently . In fact the influence of flux creep on the V-I curve and transport property of Bi2223 have been confirmed by earlier experimental and numerical results . Second, the quenched disorders or flux pinnings are inhomogeneous in polycrystalline high temperature superconductors (HTS). There are weak links and grains in sintered samples such as polycrystalline Bi2223. And the flux pinning strength of the weak links at grain boundaries is very different from that of the grains . In view of these two points, any models concerning the system with homogenous quenched disorders without flux creep could not be used to understand this novel hysteresis behavior in V-I curve.

A theoretical model describing irreversible electric-magnetic behaviors of a system with inhomogeneous quenched disorders has already been proposed . Unfortunately, flux creep was omitted there. Moreover, this model has not been used to explain any HE, including the clockwise V-I loop. Experimentally, A.M. Campbell’s group did many important works on inhomogeneous superconductors. They did observe the hysteresis V-I curves before. However, only a qualitative explanation was given. And what interests us now is how to simulated the clockwise V-I loop.

Considering flux creep and inhomogeneous flux pinning, we proposed a new model to explain this new kind of HE. We carried out numerical simulation with this model. Upon explaining the clockwise loop of V-I curve reported in reference 12, we further predicted some new phenomena. To test the predictions, electric transport measurement was conducted for characteristic V-I curve of silver-sheathed Bi2-xPbxSr2Ca2Cu3Oy tapes (Ag-Bi2223) immersed in liquid nitrogen with or without applied fields Ba. Clockwise loop of V-I curve at certain dI/dt and Ba was observed and is dependent on dI/dt and Ba, confirming the numerical result.

## II Model

The model consists of two basic assumptions: there exists flux creep in a sample and the flux pinning is inhomogeneous. The complex inhomogeneity of flux pinning is simplified by periodic strong pinning region (S) with large jc (jcS) and weak pinning region (W) with small jc (jcW), see Fig.1. The flux diffusion is characterized by collective creep of vortex glass, which has been discussed widely . To simplify the calculation, we use logarithmic U(j) relationship, which is the special case of vortex glass model with very small glass exponent $`\mu `$.

We should emphasize that our model is not a simplified version of the ”brickwall” model, which has been widely discussed before . There are some resemblances between them. However, there are at least two points that distinguish our model from the previous one. First, the ”brickwall” model is a static critical state model without flux creep. Obviously, in any static models, the transport properties studied here will be independent of the sweeping rate of applied current dI/dt, contrasting with the experimental observation. Second, in the ”brickwall” model the S and W regions are connected in series along the direction of current. While in our model, they are connected in parallel.

In our one-dimensional model, flux lines firstly enter the system from the surfaces into the weak pinning channel (W) and then diffuse into the strong pinning region (S). So W surrounds S. In other words, S is inside the slablike sample.

As an example, the model may be proper to polycrystalline samples of HTS, such as silver-sheathed Bi2-xPbxSr2Ca2Cu3Oy (Ag-Bi2223). The grain boundaries in polycrystalline HTS form weak link network of flux pinning and thus can be considered as W, whereas the intra-grain region is S. It is well known that the weak link network suppresses substantially the critical current density of HTS such as Ag-Bi2223. When current or magnetic field is applied, flux density will enter inter-grain (W) at first and then penetrate gradually into the intra-grain (S) by means of flux creep. Another example the model may be applicable to is the low temperature superconductors, such as Nb3Sn metallic compounds, where flux lines may be pinned by grain boundary (S) whereas the intra-grain is weaker pinning region (W) .

Certainly, one cannot expect that numerical result by such a simplified phenomenological model can quantitatively represent experimental data of a sample. For example, the one-dimensional geometry assumption in the model is not fulfilled for a real sample. And, the flux pinning is more complex by far than the two model parameters, jcS and jcW. Nevertheless, our model is a new one accounting for the experimental clockwise V-I loop, which was observed very recently. The effectiveness of this model was demonstrated by the electric transport measurement of hysteresis V-I loop.

## III Simulation

### A Basic Equation

Consider a slab with infinite length along y-axis, thickness d along the x-axis and width w along the z-axis. The current is applied along the y-axis. In view of w$`>`$$`>`$d, we focus on the one-dimensional case for simplicity. The non-linear diffusion equation describing the macroscopic B or j is well known . With a logarithmic dependence of U on j: $`U(j)=U_0\mathrm{ln}\left|j_c/j\right|`$, the diffusion equation of flux can be written as:

$$\frac{B}{t}=\frac{v_0}{(\mu _0j_c)^{n+1}}\frac{}{x}\left[\left|\frac{B}{x}\right|^n\left(\frac{B}{x}\right)B\right]$$
(1)
where $`n=U_0/kT`$ and $`v_0=u\omega _m`$. $`v_0`$ is the attempt velocity of the thermal-activated vortex motion, $`u`$ is the hopping distance and $`\omega _m`$ is the microscopic attempt frequency.

### B Boundary and initial conditions

In the V-I curve measurement, the current is always applied with a certain sweeping rate dI/dt. The boundary condition of Eq.1 can be obtained from the current conservation equation: $`\frac{}{t}_0^w_0^dj𝑑x𝑑z=\frac{dI}{dt}`$. Substituting $`B/x=\mu _0j`$ into it, one comes $`B(d,t)B(0,t)=\mu _0I/w.`$ Due to antisymmetry, one has $`B(d,t)=B(0,t)`$, thus the boundary condition is

$$B(0,t)=\frac{\mu _0I}{2w}=\left(\frac{\mu _0}{2w}\frac{dI}{dt}\right)t$$
(2)
As for the initial condition, there is no current at t=0 in the sample, hence

$$B(x,0)=0$$
(3)
### C Numerical method and the choice of parameters

Such kind of non-linear diffusion equation can be numerically solved with the finite difference method. First, the temporal and spatial variances are discretized with different steps respectively. Since the temporal and spatial steps are very important for the stability of the calculation, we take pains to choose suitable steps to avoid a divergence problem. Second, the equation is discretized and rewritten in an implicit difference scheme, which is better than the explicit difference scheme in the stability of calculation. Finally, based on the initial condition, the calculation by iteration will obtain the temporal and spatial magnitude of physical quantity, namely B(x,t), j(x,t) and E(x,t). The voltage can be obtained by integration of E(x,t). Note that the boundary condition plays a role of constraint in every time step during current ascending and descending, see Eq.2.

It is easy to see that $`U(j)=U_0\mathrm{ln}\left|j_c/j\right|`$ combined with $`E=E_0\mathrm{exp}(U/kT)`$ leads to the power dependence of the E-j curve $`E=E_0(j/j_c)^n`$ with $`n=U_0/kT`$. This power law E-j curves are observed frequently in transport measurements and can be used to determine the parameter $`n`$ . Certainly, $`n`$ is dependent on flux pinning strength, magnetic field and temperature . The weaker flux pinning and higher temperature bring about smaller n. In fact, many experiments including the present work show that for the Ag-Bi2223 in liquid nitrogen temperature, $`n=6`$ is a typical value.

As for the velocity of the thermal-activated vortex motion and the critical current density, typical and reasonable magnitudes were employed based on the earlier theoretical and experimental work, for example we took $`v_0=u\omega _m=1m/s`$ ( if $`u`$ $`\mathrm{~}10^6m`$, $`\omega _m\mathrm{~}10^6s^1`$) and $`j_{cW}=2\times 10^8A/m^2`$. It should be pointed out that our numerical result is not sensitive to the choice of $`n`$, $`v_0`$ and jcW.

### D Numerical results and discussions

#### 1 Effect of jcS/jcW ratio on the hysteresis loop

Since jc represents the pinning strength, the ratio jcS/jcW can be considered as a parameter reflecting the ratio of the two pinning strengths. Shown in Fig.2 is our numerical result of the dependence of HE on jcS/jcW at fixed jcW, dS/dW and dI/dt. It is very clear that the V-I loop is clockwise as reported in reference 12. To understand how jcS/jcW affects the V-I loop, the corresponding current distributions in sample were also calculated.

Note that the velocity of flux diffusion, namely the velocity of current diffusion in transport measurements, can be written as: $`v=v_0\mathrm{exp}(U/kT)`$ $`(j/j_c)^n`$ for logarithmic U(j). Hence, for the same $`j`$, the larger the jc, the smaller the $`v`$. Generally speaking, the average speed of current diffusion in S region ( $`\overline{v_s}`$ ) is smaller than that of W region ( $`\overline{v_w}`$ ). In our calculation, jcW is fixed, i.e. $`\overline{v_w}`$ is fixed. By changing jcS/jcW ratio, namely jcS, we adjust $`\overline{v_s}`$ only.

For a large ratio (jcS/jcW=10), i.e. high inhomogeneity of flux pinning, $`\overline{v_s}`$ is so small comparing with $`\overline{v_w}`$. And in such a case, there is no enough time for current diffusion in response to both the ascending branch of V-I curve (Fig.2.1(b)) and the descending one (Fig.2.1(c)). Consequently, no obvious V-I loop can be detected (Fig.2.1(a)).

For a medium ratio (jcS/jcW=4), i.e. a medium inhomogeneity of flux pinning and a larger $`\overline{v_s}`$, it is shown that the area of the loop is very large (Fig.2.2(a)). In this case, although current has already penetrated into S region during the ascending branch (Fig.2.2 (b)), much more current penetrates into S during the descending branch (Fig.2.2(c)). That is to say, the current diffusion of S region has more time in response to the descending branch than to the ascending one. As a result, the V-I loop is obvious.

For a small jcS/jcW ratio, low inhomogeneity of flux pinning, $`\overline{v_s}`$ is very large and approaches $`\overline{v_w}`$. And there is little difference between the distributions of current during the two branches (Fig.2.3(b), (c)). The minor difference is only seen when the current is small and the corresponding voltage is too small to be detected. Hence, there is almost no hysteresis loop (Fig.2.3 (a)). Especially, at the uniform case (jcS/jcW =1, $`\overline{v_s}=\overline{v_w}`$), there will be no loop at all.

#### 2 Effect of dS/dW ratio on the hysteresis loop

With a fixed total thickness, the spatial distribution of S and W region, i.e. the thickness ratio of them is also important for the V-I loop. The numerical V-I loops with different dS/dW at fixed dI/dt and jcS/jcW were displayed in Fig.3. It is found that the smaller the dS/dW, the smaller the area of the loop.

For dS/dW =0.857 and other parameters given here, the thickness of S region i.e. the distance for current diffusion in S region is so long that current diffusion in S region has less time to respond the ascending branch (Fig.3.1(b)) than to respond the descending one (Fig.3.1(c)). So the area of the V-I loop is large (Fig.3.1(a)).

For dS/dW =0.222, S region is thinner, i.e. the distance of current diffusion in S region is shorter, resulting that comparing with the dS/dW =0.857 case, the current diffusion in S region has more time in response to the ascending current (Fig.3.2(b)). In other words, for such parameters given here, the response time of current diffusion in S region to the two branches of the V-I loop shows little difference. As a result, the area of the V-I loop becomes smaller (Fig.3.2(a)).

As for dS/dW =0.105, S region is too thin and the distance of current diffusion in S region is too short. Hence there is enough time for current diffusion in response to both current ascending and descending at the given parameters (Fig.3.2(b), (c)). That is to say, there is no obvious V-I loop (Fig.3.3(a)). As dS/dW approaches 0, namely the uniform flux pinning case, the hysteresis loop will disappear thoroughly.

#### 3 Current sweeping rate (dI/dt) dependent V-I loop

The numerical result of current sweeping rate dependent V-I loops at fixed jcS/jcW and dS/dW was illustrated in Fig.4.

It is seen that if dI/dt is very small, such as 0.01A/s, flux lines will have enough time to respond the changing applied current and thus can penetrate into both W and S regions. The difference of current distribution during the ascending (Fig.4.1(b)) and descending branch (Fig.4.1(c)) is very small. As a result, the area of the V-I loop is also small (Fig.4.1(a)). As dI/dt reaches a moderate value, say dI/dt=2A/s, the difference between the current distributions during the ascending and descending branch is larger. And the area of the V-I loop becomes larger too(Fig.4.2(a)). Nevertheless, when dI/dt is too large, say dI/dt=10A/s, the area of the V-I loop decreases dramatically (Fig.4.3(a)). It is easy to understand that in such a case, comparing with dI/dt, $`\overline{v_s}`$ is so small that the current cannot penetrate into the S region during both the ascending and descending branch. And the corresponding current distributions during the two branches do not show visible difference (Fig.4.3(b), (c)). So there is no obvious V-I loop.

Now, it is clear what the effect of dI/dt on the V-I loop means. With fixed distance and average speed of the current diffusion, to change dI/dt is equivalent to shift our observing time window. When we choose time window by using dI/dt as 10A/s, the observing time is so early that there is no enough current inside S region by means of diffusion even the applied current sweeping has been finished. Of course there is no detectable V-I loop. If we use dI/dt=0.01A/s, we shift our time window so late that there are always enough current diffused into S region during both the ascending and descending branch. It is natural that no obvious V-I loop exists at all. Only when we choose a proper observing window, i.e. a proper dI/dt such as 2A/s, is there much more current in S region during the descending branch than the ascending one. In other words, only in such case, can the current diffusion in S region has more time to respond the descending branch than the ascending one. And the V-I loop will be obvious.

In conclusion, if there is more time for current diffusion of S region to respond the descending branch than to respond the ascending one, hysteresis loop in V-I curve, i.e. the HE will be obvious. On the other hand, if the current diffusion time of S region is similar or the same for the ascending and descending branch, there will be no obvious V-I loop, namely no HE. The observation of HE can be controlled by adjusting dI/dt ( the observing time window ) or dS/dW ( the distance of current diffusion ) or jcS/jcW ( the average speed of current diffusion ).

#### 4 Magnetic field dependent V-I loop

It is well known that jc is dependent on the magnetic field Ba. Commonly, jc decreases with increasing Ba. Therefore, to study the effect of jc on the V-I loop is equivalent to study the effect of Ba on it. For the inhomogeneous pinning case, jcW is more sensitive to Ba than jcS, especially when Ba is not very high . Thus, we simulated the V-I loops in different applied fields Ba, namely different jcW at fixed jcS, dI/dt and dS/dW. Shown in Fig.5 are the numerical results. It is clearly seen that the loop shifts towards smaller current with decreasing jcW while the area of the loop is insensitive to jcW. In other words, the area of the loop is weakly affected by magnetic field.

As mentioned above, the existence of HE depends on that the current diffusion of S region has different responding time for the two sweeping branches of V-I loop. Only changing jcW, i.e. Ba, would not affect the responding time of S region. As a result, the area of V-I loop will not change acutely.

## IV Experimental

Because our numerical simulation has obtained more results than the experimental one reported in reference 12, such as the effect of dI/dt on the area of loops, we measured V-I curves of Ag-Bi2223 samples in various applied magnetic fields Ba with different dI/dt to confirm the numerical prediction. It has been pointed out above that the Ag-Bi2223 sample may be proper to test our model.

The Ag-Bi2223 samples cut from long silver sheathed tapes come from the National Center for R&D on Superconductivity of China. All the samples are the c-axis textured. The c-axis is perpendicular to the wider surfaces of the samples. Including the outer silver sheath, the average size of the samples is approximately: $`4mm\times 1mm\times 4cm`$. Between each two neighbor points of four (two voltage contacts and two current contacts) is 1cm in length. The transport measurement was carried out for the sample immersed in liquid nitrogen. The external magnetic field with range 0 ~ 200 Gs was applied in the superconducting state. The current range of the supply source is 0 ~ 40 A and sweeping rate range 0.02 ~ 1 A/s. The measurements were made in the case of Ba parallel to the c-axis of the samples. More details of the sample preparation and the measurements can be found in reference 18.

The V-I loops were measured by ascending current first and then descending it. Displayed in Fig.6 are the experimental V-I loops with different dI/dt and fixed Ba. The spottiness of the experimental data came from the limited picking rate of our measurement system. As reported in reference 12, clockwise loop of V-I curve was observed. It is clearly seen that dI/dt affected not only the appearance of V-I curve but also the area of the loop. For example, the area of loop with dI/dt = 1 A/s is larger than the area of loop with dI/dt = 0.5 A/s, which is in agreement well with the numerical result.

Shown in Fig.7 are the experimental V-I loops in different applied fields Ba at a fixed dI/dt. It is found that the area of the loop is insensitive to magnetic field Ba whereas the transport critical current density jc decreases with increasing applied field very sensitively. This experimental dependence confirms the above numerical results and thus supports strongly our model. Furthermore, the agreement between the numerical and experimental V-I shifting by applied field implies that the critical current density jc of a sample with spatially non-uniform flux pinning strength, such as the silver sheathed Bi2223 tapes, is governed mainly by the weaker flux pinning region (the region with smaller jc). Hence, for these HTS materials, which have potential for technical applications, the point to enhance their critical current density is to improve the flux pinning strength of weak pinning region.

We are unable to measure more V-I curves to test all the numerical curves for the limited ability of our instrument. For example, dI/dt cannot be too small because Ohm heat at the current contacts will destroy the sample.

## V Summary

We have proposed a model to account for a new kind of history effect in an inhomogeneous flux pinning superconductor at high temperatures. The model consists of alternating weak and strong flux pinning regions whose strength was depicted by different critical current densities jcW and jcS, respectively. Based on the model, our numerical simulation successfully observed a hysteresis loop in the characteristic V-I curve as reported in references. Our numerical results also predicted some new characteristics of the V-I loop at different sweeping rates of applied current dI/dt, applied magnetic fields Ba, the ratio of the two regions’ pinning strength and their thickness.

Physically, due to the difference of flux pinning strength in S and W regions, the average speeds of flux (or current) diffusion in the two regions are different too. If dI/dt (or dS/dW, jcS/jcW) is proper, the V-I loop, i.e. the HE is obvious. Otherwise, no obvious V-I loop can be observed.

To confirm the numerical results, electric transport measurement was conducted to get the characteristic V-I curve of Ag-sheathed Bi2-xPbxSr2Ca2Cu3Oy tapes immersed in liquid nitrogen with or without applied fields. Clockwise V-I loop was observed to be dependent on dI/dt and Ba. Hence, we presented a new possible mechanism accounting for the new kind of HE in an inhomogeneous flux pinning superconductor.

Acknowledgment

The work is supported by the Ministry of Science and Technology of China (NKBRSF-G1999-0646) and NNSFC (No.19994016).

Figure Captions:

Fig.1. A schematic sketch of the one-dimensional model: a superconducting slab consisting of periodic strong (S) and weak (W) pinning regions. Vortex lines enter the system firstly from the surfaces into W region, and then diffuse from the W into S regions. So the W region surrounds the S region.

Fig.2.Numerical results, showing the effect of jcS/jcW ratio on the V-I loop (a) and the distribution of current during the ascending (b) and descending branch (c) of the loop. The arrows indicate the increasing and decreasing of current, i.e. the time evolution. The dot lines represent the boundary between S and W regions. dS/dW=1/2, dI/dt=1A/s, jcW=2$`\times `$ 108A/m2.

(2.1) jcS/jcW =10, high inhomogeneity of flux pinning. Average speed of current diffusion in S region is too small comparing with that of W region, resulting that in S region there is no enough time for current diffusion in response to both the ascending sweeping and the descending one at the given parameters. No obvious V-I loop.

(2.2) jcS/jcW =4, medium inhomogeneity of flux pinning and lager average speed of current diffusion in S region, resulting that more current has penetrated into the S region during the descending branch than the ascending one at the given parameters. Obvious V-I loop is observed.

(2.3) jcS/jcW =1, uniform flux pinning, equal average speed of current diffusion in the two regions. There is little difference between the distributions of current during the two branches. No obvious V-I loop.

Fig.3. Numerical results, showing the effect of dS/dW ratio on the V-I loop (a) and the distribution of current during the ascending (b) and descending branch (c) of the loop. jcS = 4 jcW, dI/dt=1A/s, jcW=2$`\times `$108A/m2.

(3.1) dS/dW =0.857, the thickness of S region is proper, resulting that in S region there is enough time for current diffusion in response to current descending, but the current diffusion is too busy to respond to current ascending for the given parameters, and area of the V-I loop is large.

(3.2) dS/dW =0.222, S region is thinner and the distance of current diffusion in S region is shorter, resulting that in S region there is much more time for current diffusion in response to current ascending at the given parameters, and the area of V-I loop becomes smaller.

(3.3) dS/dW =0.105, S region is too thin and the distance of current diffusion in S region is too short, resulting that in S region there is enough time for current diffusion in response to both the current ascending and descending at the given parameters, no obvious V-I loop.

Fig.4. Numerical results. Current sweeping rate dependent V-I loop (a). The corresponding distributions of current during the ascending (b) and descending branch (c) were also shown. dS/dW=1/2, jcS/jcW =4, jcW=2$`\times `$108A/m2.

(4.1) dI/dt=0.01A/s, the current sweeping rate is so slow that in S region there are enough time for current diffusion in response to both the ascending and descending current sweeping at the given parameters, no obvious V-I loop.

(4.2) dI/dt=2A/s, the current sweeping is faster, resulting that in S region there is no enough time for current diffusion in response to current ascending at the given parameters, and the V-I loop is larger.

(4.3) dI/dt=10A/s, the current sweeping is so fast that in S region there is no enough time for current diffusion in response to both the ascending and descending current sweeping at the given parameters. No obvious V-I loop.

Fig.5. Numerical results. Magnetic field dependent V-I loop, where the area of the loop is insensitive to jcW while the transport jc is strongly affected by magnetic field. dS/dW=1/2, jcS =8$`\times `$108A/m2, dI/dt=1A/s.

Fig.6. Experimental V-I loops of the Ag-Bi2223 tape with different current sweeping rates dI/dt, which is qualitatively with the numerical curves shown in Fig.4 and is a conformation of our model.

Fig.7. Experimental V-I loops of the Ag-Bi2223 tape at different applied magnetic fields, which is qualitatively with the numerical curves shown in Fig.5. It is noted that the area of the loop is insensitive to Ba whereas the value of the jc is strongly affected by magnetic field, which supports the present model further.


format: text

text:
# Abstract

## Abstract

The model of the Big Bang is an integral part of the national curriculum for England. Previous work (e.g. Baxter 1989) has shown that pupils often come into education with many and varied prior misconceptions emanating from both internal and external sources. Whilst virtually all of these misconceptions can be remedied, there will remain (by its very nature) the obstacle of ex-nihilo, as characterised by the question ‘how do you get something from nothing?’ There are two origins of this obstacle: conceptual (i.e. knowledge-based) and cultural (e.g. deeply held religious viewpoints). The article shows how the citizenship section of the national curriculum, coming ‘online’ in England from September 2002, presents a new opportunity for exploiting these.

Keywords: High School, Astronomy/Cosmology, Social/Cultural

## 1 Introduction

The model of the Big Bang (The Standard Model herein; e.g. Peebles 1993) at Key Stage 4 is a requisite part of science and physics teaching in England (i.e. contained within the framework of the National Curriculum; 3d: ‘about some ideas used to explain the origin and evolution of the Universe’). It is common for trainee teachers to be given background information (typically based upon lectures courses and planetarium visits) on the teaching of astrophysics and cosmology.

Even if a teacher is a consummate professional who has been trained in the teaching of astrophysics, however, one finds that there is a miasma of obstacles in the way of successfully teaching the Standard Model. Indeed, Comins (1993) prepares a list of over 550 individual misconceptions related to the teaching of astrophysics. These misconceptions fall roughly into three camps: those caused by incorrect internal (mental) processes, those caused by incorrect external processes (inaccurate information from teachers, parents and peers), and those caused by a combination of both internal and external factors. Whilst there exists plenty of evidence in the literature for the causes of such misconceptions, there is little in the way of strategies to deal with them.

## 2 Obstacles

Research by a number of different authors (e.g. Baxter 1989; Sneider & Pulos 1983; Nussbaum, 1979) suggests that children and juveniles will start to formulate their own ideas about the natural world long before any formal education is given. Inevitably, such self generated formulations will be riddled with both internal and external misconceptions. As their formal education progresses, these scientific misconceptions should be systematically dealt with during the course of time. For example, Baxter (1989) cites the case of a pupil persuaded towards an Aristotelian point of view.

Such misconceptions can readily be corrected. To this end, Bennett (2000) advocates the use of a five step plan based upon the un-learning of these misconceptions. By establishing a contextual framework, Bennett (2000) argues that misconceptions about astronomy can expediently be corrected. Consider the case of a pupil harbouring an Aristotelian viewpoint. Aristotle reasoned that the Earth was stationary in space and the Sun, Moon, other planets and stars orbit around it in circles. Later, Ptolomy refined this view with his theory of celestial spheres: the Earth was at the centre of the Universe and each extra-terrestrial body moved its own sphere with the outermost sphere consisting of all of the (fixed) stars. Knowing that this left room for Heaven and Hell outside of these spheres, the church adopted this theory. Observations consistently disagreed with this view, the primary one being retrograde (backwards) movement of planets across the sky. Although attempts were made to understand this movement (e.g. Ptolomy proposed that planets may execute epicycles) it was not until the seventeenth century that Galileo, having learnt of the invention of and subsequently built his own telescope, provided conclusive observations against the geocentric model of Aristotle. He showed that the satellites of Jupiter (later known as the Galilean satellites) do not orbit Earth but orbit Jupiter instead. Earth need not be at the centre of everything. Indeed, Galileo showed too that Venus exhibits phases, which can be explained if both Venus and the Earth orbit the Sun. Further, Kepler’s development from the heliocentric model of Copernicus demonstrated elegantly how retrograde motion could be explained in a simple manner. Further problematic observations for the geocentric model in the seventeenth century included Galileo demonstrating that sunspots were not shadows of planets passing in front of the Sun, rather actual spots on its surface that rotated; and the supernova of 1604 whose parallax demonstrated that it was located beyond the planets and thus the heavens were subject to change. Presentation of these observational arguments would readily challenge and correct an Aristotelian viewpoint.

There is, however, only so far that this technique can be applied before a terminal problem arises. This obstacle can be characterised by the question ‘how do you get something from nothing?’ (the ex-nihilo obstacle, if you will). This question can be divided into two aspects: a conceptual problem and a cultural problem.

Scientific understanding and knowledge of how to tackle the ex-nihilo obstacle is, perhaps, at a level above and beyond what an average physics teacher is capable of addressing. Indeed, it is probable that said teacher has had little theoretical grounding in cosmology and the Standard Model; depending, of course, where they were trained and what courses they undertook as an undergraduate (Baxter 1991; Baxter et al. 1991). Even if they do possess such knowledge, answering the characteristic question given above would require great expertise and teaching to effectively give an answer appropriate at Key Stage 4. The crucial problem is that the theory of general relativity (e.g. Einstein 1950) requires modification when density tends to infinity (as happens when we reverse time towards the Big Bang; Peebles 1993) to take account of quantum effects. Since we do not possess a complete theory of quantum gravity, we can presume that our theories are only valid for times after quantum effects became small. This limit of general relativity is the Planck time, about $`10^{43}s`$ after the Big Bang. We can only make extrapolations to the Planckian era based upon theory and conjecture, which are typically only accessible at the post-doctoral level.

Inevitably, by the very nature of cosmology, the answering of such a characteristic question is likely to come into conflict with pupils who possess strongly held religious beliefs (e.g. creationism). This conflict cannot be simply resolved by citing scientific evidence: these pupils’ views amount to belief, and are based on faith.

## 3 Discussion

In attempting to answer the ex-nihilo obstacle, there are several tactics that a teacher may employ. To confront the conceptual problem requires an acknowledgement that scientific understanding and knowledge is still being developed to this day. For example, since the Hubble Space Telescope (HST) was launched 1990 it has provided astronomers with a plethora of new data to characterize and explain. Indeed, one of the HST key projects were observations of distant supernovae to determine how the Universe is expanding. This led to the startling discovery that the Universe is not only expanding, but accelerating (e.g. Perlmutter et al. 1999). These observations are forcing astronomers to reconsider their cosmological models. We can only guess at what evidence the next generation space telescope may reveal. Science is, after all, characterised by evidence-based argument. Accordingly, there is no one theory that describes everything in the Universe without flaw.

Teachers are therefore advised to make appropriate citations to evidence for the Big Bang. Starting simply, the expansion of the Universe explains why light captured from distant galaxies displays similar line spectra to light in a laboratory but is shifted to longer wavelengths (e.g. Hubble & Humason 1931). At a more sophisticated level, the cosmological microwave background radiation (CMBR; Penzias & Wilson 1965) provides more evidence for the Big Bang.

If the Universe came from a previous epoch that was highly compressed and hot, the radiation from that epoch must still be with us today. When the temperature of the Universe was billions of Kelvin, such radiation would have been in the form of $`\gamma `$-rays. As the Universe expanded, however, it cooled and the radiation’s wavelength increased (recall Wein’s law). Gamow and others calculated that this radiation at the present epoch would be in the microwave part of the electromagnetic spectrum (e.g. if the temperature is about 5 K, then the wavelength of the radiation will be $`\lambda 0.06cm`$). Therefore, interpreted as highly redshifted radiation from a primeval fireball, the ubiquity of the CMBR gives excellent evidence of a previous state of the Universe that was extremely hot. Both the redshifted Universe and the CMBR can be explained at Key Stage 4 level. Whilst the Big Bang model explains these and many other observations, historically it has not been without criticism. For example, one may enquire, why is the Universe so close to the critical density? For the moment, scientists accept the Big Bang model as it ties together many threads of evidence. Yet, they must also be sceptical: the model has had problems in the past and will no doubt have problems in the future too.

The cultural problem is a much thornier issue. Contained within the declaration of human rights is the right to freedom of religion. The teacher, being in a responsible role, must ensure that they do not force their own (scientific) viewpoints upon such a pupil. Rather, they must allow the pupil to come to their own conclusions after presentation of the scientific evidence that supports the Standard Model. To force a viewpoint upon a pupil would be more than counter-productive. Such force could generate resentment towards the teacher, and lead to rejection of the subject material without consideration.

Since these dichotic opinions potentially lead to conflict, I advocate the use of a full classroom discussion to debate the issues surrounding the ex-nihilo obstacle. A discussion provides excellent opportunities for cross-curricular activity whilst allowing the pupils to express their own opinions. Research skills (to back up any scientific argument; for example, see Stephen Hawking’s Universe series at www.pbs.org/wnet/hawking), information technology skills (to present ideas with) and literacy can all play a significant role in such a discussion. Other subject areas could readily be drawn upon as well (e.g. drama).

More significantly, with the subject of citizenship coming ‘online’ in England from September 2002, an already squashed timetable could readily accommodate aspects of citizenship within an astronomical origins debate. Particular items in the citizenship national curriculum programme of study include (at Key Stage 4) 1b (‘the origins and implications of the diverse national, regional, religious and ethnic identities in the United Kingdom and the need for mutual respect and understanding’); 2 (all points under ‘developing skills of enquiry and communication’ relevant). Not only will such a discussion enhance pupil development, it would also actively promote citizenship though astronomy. It must be emphasized, however, that pupils need to be carefully guided to preserve mutual respect.

## 4 Conclusions

This work has reviewed some research on the misconceptions encountered during the teaching of astrophysics. When corrected, these misconceptions will generally be dealt with effectively, but there will often remain an ex-nihilo obstacle.

Such an obstacle can be divided into two separate issues: a conceptual (knowledge-based) one and a cultural one. In tackling these issues, one must address:
1. Science is still developing to this day and we do not presently possess a theory that describes everything without flaw.
2. Some pupils will hold religious values that may appear contradictory to the Standard Model.
Since both of these points lead to a potential conflict and debate, a full discussion is therefore advocated, simultaneously developing aspects of the national curriculum for citizenship.

## Acknowledgements

The preparation of this manuscript has made use of resources at both the University of Durham and St.Leonard’s R.C.V.A. School, Durham. KAP thanks Mary C. Hawkrigg for her reading through of the text and for providing useful suggestions that have improved this work. Finally, thanks to John Newman, Marion E. Jones, Peter Little and Simon B. Campbell for most unknowingly inspiring this work and for providing copious amounts of feedback on my teaching.

## References

Baxter, J., 1989, International Journal of Science Education, 11, 502
Baxter, J., 1991, Quart. Journal of the Royal Astron. Soc., 32, 147
Baxter, J., et al., 1991, Quart. Journal of the Royal Astron. Soc., 32, 159
Bennett, J., 2000, American Astron. Soc. Meeting, 197, 86.01
Comins, N.F., 1993, Bulletin of the American Astron. Soc., 25, 1430
Einstein, A., 1950, ‘The Principle of Relativity’, Metheun, London
Hubble, E. & Humason, M. L. 1931, Astrophysical Journal, 74, 43
Nussbaum, J., 1979, Science Education, 63, 83
Peebles, P.J.E., 1993, ‘Principals of Physical Cosmology’, Princeton Series in Physics, Princeton University Press, NJ
Penzias, A. A. & Wilson, R. W. 1965, Astrophysical Journal, 142, 419
Perlmutter, S. et al., 1999, Astrophysical Journal, 517, 565
Sneider, C. & Pulos, S., 1983, Science Education, 67, 205


text:
# Influence of the Potential Energy Landscape on the Equilibration and Specific Heat of Glass Forming Liquids

## Abstract

We show that a glass transition, signaled by a peak in the specific heat vs. temperature, can occur because a glassy system that shows no signs of aging progresses so slowly through the energy landscape that the time needed to obtain an accurate estimate of the thermodynamic averages exceeds the observation time. We find that below the glass transition temperature of a three dimensional binary mixture of soft spheres, the specific heat increases with measurement time spans orders of magnitude longer than previously recognized equilibration times.

As a glass forming liquid is cooled or compressed, a glass transition occurs when the system falls out of equilibrium, i.e., when the time scale for reaching equilibrium exceeds the observation time Ediger96 . Like many complex systems, such as proteins and neural networks, the dynamics of such a system is strongly influenced by the potential energy landscape where each point corresponds to a particular configuration and energy of the system Goldstein69 ; Sastry98 . The energy landscape can be used to describe the three ways in which a system can fall out of equilibrium. First a system can become trapped in a metastable minimum where it stays for the duration of the observation time. Second a system can be in an energetically unlikely part of phase space and proceed slowly to a region of the energy landscape where its configurations obey a Boltzmann distribution. As such a nonequilibrium system evolves toward more probable regions of phase space, it exhibits aging which means its properties systematically change with time and do not obey stationarity Kob97b . The aging time, after which aging stops, is equal to the $`\alpha `$ relaxation time which is the charactistic time for the system to forget its initial configuration.

The third way is not widely appreciated and is the subject of this paper. Namely, even after a glassy system no longer ages and has reached basins with appropriate energies, the system proceeds so slowly through the energy landscape that it takes a long time to accumulate the large number of statistically independent measurements needed to accurately determine a thermodynamic average. We find from molecular dynamics simulations that a glass forming liquid can undergo a glass transition, as signaled by a peak in the specific heat $`C_V`$ versus temperature, that is due to insufficient averaging by a system that shows no signs of aging. The distribution of energies that a system samples in a basin of the energy landscape is a subset of the full distribution of energies available to the system. Since this subset has a smaller variance than the full distribution, the resulting specific heat, which is proportional to the variance of the energy, will be smaller when calculated from short time spans than from long time spans. These smaller values account for the values below the peak in $`C_V`$ on the low temperature side. Going to longer time spans eliminates the peak, though at temperatures below the peak temperature $`T_p`$, these longer time spans can be orders of magnitude longer than previously recognized equilibration times such as the $`\alpha `$ relaxation time, the energy correlation time, and the aging time.

We have performed a molecular dynamics simulation on a three dimensional glass forming liquid Weber85 ; Kob94  consisting of a 50:50 mixture of two types of soft spheres, labelled A and B, which differ only in their sizes. The interaction between two particles a distance $`r`$ apart is given by $`V_{\alpha \beta }(r)=ϵ[(\sigma _{\alpha \beta }/r)^{12}+X_{\alpha \beta }(r)]`$ where the interaction length $`\sigma _{\alpha \beta }=(\sigma _\alpha +\sigma _\beta )/2`$, and the ratio of the ratio of the diameters $`\sigma _B/\sigma _A=1.4`$ ($`\alpha `$, $`\beta =`$ A, B). The cutoff function $`X_{\alpha \beta }(r)=r/\sigma _{\alpha \beta }\lambda `$ with $`\lambda =13/12^{12/13}`$. The interaction is cutoff at the minimum of the potential $`V_{\alpha \beta }(r)`$. Energy and length are measured in units of $`ϵ`$ and $`\sigma _A`$, respectively. Temperature is in units of $`ϵ/k_B`$, and time is in units of $`\sigma _A\sqrt{m/ϵ}`$ where $`m`$, the mass of the particles, is set to unity. The equations of motion were integrated using the leapfrog method Rapaport95  with a time step of 0.005. During each run the density $`\rho _o=N/L^3`$ was kept constant. $`N=N_A+N_B`$ is the total number of particles. The system occupies a cube with volume $`L^3`$ and periodic boundary conditions. According to the ideal mode coupling theory, the relaxation time diverges at a temperature $`T_C`$ Gotze92 . For our system $`T_C=0.303`$ Carruzzo01 .

We cool the system from a high temperature (T=1.5) by lowering the temperature in steps of $`\mathrm{\Delta }T=0.05`$. At each temperature we equilibrate for $`10^4`$ time steps and then measure the quantities of interest for $`10^7`$ additional steps. The results are then averaged over different runs. At the glass transition the system falls out of equilibrium and becomes trapped in a basin of the energy landscape. In order to avoid this, we have used parallel tempering together with molecular dynamics. We implement parallel tempering (PT) Kob01 ; Carruzzo01  by running molecular dynamics simulations in parallel at chosen temperatures using a temperature constraint algorithm Rapaport95  to keep the temperature of each simulation constant. At 100 time step intervals we attempt to switch the configurations of two neighboring temperatures using a Metropolis test which ensures that the energies of the configurations sampled at any given temperature have a Boltzmann distribution. Let $`\beta _1`$ and $`\beta _2`$ be two neigboring inverse temperatures, and let $`U_1`$ and $`U_2`$ be the corresponding potential energies of the configurations at these temperatures. If $`\mathrm{\Delta }=(\beta _1\beta _2)(U_2U_1)`$, then the switch is accepted with probability unity if $`\mathrm{\Delta }0`$ and with probability $`\mathrm{exp}(\mathrm{\Delta })`$ if $`\mathrm{\Delta }>0`$. The acceptance ratio is between 30% and 75%. Near $`T_p`$, the acceptance ratio was above 60%. After a swap is accepted, the velocities of the particles in each configuration are rescaled to suit their new temperature. Switching configurations allows a given simulation to do a random walk in temperature space in which it visits both low temperatures and high temperatures. This helps to prevent it from becoming trapped in a valley of the energy landscape at low temperatures. Typically we equilibrate for 2 million time steps and then do measurements for an additional 4 to 10 million molecular dynamics (md) steps. We average over different runs which have different initial positions and velocities of the particles at each temperature.

We calculate the specific heat $`C_V`$ per particle at constant volume in two ways. The first uses the fluctuations in the potential energy $`U`$ per particle: $`C_V=(3k_B/2)+Nk_B\beta ^2\left(U^2U^2\right)`$ where the first term is the kinetic energy. The second way uses $`C_V=d<U>/dT(<U(T_{i+1})><U(T_i)>)/(T_{i+1}T_i)`$. The results are shown in Figure 1. Notice that there is a sharp asymmetric peak centered at $`T_p=0.305\pm 0.003`$. The low temperature side of the peak drops steeply. The curves for 216 and 512 particles coincide, indicating that there is no size dependence. The discrete points are calculated from the energy fluctuations. The solid line is calculated from the derivative of the energy. The fact that the two coincide indicates that the system was equilibrated in all the basins of the energy landscape that were visited. For comparison we also show the result of cooling through the transition ($``$). The peak in $`C_V`$ found by cooling coincides with the peak found in parallel tempering.

We do not believe that the parallel tempering peak in the specific heat is due to the system becoming trapped in a metastable minimum in the energy landscape for the following reason. The configurations corresponding to the bottom of the basins of the energy landscape are called inherent structures Stillinger82 . We sampled the configurations that were visited during the parallel tempering runs and found the corresponding inherent structure energy $`e_{IS}`$ per particle by minimizing the potential energy locally using the method of conjugate gradients Press92 . In the inset of Figure 1 we plot the average inherent structure energies versus the temperature of the configuration that was originally saved. As the temperature decreases below 0.5, $`e_{IS}`$ decreases rather steeply Sastry98 . The temperature of the inflection point of this decrease coincides with the temperature $`T_p`$ of the peak in the specific heat. We also show $`e_{IS}`$ for a system of 512 particles that was cooled from $`T=1.5`$. At low temperatures $`e_{IS}`$ is rather independent of temperature for the cooled system, indicating that the system is trapped in an energy basin. Such a flattening off at low temperatures is not observed when parallel tempering is used, indicating that the system is able to visit deeper energy basins as the temperature decreases.

However, this does not mean that the peak in the specific heat is an equilibrium feature. We will now show that the peak is the result of not sampling enough of phase space below $`T_p`$. None of the runs we have done at $`T<T_p`$ have reached the time which is necessary in order to achieve an accurate thermodynamic average. To show this, we took a configuration of 512 particles generated by parallel tempering at $`T=0.289855`$, equilibrated for $`5\times 10^7`$ time steps, and then ran for an additional $`10^8`$ time steps during which we recorded the energy at every time step. Then we did block averaging in which we divided our 100 million time steps into equal segments, each of length $`\mathrm{\Delta }t_b`$, and calculated the specific heat from energy fluctuations for each segment Ferrenberg91 . The block averaged specific heat versus time for 2 different time spans $`\mathrm{\Delta }t_b`$ is shown in the inset of Figure 2. Notice that for any given time span, there is no sign of systematic aging. However, the specific heat averaged over time increases with $`\mathrm{\Delta }t_b`$.

The specific heat, averaged over time spans of a given size and over different runs, versus time span size $`\mathrm{\Delta }t_b`$ at several temperatures is shown in Figure 2 by the solid lines. To obtain time spans that are longer than any given run, we concatenated the energies from the runs done at a given temperature to make one huge run, and then did block averaging on the huge run. The results are shown as open symbols in Figure 2. There is good agreement between the solid lines and symbols. Notice that the specific heat initially increases with time span but then levels off when the time span is long enough to exceed the time needed to achieve the thermodynamic value of the specific heat. This time increases with decreasing temperature. Thus at $`T=0.289855<T_p`$, the specific heat continues to increase with time span up to $`\mathrm{\Delta }t_b=10^8`$ time steps which is 100 times longer than the $`\alpha `$ relaxation time $`\tau `$. Equilibrated values of $`C_V`$ near $`T_p`$ are plotted in Figure 1 and lie above the peak in $`C_V`$ found with parallel tempering. Thus the specific heat peak found with parallel tempering is the result of not sampling enough of phase space at $`T<T_p`$ to obtain the true thermodynamic value of the specific heat $`C_V^{\mathrm{true}}`$. The exploration of the energy distribution below $`T_p`$ is slow even with parallel tempering because the probability of sampling large increases in $`U`$ are exponentially small.

Note that $`C_V^{\mathrm{true}}`$ is proportional to the variance $`\sigma _U^2`$ of the distribution $`P(U)`$, i.e., $`C_V^{\mathrm{true}}=N\sigma _U^2/(k_BT^2)`$. If $`\sigma _U^2`$ is finite and if $`n`$ sample values of $`U`$ are statistically independent and identically distributed, then basic statistics dictates that the measured $`C_V`$, which is proportional to the sample variance $`S_n^2`$ of $`U`$, has an expectation value of $`C_V=C_V^{\mathrm{true}}(11/n)`$ Ferrenberg91 . Since the potential energies in a time series can be correlated, the number $`n`$ of statistically independent energies is given by $`n=\mathrm{\Delta }t_b/\tau _U`$ where $`\tau _U`$ is the energy correlation time. Fitting the data that is within 5 to 10% of $`C_V^{\mathrm{true}}`$ in Figure 2 to $`C_V=C_V^{\mathrm{true}}(11/n)`$ yields $`\tau _U3\times 10^6`$, $`1\times 10^5`$, and $`5\times 10^3`$ time steps at $`T=0.2898551`$, 0.308642, and 0.342466, respectively. These values are comparable to the $`\alpha `$ relaxation times $`\tau `$ of $`1\times 10^6`$, $`6\times 10^4`$, and $`2\times 10^4`$ time steps, respectively (see below). They imply that the energies are correlated and that $`n`$ is substantially smaller than the total number of energies.

Note that at shorter time spans and lower temperatures Figure 2 shows that $`C_V\mathrm{ln}(\mathrm{\Delta }t_b)`$. This occurs because the system does not uniformly sample $`P(U)`$ during these shorter time spans. As the system travels through the energy landscape, it samples the energies of each basin that it visits. As shown in Figure 3, we find that the distribution of energies sampled from the basins visited during a shorter time span has a smaller variance than the total distribution $`P(U)`$. Furthermore, the centers of the smaller distributions do not necessarily coincide with the center of the total distribution. Rather the smaller distributions are centered at the inherent structure energy plus the energy of vibrations around $`e_{IS}`$ Schroder00 . In support of this, we show in the inset of Figure 3 that the block averaged energy versus time approximately coincides with $`(e_{IS}(t)+3k_BT/2)`$ vs. time. As more basins are visited, the sample average $`U`$ moves towards the average of the full distribution and the variance grows. This corresponds to $`C_V`$ increasing with time span.

We now cite evidence that systems at temperatures just below $`T_p`$ have equilibrated in the sense of showing no signs of aging. First the inset of Fig. 2 shows the lack of aging in the specific heat versus time for a given value of $`\mathrm{\Delta }t_b`$. Second is the absence of aging in the $`\alpha `$ relaxation time $`\tau `$. We have calculated $`\tau `$ using the full intermediate scattering function $`F_{BB}(\stackrel{}{k},t)=(1/N_B)\rho _\stackrel{}{k}(t)\rho _\stackrel{}{k}(0)`$ for the B particles where $`N_B`$ is the number of B particles and $`\rho _\stackrel{}{k}(t)=_{i=1}^{N_B}\mathrm{exp}(i\stackrel{}{k}\stackrel{}{r}_i(t))`$. We choose $`k=k_{\mathrm{max}}`$, the wavevector of the maximum in the partial static structure factor $`S_{BB}(k)`$ for the B particles, because $`F_{BB}(k,t)`$ relaxes slowest at $`k_{\mathrm{max}}`$ Rinaldi01 . The relaxation time $`\tau `$ is defined by $`F(k,\tau )/F(k,t=0)=1/e`$. We have averaged $`F(k,t)/F(k,t=0)`$ over 40 runs after waiting times $`t_W`$ of $`5\times 10^7`$, 108, and $`1.5\times 10^8`$ time steps for a system with 512 particles at $`T=0.289855<T_p`$ Carruzzo01 . We find that the $`\alpha `$ relaxation time $`\tau =(1.0\pm 0.1)\times 10^6`$ md time steps. This value of $`\tau `$ shows no signs of aging Kob97b  in the sense that there is no systematic variation with waiting time. The lack of aging is to be expected since the aging time is equal to $`\tau `$ which is much less than $`t_W`$. We have confirmed that the aging time is the same as the equilibrium value of $`\tau `$ by starting from 11 different equilibrium configurations at $`T=1.5`$, quenching to $`T=0.289855`$, and measuring $`\tau `$ after waiting times of $`t_W=0`$, $`5\times 10^4`$, $`5\times 10^5`$, $`10^6`$, $`5\times 10^6`$, and $`2.5\times 10^7`$ md steps. For $`t_W<10^6`$ md steps, $`\tau `$ increases with $`t_W`$ which indicates aging Kob97b . However, for $`t_W10^6`$ md steps, there is no aging and $`\tau `$ equals its equilibrium value of $`10^6`$ md steps.

Third we have looked for signs of aging in the inherent structure energy versus time in our single temperature molecular dynamics runs of $`10^8`$ time steps at $`T=0.289855<T_p`$. During the run, configurations were recorded every so often. Averaging over 40 runs, we find no evidence that the inherent structure energy decreases systematically with time, though the noise in the data prevents us from seeing changes smaller than 1%.

We have also examined the root mean square displacement $`\mathrm{\Delta }r^2(t)^{1/2}`$ where $`\mathrm{\Delta }r^2(t)=(1/N)_{i=1}^N(𝐫_i(t)𝐫_i(0))^2`$. We found that $`\mathrm{\Delta }r^2(t)^{1/2}\stackrel{>}{}10\sigma _A`$ in each run examined after $`10^8`$ time steps at $`T=0.289855`$. This is comparable to the box size L = 9.48$`\sigma _A`$ for a system of 512 particles. So the system does not appear to be getting stuck in a metastable minimum of the energy landscape.

Our work is a cautionary tale for those who perform numerical simulations on slowly relaxing systems. It indicates that to obtain accurate thermodynamic averages, one must not only check that the system shows no signs of aging, but one must also check that the quantity to be measured has sampled enough of phase space to obtain a large number of statistically independent values. This sampling time can be orders of magnitude longer than previously recognized time scales such as the aging time and the $`\alpha `$ relaxation time.

We thank Francesco Sciortino, Bulbul Chakraborty, Jon Wellner and Sue Coppersmith for helpful discussions. This work was supported in part by CULAR funds provided by the University of California for the conduct of discretionary research by Los Alamos National Laboratory and by DOE grant DE-FG03-00ER45843.

Present address: Internap, Atlanta, GA 30309.


text:
# Changing Face of the Extrasolar Giant Planet, HD 209458 b

## 1. Introduction

Roughly 100 gaseous giant planets are currently known to orbit nearby sun-like stars999see, e.g., http://www.obspm.fr/encycl/encycl.html and http://exoplanets.org/almanacframe.html. Many of those planets are located at very small orbital distances from their parent stars, where tidal forces are thought to maintain a rotation rate synchronous with the orbit—thus producing permanent day and night sides on the planet. This situation presents a new régime of atmospheric circulation, not encountered in our solar system: slowly-rotating giant planets, which are continuously exposed to intense stellar heating on the same side. Measuring the resulting temperature structure on these planets is a major goal of current and future observational programs. In this Letter, we report on high-resolution, fully-turbulent global simulations of the atmospheric flow on HD 209458 b, presently the only close-in extrasolar giant planet (CEGP) with a measured mass ($`M_p`$) and radius ($`R_p`$).

The parent star of HD 209458 b shows discernible brightness decrements every 3.5 days, due to occultations by the planet as it transits across the star. This property has recently led to precise measurements of $`M_p`$ and $`R_p`$ (Charbonneau et al. 2000; Henry et al. 2000; Mazeh et al. 2000), confirming its giant nature. It has also allowed the detection of sodium absorption, providing the first probe of the planet’s atmosphere (Charbonneau et al. 2002). According to the standard planetary formation picture, HD 209458 b is expected to have formed at a large distance ($`>`$ 1 AU) from its parent star (see, e.g., Boss 1996) and migrated inward (Goldreich & Tremaine 1980; Lin et al. 1996; Murray et al. 1998), quickly ($`<`$ 10 Myr) reaching its present distance of only 0.046 AU from the star (see e.g., Burrows et al. 2000). There, it was forced by tidal effects to permanently present the same face to its star (see, e.g., Goldreich & Soter 1966)—as the Moon does to the Earth. From this synchronization, the rotation period of the CEGP is known (same as its orbital period of 3.5 days). However, unlike our Moon with its insignificant atmosphere, HD 209458 b is expected to possess vigorous meteorology and associated horizontal transport of heat and chemical species, due to the presence of a thin, stable (radiative) atmospheric region above the convective interior (Guillot et al. 1996; Seager & Sasselov 1998).

## 2. Model

We model the stable region of HD 209458 b’s atmosphere as a shallow layer of hydrostatically balanced, frictionless gas—moving under the influence of gravitational and Coriolis accelerations. The motion of such a layer enveloping a planet is governed by the equivalent barotropic formulation of the shallow-water equations on a rotating sphere (Salby 1989):

$`{\displaystyle \frac{𝐯}{t}}+𝐯𝐯`$ $`=`$ $`ghf𝐤\times 𝐯+_a,`$ (1)
$`{\displaystyle \frac{h}{t}}+𝐯h`$ $`=`$ $`𝒦h𝐯+_d,`$ (2)
where $`𝐯=𝐯(\lambda ,\phi ,t)`$ is the horizontal velocity with $`\lambda `$ and $`\phi `$ the longitude and latitude, respectively; $`h=h(\lambda ,\phi ,t)`$ is the thickness of the modeled layer, proportional to the temperature; $`f=2\mathrm{\Omega }\mathrm{sin}\phi `$ is the Coriolis parameter, where $`\mathrm{\Omega }`$ is the rotation rate of the planet; $`𝐤`$ is the local unit vector normal to the planetary surface; $`g`$ is the gravitational acceleration; $`_a=g\eta (1\mathrm{exp}\{t/\tau _a\})\mathrm{cos}(\phi )\mathrm{cos}(\lambda )`$ and $`_d=(hh_E)/\tau _d`$ are, respectively, the adiabatic and diabatic thermal forcing, representing uneven hemispheric heating and cooling due to synchronization; $`\tau _a`$ and $`\tau _d`$ are characteristic $`e`$-folding times; $`h_E=H+\eta \mathrm{cos}(\phi )\mathrm{cos}(\lambda )`$ is the equilibrium thickness; $`\eta `$, and $`H`$ are constants; and, $`𝒦=R/c_p`$, where $`R`$ is the specific gas constant and $`c_p`$ is the specific heat at constant pressure.

The nonlinear equations (1) and (2) describe circulations in which the motion field is vertically aligned over one or more pressure scale heights, $`H_p`$. They are column-integrated representation of the equations used in general circulation models and are solved numerically using the highly-accurate pseudospectral algorithm (Eliasen, Mechenhauer & Rasmussen 1970; Orszag 1970). Approximately 150 simulation runs with up to T341 (1024$`\times `$512 grid) resolution have been performed to explore the full physical and numerical parameter-space available for the HD 209458 b atmosphere. The details of the full exploration are described elsewhere (Cho et al. 2003). Here, we emphasize only the main robust features from the large number of simulations—namely, the emergence of long timescale variability due to moving polar vortices and “thermal” spots of high contrast from the background temperature field.

High-resolution, shallow-layer barotropic models in spherical geometry have been successfully used to model the atmospheric dynamics of the solar system gaseous giant planets (Williams 1978; Cho & Polvani 1996a), as well as stratospheric (stable region) phenomena on the Earth (see, e.g., Juckes & McIntyre 1987). An alternate model of observed cloud level circulation on Jupiter as a surface expression of deep convective columns, oriented parallel to the planetary rotation axis, also exists (Sun et al. 1993; Schubert & Zhang 2000). However, only the shallow-layer model has thus far qualitatively reproduced the essential features (bands, zonal winds, persistent spots, and anticyclonic dominance) on all of the giant planets in our solar system using only the observed values of physical parameters (Cho & Polvani 1996a).

For HD 209458 b, we adopt parameter values where known (Charbonneau et al. 2000; Henry et al. 2000; Mazeh et al. 2000): $`𝒦=0.29`$, $`g=8`$ m s-2, $`R_p=10^8`$ m, $`\mathrm{\Omega }=2.1\times 10^5`$ s-1, and $`H=2\times 10^6`$ m. The last (= $`H_p/𝒦`$) represents the global mean layer thickness, whose value corresponds to the equivalent blackbody temperature ($`T_{\mathrm{atm}}`$ $`1400`$ K for a Bond albedo, $`A=0.2`$; the results presented here, however, are insensitive to the precise choice of $`T_{\mathrm{atm}}`$ and $`A`$). The values for two required parameters—$`\eta `$, the amplitude of hemispheric thermal forcing, and $`U`$, the characteristic wind speed related to the global mean kinetic energy—are not known and hence varied. From run to run, $`\eta `$ is varied from 0 to $`0.4H`$, corresponding respectively to no heating and maximum substellar-point heating in the absence of atmospheric motion; $`U`$ is varied from 50 m s-1 to 1000 m s-1, roughly the observed value for Jupiter (Ingersoll 1990) and the value expected from a simple thermal wind balance (Showman & Guillot 2002), respectively. Note that the latter value leads to a nonlinearly-balanced wind field which can locally approach the sound speed ($``$2700 m s-1); on all observed atmospheres in our solar system, $`U`$ does not exceed $``$400 m s-1 (Beatty & Chaiken 1990). We have carefully studied the effects of $`U`$ on the variability and thermal contrast and find that, while a variability in the flow is a generic feature for HD 209458 b, a significant thermal contrast generally requires a small $`\eta `$ or a large $`U`$. In the following, we give an example of the expected generic flow pattern and then describe a typical flow-induced thermal variability which may be observationally relevant.

## 3. Results

Fig. 1 shows an instantaneous atmospheric flow field on HD 209458 b at dynamic equilibrium from one of our simulations with $`U=400`$ m s-1 and $`\eta `$ corresponding to a maximum day/night temperature difference of 180 K. The flow exhibits considerable complexity and variability: it is characterized by a broad equatorial band containing thin, undulating filamentary structures at low and middle latitudes and displaced, translating circumpolar vortices at the rotation poles. The large undulations at low latitudes, which can roll up into small vortices, are due to high-amplitude Rossby (planetary) waves. In addition to being responsible for weather, the planetary waves efficiently stir and mix tracers (e.g., clouds, chemical species, potential vorticity) when they break, homogenizing the region in which the breakings occur (Juckes & McIntyre 1987). In contrast, polar vortices sequester air inside, their sharp boundaries serving as robust barriers to material transport between their interiors and the well-mixed interstitial regions. The overall flow is similar to that of Earth’s winter stratosphere.

The complex flow of Fig. 1 corresponds to a global circulation pattern of 3 zonal jets, shown in Fig. 2. Jets are atmospheric flow structures, which are obtained by longitudinally averaging the eastward wind at each latitude. The pattern in Fig. 2 naturally evolves from an initially random flow without jets, with the basic profile forming approximately at day 20. The emergence of jets from a random field is a generic feature of turbulent flows in a differentially rotating shallow layer (Rhines 1975; Williams 1978; Cho & Polvani 1996b). For the extreme values of $`U`$ considered with all other parameters held fixed, 2 or 4 jets are also possible. However, once formed, the jets do not change qualitatively—even under applied thermal forcing—for up to several hundred HD 209458 b days, the duration of the simulation.

The presence of a low number of broad jets, along with polar vortices, leads to a circulation pattern which is markedly different from that of cloud-top Jupiter. Although CEGPs are sometimes dubbed “hot Jupiters”, the circulation pattern is more like that of Uranus or Neptune, or even Earth, than that of Jupiter or Saturn, which have $``$10 narrow jets (Ingersoll 1990). Dynamically, this is expected given roughly the factors of ten greater $`T_{\mathrm{atm}}`$ and ten slower rotation rate of HD 209458 b compared to Jupiter. The combined effect, giving the flows in Figs. 1 and 2, is characterized by the Rossby deformation scale, $`L_D`$$`=(𝒦gH)^{1/2}/f`$, which is to be compared with $`R_p`$. Generally, $`L_D`$/$`R_p1`$ is required for a Jupiter-like multiple-jet circulation pattern without pronounced polar vortices at the cloud-top. For HD 209458 b, $`L_D`$/$`R_p1`$.

Jets and vortices/eddies can strongly influence the overall thermal structure by transporting heat. Since the state of the HD 209458 b’s cloud-top is not likely to be at rest and driven solely by stellar heating, we delay the application of thermal forcing (in contrast to Showman & Guillot 2002) until the flow field self-organizes into the pattern in Fig. 2. Dynamically, this is a likely flow state of HD 209458 b in the presence of weak or uniform heating (i.e., at pre-synchronization). The forcing is slowly increased to the specified value of $`\eta `$ over a characteristic radiative equilibrium time, $`\tau _r=10`$ days, appropriate for $``$1 bar pressure level on HD 209458 b. At this depth, $`\tau _r\tau _c`$ (where $`\tau _c=R_p/U`$ is the characteristic circulation time), signifying a large distance advection of airmass before it warms or cools appreciably. The results presented here do not change qualitatively with the value of $`\tau _r`$ ($``$ 0–100 days), using either an adiabatic ($`_a,\tau _a=\tau _r`$) or a diabatic ($`_d,\tau _d=\tau _r`$) representation of thermal forcing.

At each point on the planet, the layer thickness is related to the temperature field, $`T(\lambda ,\phi ,t)=hg/c_p`$, via the hypsometric relation (Holton 1992; Salby 1989). In our simulations with small $`\eta `$ or large $`U`$, the $`T`$ field exhibits a generic feature which is in direct contrast to the simple permanent day/night picture. Fig. 3 shows the $`T`$ distribution (from the same run of Figs. 1 and 2) at several different times. Note that the warmest and the coldest regions are not located at the sub-stellar and the anti-stellar “thermal poles” (at the equator), respectively. Instead, the temperature extrema are located near the rotation poles, inside a coherent hot/cold spot pair at each pole. Nor are the extrema stationary: they revolve around the poles with a period of $``$25 days in this case. The movement leads to a distribution that for a time has the temperature minimum actually on the “hot” day side and the maximum on the “cold” night side. In addition, the temperature difference between the two extrema is large and strongly asymmetric; here, the difference is $``$300 K, with the diffuse hot region $``$70 K above and the sharp cold region $``$230 K below $`T_{\mathrm{atm}}`$. For the case with $`U=1000`$ m s-1, the difference reaches $``$1000 K, with a temperature minimum of $``$800 K below $`T_{\mathrm{atm}}`$ at the core of the cold spots.

## 4. Discussion

The presence of high-contrast hot and cold spots on HD 209458 b induces spatio-temporal variability, which may produce detectable fluctuations in observational signatures. Sensitive enough infrared flux measurements (e.g., such as those possible with SIRTF) could reveal variability in time during an orbit of the planet as different faces are seen from Earth, as well as from orbit to orbit as the spots (polar vortices) revolve about the rotation poles. The lack of detectable variability would thus point toward either an obscuring uniform haze overlying the modeled region or an inefficient conversion of stellar irradiation to atmospheric kinetic energy. In the latter case, the formed spots are weak (i.e., small temperature/thickness perturbation) and thermal forcing which produces bulging in the modeled layer of more than several percent overwhelms any temperature variability due to atmospheric motion. The equilibrium day/night temperature is then robustly maintained.

The extreme conditions inside the spots in our calculations suggest several additional potential observables since absorption levels, albedo, intrinsic thermal emission, and presence, type, or height of clouds could all be different within the different spots. The spatially integrated spectrum, therefore, can be different from the uniform planet case. For example, enhanced levels of CH4/CO abundance inside the cold spots are possible, given the value of $`T_{\mathrm{atm}}`$ used (Seager et al. 2000). Similarly, condensates, such as MgSiO3 (enstatite), may also be found inside the diffuse hot spots, where temperatures may be high enough (Sudarsky et al. 2000; Seager et al. 2000). If the temperature at radiative equilibrium is actually higher than assumed (e.g., $``$2000 K), the opposite situation may occur—MgSiO3 may be found in the cold spots. In addition, by sequestering chemically active species and periodically exposing them to the stellar irradiation (as the spots revolve around the poles), the spots could also affect atomic number densities and could be part of the explanation for the recently observed low abundance of Na I on HD 209458 b (Charbonneau et al. 2002).

## Acknowledgments

J.Y-K.C. acknowledges the hospitality of the Institute for Advanced Study, where the initial part of this work was completed, and R. Levine for helpful discussions. K.M. is supported by NASA under Chandra Fellowship grant PF9-10006, awarded by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-39073. B.M.S.H. is supported by the Hubble Fellowship grant HF-01120.01-99A, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. for NASA under contract NAS 5-26555. S.S. is supported by the W.M. Keck foundation. We thank the referee, W. Hubbard, for comments.


text:
# Nonlinear Behavior in Ferromagnetism: Simple Example and Possible Implications

## Abstract

Two cases of a phenomenological model for ferromagnetism are considered, discrete and continuous. And the relationship, in general, between discrete and continuous models explored. In a similar way to the logistic map behavior, the continuous case is exactly solvable while the discrete one contains the bifurcation route to chaos. Through the ferromagnetic interpretation I comment on the relevance of this to understanding evolution of systems in time, the role of the configuration space in chaotic behavior, and how this understanding may lead to new exotic magnetic phenomena.

It has been known for some time in chaos theory that models behave in different ways when the dynamics is performed using different types of time evolution D89 ; A2 . A simple example of this is the logistic map, for which the continues time case is fully solvable for all finite parameter values while for discrete time the solution ranges from solvable to chaotic in a route known as the route to chaos. In this paper I analyze another simple example of this dual behavior, for a model based on two essential characteristics of ferromagnetism: symmetry breaking of the magnetization at the critical point and dissipating massless dynamics. The model considered is related to the Landau Ginzburg model of magnetism. By keeping the model simple the limit of the continuous case can be solved exactly and since the model contains the essential characteristics of ferromagnetism we have some tools that we can use to interpret the results.

The outline of the article is as follows. The model is introduced and then analyzed with the assumption of continuous time evolution. For this case one finds an initial-condition dependent, predictable relaxation toward the fixed points of the potential energy, for all range of parameters. Next, the basis for discrete time evolution and the relationship between discrete and continuous time evolution are discussed. Finally the discrete case is explored and the same relaxation behavior is found for some parameter regime but new types of behavior (leading to chaos) for other parameter regimes. I finish by interpreting this difference, that results from a different treatment of time, and comment on the importance of this to our understanding of the way processes evolve with time, and indicate systems where such behavior may be found.

The physical picture is as follows. Assume a bulk ferromagnet and let $`(\mathrm{\Phi })`$ be its total free energy, with $`\mathrm{\Phi }`$ the total scalar magnetization in a given direction. If locally $`(\mathrm{\Phi })`$ is a proper physical function, it must have an appropriate expansion as a power series. The sign of $`\mathrm{\Phi }`$ is an artifact of the choice of south magnetic pole, so $`(\mathrm{\Phi })=(\mathrm{\Phi })`$ since reversing this choice should have no physical relevance. If we now take the fourth order approximation and neglect an additive term in the free energy we find,

$$=\frac{a}{2}\mathrm{\Phi }^2+\frac{b}{4}\mathrm{\Phi }^4$$
(1)
with $`a,b`$ the expansion coefficients. To ensure that infinite magnetization is excluded we must have $`b>0`$. Since $`b`$ is a second order parameter it is assumed for the rest that it is a constant. Equation (1) has two characteristic forms as we vary $`a`$, for $`a<0`$ there is a single potential well around $`0`$, we can identify this regime as paramagnetism. In the second regime, for $`a>0`$, there is a double well potential. This broken symmetry regime is identified with ferromagnetism. Indeed for this model we can write $`a(T_cT)`$. Where $`T,T_c`$ are the temperature and the Curie temperature respectively.

Changing the magnetization is a strongly dissipative process that carries little to no inertia. To model this I use dissipative Newtonian dynamics for a massless particle and arrive at the following equation of motion for the magnetization $`\mathrm{\Phi }`$,

$$\dot{\mathrm{\Phi }}=M(a\mathrm{\Phi }b\mathrm{\Phi }^3),$$
(2)
where $`M`$ is the mobility and $`\frac{1}{M}`$ is related monotonically to a damping coefficient also known as the Khalatnikov damping constant. We can integrate this equation immediately to find the following solution for positive initial conditions,

$$\mathrm{\Phi }(t)=\sqrt{\frac{e^{2Mat}}{\frac{b}{a}e^{2Mat}+|\frac{1\frac{b}{a}\mathrm{\Phi }_0^2}{\mathrm{\Phi }_0^2}|}}.$$
(3)
For negative initial conditions we find the solution to be $`\mathrm{\Phi }(t)`$. For the trivial initial condition we find the trivial solution. The non-trivial long time behavior exhibits an asymptote located at the fixed point of the free-energy, $`lim_t\mathrm{}\mathrm{\Phi }(t)=Re(\pm \sqrt{\frac{a}{b}})`$. This solution is completely deterministic for all finite values of parameters and initial conditions, by this it is meant that for any desired accuracy of prediction past some time $`t_0`$, one can give a requirement on the accuracy of the initial conditions at $`t_0`$ which is good for all times $`t>t_0`$. As we shall see shortly this situation is rather different in the discrete case.

Before going into the discrete case we should explore the relationship between it and the continuous case. I claim there is a physical correspondence between the discrete and the continuous models. Usually the continuous solution would be called analytic and appear more pleasing to the trained physical eye. But there is no a-priori reason to assume that systems evolve continuously with time, passing through a continuous number of states. Indeed one can imagine how discrete processes are possible, where systems make jumps from state to state and between jumps remain in their state for some time interval $`\mathrm{\Delta }t`$.

To relate continues and discrete models one needs a mesh. For simplicity I assume a mesh with constant $`\mathrm{\Delta }t`$, i.e. we have a discrete linear relationship for the magnetization,

$$\mathrm{\Phi }_n=\mathrm{\Phi }(t_0+n\mathrm{\Delta }t).$$
(4)
In general an $`m^{\text{th}}`$ order differential model will be given by the condition $`f(\mathrm{\Phi },\dot{\mathrm{\Phi }},\mathrm{},\mathrm{\Phi }^{(m)}(t))=0`$, with $`m`$ initial conditions. Here, $`\mathrm{\Phi }^{(n)}(t)=\frac{d^n\mathrm{\Phi }}{dt^n}`$. A general $`m^{\text{th}}`$ order difference model is given by, $`g_{\mathrm{\Delta }t}(\mathrm{\Phi }_n,\mathrm{\Phi }_{n+1},\mathrm{},\mathrm{\Phi }_{n+m})=0`$, with $`m`$ initial conditions.

Let $`g_{\mathrm{\Delta }t}=0`$ be the difference model that one gets when discretizing a differential model $`f=0`$, by approximating the derivative $`\dot{\varphi }\frac{\mathrm{\Phi }_{n+1}\mathrm{\Phi }_n}{\mathrm{\Delta }t}`$ and similarly 111There is some ambiguity in the choice of the discretization of the differentials. Since this ambiguity involves only translation of the order of $`\mathrm{\Delta }t`$ in time I assume that we can safely neglect it. for higher derivatives. And let $`C_{\mathrm{\Delta }t}`$ the discretizing transformation. Symbolically this can be writen as $`C_{\mathrm{\Delta }t}(f=0)=(g_{\mathrm{\Delta }t}=0)`$. From analysis we know that $`C_{\mathrm{\Delta }t}Id`$ when $`\mathrm{\Delta }t0`$, where $`Id`$ is the identity. From this general property one can show that for small enough $`\mathrm{\Delta }t`$, the $`C_{\mathrm{\Delta }t}`$ mapping preserves the physical essence of the model. Stated more precisely, for a given mesh of time there is a unique correspondence between the differential and the difference models. Because we know that going from the difference model by taking the limit $`\mathrm{\Delta }t0`$ one gets the corresponding differential model. If we show that this process is invertible, we can claim that the differential and difference equations of motion correspond to the same physical system with continuous or discrete evolution respectively. This invertiblity for $`1^{st}`$ order systems is now shown.

Assume for simplicity that a differential and a difference model are given by $`\dot{\mathrm{\Phi }}=f(\mathrm{\Phi }),\mathrm{\Phi }_{n+1}=g(\mathrm{\Phi }_n)`$ respectively, where $`f,g`$ are well behaved functions. For a given mesh size $`\mathrm{\Delta }t>0`$, let $`\mathrm{\Phi }_n`$ be given by equation (4) and let $`C_{\mathrm{\Delta }t}(\dot{\mathrm{\Phi }}=f(\mathrm{\Phi }))=(\frac{\mathrm{\Phi }_{n+1}\mathrm{\Phi }_n}{\mathrm{\Delta }t}=f(\mathrm{\Phi }_n))`$ be the correspondence between differential and difference equations. Now if,

$$C_{\mathrm{\Delta }t}(\dot{\mathrm{\Phi }}=f_1(\mathrm{\Phi }))=C_{\mathrm{\Delta }t}(\dot{\mathrm{\Phi }}=f_2(\mathrm{\Phi }))$$

it is easy to see that $`f_1=f_2`$, just by allowing the initial conditions to vary continuously, i.e. $`C_{\mathrm{\Delta }t}`$ is 1-1. To show onto, we can take $`f(\mathrm{\Phi })=\frac{1}{\mathrm{\Delta }t}[g(\mathrm{\Phi })\mathrm{\Phi }]`$ where $`g`$ is the difference model. Then,

$$C_{\mathrm{\Delta }t}(\dot{\mathrm{\Phi }}=f(\mathrm{\Phi }))=(\mathrm{\Phi }_{n+1}=g(\mathrm{\Phi }_n)).$$

This shows that, for a given mesh of time, for every difference model there exists a unique differential model and vice versa. And thus, any changes that a given model will exhibit when going from continuous to discrete time will result directly from this discretization and not from any other change to the underlying physical dynamics.

I now show the discrete case for this model of magnetism. If we discretize equation (2) we find the following difference equation for the magnetization,

$$\mathrm{\Phi }_{n+1}=g(\mathrm{\Phi }_n)=\mathrm{\Phi }_n+c_1[\mathrm{\Phi }_nc_2(\mathrm{\Phi }_n)^3]$$
(5)
where $`c_1=Ma\mathrm{\Delta }t`$ and $`c_2=\frac{b}{a}`$. $`\mathrm{\Delta }t`$ is a finite time element.

The fixed points for this map, given by the condition $`\mathrm{\Phi }_{n+1}=\mathrm{\Phi }_n`$ , are $`\mathrm{\Phi }^{}=0`$ or $`\mathrm{\Phi }^{}=\pm \frac{1}{\sqrt{c_2}}`$. To understand the nature of these fixed points their stability is analyzed,

$$g^{}(\mathrm{\Phi }^{})=1+c_1[13c_2(\mathrm{\Phi }^{})^2].$$

Or,

$$g^{}(0)=1+c_1,$$

$$g^{}(\pm \sqrt{\frac{1}{c_2}})=12c_1.$$

Thus for $`2<c_1<0`$ the $`\varphi ^{}=0`$ fixed point is stable which corresponds to the paramagnetic state. For $`0<c_1<1`$, the $`\mathrm{\Phi }^{}=\pm \sqrt{\frac{1}{c_2}}`$ is a stable fixed point for positive and negative initial conditions respectively. We recognize this as the ferromagnetic state. For larger values of $`c_1`$ the fixes points become unstable and we find higher period fixed points. If we set $`c_1`$ even larger we find a dense set of unstable periodic orbits with total measure zero D89 . This last scenario is known as chaos and the route we have taken is the period doubling route to chaos.

Remembering that $`c_1=Ma\mathrm{\Delta }t,M>0`$ and that for a forward evolving system $`\mathrm{\Delta }t>0`$, we can conclude again that $`a`$ is the parameter that determines whether the system is a paramagnet or a ferromagnet. Note that for a backward evolving system $`(\mathrm{\Delta }t<0)`$ the roles of paramagnetism and ferromagnetism change and we get the transition from ferromagnetism to paramagnetism as we lower the temperature.

The standard conclusion derived from such a treatment is that as long as we keep $`\mathrm{\Delta }t`$ small enough such that $`c_1`$ is in the stable region $`(2,1)`$ we can analyze equation (2) with equation (5). But there is some thing more fundamental that we can conclude: the nature of time evolution - whether discrete or continuous, leads to different observation in physical systems.

With this last conclusion in mind we are led quite directly to inferring from observations on the nature of ferromagnetism under discrete time evolution. For the system presented we consider two control parameters, the temperature and the mobility. Figure 1 is the bifurcation diagram as we vary the temperature, i.e. for a given temperature we iterate equation (5) many time steps (200) and we plot the final states it settles to (100 points). We start in the paramagnetic state with zero magnetization then at $`T_c`$ ($`a=0`$) we have a transition to the ferromagnetic state. As we keep increasing $`a`$ we have more transitions until we reach the chaotic state. So, if time is discrete then we might expect that if we lower the temperature enough we will get a transition from the ferromagnetic state to a state for which the magnetization oscillates between two adjacent magnetization values. If we lower the temperature even further the magnetization difference between these two states will increase. However, it could turn out that time evolution is discrete and we would not observe such exotic magnetic states. The reason is that there is a physical bound on the value of $`a`$, set by the minimum possible temperature, $`T=0K`$. It could turn out, especially if we expect a small $`\mathrm{\Delta }t`$, that in most cases the corresponding value of $`c_1`$ is well within the familiar ferromagnetic region.

The second control parameter is the mobility. Figure 2 is the bifurcation diagram keeping $`a,b`$ fixed, positive and varying the mobility. I remark that the part of the diagram for negative values of the mobility was left in mainly for esthetic reasons. As we increase the mobility we observe much the same bifurcation diagram as before. Indeed it is seen that there is a value of $`M`$ for which the familiar ferromagnetic state ends and we start a new phase for which the magnetization oscillates between nearby states. Moreover there seems to be no natural cutoff for the mobility and detection of these exotic magnetic states seems more plausible A1 . Another system where detection of exotic states may be possible is discussed in A2 .

In the frame-work of magnetization these nonlinearities are intriguing and could also provide evidence about the nature of time evolution. Indeed as this and other work suggest, bifurcation routes and chaos rise naturally from discrete dynamics, and disappear once the continuous limit is taken. Interestingly enough, experimental and numerical work on a driven magnetic system do exhibit a period-doubling route to chaos A1 ; A3 . Another possible system where the exotic nature of the magnetization may be found is in a super-critical ferromagnet where the critical temperature (and thus $`a`$) is very large compared to everyday ferromagnets. In such systems the distance between nearest neighbor magnetic moments is small such that the overlap between adjacent wave functions and hence the exchange integral are greatly increased. An example of such a system is a neutron star, where the distances between adjacent magnetic moments are of nuclear order. Such stars exhibit a large-scale ferromagnetic ordering at very high temperatures which gives a very high critical temperature. When such a star cools down $`a`$ becomes large and the possibility of exotic states arises in the realm of discrete time. Such exotic states will oscillate between different magnetization states with period of the order $`\mathrm{\Delta }t`$. If our instruments make observations over a time much larger than $`\mathrm{\Delta }t`$ we would observe only the average magnetization.

In conclusion, a simple model for the bulk magnetic state has been studied. It is found that under continuous time evolution assumption this model admits only a steady-state analytical solution that is fully deterministic for all model-parameters values. Under discrete time evolution assumption the model exhibits the same analytical solution for some parameter values but as part of the period doubling route to chaos. This different behavior resulted from the different assumption with regard to time evolution and thus such observation can aid in understanding the nature of time and of time evolution.

###### Acknowledgements.

The author would like to acknowledge the support generously provided by Natali J. Teszler.