Ultra-slow dynamics in a translationally invariant spin model for multiplication and factorization

We construct a model of short-range interacting Ising spins on a translationally invariant two-dimensional lattice that mimics a reversible circuit that multiplies or factorizes integers, depending on the choice of boundary conditions. We prove that, for open boundary conditions, the model exhibits no finite-temperature phase transition. Yet we find that it displays glassy dynamics with astronomically slow relaxation times, numerically consistent with a double exponential dependence on the inverse temperature. The slowness of the dynamics arises due to errors that occur during thermal annealing that cost little energy but flip an extensive number of spins. We argue that the energy barrier that needs to be overcome in order to heal such defects scales linearly with the correlation length, which diverges exponentially with inverse temperature, thus yielding the double exponential behavior of the relaxation time.


I. INTRODUCTION
The nature of the slow relaxation of glassy systems in the absence of disorder remains one of the most intriguing problems in condensed matter physics. There are two classes of theories of glasses: one that argues for a bona fide thermodynamic transition into a glass phase at finite temperature, such as in the Random First Order Transition theory [1][2][3]; and the other that views glassy behavior as a purely dynamical phenomenon, in the absence of any finite temperature thermodynamic transition (see Ref. 4 for a review). A notable example of this latter scenario is provided by plaquette spin models [5][6][7], which are known to display super-Arrhenius equilibration rates that scale as an exponential of a power of the inverse temperature. It is then natural to ask if there is a principled way of designing models -without disorder, frustration, or phase transitions -that display ultra-slow dynamics.
One common expectation is that the above questions can be addressed by appealing to the correspondence between some statistical mechanics models with translationally invariant (short-range) interactions and hard computational problems, such as the tiling of the plane [8,9]. It may seem natural to speculate that, when regarded as descriptions of physical systems in contact with thermal reservoirs, these models would display long glassy relaxation rates that reflect the complexity of the computation. Below we produce a counterexample which suggests that the complexity class of a computation alone does not necessarily determine the low temperature behavior of relaxation rates into the ground state of the model representing the computation.
In this paper, we argue for the possibility of novel anomalously slow relaxation rates that scale as a double exponential of the inverse temperature for a disorder-free spin model, inspired by statistical mechanics representations of reversible classical computational problems [10]. This work was motivated by our attempts to use thermal annealing in finding the ground state of a model that, for appropriate boundary conditions, encodes either multiplication of two integers or factorization of semiprimes. Our initial simulations showed that in both the factoring problem and even in the "easy" case of multiplication, a problem in complexity class P, thermal annealing leads to a remarkably long time-to-solution, consistent with a double exponential scaling with the inverse temperature. (We note that another unsuccessful attempt at using thermal annealing of a different statistical mechanics representation of the factoring problem was described in Ref. [11].) The principal result of the current paper is that the doubleexponential behavior of the relaxation time survives for open boundary conditions, in which case we prove rigorously that the model representing the multiplication circuit displays no thermodynamic phase transition. While this ultra-slow relaxation rate represents a negative result in the context of solutions of computational problems, it points to a remarkable example of a translationally invariant, short-range interacting model that displays astronomically slow glassy dynamics and no thermodynamic transitions at any finite temperature. The origin of the extremely slow relaxation rate can be traced to the propagation of errors initiated at a single defect, i.e., a bit mismatch, created during the annealing process. Such defects cost little energy but result in an extensive region of errors that cannot be removed except through healing from the boundaries. We speculate that the activation barrier for this healing process scales linearly with the correlation length, ξ ∼ exp(1/T ), the temperature dependence of which reflects the absence of a finite temperature phase transition.
The rest of this paper is organized as follows: In Sec. II, we detail the construction of a translationally invariant twodimensional spin model with local interactions, starting from a multiplication circuit. In Sec. III, we compute the partition function of this model and prove the absence of finite temperature thermodynamic phase transitions. In Sec. IV, we introduce two other computational models (which also display no finite temperature transitions) that help us highlight the special features of our principal model described in Sec. II; and we compare the resulting relaxation rates computed through thermal annealing. The different processes required in healing thermally induced defects, which are responsible for the qualitatively different behaviors of the relaxation rates of the three models, are discussed in Sec. V, where we also speculate on the origin of the double-exponentially slow relaxation. Our conclusions are presented in Sec. VI. arXiv:1812.01621v1 [cond-mat.stat-mech] 4 Dec 2018 ′ FIG. 1. Schematic of gates placed at the sites of a regular lattice with each gate connected to its neighbors. The inset depicts the details of the connection between neighboring gates: each link connects two spins (bits) representing one of the outputs of one gate and one of the inputs of its neighboring gate, respectively, The ferromagnetic interaction J ensures that the output on one gate is compatible with the input to its neighbor. The q and q label the input/output configurations of the spins (bits) that satisfy the truth table for each of the two gates, respectively. Red dashed lines demonstrate the slicing discussed in Sec. III.

II. MODEL
In this Section, we present a classical statistical mechanics spin model on a two-dimensional lattice with short-range interactions that encodes the functionality of a reversible multiplication circuit. The two elements of the model are: (a) a regular lattice covered by logic gates, the states of which are dictated by truth tables of allowed configurations of input and output bits (or Ising spins) for each gate; and (b) an interaction energy that penalizes neighboring gates for which input/output states are incompatible with one another.
We illustrate the spin model in Fig. 1, where the gates (boxes) are connected so as to form a regular lattice. The inset highlights two neighboring gates, with their connections to their neighbors, where bit lines connect outputs to inputs. In the case depicted (that is relevant to the multiplication circuit we deploy), each gate has 10 line-connections, corresponding to 5 input and 5 output bits. The states q and q of the gates take values 0, 1, . . . , 2 5 −1, which label all possible 5-bit input states; here we utilize reversible gates, so that the (5) output bits are uniquely determined by the (5) input bits, and viceversa. The input/output relation (or truth table) for each gate is energetically enforced, and we also add interactions that penalize configurations of states where the inputs of one gate do not match the outputs of its neighbors. This penalty enters in the form of a ferromagnetic interaction, depicted in the inset of Fig. 1, which shows examples of parallel and anti-parallel spins corresponding, respectively, to matched and unmatched bits between the two neighboring gates.
In Fig. 2 we show the specific circuit that implements multiplication for 4-bit numbers. (Multipliers of L-bit numbers are implemented in the same fashion.) This circuit is a reversible logic version of the standard array multiplier. It utilizes 3 types of gates: LCARRY, RCARRY and CORNER. Each of these gates has 5 inputs and 5 outputs, with truth tables shown in panel (b) of Fig. 2. It is at the boundaries of the system where the inputs and outputs of the whole multiplication circuit are written to or read from. The result of a multiplication of two L-bit numbers, a = 2 L−1 a L−1 + · · · + 2 a 1 + a 0 and b = 2 L−1 b L−1 + · · · + 2 b 1 + b 0 , namely S = a × b = 2 2L−1 S 2L−1 + · · · + 2 S 1 + S 0 , is placed in the 2L-bit registers, S 0 , S 1 , ..., S (2L−1) . The bits in the input register S, as well as in the carrier input c, are all set to 0. (The output bits of the carrier c are also 0 at exit.) Note that, in this design, the total number of gates required for the multiplication of two L-bit numbers is N gates = 2L 2 .
In order to construct the Hamiltonian of the spin model, which formalizes the above discussion, we introduce the following notation. We start by labeling the gates by an index g, and the links or wires between gates by . It is useful to define the sets of input and output wires out of gate g, denoted by w in (g) and w out (g). Notice that if an output of gate g is connected to an input of gate g by a wire , then ∈ w out (g) ∩ w in (g ). On each wire we define two classical spins σ in and σ out , which we refer to as "twin" spins. These spins, which represent the Boolean variables via the re-lation x in,out = (1 + σ in,out )/2, must align ferromagnetically if the output of one gate is to match the input of its neighbor. We shall use the shorthand notation σ in (g) for the 5 variables σ in with ∈ w in (g), and σ out (g) for the 5 variables σ out with ∈ w out (g). Finally, we denote by L ∂ the set of links at the boundary, which we use to fix inputs/outputs to the circuit (these are wires at the boundary of the circuit).
With this notation, the Hamiltonian of the model is given by The function T g is the negation of the function T g that encodes the truth table of gate g: if the input and output spins satisfy the gate, T g = 1 (T g = 0), and if they do not, T g = 0 (T g = 1). H gates penalizes configurations (i.e., costs energy ∆ g > 0) if the input and output bits (or their spin equivalents) violate the truth table of gate g. The neighbor-gate compatibility is implemented through H links , containing ferromagnetic terms of strength J > 0 that enforce alignment between connected inputs and outputs of neighboring gates (sharing a link ). (Notice that we allow for the most general form of the Hamiltonian that involves gate-and link-dependent couplings, a choice that is relevant to the discussion of Sec. III below.) Finally, H boundary describes local bias fields h that fix the states of a subset of twin spins on a boundary link ( ∈ L ∂ ) when h → ∞. These biases fix the boundary conditions of inputs or outputs appropriate for the specific computation (multiplication or factoring) implemented by the system.
To be more precise, by fixing S = 0, the bits of a and b, and ancillary inputs c = 0 and outputs c = 0 using external fields h in Eq. (1), the model performs the multiplication S = a × b upon reaching the minimum energy state. Alternatively, factorization of the 2L-bit integer S is implemented by appropriate selection of external fields h that fix the bits of S and c = c = 0, while leaving the spins belonging to a and b free. In the case of multiplication, the ground state is uniquely determined by the boundary, since there is only one product S for given factors a and b. On the other hand, depending on the number S to be factorized, there may be multiple solutions. For a semi-prime 2L-bit number that is the product of two L-bit prime numbers, there are two solutions (differing by where the numbers a and b settle along the boundaries, since As already mentioned in the introduction, this paper was motivated by the astronomically slow (double exponential) time-to-solution found by our earlier annealing studies of the the relaxation into the ground state of the statistical mechanics model encoding the multiplication circuit presented above. The initial goal of those studies was to use this model to factorize semi-prime numbers. More interesting for applications to physical systems is that the above computationally inspired model also displays double exponentially slow relaxation rates for open boundary conditions (h = 0), a result we discuss in detail in Sec. IV. Hereafter we will concentrate on this free-boundary case, for which we prove that no thermodynamic phase transition occurs down to zero temperature (see below).

III. THERMODYNAMICS
In this Section, we derive the exact expression for the partition function of the model defined in Sec. II and we show that the system lacks a finite-temperature thermodynamic phase transition. In what follows, we work with open boundary conditions (h = 0). The partition function (omitting one layer of ferromagnetic bonds in the wires at the outputs of the whole circuit) is given by This expression can be computed via transfer matrices, if we properly slice the lattice in a sequence of layers. These layers contain a single gate g, as we depict in Fig. 1. An L-bit multiplier can be thus sliced into 2L 2 layers, one for each gate. The layer has two boundaries, which enclose the gate g. For all spins not attached to links in w in (g) ∪ w out (g), the transfer matrix acts trivially; we denote the spins on this trivial line byσ (g). Therefore the transfer operator can be written as T g = t g ⊗1 1σ (g) , where t g acts only on the spins on the links attached to the gate. Specifically, the layer functions as a transfer matrix between the spins σ out , ∈ w in (g) and the spins σ out , ∈ w out (g). Explicitly, denoting { 1 , 2 , 3 , 4 , 5 } = w in (g) and { 1 , 2 , 3 , 4 , 5 } = w out (g), t g can be written as The partition function for open boundary conditions can then be expressed as where the state corresponds to the sum over all possible configurations on the boundary. The state |Σ is an eigenstate of T g with eigenvalue λ g given by as follows from the non-trivial part t g of the transfer operator in Eq. Ngates g=1 λ g .
Recall that we omitted the 5L ferromagnetic bonds at the outputs of the circuit; accounting for those links yields or equivalently, factoring out the contributions from all the gates and the links, The resulting free energy of the model then reads, Notice that this free energy has no singularities as a function of β and hence, the spin model with open boundary conditions (h = 0) displays no finite temperature thermodynamic transition, independent of the values of the couplings ∆ g , J . Moreover, the free energy is independent of the specific truth table assigned to each of the gates. In particular, the lack of a thermodynamic transition survives for a model in which the row of RCARRY − CORNER gates is replaced by a row of RCARRY − LCARRY, and the single (double) U-turns at the top (bottom) of the multiplication circuit are removed (i.e., the corresponding J couplings set to zero.) While the resulting circuit no longer implements multiplication, it is this "multiplication-circuit inspired" (MCI) model that can be regarded as a translationally invariant physical system with short ranged interactions. Below, we show that both the multiplication circuit and the MCI model with open boundary conditions display remarkably slow relaxation into their ground states.

IV. THERMAL ANNEALING
In this Section we use classical thermal annealing to study the relaxation dynamics of the MCI model with open boundary conditions. We shall focus on the case of a uniform ferromagnetic interaction J = J and ∆ g → ∞, for which all gates satisfy their truth tables. In this case, the dynamics is carried out via a Metropolis algorithm that flips a gate state q, satisfying the truth table, into another satisfying state q .
We start from a random initial spin configuration and lower the bath temperature from T 0 to T = 0 during a time τ according to: We define E τ (T ) to be the value of the energy per link when the time-dependent temperature T (t) reaches a given value T .
In the limit of infinitely slow annealing, i.e., τ → ∞, one has E ∞ (T ) = E thermal (T ). Using Eq. (10), the total thermal energy per link is where we shifted the ground state energy to 0, which is convenient for the annealing studies below. Figure 3 shows E τ (T ) for the L = 512 multiplier and different values of τ . The solid black line corresponds to the equilibrium result in Eq. (12). Red lines show the numerical results for E τ (T ) for values of τ in the range from 2 9 to 2 21 . Notice that E τ (T ) monotonically decreases, yet the system is unable to reach its ground state for finite annealing times, falling out of equilibrium. The minimum or residual energy E τ (0) ∼ 1/ 2 τ defines a length scale τ that measures the typical distance between defects in the frozen state.
To analyze the data we consider fits to three functional forms of E τ (0) vs. annealing time, τ : ln τ , and (iii) ln E τ (0) ∼ − ln ln τ . These arise from three models of relaxation via activation over a barrier corresponding, respectively, to different dependences of the barrier ∆ B ( τ ) on the length scale τ : (i) ∆ B ( τ ) ∼ constant, (ii) ∆ B ( τ ) ∼ ln τ , and (iii) ∆ B ( τ ) ∼ a τ , a > 0. In Figure 4 we plot E τ (0) as a function of τ for both the multiplication circuit and for the MCI model for L = 512. The fits to the three types of annealing behavior along with the corresponding reduced chi-squared statistic χ 2 , which we use as a figure of merit to evaluate the quality of each fit, strongly support model (iii), namely ln E τ (0) ∼ − ln ln τ , which follows from an energy barrier that scales as ∆ B ( τ ) ∼ a τ , with a ∼ 1.
The energy E τ (0) reached at the end of the protocol for a given τ can be matched to the thermodynamic value of E thermal (T ) for some T > 0; we define τ (T ) as the time scale needed in order that E τ (0) reaches E τ (0) = E thermal (T ). In that case, τ = ξ(T ), the thermal correlation length, which, given that the model displays no finite temperature transition, only diverges at zero temperature as ξ ∼ e −J/T . This translates into a relaxation time for equilibration of the double exponential form: The remarkable appearance of such a slow dynamical relaxation time in a system for which we prove that there is no finite-temperature thermodynamic phase transition is the main result of our work. While we do not have a detailed analytical understanding of the origin of this behavior, in Sec. V we will argue that its origin can be traced back to the amplification of errors induced downstream of a single-bit defect in the course of the computation.
For comparison, and in order to sharpen the discussion of Sec. V, we shall consider two additional translationally invariant spin systems derived from computational models that display the behaviors of cases (i) and (ii) above. These simpler computational models, shown in Fig. 6, represent two different ways of wiring CNOT gates in a square lattice array. The first model, referred to as CNOT 1 , is illustrated in Fig. 6(a). In CNOT 1 the output control (C) and output target (T ) bits of a given gate are connected to the input control and the input target bits of a neighboring gate, respectively. In the second CNOT 2 model, shown in Fig. 6(b), the wires are switched: the output control (C) and output target (T ) bits of a given gate are connected to the input target and input control bits of a neighboring gate, respectively. While, as in the case of the MCI model, one can prove that neither of the resulting statistical mechanics models displays a finite temperature thermodynamic transition, the simple change in the wiring drastically changes the relaxational dynamics of the corresponding models. In particular, the freezing energy E τ (0) vs. annealing time behaves as cases (i) and (ii) for CNOT 1 and CNOT 2 , namely, ln E τ (0) ∼ − ln τ , and ln E τ (0) ∼ − √ ln τ , respectively.
Indeed, Figure 5 shows E τ (0) as a function of τ for a lattice of size L = 512 for the CNOT 2 model. As in the case of the multiplication and MCI models, we show the three corresponding fits to the three kinds of annealing behaviors, along with the χ 2 for each of these fits. As advertised, the best fit to the data corresponds to case (ii), namely ln E τ (0) ∼ − √ ln τ .

V. DEFECTS AND PROPAGATION OF ERRORS
The three classes of dynamics encountered above can be rationalized by considering the nature of the defects -links where input and output spins are mismatched -and the energy barriers associated with the process of healing them. Even though a single defect only costs energy 2J, it affects an extended region of downstream bits as the computation proceeds. Since the truth tables of all gates are always satisfied, another way of visualizing errors generated by a single defect is to follow the changes induced by the computation in the state, q = 0, 1, ...2 5 − 1, of each of the gates from those in a ground state in the absence of defects (the "parent state").
More precisely, we start with a zero-energy (i.e., ground) state and create an excited state with minimum energy 2J by flipping one of the input spins of a gate at the center of the system and realigning the rest of the spins downstream of the error, so that all bonds but the original defective one are satisfied. To visualize the buildup of errors in this excited state, we compute its overlap with the "parent" ground state. Figure 7 exemplifies single defects and their downstream bit flips induced by the computation for the four circuits CNOT 1 , CNOT 2 , multiplication, and MCI. The single error leads to the presence of two distinct regions A (white region in Fig. 7) and B (colored region in Fig. 7) that, while containing no error themselves, correspond to different ground states. We remark that there is no interfacial energy penalty -the cost is again only the 2J paid at the defective link. To relax the defect and √ ln τ + γ (ii) (green) with α (ii) = 1.37 , γ (ii) = 0.51; and (iii) ln Eτ (0) = −α (iii) ln ln τ + γ (iii) (red) with α (iii) = 2.14 , γ (iii) = 1.08. The χ 2 for the three fits are shown in the legend. The best fit is to form (iii), with a χ 2 close to 1. We note that the simulation results for the multiplication circuit and the MCI models fall on top of one another within the error bars of the data. √ ln τ + γ (ii) (green) with α (ii) = 1.08 , γ (ii) = 0.29; and (iii) ln Eτ (0) = −α (iii) ln ln τ + γ (iii) (red) with α (iii) = 1.77 , γ (iii) = 0.93. The χ 2 for the three fits are shown in the legend. The best fit is to form (ii), with a χ 2 close to 1. reach a global ground state one must flip the entire region A or B.
In Fig. 7(a) we consider the case of the CNOT 1 circuit. We note that when the defect is induced in the input target (T ) spin (bit) the line of downstream defects (dashed black line in panel (c.1)) can be healed by moving the origin of the defect line downward to the lower boundary of the circuit with no energy cost, as in the case of a domain wall in a one-dimensional Ising model. When defects are initiated from a flip in the input control (C) spin the region B is built out of disconnected single lines of flipped bits, which can be healed sequentially. In this case, the defect can be moved vertically towards the boundary by generating one additional defect at cost 2J. The additional defect can then be moved horizontally all the way to the boundary at no extra cost. Therefore, the defects can effectively propagate vertically by paying a cost 2J. This is therefore an example of class (i), with a constant energy barrier 2J.
In Fig. 7(b) we show the corresponding downstream effect of a single defect initiated at the C-input in the CNOT 2 circuit for which region B is a Sierpinski gasket (a fractal), identical to that found in the triangular plaquette model of Ref. 6. (The same pattern also occurs when the errors are initiated from a defect in the T -input.) We remark that the CNOT 2 circuit provides a new realization of this type of fractal structure on a square, instead of triangular, lattice. As in the triangular plaquette model [6,7], defects can be healed in a hierarchical manner by overcoming energy barriers that are logarithmic in the linear size of the triangular region, thus providing an example of class (ii). This hierarchical relaxation leads to a qualitatively slower super-Arrhenius characteristic relaxation time as function of temperature of the form τ ∼ exp 1/T 2 .
Finally, the downstream effects of a single defect in the multiplication circuit and MCI model (at the input of a LCARRY gate at the center of the system) are shown in Fig. 7(c). In this case there are three kinds of defects, depend- ing on which of the inputs (a; b or c; or S) to a gate is flipped to initiate the error. Fig. 7(c) shows the different shapes of the regions that are misaligned with the parent state for the different types of defects. For all but S defects, errors spread into a two-dimensional wedge. As detailed in the figure, the multiplication circuit and MCI model show the same behavior when the errors are initiated from the b-, cor S-inputs, but display error regions of the same size but of different orientation for defects originated from the a-input. We note that, unlike the CNOT models, here the overlaps do depend explicitly on the parent state. Panel (c.5) of Fig. 7(c) shows the distribution of overlaps for each type of defect computed for 1024 randomly chosen parent states.
Qualitatively we expect that healing for the case in Fig. 7(c) should proceed as depicted in the cartoon in Fig. 8. In order to flip region B, one progressively increases the bulge of region A invading B. The boundary or interface of the bulge must contain additional defects and hence the barrier scales as the length scale of the bulge, where defects are created to produce the curvature. It is this scaling of the energy barrier with the length of the buldge that makes this model qualitatively different from classes (i) and (ii), and consistent with the scaling of the barrier within class (iii). This cartoon rationalizes the sur- prisingly slow relaxation times observed in the data presented in Fig. 4, which are consistent with the double exponential behavior of the rate with the inverse temperature in Eq. (13). We summarize the results of this Section in Table I.

VI. CONCLUSION
We have introduced a translationally invariant twodimensional spin model, inspired by a reversible multiplication circuit, and studied its thermodynamic and dynamic prop-  erties. We have proven that this model displays no phase transition at any finite temperature. Even though the thermodynamic behavior is trivial, healing thermally excited defects is extremely difficult. Single defects cost little energy but healing them requires flipping an extensive number of spins (bits). In turn this leads to ultraslow dynamics with a remarkably long relaxation time that scales as a double exponential of the inverse temperature.
It is important to stress that, since for different boundary conditions the same statistical mechanics model represents computations of different complexity -for example, multiplication vs. factoring -an efficient algorithm for reaching the ground states and thus solutions to these problems should generically lead to qualitatively different times-to-solution. However, our results show that thermal annealing leads to astronomically long times-to-solution, independent of boundary conditions, and thus independent of the complexity of the computation. Thermal annealing is therefore ineffective in solving even P-class computational problems, such as multiplication. On the other hand, thermal annealing is the generic path to relaxation for physical systems in contact with a thermal reservoir. From this point of view, some computation-inspired statistical mechanics models may provide a framework for identifying disorder-free systems with ultra-slow dynamics despite of the absence of a finite temperature phase transition.