L D C Jaubert

M A M De Aguiar

Semiclassical Tunneling of Wavepackets with Real Trajectories

Keywords: numbers: 03, 65, Sq,31, 15, Gy

Semiclassical approximations for tunneling processes usually involve complex trajectories or complex times. In this paper we use a previously derived approximation involving only real trajectories propagating in real time to describe the scattering of a Gaussian wavepacket by a finite square potential barrier. We show that the approximation describes both tunneling and interferences very accurately in the limit of small Planck's constant. We use these results to estimate the tunneling time of the wavepacket and find that, for high energies, the barrier slows down the wavepacket but that it speeds it up at energies comparable to the barrier height.

I. INTRODUCTION

The success of semiclassical approximations in molecular and atomic physics or theoretical chemistry is largely due to its capacity to reconcile the advantages of classical physics and quantum mechanics. It manages to retain important features which escape the classical methods, such as interference and tunneling, while providing an intuitive approach to quantum mechanical problems whose exact solution could be very difficult to find. Moreover, the study of semiclassical limit of quantum mechanics has a theoretical interest of its own, shedding light into the fuzzy boundary between the classical and quantum perspectives.

In this paper we will apply the semiclassical formalism to study the scattering of a one dimensional wavepacket by a finite potential barrier. In the case of plane waves, the tunneling and reflection coefficients can be easily calculated in the semiclassical limit, giving the well known WKB expressions [START_REF] Merzbacher | Quantum Mechanics, cap. 7[END_REF]. For wavepackets, however, the problem is more complicated and few works have addressed the question from a dynamical point of view [START_REF] Xavier | [END_REF]3,4]. The time evolution of a general wavefunction with initial condition ψ(x, 0) = ψ 0 (x) can be written as

ψ(x, T ) =< x|K(T )|ψ 0 >= < x|K(T )|x i > dx i < x i |ψ 0 >, (1)
where K(T ) = e -iHT / is the evolution operator and H is the (time independent) hamiltonian. The extra integration on the second equality reveals the Feynman propagator < x|K(T )|x i >, whose semiclassical limit is known as the Van-Vleck formula [START_REF] Van Vleck | Proc. Nat. Acad. Sci. USA[END_REF] (see next section). When the Van-Vleck propagator is inserted in Eq.( 1) we obtain a general semiclassical formula which involves the integration over the 'initial points' x i :

ψ sc (x, T ) = < x|K(T )|x i > Van Vleck dx i < x i |ψ 0 > . (2)
If this integral is performed numerically one obtains very good results, specially as goes to zero. However, doing the integral is more complicated than it might look, because for each x i one has to compute a full classical trajectory that starts at x i and ends at x after a time T , which may not be simple task. Alternative methods involving integrals over initial conditions (instead of initial and final coordinates) in phase space have also been developed and shown to be very accurate [START_REF] Kenneth | [END_REF]7,8]. All these approaches sum an infinite number of contributions and hide the important information of what classical trajectories really matter for the process.

In a previous paper [9] several further approximations for this integral were derived and applied to a number of problems such as the free particle, the hard wall, the quartic oscillator and the scattering by an attractive potential. The most accurate (and also the most complicated) of these approximations involves complex trajectories and was first obtained by Heller and collaborators [10,11]. The least accurate (and the simplest to implement) is known as the Frozen Gaussian Approximation (FGA), and was also obtained by Heller [12]. It involves a single classical trajectory starting from the center of the wavepacket.

However, other approximations involving real trajectories can be obtained [4,9,13]. These are usually not as accurate as the complex trajectory formula, but are much better than the FGA and can be very good in several situations. Moreover, it singles out real classical trajectories from the infinite set in Eq.( 2) that can be directly interpreted as contributing to the propagation.

In this paper we apply these real trajectory approximations to study the tunnel effect.

Since this is a purely quantum phenomena, it is a very interesting case to test the semiclassical approximation and to understand what are the real trajectories that contribute when the wavepacket is moving 'inside' the barrier. More specifically, we will consider the propagation of a Gaussian wavepacket through a finite square barrier. We shall see that the semiclassical results are very accurate, although some important features of the wavepacket propagation cannot be completely described.

This paper is organized as follows: in the next section we review the semiclassical results derived in [9], which are the starting point of this work. Next we describe the evolution of a Gaussian through a square potential barrier in its three separate regions: before, inside and after the barrier. Finally in section IV we discuss the calculation of tunneling times, as proposed in [START_REF] Xavier | [END_REF]. We find that the barrier slows down the wavepacket at high energies, but that it speeds it up at energies comparable to the barrier height. Finally, in section V we present our conclusions.

II. APPROXIMATION WITH COMPLEX AND REAL TRAJECTORIES

One important class of initial wavefunctions is that of coherent states, which are minimum uncertainty Gaussian wavepackets. In this paper we shall consider the initial wavepacket |ψ 0 > as the coherent state of a harmonic oscillator of mass m and frequency ω defined by

|z = e -1 2 |z| 2 e zâ † |0 , (3)
where |0 is the harmonic oscillator ground state, â † is the creation operator and z is the complex eigenvalue of the annihilation operator â with respect to the eigenfunction |z .

Using the position and momentum operators, q and p respectively, we can write

â † = 1 √ 2 q b -i p c z = 1 √ 2 q b + i p c , (4)
where q and p are real numbers. The parameters b = ( /mω)

1 2 and c = ( mω) 1 2
are the position and momentum scales respectively, and their product is .

In order to write the Van Vleck formula of the Feynman propagator, we need to introduce the tangent matrix. Let S ≡ S(x f , T ; x i , 0) be the action of a classical trajectory in the phase space (X, P ), with x i = X(0) and x f = X(T ). A small initial displacement (δx i , δp i ) modifies the whole trajectory and leads to another displacement (δx f , δp f ) at time T . In the linearized approximation, the tangent matrix M connects these two vectors of the phase space

      δx f b δp f c       =        - S ii S if - c b 1 S if b c S if -S f f S ii S if - S f f S if              δx i b δp i c       ≡      m qq m qp m pq m pp            δx i b δp i c       (5)
where

S ii ≡ ∂ 2 S/∂x 2 i , S if ≡ ∂ 2 S/∂x i ∂x f ≡ S f i and S f f ≡ ∂ 2 S/∂x 2 f .
In terms of the coefficients of the tangent matrix, the Van Vleck propagator is [START_REF] Van Vleck | Proc. Nat. Acad. Sci. USA[END_REF] x

f |K(T )|x i V anV leck = 1 b 2πm qp exp i S(x f , T ; x i , 0) -i π 4 . (6)
For short times m qp is positive and the square root is well defined. For longer times m qp may become negative by going through zero. At these 'focal points' the Van Vleck formula diverges. However, sufficiently away from these points the approximation becomes good again, as long as one replaces m qp by its modulus and subtracts a phase π/2 for every focus encountered along the trajectory. We shall not write these so-called Morse phases explicitly.

Assuming some converging conditions, the stationary phase approximation allows us to perform the integral over x i in Eq. ( 2) (for more details, see [9]). We obtain

ψ(z, x f , T ) sc = b -1/2 π -1/4 m qq + im qp exp i S(x f , T ; x 0 , 0) + i p(x 0 -q/2) - (x 0 -q) 2 2b 2 , (7)
where x 0 is the value of the initial coordinate x i when the phase of the propagator is stationary. It is given by the relation

x 0 b + i p 0 c = q b + i p c where p 0 = - ∂S ∂x i x 0 . (8)
The end point of the trajectory is still given by X(T ) = x f . In spite of q and p being real, x 0 and p 0 are usually complex. This implies that the classical trajectories with initial position

x 0 and momentum p 0 are complex as well, even with x f ∈ R. Eq. ( 7) was first obtained by Heller [10,11] and it is not an initial value representation (IVR). There are a priori several complex trajectories satisfying the boundary conditions. Thanks to the stationary phase approximation, we were able to replace an integral over a continuum of real trajectories (2) by a finite number of complex ones (7). The problem is now solvable, but still quite difficult to compute. However, it turns out that, in many situations, these complex trajectories can be replaced by real ones, which are much easier to calculate [4,9].

Therefore, we look for real trajectories that are as close as possible to the complex ones.

Let (X(t), P (t)) ∈ C × C be the coordinates of a complex trajectory, and (u(t), v(t)) a new set of variables defined by

u = 1 √ 2 X b + i P c , v = 1 √ 2 X b -i P c . (9)
According to Eq. ( 8), the boundary conditions become

u(0) = 1 √ 2 x 0 b + i p 0 c = 1 √ 2 q b + i p c = z and X(T ) = x f . (10)
The initial condition is then the complex coordinate z and the final condition is the real position x f . The real and imaginary parts of z are related to the central position q and the central momentum p of the initial wavepacket respectively. This gives us three real parameters that we may use as boundary conditions to determine the real trajectory. But a particle whose initial conditions are q and p will not a priori reach x f after a time T .

Although it is possible to satisfy such final condition, it will not usually happen because X(T ) is imposed by q and p. Likewise, fixing the initial and final positions q and x f will not generally lead to P (0) = p. Therefore we need to choose only two boundary conditions among the three parameters, and use the hamiltonian of the system to calculate analytically or numerically the third one. This means that the relation (8) will not be generally fulfilled and the hope is that it will be fulfilled approximately. For a discussion about the validity of this approximation, see the beginning of the third section in [9]. If we fix (q, p), we obtain the Frozen Gaussian Approximation of Heller [12]. This is an initial value representation that involves a single trajectory and is unable to describe interferences or tunneling, which are the aim of this paper. However, we can fix X(0) = q , X(T ) = x f and calculate P (0) = p i .

When the complex quantities in Eq.( 7) are expanded about this real trajectory we obtain [9] ψ(z, x f , T

) sc = b -1/2 π -1/4 m qq + im qp exp i S(x f , T ; q, 0) + i 2 pq - 1 2
im qp m qq + im qp p -p i c 2 . ( 11
)
Eq. ( 11) is the semiclassical formula we are going to use in this paper. We shall show that, although still very simple, it can describe tunneling and interferences quite well.

III. THE 1-D SQUARE BARRIER

Consider the specific case of a particle of unit mass scattered by the 1-D square barrier defined by (see fig. 1)

V (x) =      V 0 if x ∈ [-a, a] where a ∈ R + 0 otherwise . ( 12
)
The initial state of the particle is a coherent state ψ(z, x, 0) = x|z with average position q < -a and average momentum p > 0, i.e., the wavepacket is at the left of the barrier and moves to the right. In all our numerical calculations we have fixed V 0 = 0.5 and defined the

critical momentum p = √ 2V 0 = 1.
The application of the semiclassical formula Eq.11 requires the calculation of classical trajectories from q to x f in the time T . For the case of a potential barrier, the number of such trajectories depends on the final position x f . This dependence, in turn, causes certain discontinuities in the semiclassical wavefunction.

Since the initial wavepacket starts from q < -a, it is clear that for x f > a (at the right side of the barrier) there is only one trajectory satisfying x(0) = q and x(T ) = x f . This 'direct trajectory' has p i > √ 2V 0 and x(t) increases monotonically from q to x f . For x f < -a, on the other hand, in addition to the direct trajectory there might also be a reflected trajectory, that passes through x f , bounces off the barrier and returns to x f in the time T . The initial momentum of such a reflected trajectory must be greater than that of the direct one, since it has to travel a larger distance. However, if this distance is too big, i.e., if x f << -a, the initial momentum needed to traverse the distance in the fixed time T becomes larger than √ 2V 0 and the reflected trajectory suddenly ceases to exist (see next subsection for explicit details for the case of the square barrier and figure 2 for examples).

This qualitative discussion shows that reflected trajectories exist only if x f is sufficiently close to the barrier. The points where these trajectories suddenly disappear represent discontinuities of the semiclassical calculation. Fortunatelly, this drawback of the approximation becomes less critical as goes to zero, since the contribution of the reflected trajectory at those points become exponentially small as compared to the direct one (see for instance figure 2(g)).

In the remaining of this paper we are going to obtain explicit expressions for ψ(z, x f , T ) sc before, inside and after the barrier. For fixed q we will calculate the classical trajectories for each x f , extracting the initial momentum p i , the action S ≡ S(x f , T ; q, 0) and its derivatives (in order to obtain m qq and m qp ).

A. Before the barrier:

x f < -a
The specificity of this region is that there may exist two different paths connecting q to

x during the time T : a direct trajectory and a reflected one (fig. 1) whose initial momenta, action and tangent matrix elements are given by

p i d = x -q T ; S d = (x -q) 2 2T ; m qq d = 1; m qp d = T λ , (13)
p i r = - x + q + 2a T ; S r = (x + q + 2a) 2 2T ; m qq r = -1; m qp r = - T λ , (14)
where λ = b/c. The contribution of each trajectory to the wavefunction at x f , ψ d and ψ r , is

ψ d = b -1/2 π -1/4 √ 1+i T λ exp i (x f -q) 2 2T + i 2 pq - 1 2 iT λ + iT pT -x f + q cT 2 , ψ r = b -1/2 π -1/4 √ 1+i T λ exp iθ + i (x f + q + 2a) 2 2T + i 2 pq - 1 2 iT λ + iT pT + x f + q + 2a cT 2 . ( 15
)
Notice that we have added an extra phase θ in ψ r . Without this extra phase (that includes the minus sign coming from the tangent matrix elements in Eq.( 14)), the wavepacket would not be continuous as it goes through the barrier. For a hard wall, for instance, we impose θ = π to guarantee that the wavefunction is zero at the wall. For smooth barriers this phase would come out of the approximation automatically, but for discontinuous potentials we need to add it by hand. To calculate θ we rewrite the previous expressions in complex polar representation,

ψ d = D(x f )e iϕ d (x f ) , ψ r = R(x f )e iϕr(x f )+iθ
|ψ(z, x f , T ) sc | 2 = 1 b √ π 1 1 + T 2 λ 2 exp - λ 2 λ 2 + T 2 x f -q -pT b 2 + exp - λ 2 λ 2 + T 2 x f + q + pT + 2a b 2 + 2 cos 2(x f + a) (λ 2 + T 2 ) λ 2 p -(q + a)T -θ exp - λ 2 λ 2 + T 2 (pT + q + a) 2 + (x f + a) 2 b 2 . ( 17
)
This is the same result as obtained in [9] for a completely repulsive barrier (V 0 → ∞), except for the phase, because of the different boundary condition at

x f = -a (|ψ(-a) sc | 2 = 0
for the hard wall). However, as discussed in the beginning of this section, an additional difficulty appears when the wall is finite: the reflected trajectory does not always exist.

From the classical point of view, there is no reflected part if the energy E = p 2 i r /2 > V 0 . The maximum initial momentum allowed is then √ 2V 0 and a particle with such momentum takes the time T c = -a + q √ 2V 0 to reach the barrier. Furthermore, for T > T c the reflected trajectory only exists if p i r = -

x f + q + 2a T √ 2V 0 i.e. if |x f | = -x f x c = q + 2a + √ 2V 0 T .
Therefore, if T T c and |x f | x c , the probability density is given by Eq.( 17), otherwise we only have the contribution of the direct ψ d and

|ψ(z, x f , T ) sc | 2 = 1 b √ π 1 1 + T 2 λ 2 exp - λ 2 λ 2 + T 2 x f -q -pT b 2 . ( 18
)
As a final remark we note that the calculation of θ might involve a technical difficulty depending on the values of , p and T . For some values of these parameters the contribution of the direct and reflected trajectories might become very small at x=-a (see for instance fig. 2(f), which shows the reflected wavepacket in a case of large transmission). In these cases the probability density becomes very small at x=-a and the value of the phase θ is irrelevant. In some of these situations, where the value of θ does not affect the results, we actually found that cos(theta) = A 2 /2 -1 > 1, which cannot be solved for real θ. For the sake of numerical calculations we have set θ = 0 in these cases.

The semiclassical wavepacket is now completely described for x f < -a. The probability density |ψ sc | 2 is a function of q, p, x f , T and depends on several parameters, a, b, and V 0 .

In our numerical calculations we fixed a = 50. This makes the barrier large enough so whereas the exact one is decreasing continuously. On the one hand this means that the approximation is not perfect but, on the other hand, the semiclassical approximation explains that the fast rundown of the exact quantum wavepacket comes from the progressive disappearance of the reflected classical trajectory due to the finite size of the barrier. Finally, for = 0.1, Fig. 2(g), the approximation becomes nearly perfect. As expected, the semiclassical approximation works better and better when is decreasing, i.e. when the quantum rules give way to the classical ones.

To end this subsection, we mention that the quality of the approximation is independent From the classical point of view there is only the direct trajectory in this region (see Fig. 1), since a reflection on the other side of the barrier (at x = a) can not be considered without quantum mechanics. Calling p 1 = p i > √ 2V 0 and p 2 the momentum of this trajectory before and inside the barrier respectively, energy conservation gives p 2 1 /2 = p 2 2 /2 + V 0 . This is the first equation connecting p 1 to p 2 , but we need a second one which is imposed by the propagation time T = t 1 + t 2 where: t 1 = -a+q p 1 is the time to go from q to -a with momentum p 1 ; t 2 = x f +a p 2 is the time to go from -a to x f with momentum p 2 .

The combination of these two equations gives

T = - a + q p 1 + x f + a p 2 1 -2V 0 (19)
which can be rewritten as

(p 2 1 -2V 0 )(p 1 T + a + q) 2 = (x f + a) 2 p 2 1 . (20)
This is a quartic polynomial, which we solve numerically. We obtain four solutions: one is always negative, which we discard since we fixed the initial position q on the left side of the barrier; two are sometimes complex and, when real, have p 1 < 1; finally, one of the roots is always real, larger than 1 and tends to

x f -q T when V 0 is negligible (the limit of a free particle). We take this last root as the initial momentum p 1 .

The action S is also a function of p 1 given by

S(z, x f , T ) = t 1 0 p 2 1 2 dt + T t 1 p 2 2 2 -V 0 dt = p 2 1 2 t 1 + p 2 2 2 -V 0 t 2 = - 1 2 (a + q)p 1 + 1 2 (x f + a) p 2 1 -2V 0 - V 0 (x f + a) p 2 1 -2V 0 . (21)
We calculate the derivatives of S numerically by computing p 1 and S for the initial conditions (q, x f ), (q + dq, x f ), (q, x f + dx f ) . . . and approximate ∂S ∂x f (z, x f , T ) by [S(z, x f + dx f , T ) -S(z, x f , T )]/dx f , etc. Finally, the propagator inside the barrier is given by Eq.( 11) plus the phase correction ξ calculated in the previous subsection. The probability density, which in independent of ξ, becomes

|ψ(z, x f , T ) sc | 2 = 1 b √ π 1 m 2 qq + m 2 qp exp - m 2 qp m 2 qq + m 2 qp p -p 1 c 2 . ( 22
)
Figure 3 shows |ψ sc | 2 as a function of x f for the same parameters as in subsection III A.

Although the semiclassical approximation also improves for small , here we shall fix = 1. This is because the behavior of the propagator becomes trivial for small : if p < 1 the wavepacket bounces off the barrier almost completely, and otherwise it simply passes over the barrier barely noticing the presence of the potential.

The first remark is that the wavepacket is continuous at x f = -a: the extra phase θ does play its role correctly. As in the case before the barrier, the comparison between exact and semiclassical calculations is always at least qualitatively good, and sometimes even quantitatively so. However, there are two main effects that the semiclassical approximation cannot take into account.

1. there is a gap between the exact and semiclassical curves, which decreases progressively as x f increases, and is bridged near the local maximum of the probability density.

The reason may come from the fact it is not possible to impose the continuity of the derivative of ψ sc with respect to x f at -a.

2. there are oscillations on the exact curve (especially for p = 2 and T = 50) close to the right side of the barrier, that are not present in the semiclassical approximation. This is a purely quantum effect, because classical mechanics can not account for a reflected trajectory which would interfere with the direct one in this case. |ψ sc | 2 is in fact the mean-value of the oscillations, and that is why there is a discontinuity of the wavepacket at x f = a, since the exact curve is beginning at the bottom of an oscillation.

If we want to stay strictly in the semiclassical limit, there is nothing we can do about the lack of interferences in the barrier region: this is the limit of our approximation. But if we want to use the semiclassical point of view in order to provide a more intuitive picture of the quantum world, we can add a 'ghost' trajectory that reflects at x f = a and see if it can account for the interferences. Similar ideas have been applied to the frequency spectrum of microwave cavities with sharp dielectric interfaces [14] and, more recently, to the spectrum of step potentials confined by hard walls [15]. The argument will be the same as in subsection III A, except of course that the reflected trajectory will now bounce on the right side of the barrier. The equation for p i = p 1 is again a quartic polynomial given by

(p 2 1 -2V 0 )(p 1 T + a + q) 2 = (3a -x f ) 2 p 2 1 . ( 23
)
We know that p 1 direct is the same as p 1 ref lected at x f = a and we choose the only solution of (23) which satifies this condition. The expression of the new action is:

S r (z, x f , T ) = - 1 2 (a + q)p 1 + p 2 1 2 -2V 0 3a -x f p 2 1 -2V 0 . ( 24
)
The expressions of ψ d and ψ r are the same as eq. ( 11) but with p i , S, m qq and m qp indexed by d or r. After some calculations, the new expression of the probability density inside the barrier becomes

|ψ(z, x f , T ) sc | 2 = 1 b √ π 1 m 2 qq d + m 2 qp d exp - m 2 qq d m 2 qq d + m 2 qp d p -p 1 d c 2 + 1 b √ π 1 m 2 qq r + m 2 qp r exp - m 2 qq r m 2 qq r + m 2 qp r p -p 1 r c 2 + 2 b √ π cos (ϕ r -ϕ d + θ ′ ) 1 4 m 2 qq r + m 2 qp r m 2 qq d + m 2 qp d × exp - 1 2 m 2 qq d m 2 qq d + m 2 qp d p -p 1 d c 2 - 1 2 m 2 qq r m 2 qq r + m 2 qp r p -p 1 r c 2 , ( 25
)
where θ ′ is the new extra phase (that absorbs the previously computed ξ) and

ϕ r -ϕ d = S r -S d + 1 2 arctan m qp d m qq d - 1 2 arctan m qp r m qq r + 1 2
m qq d m qp d m 2 qq d + m 2 qp d p -p 1 d c 2 - 1 2
m qq r m qp r m 2 qq r + m 2 qp r p -p 1 r c 2 . ( 26
)
The results of such an expression, however, are not good: the oscillations become too big, which means that the reflected trajectory needs to be attenuated by a reflection coefficient ρ. To calculate ρ we use the following reasoning: for each point x f inside the barrier there corresponds a reflected trajectory from q to x f with a certain value of p 1 > 1 computed with Eq. ( 23). We take for ρ the same attenuation coefficient a plane wave with momentum p 1 would have. Let (F e iκx f + G e -iκx f ) and C e ikx f be such a plane wave inside and after the barrier respectively, where κ = 2(E -V 0 )/ = p 2 1 -p2 / and k = √ 2E/ = p 1 / . The continuity of this function and its derivative at x f = a give us the relative weight of the reflected trajectory with respect to the direct one:

ρ(p 1 , V 0 ) = G F = 1 -κ/k 1 + κ/k = 1 -1 -p2 /p 2 1 1 + 1 -p2 /p 2 1 . ( 27
)
The expression for the total propagator becomes ψ sc d + ρ ψ sc r e iθ ′ . We use the same argument as in subsection III A to compute the extra phase θ ′ , adding another correction ξ ′ to the wavefunction on the right side of the barrier. Because there is always a single trajectory on the right side, ξ ′ does not affect the probability density there. We find that cos θ ′ = A 2 + /2 -1 where A + = W (a)/D(a). The new results are displayed in figure 4. The gap is still present, but the agreement between exact and semiclassical on the right side is nearly perfect. The interferences are indeed coming from a real 'ghost' trajectory that bounces off at the end of the barrier like a quantum plane wave. Since the left side of the figure has not changed much, the reflected trajectory has no effect on this part of the wavepacket and we don't need to consider additional reflections. Furthermore, we don't have to take ρ ψ sc r into account when we calculate θ in subsection III A. We finish this subsection with two comments: first, the approximation with the ghost trajectory is accurate even for larger values of . Second, the wavepacket becomes continuous at x f = a. That is very interesting because continuity comes only when we include the reflected trajectory, whereas the part of the wavepacket which goes through the barrier is calculated independently with a single direct trajectory (see next subsection). This means that the semiclassical propagator after the barrier somehow knows there is a reflected part.

In the next subsection, we will briefly present the computation of the wavefunction at the right side of the barrier.

C. After the barrier: a < x Following the same arguments as in subsection III B, we use the energy conservation

p 2 1 /2 = p 2 2 /2 + V 0 = p 2 3
/2 (the index 3 refers to the right of the barrier) and the different times t 1 = -a+q p 1 , t 2 = 2a p 2 and t 3 = x-a p 3 to calculate the initial momentum of the direct trajectory. We obtain

(p 2 1 -2V 0 )(p 1 T + 2a + q -x) 2 = (2a) 2 p 2 1 , (28)
whereas the action becomes

S(z, x, T ) = t 1 0 p 2 1 2 dt + t 1 +t 2 t 1 p 2 2 2 -V 0 dt + T t 1 +t 2 p 2 3 2 dt = 1 2 (x -q -2a)p 1 + a p 2 1 -2V 0 - 2a V 0 p 2 1 -2V 0 . (29)
In this region, no reflection is possible and the probability density |ψ sc | 2 is simply given by Eq. ( 22). The results are presented in figure 5. For any values of p, T or , there is still a very small difference between the exact and semiclassical curves for the ascending part of the wavepacket, whereas the agreement is perfect when the function is decreasing.

The conclusion of this section is that the semiclassical approximation with real trajectories gives very good results and is indeed able to describe some important quantum effects.

Interference on the left side of the barrier appears naturally when the wavepacket hits the barrier and the comparison with the exact solution gets better as gets smaller. However, these interferences cannot be obtained in the barrier region, since there are no reflected trajectories in the classical dynamics. We showed that these interferences can be recovered if a 'ghost' trajectory that reflects at x = a is added and assumed to contribute with the same coefficient of a plane wave of initial momentum p i . With this addition the semiclassical approximation becomes again very accurate inside the barrier. In the next section we shall briefly discuss the possibility of using our results to calculate the tunneling time as defined in [START_REF] Xavier | [END_REF].

IV. SEMICLASSICAL TUNNELING TIMES

The question of how much time a particle spends in the classically forbidden region during the tunneling process has been attracting the attention of physicists for a long time [START_REF] Xavier | [END_REF]16,17,18,19,20,21,22]. The very concept of a 'tunneling time' is, however, debatable [18]. Nevertheless, in a semiclassical formulation where real trajectories play crucial roles in the tunneling process, the temptation to estimate such a time is irresistible.

Since we are considering a wavepacket, and not a classical state localized by a point in the phase space, we can only define a mean value of the tunneling time. Let us fix the initial conditions q, p (such that p < 1) and x f > a. The probability of finding the initial Gaussian state at x f after a time T is given by | < x f |K(T )|z > | 2 . Therefore, the particle can reach x f from (q, p) in several different time intervals T . For each value of the time T there corresponds a single real trajectory whose initial momentum p 1 (T ) > p = √ 2V 0 is given by Eq. ( 28). This trajectory spends a time

τ (T ) = 2a p 2 (T ) = 2a √ p 2
1 (T )-p 2 in the region -a < x < a. Notice that the average energy of the wavepacket is below the barrier but the contributing trajectory always has energy above the barrier. Therefore, for fixed q, p, x f , the probability that the wavepacket crosses the barrier in a time τ (T ) is proportional to

| < x f |K(T )|z > | 2 .
Following ref. [START_REF] Xavier | [END_REF], we can define the mean value of the tunneling time as

τ = N -1 +∞ 0 τ (T ) | < x|K(T )|z > | 2 dT ( 30
)
where

N = +∞ 0 | < x|K(T )|z > | 2 dT (31)
is the normalization factor. It is not equal to 1 because only the part of the wavepacket which goes through the barrier is considered. This is important in our case, since the semiclassical approximation is better for x f > a.

We calculated these integrals numerically, performing a discrete sum over T n = nδT , with n = 1, 2, . . . , N and δT = T max /N. If an observer is placed at a fixed position x f > a, as the time T slips by, he/she sees the wavepacket arriving from the barrier, becoming bigger and bigger, reaching a maximum and then decreasing and disappearing. We ended the sum

at T max such that | < x|K(T )|z > | 2 < 10 -4 ∀ T > T max .
An important remark is that τ is independent of the observer's position x f (except for small fluctuations due to the numerical computation), since Eq.( 30) measures only the time inside the barrier. The three different times we are going to use for comparison are:

• τ barr is the tunneling time computed according to Eq. (30)

• τ f ree is obtained from the same way as τ barr , but in a system without barrier; τ f ree is simply the time for a free wavepacket to go from -a to a.

• τ class = 2a √ p 2 -p 2 is the time required by a classical particle to cross the barrier. Fig. 6 shows the dependence of these functions with respect to p. The curves become very similar as p increases, because the barrier becomes more and more negligible. The wavepacket spreads but stays centered around p, which explains why it behaves like a particle of momentum p. When the influence of the barrier is more important, the wavepacket gets trapped by the barrier and slows down ( τ barr is above τ f ree ), but for p < 1.8, τ f ree and τ class start to increase very fast ( τ class actually diverges at p = 1), whereas τ barr stays finite until p is very close to 0: thanks to the tunnel effect the wavepacket is accelerated by the barrier, which acts like a filter for the wavepacket and cuts off the contributions of its slowest components (see Fig.

V. CONCLUSION

In this paper we used the semiclassical approximation Eq. ( 11), derived in [9], to study the propagation of a wavepacket through a finite square potential barrier. One of the main purposes of this work was to test the validity and accuracy of the approximation, which involves only real trajectories, in the description of tunneling. Surprisingly, we have shown that the approximation is very good to describe the wavepacket after the barrier, even when the average energy of the wavepacket is below the barrier height. The region before the barrier is also well described by the approximation, although discontinuities are always observed because of the sudden disappearance of the reflected trajectory. The continuity of the wavefunction between this region and the region inside the barrier also depends on the calculation of an extra phase θ. Finally, inside the barrier the semiclassical formula is not able to describe interferences. These, however, can be recovered when a ghost trajectory, that reflects on the right side of the barrier, is included and attenuated with the proper coefficient. In all regions the approximation becomes more accurate as becomes smaller.

The semiclassical approximation ( 11) is particularly relevant because the propagated wavepacket is not constrained to remain Gaussian at all times, as in the case of Heller's Thawed Gaussian Approximation [12], and also because it uses only a small number of real trajectories. These are much easier to calculate than complex ones, especially in multidimensional problems. The demonstration of its ability to describe tunneling and interferences is important to establish its generality and also to provide a more intuitive understanding the processes themselves. In particular, using the underlying classical picture, we have computed a tunneling time which shows that the wavepacket can be accelerated or restrained by the barrier depending on the value of the initial central momentum p.

Some interesting perspectives of this semiclassical theory are the study of propagations through smooth potential barriers (which are more realist and more adapted to semiclassical calculations), the study of time dependent barriers and the extension of the method to higher dimensions and to chaotic systems.

  , and let W (x f )e iϕw(x f )+iξ be the wavefunction inside the barrier, where ξ is the corresponding phase correction. The continuity of the wavefunction at x f = -a imposes D(-a)e iϕ d (-a) + R(-a)e iϕr (-a)+iθ = W (-a)e iϕw(-a)+iξ . (16) Eqs. (15) show that R(-a) = D(-a) and ϕ d (-a) = ϕ r (-a). Denoting ϕ = ϕ w (-a)ϕ d (-a), Eq. (16) becomes 1 + e iθ = Ae i(ϕ+ξ) where A = [W (-a)/D(-a)]. This complex relation represents in fact two real equations for the unknown variables θ and ξ. The solutions consistent with the boundary conditions are cos(θ) = A 2 /2 -1 and cos(φ + ξ) = A/2. In the limit where p goes to zero (or the potential height V 0 goes to infinity) we obtain θ = π as expected. Finally, the full wavefunction before the barrier is ψ(z, x f , T ) sc = ψ d + ψ r and the probability density can be written as

  Fig.2shows snapshots of the wavepacket as a function of x f at time T = 50. Consider first the panels (a)-(c) with = 1. The agreement between the exact and the semiclassical curves is qualitatively good for p √ 2V 0 = 1. The interference peaks occurs at about the same positions, but the height of the peaks are not exactly the same. Also the intervals between peaks are a little bigger for the semiclassical curve than for the exact one. On the other hand, when p is increased, the comparison gets worst and the approximation is not really accurate for p = 2. However, we see that the value of |ψ sc | 2 at p = 2 is only a tenth of its value at p = 0.5: the most important part of the wavepacket is in fact inside and after

  of the time T , except for times slightly smaller than T c . In this time interval only the direct trajectory contributes but the exact wavepacket already shows interferences that can not be described by |ψ sc | 2 (fig. 2, T = 10). We now enter the heart of the matter, and consider what's happening inside and after the barrier. B. Inside the barrier: -a x f a

  7.(a)). On the other hand when p increases, the fraction of the trajectories with p < √ 2V 0 = 1 becomes negligible and the barrier simply restrains the propagation of the wavepacket (fig. 7.(b)).

FIG. 1 : 10 FIG. 2 :

 1102 FIG.1: Direct and reflected trajectories from q to x 1 < -a. For -a < x 2 < a or x 3 > a only the direct trajectory exists.

FIG. 3 :FIG. 4 :FIG. 5 :

 345 FIG. 3: (Color online) Exact (blue thin lines) and semiclassical (red thick lines) wavepacket inside the barrier for various values of p and T and = 1.

FIG. 6 :FIG. 7 : 1 -

 671 FIG. 6: (Color online) Tunneling time as a function of p. The red thick line and the green curveshow the semiclassical result according to Eq.(30) for the square barrier and the free particle respectively. The thin black line is the classical time for the square barrier potential.

Acknowledgements

It is a pleasure to thank F. Parisio, A.D. Ribeiro and M. Novaes for many interesting comments and suggestions. This work was partly supported by the Brazilian agencies FAPESP, CNPq and FAEPEX.


L D C Jaubert

P C W Holdsworth

Signature of Magnetic Monopole and Dirac string Dynamics in Spin Ice

Magnetic monopoles have eluded experimental detection since their prediction nearly a century ago by Dirac [1]. Recently it has been shown that classical analogues of these enigmatic particles occur as excitations out of the topological ground state of a model magnetic system, dipolar spin ice [2].

These quasi-particle excitations do not require a modification of Maxwell's equations, but they do interact via Coulombs law and are of magnetic origin. In this paper we present an experimentally measurable signature of monopole dynamics and show that magnetic relaxation measurements in the spin ice material Dy 2 T i 2 O 7 [START_REF] Snyder | Low-temperature spin freezing in the Dy 2 T i 2 O 7 spin ice[END_REF] can be interpreted entirely in terms of the diffusive motion of monopoles in the grand canonical ensemble, constrained by a network of "Dirac strings" filling the quasi-particle vacuum. In a magnetic field the topology of the network prevents charge flow in the steady state, but there is a monopole density gradient near the surface of an open system.

Spin ice systems [START_REF] Harris | Geometrical Frustration in the Ferromagnetic Pyrochlore Ho 2 T i 2 O 7[END_REF][START_REF] Bramwell | Spin Ice State in Frustrated Magnetic Pyrochlore Materials[END_REF][START_REF] Bramwell | Frustrated Spin Systems[END_REF] such as Dy 2 T i 2 O 7 and Ho 2 T i 2 O 7 , can be described by a corner sharing network of tetrahedra forming a pyrochlore lattice of localized magnetic moments, as shown in figure 1a. The pairwise interaction is made up of both exchange and dipolar terms

H = Jm 2 i,j S i • S j + Dm 2 i,j S i • S j |r ij | 3 - 3 (S i • r ij ) (S j • r ij ) |r ij | 5 (1)
where the rare earth ions carry a moment of ten Bohr magnetons, m ≈ 10µ B and where S i is a spin of unit length. The coupling constants are on the 1 K energy scale ;

for example for Dy 2 T i 2 0 7 |J|m 2 ≈ 3.72 K and Dm 2 ≈ 1.41 K [START_REF] Den Hertog | Dipolar Interactions and Origin of Spin Ice in Ising Pyrochlore Magnets[END_REF]. These energy scales are 100 times smaller than the crystal field terms [START_REF] Ehlers | Dynamical crossover in 'hot' spin ice[END_REF] that confine the spins along the axis joining the centres of two adjoining tetrahedra. As a result, on the 1 K energy scale the moments behave as Ising spins along this axis. Remarkably, within this Ising description the long ranged dipolar interactions are almost perfectly screened [START_REF] Den Hertog | Dipolar Interactions and Origin of Spin Ice in Ising Pyrochlore Magnets[END_REF][START_REF] Isakov | Why spin ice obeys the ice rules[END_REF] at low temperature, with the result that the low energy properties are almost identical to those of an effective frustrated nearest neighbour model with antiferromagnetic interactions [START_REF] Anderson | Ordering and antiferromagnetism in ferrites[END_REF] of strength J ef f = (5D -J)m 2 /3 . This is equivalent to Pauling's model for proton disorder in the cubic phase of ice [START_REF] Pauling | The Structure and Entropy of Ice and of Other Crystals with Some Randomness of Atomic Arrangement[END_REF], which has extensive ground state entropy and violates the third law of thermodynamics [START_REF] Bramwell | Frustrated Spin Systems[END_REF]. It successfully reproduces the thermodynamic behaviour of both ice [START_REF] Giauque | The entropy of water and the third law of thermodynamics. The heat capacity of ice from 15 to 273 K[END_REF] and spin ice [START_REF] Ramirez | Zero-point entropy in spin ice[END_REF] and describes the microscopic properties of the latter to a good approximation. The extensive set of spin ice states satisfy the Bernal -Fowler ice rules [START_REF] Bernal | A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions[END_REF] ; a 3d analogue of the 6 vertex model with topological constraint consisting of two spins pointing into and two out of each tetrahedron (2 in -2 out), as shown in figure 1a. Flipping one spin breaks the constraint leaving neighbouring tetrahedra with 3 in -1 out and with 3 out -1 in, which constitute a pair of topological defects. Within the nearest neighbour model, creation of the defect pair costs energy 4 J ef f , while further spin flips can move the defects at zero energy cost. It has recently been shown [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF] that including the full dipolar Hamiltonian of equation (1) leads to an effective Coulombic interaction between the topological defects separated by distance r , µ 0 q i q j /4πr, where µ 0 is the permeability of free space, q i = ±q = ±2m/a, and a is the distance between two vertices of the diamond lattice (see figure 1) ; that is, to a Coulomb gas of magnetic monopoles.

Standard electromagnetic theory does allow for such excitations [START_REF] Jackson | Classical Electrodynamics Ch[END_REF], which correspond to divergences in the magnetic intensity H, or magnetic moment M, rather than in the magnetic induction : ∇ • B = ∇ • (H + M) = 0. On all length scales above the atomic scale, a 3 in -1 out defect appears to be a local sink in the magnetic moment and therefore as a source of field lines in H . It can lower its energy by moving in the direction of an external field and therefore carries a positive magnetic charge [START_REF] Ryzhkin | Magnetic relaxation in rare-earth oxide pyrochlores[END_REF]. What is remarkable about spin ice is that it allows for the deconfinement of these effective magnetic charges so that they occur in the bulk of the material on all scales, rather than just at the surfaces within a coarse grained description [START_REF] Jackson | Classical Electrodynamics Ch[END_REF]. A two dimensional equivalent may exist in artificial spin ice, constituting arrays of nanoscale magnets [START_REF] Wang | Artificial spin ice in a geometrically frustrated lattice of nanoscale ferromagnetic islands[END_REF].

Given the accessibility of these magnetic quasi-particles, the development of an experimental signature is of vital importance and interest. The "Stanford" superconducting coil experiment [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF][START_REF] Cabrera | First results from a superconductive detector for moving magnetic monopoles[END_REF] could in principle detect the passage of a single magnetic quasi-particle, but this seems highly unlikely given that the charges have no mass and therefore have diffusive, rather than Newtonian dynamics. A more promising starting point is therefore to look for a monopole signal from magnetic relaxation of a macroscopic sample [START_REF] Snyder | Low-temperature spin freezing in the Dy 2 T i 2 O 7 spin ice[END_REF][START_REF] Ehlers | Dynamical crossover in 'hot' spin ice[END_REF][START_REF] Matsuhira | Low temperature magnetic properties of frustrated pyrochlore ferromagnets Ho 2 Sn 2 O 7 and Ho 2 T i 2 O 7[END_REF]. The general dynamic behaviour of spin ice is illustrated in figure 2 by the magnetic relaxation time, as a function of temperature for Dy We have tested this idea by directly simulating a Coulomb gas of magnetically charged particles (monopoles), in the grand canonical ensemble, occupying the sites of the diamond lattice. The magnetic charge is taken as q i = ±q. In the grand canonical ensemble the chemical potential is an independent variable, whose value in the corresponding magnetic experiment is unknown. In a first series of simulations we have estimated it numerically by calculating the difference between the Coulomb energy gained by creating a pair of neighbouring magnetic monopoles and that required to produce a pair of topological defects in the dipolar spin ice model, with parameters taken from reference [START_REF] Den Hertog | Dipolar Interactions and Origin of Spin Ice in Ising Pyrochlore Magnets[END_REF], giving a configurationally averaged estimate µ/k B = 8.92 K. In a second series of simulations µ was taken as the value required to reproduce the same defect concentration as in a simulation of dipolar spin ice at temperature T . Here µ varied by 3% only, with the same mean value as in the first series, showing that our procedure is consistent.

The chemical potential used is thus not a free parameter. As the Coulomb interaction is long ranged, we treat a finite system using the Ewald summation method [START_REF] De Leeuw | Simulation of electrostatic systems in periodic boundary conditions. I. Lattice sums and dielectric constants[END_REF][START_REF] Frenkel | Understanding Molecular Simulation : from Algorithms to Applications Ch[END_REF].

The monopoles hop between nearest neighbour sites via the Metropolis Monte Carlo algorithm, giving diffusive dynamics, but with a further local constraint : in the spin model a 3 in -1 out topological defect can move at low energy cost by flipping one of the three in spins, the direction of the out spin being barred by an energy barrier of 8 J ef f .

An isolated monopole can therefore hop to 3 out of 4 of its nearest neighbour sites only, dictated by an oriented network of constrained trajectories similar to the ensemble of classical "Dirac string" [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF] of overturned dipoles [START_REF] Jackson | Classical Electrodynamics Ch[END_REF]. The positively charged monopoles move in one sense along the network while the negative charges move in the opposite direction (see figure 1b). The network is dynamically re-arranged through the evolution of the monopole configuration. The vacuum for monopoles in spin ice thus has an internal structure ; the Dirac strings which, in the absence of monopoles, satisfy the ice rules at each vertex. This structure is manifest in the dynamics and influences the resulting time scales. In fact the characteristic time scale that we compare with experiment comes from the evolution of the network of Dirac strings rather than from the monopoles themselves.

Indeed, the monopole autocorrelation time, as extracted from the monopole densitydensity correlation function [START_REF] Barrat | Basic Concepts for Simple Liquids Ch[END_REF] turns out to be small for this range of temperature.

We locally define the string network by an integer σ = ±1, giving the orientation of the Dirac string along each bond of the diamond lattice and define the autocorrelation function

C(t) = 1 N i σ i (t)σ i (0), ( 2
)
where t is the Metropolis time and N is the number of bonds (up to N ≈ 25000). For the initial conditions we take an ordered network with no monopoles, which we let evolve at temperature T until an equilibrium configuration is attained. This defines Agreement between the experimental and numerical data now looks excellent, showing clearly that the experimental relaxation is due to the creation and proliferation of quasi-particle excitations that resemble classical monopoles in the magnetic intensity H.

t = 0. C(t)
As the temperature increases, towards the end of the plateau region, a small systematic difference occurs. This is because the spin system can access double defects at finite energy cost, while this state corresponds to two like charges superimposed on the same site, which is excluded by the Coulomb interaction. The inset of figure 2 shows results at low temperature illustrating in detail the extent of the improvement in comparison with experiment. Allowing the variation of µ provides a further evolution towards the experimental data, compared with that for fixed µ and the blue circles represent our best numerical results. We now have quantitative agreement between experiment and theory down to low temperature, showing that the Coulomb interactions are responsible for the non-Arrhenius temperature dependence of the relaxation time scales. Differences remain at this level of comparison below 1 K but to go further would require an even more detailed modelling of spin ice [START_REF] Yavors'kii | 2 T i 2 O 7 Spin Ice : A Test Case for Emergent Clusters in a Frustrated Magnet[END_REF] as well as complementary experimental measurements.

Finally we consider the response of monopoles to an external magnetic field, h, placed along one of the [100] directions. Applying such a field to a system for closed circuit geometry (periodic boundaries), one might expect the development of a monopole current in the steady state [START_REF] Ryzhkin | Magnetic relaxation in rare-earth oxide pyrochlores[END_REF]. This is not the case, at least for the nearest neighbour model, where we find that a transient current decays rapidly to zero (see figure 3a).

The passage of a positive charge in the direction of the field re-organises the network of strings, leaving a wake behind it that can be followed either by a negative charge, or by a positive one moving against the field, with the result that the current stops. This is a dynamic rather than static effect and is not related to confinement of monopole pairs by the background magnetization [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF]. Reducing the temperature at finite field, the magnetization saturates around a critical temperature ; a vestige of the Kasteleyn transition [START_REF] Jaubert | A three dimensional Kasteleyn transition : spin ice in a [100] field[END_REF], which is unique to topologically constrained systems. Confinement occurs here, as the Zeeman energy outweighs the entropy gain of free monopoles. The transient currents suggest the development of charge separation in an open system. This is indeed the case despite the fact that monopole numbers are not conserved at open boundaries.

In figure 3b we show the profile of positive charge density across a sample of size L, with open boundaries, for varying fields. There is a clear build up of charge over a band of 4-5 lattice spacings, although including long range interactions may lead to a quantitative change in this value. As the ratio T /h and the monopole density go to zero the band narrows and the system forms a conventional layer of magnetic surface charge as one expects for any magnetically ordered system [START_REF] Jackson | Classical Electrodynamics Ch[END_REF]. In the absence of topological defects the magnetization is conserved from one layer to another, so that a charge density profile manifests itself as a magnetization profile. The data here suggest charge build up in a layer several nanometres thick, making it in principle a measurable effect.

2 T i 2 O 7

 227 [START_REF] Snyder | Low-temperature spin freezing in the Dy 2 T i 2 O 7 spin ice[END_REF], taken from bulk susceptibility measurements. The energy scales discussed above give rise to different regimes : the time scale increases in the thermally activated high temperature regime, entering a quasi-plateau region below 12 K associated with quantum tunnelling processes[START_REF] Ehlers | Dynamical crossover in 'hot' spin ice[END_REF], before experiencing a sharp upturn below 2 K. The spins are Ising like below 12 K and the configuration evolves by quantum tunnelling through the crystal field barrier, while above this temperature higher crystal field levels are populated and the time scale drops dramatically. The quantum tunnelling plateau regime can therefore be well represented by an Ising system with stochastic single spin dynamics and hence should be dominated by the creation and propagation of monopole objects. This is illustrated, in a first approximation, by comparing the data with an Arrhenius law τ = τ 0 exp(2 J ef f /k B T ), as shown by the red curve in figure2. The time scale τ 0 is fixed by fitting to the experimental time at 4 K with J ef f = 1.11 K, the value estimated for Dy 2 T i 2 0 7[START_REF] Den Hertog | Dipolar Interactions and Origin of Spin Ice in Ising Pyrochlore Magnets[END_REF]. 2 J ef f is the energy cost of a single, free topological defect in the nearest neighbour approximation and is half that for a single spin flip. The calculation fits the data over the low temperature part of the quasi-plateau region, where one expects a significant defect concentration without any double defects (4-in or 4-out), and gives surprisingly good qualitative agreement at lower temperature, as the concentration decreases. Although still in the tunnelling regime, the plateau region corresponds to high temperature for the effective Ising system. Good agreement here provides a stringent test and any theory not fitting must be discarded. The above expression clearly does a good job, allowing us to equate τ 0 with the microscopic tunnelling time.This test therefore already provides very strong evidence for the fractionalization of magnetic charge[START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF] and the diffusion of unconfined particles. However, this (or any other) Arrhenius function ultimately fails, underestimating the time scale at very low temperature : while it is possible to fit the data reasonably below 2 K by a single exponential function by varying the barrier height, simultaneous agreement along the plateau and at lower temperature is impossible. The role of the missing Coulomb interaction is therefore clear : although non-confining it must considerably increase the relaxation time scale by modifying the defect concentration and slowing down diffusion through the creation of locally bound pairs.

  decays almost exponentially, with characteristic time, τ , that varies with temperature. To avoid initial transient effects we define τ such that C(τ ) = 0.8. The time is re-set to zero when C(t) decays beyond 0.01 and the process is repeated many times to give the configurationally averaged decay time. In figure2we compare our simulations with the experimental data of reference[START_REF] Snyder | Low-temperature spin freezing in the Dy 2 T i 2 O 7 spin ice[END_REF]. The Metropolis time is again scaled to the experimental time at 4 K and there is again no scale factor on the temperature axis. Data for fixed chemical potential are shown by the pink triangles, while data with µ varying are shown by the blue circles. There is a quantitative evolution of the simulation data compared to the nearest neighbour spin ice model.

Figure 1 :

 1 Figure 1: Spin Ice structure and emergence of monopoles. a) The magnetic ions (Ho 3+ or Dy 3+ ) lie on the sites of the pyrochlore lattice and are constrained to the bonds of the duel diamond lattice (dashed lines). Local topological excitations 3 in -1 out or 3 out -1 in correspond to magnetic monopoles with positive (blue sphere) or negative (red sphere) charges respectively. b) The diamond lattice provides the skeleton for the network of Dirac strings with the position of the monopole restricted to the vertices. The orientation of the Dirac strings shows the direction of the local field lines in H.

Figure 2 :

 2 Figure 2: Relaxation time scales τ in Dy 2 Ti 2 O 7 : experiment and simulation. The experimental data (×) are from Snyder & al. [3]. The Arrhenius law (red line) represents the free diffusion of topological defects in the nearest neighbour model. The relaxation time scale of the Dirac string network driven by Metropolis dynamics of magnetic monopoles has been obtained for fixed chemical potential (pink ) and with µ varying slowly to match the defect concentration in dipolar spin ice (blue •). The temperature scale is fixed without any free parameters. Inset : Same data shown in the low temperature region.

Figure 3 :

 3 Figure 3: Monopole density profile. a) A magnetic field h is applied at t = 0 along the z-axis ([100] direction). We display the transient flux of positive charges, Φ, passing through a plane perpendicular to the field, as a function of Metropolis time, t. The simulations are obtained using the nearest neighbour spin ice model with periodic boundary conditions ( ) and open boundaries, with current measured either at the surface (blue ) or in the bulk (red •). b) Density of positive defects in the horizontal planes for T = 1 K, as a function of z and h (in units of k B T /mµ 0 ).

Acknowledgements We thank S.T. Bramwell, C. Castelnovo, J.T. Chalker, M. Clusel, M.J.P. Gingras and R. Melko for useful discussions and P. Schiffer for providing the data from reference [3]. We acknowledge financial support from the European Science Foundation PESC/RNP/HFM.

Author Contributions

The authors contributed equally to the manuscript.

Author Information Correspondence should be addressed to L.J. (e-mail : ludovic.jaubert@ens-lyon.org).