MULTISTABILITY AND FRACTAL PROPERTIES OF HAMILTONIAN LATTICE MODELS

Karl E. Kürten

Institute for Experimental Physics, University of Vienna 5, Boltzmanngasse, A-1090 Vienna, Austria

Christian Krattenthaler

Institut Girard Desargues, Université Claude Bernard Lyon-I, 21, avenue Claude Bernard, F-69622 Villeurbanne Cedex, France

1. INTRODUCTION

Numerous condensed matter systems are effectively discrete by nature because the relevant length scales are of the order of the interparticle distance. Such systems are described by a Hamiltonian that is discrete in space, while their time evolution is considered as continuous. Their remarkable behavior, exemplified in models of charge-density waves, magnetic spirals, disordered crystals, magnetic multilayers or adsorbed monolayers, stems from a competition between two or more forces that leads to multistable spatially modulated structures [1–5]. One other prominent example is modeling proton dynamics in hydrogen bonded chains [6]. Spatial disorder is associated with a highly complex energy landscape in configuration space, where the number of locally stable configurations typically increases exponentially with the size of the system. The simplest model system can be envisioned as a onedimensional chain of N particles with nearest-neighbour interactions, each particle, which is considered as a classical spin, also subject to an external multi-well potential field. Specifying the displacement of particle i at time t from its lattice site by the variable x_{i}, and assuming periodic boundary conditions, the total energy of the system is given by the N-particle Hamiltonian

H(x_{1}, x_{2}, …, x_{N}) = ∑_{i=1}^{N} (\frac{1}{2}ẋ_{i}^{2} + E_{I}(x_{i−1}, x_{i}, x_{i+1}) + U(x_{i})) .  (1)

The first term is the kinetic energy, the second term defines the inter-site energy, specified by nearest-neighbour interactions, while the last term is the potential energy due to an external on-site potential.

A standard harmonic inter-site potential is given by

E_{I}(x_{i−1}, x_{i}, x_{i+1}) = \frac{1}{4}(x_{i} − x_{i−1})^{2} + \frac{1}{4}(x_{i} − x_{i+1})^{2},

while the on-site potential is specified by the Φ^{4} double-well potential

U(x_{i}) = \frac{α}{4}(1 − x_{i}^{2})^{2},  (2)

depicted in Fig. 1. It is often used for the description of displacive and order-disorder transitions. Here, the control parameter α specifies the height of the barrier U(0) = \frac{α}{4} which separates the two symmetrical sinks of the two potential wells.

Figure 1. Symmetric Φ^{4} double-well potential for α = 8.

Based on the Φ^{4} potential, the total energy Eq. (4) takes the form

E = ∑_{i=1}^{N} (\frac{1}{2}ẋ_{i}^{2} + \frac{1}{4}(x_{i} − x_{i−1})^{2} + \frac{1}{4}(x_{i} − x_{i+1})^{2} + \frac{α}{4}(1 − x_{i}^{2})^{2}) .  (3)

The ground state with zero energy is obtained, when all particles are placed in one of the two sinks of the double-well potential. According to the variational principle, stationary configurations obey the set of N non-linear coupled equations

−x_{i−1} + 2x_{i} − x_{i+1} − αx_{i}(1 − x_{i}^{2}) = 0,  i = 1, …, N.  (4)

However, solutions x^{∗} = (x_{1}^{∗}, x_{2}^{∗}, …, x_{N}^{∗}) specify a local minimum only if the eigenvalues of the triangular Hessian ∇^{2}H(x)_{|x=x^{∗}} = (\frac{∂^{2}H(x)}{∂x_{i}∂x_{j}})_{|x=x^{∗}} are all positive. Since the diagonal entries are 2 + U′′(x_{i}^{∗}) and the two off-diagonal entries take the value −1, Gershgorin’s disk theorem provides a sufficient stability condition U′′(x_{i}^{∗}) > 0, where all particles occupy the convex part of the potential well. For the double well potential this condition is satisfied for x_{i}^{2} > \frac{1}{3}.

There are two macroscopic quantities of interest. First, the total energy, which takes the simple form, reminiscent of the virial theorem in the continuous counterpart of the nonlinear Klein-Gordon equation,

E = \frac{α}{4} ∑_{i=1}^{N}(1 − x_{i}^{4}),  (5)

which is obtained after insertion of the equilibrium condition Eq. (4). Second, the average displacement

M = \frac{1}{N} ∑_{i=1}^{N} x_{i},  (6)

reminiscent of the total magnetization of a classical spin chain of interacting magnetic particles, a macroscopic variable, which is often experimentally accessible. In this case the intersite energy E_{I} specifies the exchange interaction between two neighbouring particles, while the on-site potential plays the role of the uniaxial anisotropy energy with strength parameter α [7].

2. FEW PARTICLE CASE

Let us first consider the situation for N = 2 particles, where the energy takes the symmetric form

E(x_{1}, x_{2}) = (x_{1} − x_{2})^{2} + \frac{α}{4} ((1 − x_{1}^{2})^{2} + (1 − x_{2}^{2})^{2}) .

The set of coupled equations

−x_{2} + 2x_{1} − x_{2} − αx_{1}(1 − x_{1}^{2}) = 0,
−x_{1} + 2x_{2} − x_{1} − αx_{2}(1 − x_{2}^{2}) = 0.

has 9 real valued solutions representing stationary configurations. Fig. 2 shows the energy landscape E(x_{1}, x_{2}) for α = 8 characterized by one relative maximum (0,0) and four local minima ±(1, 1), and ± (\sqrt{1 − 4/α}, −\sqrt{1 − 4/α}). The first two minima represent the two degenerate ground states with energy E = 0, while the latter represent two energetically higher lying antiferromagnetic states with energy per particle E/2 = 2 (1 − \frac{2}{α}) . Note that these four local minima are separated by four saddle points with energy per particle E/2 = (2 + α)^{2}/(8α).

The saddle points act as energy barriers, whose height increases linearly with increasing control parameter α. It is important to note that the height of the barriers is of the order of α, while the energy values of the local minima are all of the order of one. Accordingly, the energy difference between these states is much smaller than the energy required to overcome the barrier in order to switch to another state of equilibrium.

Figure 2. Energy landscape E(x_{1}, x_{2}) for α = 8.

The same arguments hold for larger values of N, where we will show in Section 3 that the energy landscape exhibits an extremely complicated multi-valley structure with an exponentially increasing number of local minima and saddle points that also allows the appearance of structural disorder. For N = 3 we have a ”soliton” solution

x^{∗} = {x_{1}, x_{2}, x_{1}},

with

x_{1} = [\frac{2}{3} − \frac{1}{α} + \frac{2}{3} (1 − \frac{3}{α})^{\frac{1}{2}} cos \frac{ϕ + π}{3} ]^{1/2} and x_{2} = x_{1}(1 − α + cx_{1}^{2}),

and with

cos ϕ = \frac{1 − \frac{9}{2α} + \frac{27}{α^{3}}}{(1 − \frac{3}{α})^{\frac{3}{2}}} .

This solution represents a local minimum only for α > 3 + \sqrt{9 + 6\sqrt{3}} ≈ 7.40367 with the expansion x_{1} = 1 − \frac{2}{α} − ⋯ and x_{2} = −1 + \frac{1}{α} − ⋯, respectively.

3. EXACT ENUMERATION OF ALL LOCALLY STABLE STRUCTURES

Let us now further concentrate on large values of α, where the on-site potential is sufficiently strong. One might ask: which kind of landscape- or domain structures can evolve due to the competition between the exchange- and the on-site energy? Depending on the number of particles N, how many distinct non-equivalent structures can there exist? When the control parameter α is varied, the energy balance is changed such that a rearrangement of the energy landscape can take place. Stationary configurations of the Hamiltonian (1) are described by the N non-linear coupled equations

−x_{i−1} + 2x_{i} − x_{i+1} − αx_{i}(1 − x_{i}^{2}) = 0,  i = 1, …, N.  (7)

In order to determine all possible locally stable spatial structures we first consider the system for large values of α. In the large α-limit, Eq. (7) reduces to

x_{i}(1 − x_{i}^{2}) = 0,  i = 1, …, N.  (8)

Ignoring the solutions x_{i} = 0, which correspond to saddle point configurations, this coupled set of equations has exactly 2^{N} distinct solutions consisting of N-dimensional binary strings

x^{∗} = (x_{1}^{∗}, x_{2}^{∗}, …, x_{N}^{∗})  (9)

with components x_{i}^{∗} ∈ {−1, 1}. Provided that α is sufficiently large, the implicit function theorem guarantees the existence of a solution of Eq. (7) in a neighbourhood of any of these binary strings, all of these solutions being distinct. (Experimentally, the critical α depends weakly on the particle number N, however, numerical calculations suggest that in fact all 2^{N} locally stable solutions of Eq. (7) exist in a neighbourhood of any of these binary strings provided that α ≈> 8.) If these are expressed in terms of a power series expansion about the asymptotic strings x^{∗}, then the first order correction term

Δx_{i}^{(1)} = \frac{x_{i−1}^{∗} − 2x_{i}^{∗} + x_{i+1}^{∗}}{2α}  (10)

is proportional to 1/α. The particular form of the first- and the higher order correction terms then implies that |x_{i}| ≤ 1 rigorously holds. It should be noted however that, depending on a specific solution, the stability range can vary widely.

Since the particles reside on a circle, we will consider configurations that only differ by a rotation (x_{i} → x_{s+i})_{1≤i≤N} (for some fixed s with 1 ≤ s ≤ N) or by a reflection (x_{i} → x_{s−i})_{1≤i≤N} (for some fixed s with 1 ≤ s ≤ N) as equivalent. Together, these rotations and reflections form the dihedral group D_{N}. Since for large α, the stable solutions of (7) can be identified with the 2^{N} strings of 1s and −1s of length N (for example, the “soliton” solution for N = 3 presented in Section 2 would be identified with the string (−1, 1, −1)), the number of inequivalent solutions will be equal to the number of orbits of these (circular) strings under the action of the dihedral group D_{N}. In more colloquial terms, we want to count the number F_{D}(N) of different necklaces one can make out of N beads, each of them either coloured white or black. It follows from Redfield-Pólya theory (see for example [8]), that the number of these orbits (or necklaces), and thus the number of inequivalent stable solutions to (7), is

F_{D}(N) = \frac{1}{2N} ∑_{d|N} φ(d) 2^{N/d} + \{ \begin{matrix}\frac{3}{4}2^{N/2} & for N even, \\ 2^{(N−1)/2} & for N odd, \end{matrix}  (11)

where φ(x) is the Euler totient function. For N = 1, 2, 3, … this is the sequence {2, 3, 4, 6, 8, 13, 18, 30, 46, 78, 126, 224, …}. Since φ(1) = 1, the number F_{D}(N) is asymptotically equal to 2^{N−1}/N as N tends to ∞.

In order to enumerate the orbits of strings that are invariant under reflections and those which are not separately, we also count the number F_{C}(N) of different necklaces made out of N beads, each of them coloured white or black, which are embedded in the plane, that is, the orbits of strings of 1s and −1s of length N under the action of the cyclic group C_{N} of rotations (x_{i} → x_{s+i})_{1≤i≤N}. Redfield–Pólya theory again shows that

F_{C}(N) = \frac{1}{N} ∑_{d|n} φ(d)2^{N/d}.

The first few terms of this sequence are {2, 3, 4, 6, 8, 14, 20, 36, 60, 108, 188, 352, …}.

As is easy to see, the number of orbits under the action of the dihedral group D_{N} which are not invariant under reflection is then F_{C}(N)−F_{D}(N), while the number of the orbits under the action of the dihedral group D_{N} which are is 2F_{D}(N) − F_{C}(N).

The spectrum of the total magnetizations M as well as the energy spectrum, consisting of 8267 distinct values, both reminiscent of Cantor-set structures, are depicted in Fig. 3 and Fig. 4, respectively.

Figure 3. Spectrum of the total magnetizations Eq. (6) for N = 15 and α = 8.

Figure 4. Energy spectrum Eq. (5) for N = 15 and α = 8.

Note that a specific domain structure corresponds to the system being trapped in one of these local energy minima. If thermal fluctuations are neglected and the energy barriers separating this minimum from neighbouring ones are large enough, the system will indefinitely remain in such a metastable state. However, a slight change of the strength of the on-site potential U(x) can easily destabilize a specific domain structure such that it is sufficient that a local minimum of the energy landscape is transformed into a saddle point and the system can evolve towards some other metastable configuration.

4. EXACT ENUMERATION OF ALL DISPLACEMENTS

In order to determine the distribution or the ”spectrum” of all possible displacements |x_{i}|, i = 1, 2, …, N, as a function of the particle number N, we exploit the fact that two-dimensional maps allow us to represent any stationary configuration of a Hamiltonian of the form (1) by a trajectory of a dynamical system. Therefore, the mathematical situation is identical to that of temporal evolution, although the problem has been formulated in terms of spatial arrangements. Introducing the auxiliary variable y_{n+1} = x_{n}, Eq. (7) takes the form of a two-dimensional iterative map in discrete time,

( \begin{matrix}x_{n+1} \\ y_{n+1} \end{matrix} ) = F(x_{n}, y_{n}) = ( \begin{matrix}−y_{n} + 2x_{n} − αx_{n}(1 − x_{n}^{2}) \\ x_{n} \end{matrix} )  (12)

Due to the periodic boundary conditions, we have to determine all F_{D}(N) inequivalent periodic orbits of the map equation. These can be easily calculated, once the fixed points

( \begin{matrix}x^{∗} \\ y^{∗} \end{matrix} ) = F^{(N)}(x^{∗}, y^{∗})  (13)

of the N-th iterate of the map are known. Since periodic orbits of the map equation can be shown to be unstable with respect to iterating the map [9], we solve the nonlinear set Eq. (13) with a standard Newton-Raphson procedure, where the asymptotic strings Eq. (9) serve as initial conditions.

We first calculate all possible displacements as the length N of the chain is varied, where we first concentrate on N odd. Fig. 5 shows the displacement spectrum for odd N = 3, 5, 7, 9. Experimentally, with increasing number of particles, the spectrum evolves by consecutive tripling and quadrupling of spectral substructures. It obviously grows in a quasi self-similar manner, where structures replicate themselves on successively smaller scales with respect to their statistical properties. Triplets split into a quadruplet surrounded by two triplets, while quadruplets split into four quadruplets. Hence, it seems that the spectrum grows according to

S(N) = 2S(N − 2) + 4^{\frac{N−3}{2}} .  (14)

Figure 5. Spectrum of the displacements for N = 3, N = 5, N = 7, and N = 9 (from above) for α = 8 .

For N = 3, 5, 7, and N = 9, we have 3, 10, 36, and 136 distinct values, respectively. An estimate for the fractal dimension can be obtained by using a box counting algorithm, and it yields the value d_{f} = 0.58.

A rigorous verification of the recurrence (14) is equivalent to rigorously determine the numbers S(N) of distinct displacements. We shall be able to do this under the hypothesis, which is confirmed by our experiments, that displacements in one orbit are distinct from the displacement in another orbit, and that the number of distinct displacements in an orbit is equal to the minimal period s under rotation of the corresponding binary string x^{∗} = (x_{1}^{∗}, x_{2}^{∗}, …, x_{N}^{∗}) if the orbit is not invariant under reflection, and equal to (s + 1)/2 otherwise. Thus, one way to determine S(N) is to calculate the sum \frac{1}{2} ∑_{x^{∗}} w(x^{∗}), where the sum is over all orbits x^{∗}, and w(x^{∗}) = s, respectively w(x^{∗}) = (s + 1)/2, depending on the two cases described above. (The factor \frac{1}{2} results from the invariance x_{i} → −x_{i} in Eq. (7)). There is a similar weighted sum for even N.

It is possible to do this calculation directly. There is, however, an elegant way around this. Consider the set of all binary strings constructed in the following manner: For any orbit, let x^{∗} = (x_{1}^{∗}, x_{2}^{∗}, …, x_{N}^{∗}) be the corresponding binary string, with weight w(x^{∗}) = w, say. List the strings

(x_{1}^{∗}, x_{2}^{∗}, …, x_{N−1}^{∗}, x_{N}^{∗})
(x_{2}^{∗}, x_{3}^{∗}, …, x_{N}^{∗}, x_{1}^{∗})
(x_{3}^{∗}, x_{4}^{∗}, …, x_{1}^{∗}, x_{2}^{∗})
. . . . . . . . . . . . . . . . . . . . . . . .
(x_{w}^{∗}, x_{w+1}^{∗}, …, x_{w−2}^{∗}, x_{w−1}^{∗}).

Finally form the union of all these strings over all orbits. The upshot is that this set of strings is a set of representatives of orbits of strings of 1s and −1s of length N under the action of the group generated by the two operations of reversal of strings r : (x_{i}^{∗}) → (x_{N−i+1}^{∗}), and the operation of simultaneous exchange of 1s and −1s e : (x_{i}^{∗}) → (−x_{i}^{∗}). This group consists only of 4 elements, the identical operation (x_{i}^{∗}) → (x_{i}^{∗}), the reversal r, the exchange e, and the composition re : (x_{i}^{∗}) → (−x_{N−i+1}^{∗}) of reversal and exchange. Thus, with the help of Redfield-Pólya theory, we arrive at the explicit expression

S(N) = 2^{N−2} + \{ \begin{matrix}2^{(N−2)/2} & for N even, \\ 2^{(N−3)/2} & for N odd, \end{matrix}  (15)

which had been empirically found earlier in [7] and confirms indeed the empirically observed recurrence Eq. (14).

4. DISCUSSION

We have shown that the energy surface of our model consists of exponentially many locally stable minima separated by large barriers specified by unstable saddle points. With the aid of Redfield-Pólya theory we give the exact number of inequivalent locally stable configurations as well as the exact number of distinct values of the displacement spectrum. Both quantities increase exponentially with increasing particle number N. As N tends to ∞ they are asymptotically 2^{N−1}/N and 2^{N−2} , respectively. The spectrum, which grows exponentially in a quasi self-similar manner is shown to be a fractal with fractal dimension d_{f} = 1.55.

The model with slight modifications can also describe the formation of magnetic structures in physical novel systems such as antiferromagneticly coupled magnetic multilayers, where modulated structures arise due to the competition between the exchange energy, the anisotropic on-site energy and an additional external magnetic field [7]. Here the fractal structure is intimately connected with the experimentally observed Barkhausen jumps which exhibit fractal properties associated with random Cantor dust of a fractal dimension d_{f} depending on the strength of the external magnetic field [10].

ACKNOWLEDGMENTS

The first author thanks S. Charpentier, M. Olivier and M. Martin for helpful computer assistance and F.V. Kusmartsev for various useful discussions.

REFERENCES

[1] W. Selke, in Phase Transitions and Critical Phenomena, Vol. 15,1, (Academic Press Limited, 1992).

[2] P. Bak, Rep. Prog. Phys. 45,587 (1982).

[3] S. Aubry, Physica D 7,240 (1983).

[4] K. E. Kürten, in Condensed Matter Theories, Vol. 13, edited by B. Malik and J. da Providencia , (Nova Science Publisher, New York, 1998), p.473.

[5] F. V. Kusmartsev and K. E. Kürten, in Lecture Notes in Physics, Vol. 284 (1997), edited by J. W. Clark and M. L. Ristig, Springer-Verlag, New York Heidelberg Berlin.

[6] R. Grauer, K.H. Spatschek and A.V. Zolotaryuk, Phys.Rev. E 47,1,236 (1993).

[7] K. E. Kürten, in Condensed Matter Theories, Vol. 15, edited by J. da Providencia , (Nova Science Publishers, New York, 2003).

[8] E. N. Gilbert and J. Riordan, Illinois J.Math. 5, 657–665 (1961).

[9] K. E. Kürten, in Condensed Matter Theories, Vol. 14, edited by D. Ernst, I. Perakis, and S. Umar (Nova Science Publishers, New York, 1999).

[10] G. Durin, A. Magni, and G. Bertotti, J. Magn. Magn. Mater. 140–144, 1835– 1836 (1995).