Shaken Dynamics on the 3-D Cubic Lattice

On the space of $\pm 1$ spin configurations on the 3$d$-square lattice, we consider the \emph{shaken dynamics}, a parallel Markovian dynamics that can be interpreted in terms of Probabilistic Cellular Automata. The transition probabilities are defined in terms of pair ferromagnetic Ising-type Hamiltonians with nearest neighbor interaction $J$, depending on an additional parameter $q$, measuring the tendency of the system to remain locally in the same state. Odd times and even times have different transition probabilities. We compute the stationary measure of the shaken dynamics and we investigate its relation with the Gibbs measure for the 3$d$ Ising model. It turns out that the two parameters $J$ and $q$ tune the geometry of the underlying lattice. We conjecture the existence of unique line of critical points in $J-q$ plane. By a judicious use of perturbative methods we delimit the region where such curve must lie and we perform numerical simulation to determine it. Our method allows us to find in a unified way the critical values of $J$ for Ising model with first neighbors interaction, defined on a whole class of lattices, intermediate between the two-dimensional hexagonal and the three-dimensional cubic one, such as, for example, the tetrahedral lattice. Finally we estimate the critical exponents of the magnetic susceptibility and show that our model captures a dimensional transition in the geometry of the system at $q = 0$.


Introduction
Probabilistic Cellular Automata (PCA) are discrete-time Markov chains on a product space S Λ (configuration space) whose transition probability is a product measure, i.e. given two generic configurations τ = (τ 1 , . . ., τ N ) and σ = (σ 1 , . . ., σ N ): so that for each time n, the components of the configuration are independently updated.From a computational point of view, the evolution of a Markov chain of this type is well suited to be simulated on parallel processors.
Recently, a class of PCA has been introduced in order to study nearest neighbors spin systems on lattices and, more generally, spin systems on arbitrary graphs G = (V, E), where the interaction Hamiltonian is given by with both J xy and λ x in R, and σ ∈ {−1, +1} V a configuration on G.In this context, the transitions probability from a configuration σ to a configuration τ is defined in terms of a pair Hamiltonian H(σ, τ ) and these transitions are such that, at each time step, the value The layers structure of the induced bipartite graph (red edges) associated to the original graph (gray edges).On the induced graph the vertical edges, corresponding to the self interaction, have weight q whereas the slanted edges have weight J. of all spins is simultaneously updated (see: [3,7,4,10,11]).In this framework, a new parallel dynamics for Ising-like models on general finite graph called shaken dynamics has been introduced in [1] and has been extensively investigated in the case of the two dimensional square lattice.The distinctive feature of the shaken dynamics is the fact that transitions between states are obtained through a combination of two half steps.In each of these half steps the value of the spin at site x is updated according to a probability distribution depending, through a self interaction parameter q > 0, on the value of the spin at site x itself and the values of the spins sitting at a suitable subset of the sites adjacent to x in such a way that all neighbors of x are considered exactly once in the whole step.It is worth noting that a shaken dynamics on a given graph structure can be naturally associated to a dynamics on an induced bipartite graph where the spins in each partition are alternatively updated.The vertex set of this bipartite graph consists of two copies of the original vertex set so that this induced graph can be thought as to have two layers (see Fig 1).
The sub-configuration on one of the layers is the "current" configuration of the shaken dynamics whereas the sub-configuration on the other layer is the "intermediate" configuration reached through the first half step.In this view, the shaken dynamics can be seen as the evolution taking place on one of the layers of the associated alternate dynamics.The geometry of the induced bipartite graph where the alternate dynamics lives varies continuously with q.For instance, in the case of the shaken dynamics on Z 2 with J xy = J for all pairs of nearest neighbors {x, y}, the bipartite graph where the associated alternate dynamics evolves is a non homogeneous hexagonal lattice.For q = J this graph becomes the homogeneous hexagonal lattice; in the limit q → ∞ the hexagonal lattice "collapses" onto the square lattice whereas for q = 0 the hexagonal lattice becomes a collection of independent one dimensional lattices.
In this work we study the shaken dynamics on the 3d cubic lattice with J x,y = J for all {x, y} that are nearest neighbors.We determine the stationary measure of the shaken dynamics and show that if the self interaction q is sufficiently large, then this equilibrium measure tends to the Gibbs measure in total variation distance.We argue that the associated alternate dynamics takes place on a suitable tetrahedral lattice that becomes homogeneous if J = q, becomes the cubic lattice in the limit q → ∞ and reduces to a collection of independent 2d hexagonal lattices if q = 0.
It is reasonable to assume that, in the case of null external magnetic field, there is a critical curve J c (q) in the J − q plane which separates the ordered phase from the disordered one.To gain some information on J c (q), we determine two curves in the J − q plane such that above the "upper curve" the system is in an ordered phase (low temperature regime), whereas below the "lower curve" the system is in a disordered phase (high temperature regime).The critical curve must lie in the region delimited by these two curves.Further we provide a numerical estimate for J c (q).We see that our estimates for J c (0), J c (J) and J c (∞) are, respectively, in good agreement with the critical temperature of the Ising model on the hexagonal lattice and the numerical estimates available for the critical temperature of the Ising model on the tetrahedral and the cubic lattice.This suggests that the numerically determined critical curve should be not too far apart from the "real" one.Moreover we study numerically the critical exponents for the magnetic susceptibility as a function of the self interaction q and provide some evidence that the system retains its three dimensional structure as long as q > 0 whereas it becomes two dimensional when q = 0.In other words our model is able to capture the dimensional transition at q = 0.
In the next section we define the lattice spin model and the shaken dynamics, we describe the first properties of the model, and highlight its relation with the alternate dynamics on the tetrahedral lattice.Further we state our main results.Section 3 is devoted to the proofs.Finally, in section 4 we present our numerical findings concerning the critical curve and discuss the behavior of the critical exponents.

Model description and main results
2.1.Definitions and first properties.Let Λ be a L × L × L square box in Z 3 and let X be the set of all possible spin configurations, i.e.X = {σ : σ = {−1, 1} |Λ| }.
We call B Λ be the set of all pairs of nearest neighbors when periodic boundary conditions are imposed on Λ.Let σ, τ ∈ X be two spin configurations and define a pair Hamiltonian in the following manner: where q J > 0 represents the ferromagnetic interaction constant, q q > 0 represents the inertial (or self-interaction) term, q λ > 0 represents the intensity of external magnetic field, q x u is the site above x at lattice distance 1 from x itself, q x r is the site on the right of x q x f is the site in front of x q x d is the site below (down) x q x l is the site on the left of x q x b is the site behind x q σ x (resp.τ x ) is the spin at site x in configuration σ (resp.τ ) q σ d x is the spin at site x d in configuration σ (σ l x , σ b x , τ u x , . . .are defined likewise) See Fig. 2.
It is straightforward to check that Hamiltonian (2.1) is linked tightly to the standard Ising one on the standard cubic lattice.In particular the following proposition holds.
Proposition 2.1. ) the standard Ising Hamiltonian with external magnetic field twice that of Hamiltonian (2.1).
In the same spirit of [1] we want to define a shaken dynamics on X .To this end, consider a Markov Chain that updates the spin configuration with transitions probability P dlb at odd times and P urf at even times where: with Then, the shaken dynamics is defined through the composition of an "odd" and an "even" step.More precisely: Though, strictly speaking, the shaken dynamics (2.6) is not a PCA in the sense of (1.1), it is the composition of two steps each having a factorized transition probability.Indeed: where ) are the local fields felt at site x at, respectively, the odd and the even "half steps".
Observe that the Hamiltonian (2.1) is not symmetric: H λ (σ, τ ) = H λ (τ, σ).This implies that a dynamics evolving solely according to P dlb or P urf is not reversible.However, when the shaken dynamics (2.6) is considered, then the following result holds: Blue and red dots represent respectively the sites in the two three-dimensional lattices V 1 and V 2 .The black segments correspond to the interaction governed by the parameter J, the green segments correspond to the self-interaction, governed by the parameter q.
Proof.The detailed balance condition is readily established, indeed: 2.2.Alternate dynamics.Let V 1 and V 2 be two copies of V the vertex set of Λ and let Λ D be a finite graph with vertex set given by V 1 ∪ V 2 .For each site x in Λ, x 1 , x 2 are the two copies of x in, respectively, V 1 and V 2 and are called corresponding sites.Consider a spin configuration σ 1 on V 1 and a spin configuration σ 2 on V 2 .Then, the Hamiltonian H(σ 1 , σ 2 ) defines the set of edges on Λ D .In particular each pair x 1 , x 2 of corresponding sites is connected by an edge with weight q.Moreover, each Note that the graph Λ D is bipartite by construction and each edge has one endpoint in V 1 and one in V 2 .A graphical representation of Λ D is given in Fig. 3.
The parameter q determines the geometry of the lattice Λ D .Indeed, thinking to the edge weight as to be proportional to the inverse of the geometrical distance between the vertices we have: q the limit q → 0 correspond to erasing the q−edges obtaining, from the lattice Λ D , "independent" copies of the two-dimensional honeycomb lattice; q when J = q, the q−edges and the J−edges become of the same length, so the lattice Λ D becomes a tetrahedral lattice, that we can imagine like a diamond structure; q the limit q → ∞ correspond to identify the two vertices linked by the q−edge, in this case the lattice Λ D degenerates into a simple cubic lattice.
Consider a dynamics (in the remainder referred to as alternate dynamics) that, alternatively, at each step updates all spins in one of the two layers.Then the shaken dynamics can be seen as the projection of the alternate dynamics onto one of the layers.
To make this statement precise, let σ = (σ 1 , σ 2 ), τ = (τ 1 , τ 2 ) ∈ X D .Then, the alternate dynamics on X D is defined by the transition probabilities and the transition probabilities of the shaken dynamics can be written in the form Let Z Λ D = σ,τ e −H(σ,τ ) .As far as the stationary measure of the alternate dynamics is concerned we have the following: Proposition 2.3.The alternate dynamics defined on X D with transition probability P alt ( σ, τ ) has the following stationary measure: (2.12) Moreover, in general, this dynamics is irreversible. Proof.
Remark 2.4.The stationary measure of the shaken dynamics is the marginal of the stationary of the alternate dynamics, that is (2.15) 2.3.Results.In the case of null external magnetic field, we identify two regions of analiticity of the partition function in the thermodynamic limit.These two regions correspond, respectively, to a low temperature and a high temperature regime for the system.The estimation of the critical curve expected to lie outside of these two regions is provided below (see Section 4).Moreover, for large q we establish a link between the equilibrium measure of the shaken dynamics and the Gibbs measure for the Ising model defined on the standard cubic lattice.This result, presented in Theorem 2.5, gives a quantitative support to the statement that, if q is sufficiently large, the lattice where the alternate dynamics takes place, tends to the simple cubic one.
The stationary measure of the Markov chain defined above, that is π Λ (σ), is linked to the Gibbs measure π G Λ (σ) by the following result: , then there exist a J and J, with J > J, such that for J > J and J < J we have:

16)
Remark: In this paper we are considering periodic boundary condition on Λ. However Theorem 2.5 holds also if "plus" or "minus" boundary conditions are imposed at the exterior of the box.Indeed, if we consider a set B of fixed spins in Λ, such that |B|/|Λ| → 0 as |Λ| → ∞, the proof of the theorem requires only minor modifications (see below).Using this approach, any kind of external boundary conditions can easily be mimicked by fixing the spins on three mutually orthogonal planes.
To identify the low temperature regime, we look at the magnetization of the origin and determine parameters such that this magnetization is positively correlated with the magnetization at the boundary even in the thermodynamic limit.
Let π + Λ (σ) (resp.π − Λ (σ)) be the equilibrium measure of the shaken dynamics when + (resp.−) boundary conditions are taken into account and let σ 0 + Λ (resp.σ 0 − Λ ) be the expected value of σ 0 with respect to this probability measure, that is σ 0 Then, for λ = 0: Theorem 2.6.In the thermodynamic limit the mean magnetization of the origin depends on the boundary conditions, that is if J and q are sufficiently large.The explicit description of the low temperature region is given in (3.63).
Conversely, when the system is at high temperature, it is possible to bound the analyticity region of the free-energy density by considering a suitable (high temperature) expansion.More precisely, for λ = 0: Theorem 2.7.The region in the q − J plane where the free-energy density is analytic in the thermodynamic limit contains a well defined region specfied by (3.101) below.

Proofs of the main results
3.1.Proof of theorem 2.5.We start proving the first part of the theorem: if lim |Λ|→∞ δ 2 |Λ| = 0, then there exists a J such that, for J > J we have: To this end, we need some preliminaries lemmas. where: Proof.It follows using the same steps of [1,Equation 25] In order to compare the stationary measure π Λ with the Gibbs measure π G Λ , it is convenient to rewrite the previous expression for Recalling the definition of the Gibbs measure: then π Λ can be written as: where: With this notation, the following lemma provides a bound for the difference between the two measures on X . with: (3.10) Hence, to prove the first part of theorem (2.5), we need to show that: Writing: then we need to show that ∃ J such that, for J > J, we have: Indeed, if (3.13) holds, the Taylor expansion of the exponential of the r.h.s. of (3.12) is well defined.(3.14) follows by noting that the first order terms of the Taylor expansions are zero.
To prove both claims, it is convenient to partition the sites of the finite cubic lattice Λ according to the value of their spin and the sum of the spins at the downwards, leftwards and backwards neighboring sites.To this purpose we define the sets Then, arguing as in [1], it is possible to rewrite f (σ) in this way: with: To bound ξ, we rewrite H 2λ (σ) in terms of 3d-Peierls contours defined in the following way: for each pair nearest neighboring sites x and y such that σ x σ y = −1, we build a square unitary plate that is orthogonal to the segment between σ x and σ y and passing through the midpoint of this segment.In this way, starting from a spin configuration σ, we can introduce a family of closed polyhedra (or 3d-Peierls contours configuration) Γ(σ) = {γ 1 , . . ., γ N } separating the regions with spin + 1 from those with spin −1. 2   Denote by B − the total number of −1 bonds in B Λ , that is the total number of edges with spins of opposite sign on its endpoints and by B T OT the total number of bonds in B Λ .
Denoting by .17) 2 The correspondence between σ and Γ(σ) is one to two for periodic boundary conditions and one to one if at least one spin is fixed, this includes the case of "plus" or "minus" boundary conditions.
the total number of plates of the contours configuration, we clearly have Performing simple algebraic calculations, we can rewrite H 2λ (σ): where V − (σ) is the number of sites with negative spin.
We have: B T OT = 3|Λ|.Moreover, using (3.18): Using this result, we can write, for k = 1, 2,: with: It is now straightforward to prove the next technical lemma:

.23)
Proof.A simple algebraic calculation shows that each factor of ξ k (σ, λ) is less or equal to the respectively factor of ξ k (σ, 0) for both k = 1 and k = 2.
As a consequence of the previous lemma: And hence: with: The first three terms of the r.h.s. of (3.27) do not depend on δ and the fourth one is analytic in δ.Therefore, to prove that the r.h.s. is analytic, it must be shown that the last term is analytic.As in [10], Ξ := Γ e −2J|Γ| ξ k (Γ, 0) can be written as the partition function of an abstract polymer gas.The analyticity of log Ξ |Λ| follows by showing that the activity of each polymer is sufficiently small.In our case a polymer γ is a single 3d-Peierls' contour (defined above) and its activity is ρ(γ) = e −2J|γ|ξ k (γ) , where ξ k (γ) is defined as in (3.26), when Γ consists of the single contour γ.Then the proof can be concluded following the same steps of the proof of [10,Lemma 2.2].This establishes the first part of Theorem 2.5.The second part of theorem says: if lim |Λ|→∞ δ 2 |Λ| = 0, then exist a J such that, for This follows applying [3, Theorem 1.1] and noting that, for J sufficiently small, Remark: If some spin in Λ is kept fixed, then, as already mentioned, the factor 2 in the last equality of (3.24) becomes a 1. Consequentely, |Λ| log Γ e −2J|Γ| ξ k (Γ, 0) still holds since the set of Peierls' countours Γ when some spin is fixed is a strict subset of the set of Peierls' contours when periodic boundary conditions are taken into account.
3.2.Proof of theorem 2.6.We want to show that the mean value of a spin, in the low temperature regime, at the centre of the lattice depends on the boundary conditions in the case of finite volume and continues to depend on the boundary even in the limit of infinite volume.We interpret this as the fact that the system is in the ordered phase.Of course if we have a external magnetic field different from zero, all spins follow the orientation of such external field.So to study the spontaneous behavior of the system, we go back to Hamiltonian (2.1), (2.2) and set λ = 0.Moreover, from now on, we fix the external spins of Λ (that we denote by ∂ ext Λ) to assume the value +1 that is we impose +1 boundary conditions.We have From Proposition 2.1 it follows: with where B + Λ is the set of all nearest neighbors pairs in Λ ∪ ∂ ext Λ: From now on, for convenience, we will omit the over-script + in the Hamiltonian, that is we write H in place of H + .
From Lemma 3.1, we have: where: We can, therefore, compute the mean value of σ 0 .We assume that the lattice goes from −L/2 to +L/2 in all the three directions, so σ 0 is the spin at the centre of the lattice, that is the furthermost from the boundary.The mean value of σ 0 with positive boundary conditions is: σ,τ e −H(σ,τ ) .
(3.35) Using (3.33) the last expression can be written as follow: Adapting the notation of the previous sections to the present one, we have: ) and, consequently, We now compute σ 0 + Λ in the low temperature regime J 1. Clearly: We can now estimate P + Λ (σ 0 = −1) using a contour representation, defining the 3d-Peierls contours as in the previous section.Let Γ(σ) = {γ 1 , . . ., γ N } be the family of 3d-Peierls contours associated to the spin configuration σ.
We now see how ω G (σ) and f (σ) can be written in terms of Peierls contours.
where |Γ| is the total length of all Peierls contours in Γ(σ): Proof.Let B + be the number of bonds in of positive sign B + Λ (the number of nearest neighbors with same sign) and with B − the number of bonds of negative sign (the number of nearest neighbors with opposite sign).Then (3.47) Now, let {γ 1 , γ 2 , . . ., γ N } be a contours configuration associate to the spin configuration σ.Then, by construction: as above.Moreover: Then: and so: Finally, is easy to show that B T OT = 3L 2 (L + 1).Then, we have: with: We can now compute P + Λ (σ 0 = −1) using equations (3.43), (3.44), (3.45) and (3.52) and noting that the terms e 3JL 2 (L+1) and 1 + δe −6J |Λ| appear both in the numerator and in the denominator as factors for all Γ: Since the denominator contains more terms than the numerator, it is a standard task (Peierls' argument) to show that (3.55) Now observe that Indeed each site in N 0 (γ 0 ) contributes to the length of γ 0 with 3 unitary plates; each site in N 1 (γ 0 ) contributes with two plates and each site in N 2 (γ 0 ) contributes with one plate.Then, a simple algebraic computation leads to: (3.57) Let us proceed now by transforming the sum on the contours of Peierls γ 0 {0} in a sum over their lengths: We indicate with η 0 (k) the number of Peierls contours of length k that surrounds the site 0. Moreover k ≥ 6 for closed contours.Then, we can write: For J large enough (low temperature regime): Finally, an estimate of the number of contours of length k surrounding the origin is given by Ruelle's lemma: Then (3.58) becomes: This expression identifies a geometric series with common ratio r = 3 e −6J + e −2q 1/3 .
Under this conditions, the series converges to the value: ∞ k=6 r k = r 6 1−r .
Then (3.61) becomes: . This expression tends toward zero for J 1 and q 1 uniformely in Λ.Therefore, for J and q large enough: This happens if the following condition holds: Moreover, in the same manner, we can show that, under the condition (3.63) we have: And this inequality holds even in the thermodynamic limit Λ → ∞: So, the 1-point correlation functions are not unique, in the low-temperature regime, with respect to boundary conditions.This suffice to asserts that the equilibrium state, that is the family of all n-points correlation functions in the thermodynamic limit (with n = 1, 2, ..., |Λ|), is not unique at low temperature, but it depends on boundary conditions.

Proof of Theorem 2.7.
Proof.We start by rewriting the Hamiltonian (3.29): where E Λ D = E J ∪ E q is the edge set of the finite tetrahedral lattice Λ D with periodic boundary conditions and e 1 , e 2 are two sites in Λ D linked by the edge e. Moreover: Then, the partition function can be rewritten as: We can naturally split a non vanishing graph into non intersecting connected components which we will call lattice animals.A lattice animal γ is thus nothing but a graph g with edge set E γ = {e 1 . . .e N } ∈ E Λ D formed by nearest neighbor links.The allowed lattice animals are only those γ with incidence number at the vertices equal to two or four.We denote with L Λ D the set of all possible lattice animals in Λ D .Two lattice animals γ and γ are non overlapping (i.e.compatible), and we write γ ∼ γ if and only if γ ∩ γ = ∅.We will denote shortly |γ| = |E γ |.In conclusion we can write (from 3.73): is the partition function of a hard core polymer gas in which polymers are lattice animals, i.e. elements of L Λ D with the incompatibility relation γ γ .Each polymer γ has an activity given by: ) where 2k and 2m denote respectively the number of J−edges and the number of q−edges of the closed polymer γ, so 2k + 2m = |γ|.Actually, we can imagine the lattice Λ D as a collection of layer of hexagonal lattice, the sites of which are linked only by J−edges, and these layers are linked to each other by edges of type q.So, a closed circuit must have an even number of q−segments and an even number of J−segments and the number of q−type segments must be less then or equal to half the number of J−type segments.To control the analyticity of Z Λ D , we can apply the Fernandez-Procacci convergence criterion (see [6]) to the Ξ Λ D (J e ).Namely, we need to find numbers µ(γ) ∈ (0, +∞) such that: where L γ = {γ ∈ L Λ D : γ γ} is the set of all polymers in Λ D incompatible with γ (i.e. the set of all polymers that intersect γ).Setting µ(γ) = ξ(γ)e a|γ| , the condition (3.79) becomes where with the sum (γ 1 ,...,γn)⊂L n γ , γ i ∼γ j (•) running over all ordered n-tuples of polymers.
Consider now the factor: we have thus to choose n lattice animals (γ 1 , . . ., γ n ) all incompatible with a given lattice animal γ and all pairwise compatible.We recall that two lattice animals are incompatible if they share a vertex in Λ D .Since γ has |γ| edges, we can find at most |γ| animals incompatibles with γ and pairwise compatible.Thus, the factor above is zero whenever n > |γ|, γ i γ.
We want to rearrange the sum over all ordered n-tuples of polymers (γ 1 , . . ., γ n ) ⊂ L n γ in a sum over all polymers {γ 1 , . . ., γ n } ⊂ L Λ D , regardless of their order.Recalling that the factor above is zero whenever n > |γ|, we have: Then, (3.81) becomes: The convergence condition (3.80) becomes: Observe finally that, due to the structure of the lattice, the function where 0 is the "origin" in Λ D .Now, recalling (3.78), the above condition becomes: We want to convert the sum over γ's passing through 0 in a sum over their lengths |γ| = 2k + 2m.To this end, we observe that in a closed circuit, the number q−segments must be less then or equal to the number of J−segments: 2m ≤ 2k; and the minimal number of edges must be 6: 2k + 2m ≥ 6.So, condition (3.90) becomes: [tanh(J)] 2k [tanh(q)] 2m e a(2k+2m) 1 ≤ e a − 1. (3.91) The sum: corresponds to the number of closed circuits of length 2k + 2m passing through 0.
We can find this number imagining to start from a certain point and doing 2k steps of J−type and 2m steps of q−type, until returning to the starting point.As long as we move on a layer we carry out all J−type steps, while when we change layer we carry out a q−type step.As long as we move on a layer, we have 2 2k−2m possible circuits (at each node, we have 2 possible directions); however, occasionally we have to insert a change of layer (a segment of type q), for a total of 2m segments of this type.So we have to insert 2m step of type q among the 2k steps of type J.We can do this in 2k 2m ways.Once we change layers, we have 3 possible directions for the first step in this new layer, so we need to consider also a factor 3 2m .So, the total number of such circuits is: (3.95) We observe that the two sums must satisfy the constraint k + m ≥ 3, that is the close circuits condition.We can then extend these sums to all value k ≥ 2, q ≥ 0 without any constrain if we subtract by hand the only term forbidden by the constraint; this term corresponds to k = 2 and m = 0. So, we finally have: We can finally perform the remaining sums over k: (3.100) Finally, recalling the form of x and y, we have: (2 tanh(J)e a ) 4 ( 3 2 tanh(q)e a − 1) 4 1 − (2 tanh(J)e a ) 2 ( 32 tanh(q)e a − 1) 2 + (2 tanh(J)e a ) 4 ( 3 2 tanh(q)e a + 1) 4 1 − (2 tanh(J)e a ) 2 ( 32 tanh(q)e a + 1) 2 + (2 tanh(J)e a ) 4 ≤ e a − 1. (3.101) This expression is the condition that J and q must satisfy to be sure that ) is an analytic function for in J and q.Numerical evaluations show that a good value of a is a = 0.15.For this value, the expression (3.101) identifies the region below the lower curve in Fig. 4. Hence, for values of J and q small enough (i.e. in the aforementioned region), f Λ D (J, q) is analytic.Thus, we have shown that in low-temperature regime the system is in the ordered phase, while in the high-temperature regime it is in the disordered one.Therefore there must be a critical line in J − q plane that separates the ordered phase from the disordered one and this curve must lie in the region between the two curves as shown in Figure 4.This fact is well supported by numerical simulations.

Numerical simulations
In this section we present our numerical results (obtained with techniques similar to those used in [5]) concerning the critical curve in the (q, J) plane and discuss the behavior of the critical exponents of the magnetic susceptibility as q varies.
To this end we consider the alternate dynamics taking place on a class of "tetrahedral" non homogeneous lattices.It is possible to give a geometric interpretation to the pair interaction thinking to J and q to be proportional to the inverse of the distance of the lattice points.In this way, the three dimensional lattice can be thought to be a collection of honeycomb layers (where the the side of each hexagon is 1 J ) at distance 1 q one from another.With this picture in mind, we observe that when q = J the dynamics lives on the diamond lattice whereas for q → ∞ the lattice becomes the simple cubic one.Moving in the other direction (towards smaller values of q) the interaction between the layers becomes weaker and weaker up to the point (q = 0) when they become independent so that the system resembles a collection of well separated graphene sheets.In this framework, estimating the critical curve in the (q, J) plane, amounts to finding the critical "size" J of the hexagons for each "distance" q between the sheets.Moreover, it is reasonable to think that the transition from the three dimensional model to a collection of independent two dimensional ones implies that the critical exponents of the magnetic susceptibility undergoes a sharp change at q = 0.
Both the critical values of J and the critical exponents can be estimated by looking at the variance of the magnetization of the system.Indeed, at the critical J this variance diverges.
To estimate the variance of the magnetization, we considered its sample variance computed over a long run for a somewhat large system, see below for more details.The critical values of J that we found are consistent with the actual critical value of J for the Ising model on the honeycomb lattice (q = 0) and with recent numerical estimates for the critical J for the diamond lattice (q = J) and the simple cubic lattice (q large).This fact gives us some confidence on our findings concerning the whole critical curve.
Our results are summarized in Figs. 4, 5 and 6.In particular Fig. 4 shows that the estimated critical curve lies in the region between the curves of equations (3.63) and (3.101) delimiting the low and the high temperature region respectively.Fig. 5 shows the values of the normalized standard deviation of the magnetization as a function of J for q = 0, q = J and q = 2 corresponding, respectively, to the collection of two dimensional honeycomb lattices, the diamond lattice and, ideally, the simple cubic lattice.Our estimate of the critical J for each value of q, denoted by Ĵc (q), is given by the value at which the variance is maximal.We have Ĵc (0) = 0.659.In this case the analytical critical value is J h.l.c ≈ 0.659 (see [1]) For q = J we obtained Ĵc (J) = 0.370 whereas in [8] the numerical estimate is J d.l.c = 0.370.Finally, setting q = 2 we estimated Ĵc (2) = 0.226.In [12] J s.c.c = 0.222 Figure 5.The standard deviation of the magnetization as a function of J for q = 2 (ideally the limit q → ∞), q = J and q = 0. Our estimates for the critical temperature are J c = 0.226 for the cubic lattice, J c = 0.370 for the tetrahedral diamond lattice and J c = 0.669 for the 2d honeycomb lattice.Approaching the critical temperature, the magnetic susceptibility χ (the variance of the magnetization) diverges with a power law with some critical exponent γ (see, e.g., [9,13]) Recalling that, in this paper we wrote J in place of the usual βJ, where β = 1 k b T , we get Figure 7.The regression lines for γ for some values of q.In each chart the red line is obtained looking at the values J < J c (high temperature) whreas the blue line is obtained looking at the values J > J c (low temperature).
for a suitable constant C.
Note that the value of γ is related to the dimension of the system and it is the same for the whole class of Ising-like lattice systems with the same dimension (see [9]).We estimated γ for several values of q ranging from 0 to 2, that is for geometries ranging from a collection of 2d honeycomb lattices to the simple cubic lattice.In this case, our estimates are not meant to determine the values of the critical exponent with high accuracy.Rather, as long as our values are consistent with those available in the literature, we want use them to support our conjecture that the systems retains a three dimensional structure for all positive qs.
Our results are summarized in Fig. 7 and Fig. 8.There it is possible to see that, for q > 0 both the high and low temperature critical exponents are quite close to the value γ ≈ 1.237 that is the critical value for three dimensional Ising systems (see [2]).On the other hand, as soon as q = 0, our estimate jumps to a value that is much closer to the critical value for two dimensional Ising system (γ = 7 4 , see [13]).These findings show that our model is able to capture through the variation of the parameter q the dimensional transition in the geometry of the system.4.1.Numerical details and heuristic discussion.We simulated the shaken dynamics, without external magnetic field, on a 96 × 96 × 96 grid on which we imposed periodic boundary conditions.
Figure 8.The low temperature (J > J c , in blue) and high temperature (J < J c , in red) coefficients for values of q between 0 and 2.0.The solid lines correspond to the "true" values of the critical exponents for the two and three dimensional systems.The numerically determined values appear to be close to the critical value of the three dimensional system for q > 0 and to the critical value of the two dimensional system for q = 0.The behavior of numerical estimate for the low temperature critical exponent in the case q = 2 is likely to be due to the limited size of the simulated system.
We considered a grid of points in the (q, J) and for each point in the grid, we started the simulation from the configuration with all spin set to −1 and let the system perform 510000 steps of the shaken dynamics (that is, 1020000 half steps).We considered the first 10000 steps as a "transient" and collected statistics on the final 500000 steps.In particular, for each pair of parameters (q, J) we computed the the average and variance (over time) of the magnetization.
Simulations have been carried over using the language "julia".
A heuristic insight on why this procedure should be useful it is possible to argue as follows.Letting the dynamics start from the configuration with all spins taking value −1, it is expected to reach very rapidly a local minimizer of the free energy and start visiting configurations that are close to this minimizer.
In the high temperature regime, the minimizers of the free energy are expected to have all zero mean magnetization and if the parameters (q, J) are in the high temperature region, the dynamics is likely to return very quickly to a state where the number of plus and minus spins is essentially the same.As a consequence, it is possible to conjecture that the average (over time) of the magnetization is very close to zero and that its variance is very small (see figure 6).
In the low temperature region, the free energy has minimizers whose mean magnetization is closer (and closer as the system freezes) to ±1.The colder the system, the higher the free energy barriers separating the "positive magnetization" minimizers from the "negative magnetization" ones.As the chain evolves, the dynamics will overcome a free energy barrier of magnitude ∆ with a probability that is exponentially small in ∆.Therefore, since the system starts from the configuration where all spins are −1, it will very likely reach the vicinity of one of these −1 minimizer and will stay, with very high probability, in the region where the minimizers of the free energy have negative mean magnetization.The typical time to observe a transition to the +1 minimizers are exponentially large in the volume, and hence they are way beyond the possibility of a numerical simulation.With probability very close to 1 the system will remain captured by the −1 minimizers.Consequently, also in the low temperature region we can expect a very small variance for the average magnetization whereas its mean is likely to be more and more negative as the system becomes colder (see figure 6).
Around the critical temperature, the free energy has minimizers with both positive and negative mean magnetization.However, the "valleys" of the free energy landscape where these minimizers sit are rather shallow and, therefore, the dynamics is expected to move between minimizers whose mean magnetization has opposite signs.An evolution of this type will produce an average magnetization that is close to zero.Nevertheless, the variance of the magnetization is expected, in this case, to increase when the temperature approaches its critical value.Note that the general theory of critical phenomena (see again [9]) shows that it is rather delicate to measure the features of the systems close to the critical temperature.The relatively good results we obtained with the simulations presented above show that the shaken (or alternate) dynamics is able to capture the features of the system also when the parameters are close-to-critical.

Figure 3 .
Figure3.The graph Λ D .Blue and red dots represent respectively the sites in the two three-dimensional lattices V 1 and V 2 .The black segments correspond to the interaction governed by the parameter J, the green segments correspond to the self-interaction, governed by the parameter q.
3.73) Developing the product e∈E Λ D [1 + b e tanh(J e )], we get terms of the type: [tanh(J e )] N b e 1 b e 2 . . .b e N , (3.74) which has a clear geometric interpretation: the set of bonds b e 1 . . .b e N form a graph (connected or not) in Λ D whose links are nearest neighbors.Performing the sum over σ we get that the only graphs which yield a non vanishing contribution to σ b e 1 . . .b e N , and hence to the partition function, are those whose vertices have incidence number two or four, while all other graphs are zero once the sum over configurations σ has been done.Graphs of this type are called non-vanishing.If the graph b e 1 . . .b e N is non-vanishing, then: σ b e 1 . . .b e N = 2 |Λ D | .(3.75) x∈Λ D x∈γ∈L Λ D |ξ(γ)|e a|γ| ≤ e a − 1.(3.87)

Figure 4 .
Figure 4.The curves delimiting the Low Temperature and High Temperature regions identified in equations (3.63) and (3.101) together with a numerical approximation (in red) of the conjectured critical curve.

Figure 6 .
Figure 6.The "critical behavior" of the mean magnetization (left) and the standard deviation (right) for the tested values of J and q