{'fraction_non_alphanumeric': 0.04917998581145933, 'fraction_numerical': 0.014188540666861411, 'mean_word_length': 4.172350528815023, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ground-state configuration space heterogeneity of random finite-connectivity spin glasses and random constraint satisfaction problems

We demonstrate through two case studies, one on the p-spin interaction model and the other on the random K-satisfiability problem, that a heterogeneity transition occurs to the ground-state configuration space of a random finite-connectivity spin glass system at certain critical value of the constraint density. At the transition point, exponentially many configuration communities emerge from the ground-state configuration space, making the entropy density s(q) of configuration-pairs a non-concave function of configuration-pair overlap q. Each configuration community is a collection of relatively similar configurations and it forms a stable thermodynamic phase in the presence of a suitable external field. We calculate s(q) by the replica-symmetric and the first-step replica-symmetry-broken cavity methods, and show by simulations that the configuration space heterogeneity leads to dynamical heterogeneity of particle diffusion processes because of the entropic trapping effect of configuration communities. This work clarifies the fine structure of the ground-state configuration space of random spin glass models, it also sheds light on the glassy behavior of hard-sphere colloidal systems at relatively high particle volume fraction.


Introduction
Spin glass models defined on finite-connectivity random graphs have two control parameters, one is the temperature T and the other is the interaction (or constraint) density, defined as the number M of interactions (constraints) versus the number N of vertices, α ≡ M/N. Theoretical work of the last ten years [1,2,3] have established the understanding that, a random spin glass system with many-body interactions will experience a clustering transition if either the temperature T is lowered below certain threshold value or the constraint density α is increased beyond certain threshold value. At the clustering transition, the configuration space of the system splits into an exponential number of Gibbs pure states and ergodicity is broken. This transition is followed by another phase transition of the configuration space, the condensation transition, as T is further lowered or α further increased. At the condensation transition, a sub-exponential number of large Gibbs states start to dominate the configuration space and hence the equilibrium property of the system [3]. These properties were first observed in fully-connected mean-field spin glass models [4]. Whether they are also valid in D-dimensional real-world spin glass systems (D = 3 or D = 2) is still a debated open issue.
Although the configuration space geometric properties of random finite-connectivity spin glass models at or after the clustering transition have been well characterized by statistical physics methods, much less is known about the configuration space structure just before the clustering transition. Why do an exponential number of Gibbs states suddenly appear at the clustering transition? Are they preceded by precursor structures in the ergodic phase of the configuration space? If yes, when do the Gibbs-state precursors start to form and how to describe their evolution? What are the impacts of these precursor structures to the equilibrium dynamical properties of the system? These questions are important for a full understanding of the structural evolution of the configuration space.
In this paper, as a continuation of our recent efforts [5,6,7,8], we study the fine structure of the configuration space of random finite-connectivity spin glass models in the vicinity of the clustering transition. We focus on the evolution of the groundstate configuration space (corresponding to T = 0) using the constraint density α as the control parameter. As demonstrated through two case studies, a heterogeneity transition occurs to the ground-state configuration space at certain critical value of α. At the transition point, exponentially many communities of configurations emerge from the ground-state configuration space. Each configuration community is a collection of relatively similar configurations and it forms a stable thermodynamic phase in the presence of a suitable external field; they are the precursors for the Gibbs states at the clustering transition. The entropy density s(q) of configuration-pairs as a function of configuration-pair overlap q is calculated by the replica-symmetric (RS) and the firststep replica-symmetry-broken (1RSB) cavity methods. Extensive numerical simulations are performed to confirm that, the entropic trapping effect of ground-state configuration communities leads to strong dynamical heterogeneity of diffusion processes within the configuration space.
Dynamical heterogeneity in the fully-connected p-spin interaction spherical model (with continuous spins) [9] was quantitatively investigated by Donati and co-workers [10] within the framework of the Franz-Parisi effective potential theory [11,12]. The results of the present paper for finite-connectivity systems with discrete spins are qualitatively similar to the results of [10]. Some aspects of the heterogeneity transition in random finite-connectivity systems are also shared by real-world D-dimensional spin glass and structural glass systems, such as hard-sphere poly-disperse colloidal systems, where the particle volume fraction plays the role of the constraint density α [10]. In real-world glass systems and supercooled liquid systems, dynamic heterogeneity occurs in real space. The present work is also related to the work of Krzakala and Zdeborova on the adiabatic evolution of a single Gibbs state of a finite-connectivity spin glass system as a function of temperature T [13,14].
Two model systems are studied in the paper. Section 2 concerns with the pbody spin glass system, which is equivalent to the random K-XORSAT (exclusiveor-satisfiability) problem of computer science. Replica-symmetric and 1RSB mean-field calculations are carried out to obtain the entropy density of ground-state configurationpairs. Section 3 focuses on the heterogeneity transition and the dynamical heterogeneity of the random K-satisfiability (K-SAT) problem. Results obtained by RS calculations and numerical simulations are reported in this section. We make further discussions in section 4.

The random p-body spin glass model
The random p-body spin glass model is defined by the energy function where σ ≡ (σ 1 , σ 2 , . . . , σ N ) denotes a spin configuration for the N vertices i ∈ [1, N], with each spin variable σ i ∈ {−1, +1}; the index a ∈ [1, M] denotes one of the M interactions of the system and J a is the quenched random coupling constant, whose value is fixed to J a = +1 or J a = −1 with equal probability; the set ∂a includes all the vertices that participate in the interaction a, its size is fixed to p, and each of its p different elements is randomly and uniformly chosen from the whole set of N vertices. Similar to ∂a, we denote by ∂i the set of interactions that involve vertex i. While each interaction a in the p-body spin glass model affects the same number p of vertices, the sets ∂i may have different sizes for different vertices i. Actually, when N is large enough, the probability that a randomly chosen vertex participates in k interactions is governed by the Poisson distribution P P (k) = e −c c k /k!, with mean value c = pM/N. In the p-body spin glass model (1), each interaction contributes either a positive energy +1 or a negative energy −1 to the total energy. If we set p = K, this model is equivalent to the random K-XORSAT problem with the energy function For the random K-XORSAT problem, each interaction a is also referred to as a constraint, whose energy (1 − J a i∈∂a σ i )/2 is either zero (constraint being satisfied) or unity (constraint violated). The energy E xorsat ( σ) counts the total number of violated constraints by the spin configuration σ. The constraint density α of the system is by α = M/N. The random K-XORSAT problem (or equivalently, the random p-body spin glass model) has been well studied in the statistical physics community. It serves as an interesting system for understanding the low-temperature equilibrium property of finiteconnectivity spin glasses [1,15,16,17], and for understanding the dynamical property of glassy systems [18,19]. This model is also closely related to error-correcting code systems of information science, such as the Sourlas code [20]. The ground-state configuration space structure of the random K-XORSAT problem has been investigated in great detail [16,17,21] and was found to depend only on K and the constraint density α in the limit of N → ∞.
For a given value of K ≥ 2 there is a satisfiability threshold α s (K). When the constraint density α is below α s (K), the ground-state energy of model (2) is zero (the system is in the SAT phase), but it becomes positive when α > α s (K) (the UNSAT phase). A zero-energy spin configuration of model (2) is referred to as a solution, and all the solutions form the solution space S of this system. The solution space of a large random K-XORSAT problem is non-empty only if its constraint density is in the range of α ≤ α s (K). We have α s (2) = 0.5, α s (3) = 0.918, and α(4) = 0.977 [16].
Before the solution space of the random K-XORSAT (K ≥ 3) becomes empty at α > α s (K), it experiences an ergodicity-breaking (clustering) transition at the threshold value α = α d (K), with α d (3) = 0.818 and α d (4) = 0.772 [16]. For α < α d (K), the whole solution space forms a single Gibbs state (the meaning of a Gibbs state is explained geometrically in the following subsection). On the other hand, for α > α d (K), the solution space is no longer ergodic but is composed of exponentially many solution clusters (Gibbs states), each of which containing exactly the same number of solutions. For α d (K) < α < α s (K), solutions from different solution clusters are separated by energy barriers that are proportional to the vertex number N [18].

Non-concavity of the entropy function and solution space heterogeneity
The similarity of two solutions σ 1 , σ 2 in the non-empty solution space S of a random K-XORSAT problem can be measured by the overlap value We denote by N (q) the total number of solution-pairs in the solution space S with an overlap value q. This number is exponential in N in the SAT phase α ≤ α s (K).  (4) and the large scale organization of the solution space of a random constraint satisfaction problem. s(q) is a concave function of the solution-pair overlap q when the constraint density α is low enough (left-most panel); at α = α cm , s(q) starts to be non-concave, marking the formation of many solution communities in the solution space (second panel from left); then as α increases to α = α d , s(q) becomes non-monotonic, the solution communities separate into different solution clusters and ergodicity of the solution space is broken; the solution clusters shrink in size as α further increases, and eventually the number of solution-pairs with intermediate overlap values q ∈ [q 1 , q 2 ] will be zero and s(q) is only positive for q > q 2 and q < q 1 .
Therefore in the limit of large N, an entropy density s(q) is defined as The entropy density function s(q) contains rich information about the structure of the solution space S (see figure 1). The shape of s(q) has several qualitative changes as the constraint density α increases. The first qualitative change occurs at α = α cm , where s(q) becomes non-concave. This concavity change corresponds to the formation of (exponentially) many solution communities in the solution space (the large-scale homogeneity of the solution space is then broken) [6,7]. We refer to α cm as the heterogeneity transition point. We introduce the following partition function Z(x) for the solution space S [7] where x is a coupling field between the solutions. When x is positive, solution-pairs with larger overlap values have larger weights in the partition sum (5). Under the field x, the mean solution-pair overlap value is q(x) = (1/N)d ln Z(x)/dx. For N → ∞, we see from (5) that, q(x) = arg max q (s(q) + xq). At x = 0, q(0) is equal to q 0 , the most probable solution-pair overlap value of the solution space. If the entropy density s(q) is concave in q ∈ [q 0 , 1], then q(x) increases continuously with x for x ≥ 0. However if s(q) is non-concave, because there exists a line of slope −x * which touches s(q) at two different points q = q a and q = q b > q a (figure 1), then the value of q(x) jumps from q a to q b at x = x * [6,7]. This discontinuity of q(x) reveals the existence of a field-induced first-order phase transition at x = x * . The solutionpairs exhibit two different levels of similarity. For x > x * , the partition function Z(x) are contributed overwhelmingly by (intra-community) solution-pairs with overlap values ≥ q b (overlap-favored phase), while for x < x * , Z(x) is dominated by (inter-community) solution-pairs with overlap values ≤ q a (entropy-favored phase). The difference q b − q a increases from zero as the constraint density α exceeds the critical value α cm . At α = α cm , the solution space is in a critical state, where the boundaries between different solution communities are elusive, and the field-induced phase transition is second-order.
The Hamming distance between two solutions σ 1 and σ 2 is defined as , which is related to the solution-pair overlap by D( σ 1 , σ 2 ) = (N/2)(1− q( σ 1 , σ 2 )). The solution space S can be represented by a graph of nodes and edges. Each node of this solution graph denotes a solution, and an edge is linked between two nodes of the graph if the corresponding two solutions has a Hamming distance not exceeding a specified value D 0 . For the random K-XORSAT problem at α < α d (K), there exists a minimum value of D 0 ≪ N such that all the nodes of the solution graph are in a single connected component [18]. We take this minimum value as our edge linking criterion. (For α > α d (K), if D 0 is not of the same order as N, the solution graph will be a collection of exponentially many disjointed components.) In the ergodic phase of α < α d (K), we may introduce two particles to the solution graph. Initially these two particles are residing on the same node, say σ 0 , of the graph. In case the particles are uncoupled, then each particle performs a random diffusion in the solution graph independent of the other: Suppose at time t the particle is at node σ, then at the next time step it will, with probability k σ /k max , make a move to a randomly chosen nearestneighbor of this node, where k σ is the number of attached edges (the degree) of node σ, and k max is the maximal node degree in the solution graph. In the case the particles are coupled by a field x, however, the particle diffusions are mutually influenced, and the visited node-pairs (say σ 1 and σ 2 ) are no longer uniformly distributed but are favored to more similar pairs by a factor of e xN q( σ 1 , σ 2 ) . In the thermodynamic limit N → ∞, when the coupling field x is larger than x * , then even at time approaching infinity, the two particles will still be diffusing in the neighborhood of each other and in the neighborhood of the initial node σ 0 . Such a strong memory effect at field value x > x * is a dynamical manifest of the existence of communities in the solution space.

Annealed approximation for s(q)
If one knows a solution σ 1 = (σ 1 1 , σ 1 2 , . . . , σ 1 N ) for the K-XORSAT system (2), it is convenient to perform a gauge transform σ i ← σ i σ 1 i to the spin value of each vertex i. Under this transform, (2) is simplified to (1/N) i σ i . Equation (6) is independent of the reference solution σ 1 . This is a wellknown property of the K-XORSAT problem, namely its solution space S has the same local and global structure when viewed from any of its solutions. Because of this nice property, instead of calculating the solution-pair number N (q), we calculate the number N 1 (q) of solutions which have an overlap value q with a reference solution σ 1 . We denote the corresponding entropy density also as s(q), i.e., s(q) = (1/N) ln N 1 (q). (With this slight abuse of notation, the solution-pair entropy density as defined by (4) is s(q)+s 0 (α), where s 0 (α) is the entropy density of the whole solution space at constraint density α.) The average value of N 1 (q) over the ensemble of random K-XORSAT problems with fixed vertex number N and constraint density α is For N → ∞ we obtain the following annealed approximation for s(q) as For K = 3, the function s ann (q) is concave in q only for α < 0.577. For 0.577 < α < 0.778, s ann (q) is non-concave in q but is still monotonic in the range q ∈ [0, 1]; the monotonicity of s ann (q) in q ∈ [0, 1] is lost when α exceeds 0.778 (figure 2). The same qualitative results are obtained for the random K-XORSAT problem with K ≥ 4.
The annealed approximation s ann (q) is an upper bound of the entropy density s(q) for a typical random K-XORSAT problem. These two quantities are identical only when N 1 (q) has the self-averaging property, i.e., the value distribution of N 1 (q) among the ensemble of random K-XORSAT problems of constraint density α approaches a delta function in the limit of N → ∞. To check whether this self-averaging property is valid, we need to calculate the mean value of N 2 1 (q). After some combinatorial analysis, we obtain where with s 2 (q, ρ) being expressed as Self-averaging of N 1 (q) requires that Numerical calculations reveal that (12) is satisfied only for values of q very close to 0 or very close to 1. In figure 2, s ann (q) is plotted as a black solid line if ∆ < 10 −3 and as a red dashed line if ∆ ≥ 10 −3 . As s ann (q) is only an upper bound to s(q), the fact that s ann (q) becomes nonconcave at α < α d (K) can not be used as a proof that s(q) also becomes non-concave at α < α d (K). We proceed to calculate s(q) by the cavity method of statistical physics.

Replica-symmetric mean-field analysis
For the gauge-transformed K-XORSAT system (6), if vertex i is involved in constraint a, we define a cavity probability p + i→a as the probability that σ i takes the value σ i = +1 when the constraint a is absent. For each constraint a we exploit the Bethe-Peierls approximation [1,22,23] and assume that, the spin states of the vertices j ∈ ∂a are mutually independent in the absence of a. Under this approximation, we obtain that if the vertex i ∈ ∂a takes the spin value σ i , the probability that constraint a being satisfied is equal to [ is discontinuous and has a hysteresis loop.
that is missing element i. Then the following belief-propagation equation can be written down for each vertex-constraint association (i, a): For a given energy function (6), there are M × K iteration equations, which form the replica-symmetric cavity theory [1,22,23]. After this set of belief-propagation iteration equations has reached a fixed point, the mean overlap with the reference solution is calculated as And the entropy density as a function of x is The mean value of s(x) as averaged over the ensemble of random K-XORSAT problems (at fixed value of α) can be calculated using the population dynamics technique [1,21]. By eliminating x from s(x) and q(x) we obtain the entropy density function s(q). The numerical results are shown in figure 2 (blue dots) for K = 3 at different values of α. The mean overlap function q(x) is shown in figure 3 for K = 3 and α ≈ 0.62. Similar results are obtained for the cases of K ≥ 4.
For the random 3-XORSAT problem, q(x) is a continuous and smooth function of field x when α < α cm (3) ≃ 0.6182. As α approaches α cm from below, however, the maximal slope of q(x) is proportional to (α cm (3) − α) −1 and diverges at α cm (3). This divergence is a consequence of the fact that the entropy density s(q) starts to be nonconcave at α = α cm (3). For α > α cm (3), q(x) as calculated by the RS cavity theory shows discontinuity and hysteresis behavior when x is close to certain threshold value x * , indicating the existence of two distinct phases of the solution space as viewed from the reference solution ( σ 1 ). One of the phases contains solution σ 1 and the other similar solutions, whose overlap with σ 1 is larger than certain characteristic value ≈ (q a + q b )/2, see figure 1. We regard this phase as the solution community of solution σ 1 . If we choose another solution outside the solution community of σ 1 , we will find that this new reference solution is also associated with a different solution community.
For α cm (K) < α < α d (K), the solution space of the random K-XORSAT problem is therefore formed by exponentially many solution communities. As the solution space is very heterogeneous at this range of α, the replica-symmetric mean-field theory probably is not sufficient to describe its statistical property. We now proceed to study the solution space heterogeneity using the 1RSB cavity theory.

First-step replica-symmetry-broken mean-field analysis
For the gauge-transformed model (6) under the coupling field x, to apply the 1RSB mean-field theory, the solution space is first divided into an exponential number of Gibbs states γ [1,21]. Each Gibbs state γ represents a subspace S γ ⊂ S, and its partition function is defined as We can then define a 'free energy' density f γ (x) as f γ = (1/N) ln Z γ (x). This free energy density is the sum of two parts, f γ = s γ + xq γ , where s γ is the entropy density of Gibbs state γ and q γ is the mean overlap level of solutions in S γ to the reference solution σ 1 . The total 1RSB partition function of the system is In the above equation, m is the Parisi parameter, and Σ(f ) is the complexity, which measures the 'entropy density' of Gibbs states at the free energy density level f [1,21]. If we set m = 1 in (17), each Gibbs state contributes a term e N fγ to Z 1RSB , and then Z 1RSB is the total partition function of the system, provided that the complexity Σ(f ) calculated at m = 1 is non-negative. The existence of a nontrivial 1RSB solution at m = 1 is also a signature of the instability of the RS theory of the previous subsection [2,3]. We therefore set m = 1 in our calculations. The details of the 1RSB cavity mean field theory are presented in Appendix A, and here we discuss some of the numerical results obtained by population dynamics on the random 3-XORSAT problem. At α < α cm (3), no nontrivial 1RSB mean field results are obtained. The fixed point of the 1RSB population dynamics reduces to that of the RS theory for all the overlap values 0 ≤ q ≤ 1. As α exceeds α cm (3), the replica-symmetric mean field theory becomes unstable for intermediate values of overlap q, and nontrivial fixed points of the 1RSB mean-field population dynamics are observed. For example, at α = 0.80 we find that the complexity Σ(q) is exactly zero for overlap value q ≤ 0.233 and q ≥ 0.675, and Σ(q) is negative for 0.233 < q < 0.675 (see figure 4). The fact that Σ(q) = 0 for q ≥ 0.675 suggests that, the solutions with overlap levels q ≥ 0.675 to σ 1 are in a single solution community. In the corresponding solution subgraph of this solution community, any two solutions (nodes) with the same Hamming distance D to σ 1 are connected by at least one path that involves only other solutions with the same Hamming distance D to σ 1 (i.e., the subspace of solutions having the same overlap q (≥ 0.675) with σ 1 is ergodic within itself). The second message we get, from the fact that Σ(q) < 0 for 0.233 < q < 0.675, is that the solution subgraph formed by all the solutions at the same overlap level q [∈ (0.233, 0.675)] to σ 1 is not ergodic within itself but is divided into greatly many disjointed connected sub-components, with a sub-exponential number of dominating ones. In this overlap range, the community entropy density and the total entropy density as obtained by the m = 1 1RSB mean-field theory can only be regarded as upper bounds for the true values. One needs to work with m < 1 to obtain better estimates for the community entropy density and the total entropy density.
The third message we get from Σ(q) = 0 for q ≤ 0.233 is that, the solutions with overlap levels to σ 1 less than 0.233 form an ergodic subspace within itself. More precisely, all the solutions at each overlap level q (≤ 0.233) to σ 1 are ergodic within themselves. The solutions in such a subspace of fixed overlap q come from different solution communities, but they are connected with each other in the solution graph even when only the edges inside the subgraph are remained.
The theoretical results as shown in figure 5   When the constraint density α exceeds the clustering transition point α d (3), the complexity Σ(q) at m = 1 becomes positive at q = 0 (see the exemplar case of α = 0.819 in figure 6). This result means that the largest subspace of solutions (with overlap level q = 0 to σ 1 ) becomes non-ergodic at α > α d (3), as expected.

Summary for K-XORSAT
The calculations of this section demonstrated that the solution space of the random K-XORSAT problem experiences a heterogeneity transition as the constraint density α approaches α cm (K), where an exponential number of solution communities start to form in the solution space. These solution communities separate into different solution clusters at a larger constraint density value α d (K), where an ergodicity-breaking transition occurs. For α cm (K) < α < α d (K), although the solution space as a whole is ergodic, the subspaces of solutions with intermediate overlap levels to a reference solution are non-ergodic within themselves.
Our theoretical results, combined with the results of [21] for α > α d (K), give a complete picture on the structural evolution of the solution space of the random K-XORSAT problem.

The random K-satisfiability problem
The random K-SAT problem is a famous model system for the study of typical-case computational complexity of NP-complete combinatorial satisfaction problems [24]. Its energy function E( σ), like the random K-XORSAT problem, is defined as a sum over M = αN constraints a: In (18), each constraint a affects a set ∂a of K randomly chosen vertices from the vertex set {1, 2, 3, . . . , N}; J i a is the preferred spin state of constraint a on the vertex i ∈ ∂a, it takes the quenched value +1 or −1 with equal probability. If at least one of the vertices i ∈ ∂a takes the spin value σ i = J i a , the energy of the constraint a is zero, otherwise its energy is unity. The solution space S of model (18) is formed by all the spin configurations σ of zero total energy (i.e., satisfying all the constraints). In model (18), each vertex i is constrained by a set (denoted as ∂i) of constraints.
After the experimental demonstration of a satisfiability phase-transition in the random 3-SAT problem by Kirkpatrick and Selman [24], studies on the solution space structure of the random K-SAT problem have been carried out through rigorous mathematical methods (see, e.g., [25,26,27]) and through statistical physics methods (see, e.g., [28,29,30,3,31,32]). The threshold constraint density α s (K) for the solution space to be empty is calculated to be 4.2667 for the random 3-SAT problem [30] and its values for K ≥ 4 are also predicted by the 1RSB zero-temperature energetic cavity method [33]; a lower bound on α s (K) is calculated by the zero-temperature long-range frustration theory [34,8]. Ergodicity of the solution space is broken at the clustering transition point α = α d (K) < α s (K) [29,30]. At α d (K) ≤ α ≤ α s (K), the solution space contains an exponential number of Gibbs states. The value of α d (K) is calculated by the 1RSB zero-temperature entropic cavity method in [3,35], reporting α d (3) ≃ 3.87 and α d (4) ≃ 9.38.
We are interested in the heterogeneity of the ergodic solution space at α < α d (K). In the following subsection we calculate the solution-pair mean overlap as defined by (3) using the replica-symmetric cavity method. The heterogeneity transition point α cm (K) is determined by this RS mean-field theory.

Replica-symmetric mean-field analysis
The partition function (5) is a weighted sum over all the solution-pairs ( σ 1 , σ 2 ). Each vertex is then associated with a spin vector-state (σ 1 i , σ 2 i ). Under the coupling field x, we define the cavity probabilityp a→i (σ, σ ′ ) as the probability that constraint a is satisfied if the vertex i ∈ ∂a takes the spin vector-state (σ, σ ′ ). Similarly, we define the cavity probability p i→a (σ, σ ′ ) as the probability that vertex i is in the spin vector-state (σ, σ ′ ) in the absence of the constraint a ∈ ∂i.
If we exploit the Bethe-Peierls approximation as mentioned in section 2.3, the following RS cavity equations can be written down forp a→i and p i→â where δ n m is the Kronecker symbol (δ n m = 1 if m = n and δ n m = 0 if m = n). There is a symmetry requirement for p i→a (σ, σ ′ ), namely that p i→a (+1, −1) = p i→a (−1, +1). This condition is a result of the fact that, the solution-pair ( σ 1 , σ 2 ) has the same contribution to the partition function (5) as the solution-pair ( σ 2 , σ 1 ).
After a fixed-point has been reached for the RS iterative equations, the mean value q(x) of the solution-pair overlap and the solution-pair entropy density s(q) as defined by (4) can then be calculated. For example, the mean solution-pair overlap is expressed as where . . . x means averaging under the coupling field x. The predicted mean overlap function q(x) for the random 3-SAT problem is shown in figure 7(a) for α in the vicinity of 3.75. Similar to the results of the random 3-XORSAT problem shown in figure 3, the continuity of q(x) changes as α exceeds a critical value α cm (3) ≈ 3.75. The jumping and hysteresis behavior of the mean overlap q(x) for α > α cm (3) indicates that solution communities start to emerge in the solution space of the random 3-SAT problem at α ≈ 3.75.
The overlap susceptibility χ(x) measures the sensitivity of the mean solution-pair overlap with the coupling field x, it is defined by χ(x) ≡ dq(x)/dx. The susceptibility χ(x) is related to the overlap fluctuation by Figure 7(b) demonstrates that, as α approaches 3.75 from below, the peak value of χ(x) becomes more and more pronounced and finally diverges. From the divergence of the peak value of χ(x), we obtain that α cm (3) = 3.7497 for the random 3-SAT problem. This value is much below the value of α d (3) = 3.87. Similar mean-field calculations are performed for the random 4-SAT problem and we obtain that α cm (4) = 8.4746. This value is again much below the clustering transition point α d (4) = 9.38.

Dynamical heterogeneity of Glauber dynamics
The solution space heterogeneity of the random K-SAT problem influences the dynamics of random walking diffusion processes [6]. Similar to the random K-XORSAT problem, we can represent the solution space of the random K-SAT problem as a solution graph of nodes and edges, with the nodes denoting individual solutions and the edges connecting pairs of solutions of unit Hamming distance. When the constraint density α of the random K-SAT problem is less than α d (K), this solution graph has a single giant connected component that includes almost all the solutions of the solution space. In this ergodic phase of α < α d (K), the structure of this huge solution graph become heterogeneous when the solutions aggregate into many different communities. There are many domains of high edge density in the solution graph, corresponding to the different solution communities. The nodes of different domains are also connected by many edges, but the density of inter-domain edges is much lower than the density of intradomain edges. As was demonstrated in [6], this heterogeneity of edge density causes an entropic trapping effect to diffusive particles on the solution graph. The dynamics of a diffusive particle can be decomposed into a trapping mode (the particle diffuses within a relatively dense-connected domain of the solution graph) and a transition mode (the particle escapes from one domain of the solution graph, wonders for a while, and then enters into another domain of the solution graph). As the trajectory of the diffusive particle oscillates between the trapping mode and the transition mode, if a clustering analysis is performed on a set of solutions sampled from this trajectory at equal time interval, a clear community structure can be observed among the sampled solutions [6,8].
The solution space diffusion process can be turned into a stochastic search algorithm. One such algorithms, the SEQSAT of [36], constructs a solution for a random K-SAT problem in a sequential manner. Constraints of the problem are added one after another in a random order, and as a new constraint is added, a random walk process of single-spin flips is performed in the solution space of the satisfied sub-problem to reach a solution that also satisfies the new constraint. It was observed [36,8] that, when the constraint density α of the satisfied sub-problem exceeds α cm (K), the SEQSAT search process becomes viscous, the mean waiting time needed to satisfy a new constraints starts to increases rapidly with α, and the sequence of waiting times starts to have large fluctuations. This dynamical behavior is easily understood in terms of the heterogeneity of the underlying solution space.
Solution space structural heterogeneity results in dynamical heterogeneity of diffusion processes. To demonstrate this point more clearly and to estimate a typical relaxation time, we study in this subsection a simple solution space Glauber dynamics. For a large random K-SAT problem with N variables and M = αN constraints (α < α d (K)), first we construct through SEQSAT a spin configuration that satisfies all the M constraints. (Of cause we can also use other heuristic algorithms to generate an initial solution. The only requirement is that this solution should be a typical solution, or in other words, it should belong to the single largest solution cluster of the solution space. This requirement is satisfied in our simulation studies.) Then we set the initial time as t = −T e and denote the initial solution as σ(−T e ) ≡ (σ 1 (−T e ), σ 2 (−T e ), . . . , σ N (−T e )). The spin configuration is then updated by single-spin flips at each elementary time step ∆t = 1/N. Let us suppose at time t the spin configuration is σ(t). Then a vertex i is chosen uniformly randomly from the whole vertex set {1, 2, . . . , N}; a candidate spin configuration σ ′ is constructed, with σ ′ j = σ j (t) if j = i and σ ′ i = −σ i (t). If σ ′ is not a solution of the K-SAT problem, then at time t ′ = t + ∆t, the old spin configuration is kept, i.e., σ(t ′ ) = σ(t). However, if σ ′ is also a solution, then with probability one-half σ(t ′ ) = σ(t) and with the remaining probability one-half σ(t ′ ) = σ ′ .
A unit time of the above-mentioned Glauber dynamics corresponds to N spin-flip attempts. The actual number of accepted spin flips in a unit time is about 0.1N for the problem instances of figure 8 and figure 9. If this random walk diffusion process is simulated for an extremely long period of time, every solution in the connected component of the solution graph which the initial solution σ(−T e ) belongs will have the same frequency of being visited. The time T e is set to be large enough (e.g., T e ∼ 10 8 ) to ensure that the diffusion process has completely forget the initial solution σ(−T e ) at time t ≥ 0.
For the trajectory of spin configurations at t ≥ 0, the quantity q(τ ) = (1/N) N i=1 σ i (t)σ i (t + τ ) measures the overlap between two spin configurations σ(t) and σ(t + τ ) that are separated by a time τ . q(τ ) is a random variable, its value fluctuates with different choices of the time t. We are mainly interested in the mean value q(τ ) and the variance ∆(τ ) of the overlap q(τ ). The mean overlap q(τ ) is calculated by where . . . t means averaging over different starting times t ≥ 0 along the diffusion trajectory. The variance ∆(τ ) is expressed as Comparing (24) with (22), a dynamical susceptibility χ 4 (τ ) can be defined as χ 4 (τ ) ≡ N∆(τ ) (such a quantity was introduced in [37,10], see also the review [38]). The upper panel of figure 8 shows the relaxation behavior of q(τ ) for a large random 3-SAT formula with N = 10 5 vertices. To study how the shape of q(τ ) changes with constraint density α, the first M = αN constraints of this same formula is used for Glauber dynamics simulation at each value of α. We focus on α values in the vicinity of α cm (3). For α = 3.70, the mean overlap q(τ ) reaches its plateau value of 0.37 at τ ≈ 10 6 . For α = 3.75 ≃ α cm (3), the plateau of q(τ ) is reached at τ ≈ 10 6.5 , and its value increases to 0.38. For α = 3.80, the plateau of q(τ ) is reached at τ ≈ 10 7.2 , with a value of about 0.39. These and our other unshown simulation results clearly confirm that, as α exceeds α cm (3), the relaxation of the solution space diffusion process is slowed down greatly.
The peak time τ dh of the variance ∆(τ ) gives a measure of the typical time scale of dynamical heterogeneity of the Glauber dynamics. For τ ≪ τ dh , two solutions σ(t) and σ(t + τ ) have a large chance of being in the same solution community, therefore the fluctuation of overlap values q(τ ) is small; a large fraction of vertices are inactive in this relatively short time window of τ , as their spin values are flipped only with low frequencies. On the other hand, for τ ≫ τ dh , the compared two solutions σ(t) and σ(t + τ ) have a large chance of belonging to two different solution communities, and therefore the fluctuation of their overlap values is again small; the spin values of most of the vertices have been flipped many times during such a large time window, and hence dynamical heterogeneity is destroyed. When the time window τ is comparable to τ dh , the solutions σ(t) and σ(t + τ ) have comparable probabilities of being in the same solution community and being in two different communities; this causes relatively large fluctuation of the overlap values. As α increases from α cm (3) further, it takes more time for the diffusion process to escape from a solution community, and the difference between intra-and inter-community overlap values become larger, these two facts make the peak value ∆(τ ) and the peak time τ dh both to increase rapidly. Figure 9 reports the simulation results for a random 4-SAT problem instance of N = 10 5 vertices and M = αN constraints. Similar to the case of random 3-SAT, we find that the viscosity of the Glauber dynamics increases rapidly with constraint density α. For α = 9.10 > α cm (4), the mean overlap q(τ ) reaches its plateau value only at a very large value of τ ≈ 10 8 . As demonstrated by the variance ∆(τ ) of solution-pair overlaps, dynamical heterogeneity become to manifest itself at α = 8.50 ≃ α cm (4) and it becomes more and more pronounced as α further increases.
Comparing figure 8 and figure 9 we can also notice an important dynamical difference between the random 3-SAT and the random 4-SAT case. For the random 4-SAT system, we observe that the relaxation curve of q(τ ) is 'Z' form in shape at α = 9.10 (i.e., 97% of α d (4)), with a plateau value of q ≈ 0.7 at intermediate time intervals τ before it finally decay to a much lower q ≈ 0.3 at τ ∼ 10 8 . The peak of the overlap variance ∆(τ ) corresponds to the time interval τ at which q starts to decay from the larger plateau. However, for the random 3-SAT system, we find that the curve of q(τ ) has a 'L' shape even at α = 3.80 (98% of α d (3)), and at the time when ∆(τ ) reaches its peak value, the mean overlap q(τ ) is already close to its asymptotic value of τ → ∞. We suggest that this dynamical difference is a strong reflect of the difference between the community structures of the random 3-SAT and the general random K-SAT problems with K ≥ 4.
For the random 3-SAT problem, in the heterogeneity phase of α cm (3) < α < α d (3), the solution space probably is dominated by a sub-exponential number of largest solution communities. Because of this predominance, the mean intra-community solution-pair overlap of each of these dominating communities is only slightly higher than the mean solution-pair overlap of the whole solution space. Then at typical time τ dh when the diffusive particle jumps between different dominating solution communities, one will not observe much drop in the mean overlap value q. Such a picture is consistent with the prediction that, at α > α d (3), only a sub-exponential number of solution Gibbs states dominate the solution space [3].
For the random K-SAT problem with K ≥ 4, however, in the range of α ∈ (α cm (K), α d (K)) the solution space probably is contributed mainly by a exponential number of median-sized solution communities. The sub-exponential number of the largest solution communities has only a negligible contribution to the whole solution space and therefore barely influence the dynamics of the diffusive particle. For such a community structure, then the mean intra-community solution-pair overlap will be much larger than the mean solution-pair overlap of the whole solution space, resulting in a change of trend of q(τ ) at τ ∼ τ dh . Such a picture for K ≥ 4 is again consistent with the prediction that, at α > α d (K) the solution space is dominated by an exponential number of median-sized Gibbs states [3].

Summary for K-SAT
In this section, we confirmed that the solution space of the random K-SAT problem becomes heterogeneous at α ≥ α cm (K) and determined by the replica-symmetric cavity method that α cm (3) ≃ 3.75 and α cm (4) ≃ 8.47. We demonstrated that the existence of many solution communities in the solution space caused heterogeneous behavior in the dynamics of a solution space diffusion process. The typical time scale of dynamical heterogeneity of this diffusion process is determined by computer simulations.

Outlook
A heterogeneity transition was found to occur in the ground-state configuration spaces of two multiple-spin interaction systems, the random K-XORSAT problem and the random K-SAT problem. We expect that this transition is a general phenomenon that occurs in many other spin glass systems before the ergodicity-breaking transition of the groundstate configuration space. Such a heterogeneity transition is unlikely to be special to the ground-state configuration space but should also be observed as the energy level (or equivalently the temperature T ) of the configuration space is lowered to certain critical level. If the configuration space of a spin glass system is ergodic but highly heterogeneous at certain temperature T , this structural heterogeneity probably will manifest itself through heterogeneous behaviors in various spin relaxation dynamical processes of the system. More deep understanding on the relationship between the phenomenon of dynamical heterogeneity and the structural heterogeneity of the configuration space is very desirable. Such efforts may bring new ways of probing configuration space heterogeneity from observing features of dynamical heterogeneity.
The solution space heterogeneity of the random K-SAT problem has been studied analytically only through the replica-symmetric cavity method. However, from the experiences gained on the random K-XORSAT problem, we believe the replicasymmetric cavity theory is not sufficient for a heterogeneous solution space. A complete study on the heterogeneity transition of the random K-SAT using the 1RSB mean-field cavity theory will be reported in a later paper.
In the case of the random K-SAT problem, we have not yet performed a systematic investigation on the scaling behaviors of the peak value of overlap variance (equivalently, the overlap susceptibility χ(τ )) and the typical time τ dh of dynamical heterogeneity. Probably both the peak value of χ(τ ) and the characteristic time τ dh diverge at the clustering transition point α = α d (K). To get unambiguous results, we need a more efficient protocol of simulating the diffusion dynamics on the solution space.
Simulations results of [6] and [5] indicated that, for the random 3-SAT and random 4-SAT problem, the single solution space Gibbs states at α > α d (K) are themselves very heterogeneous in internal structure. To study analytically the heterogeneity of single solution clusters in the ergodicity-breaking phase, however, appears to be a very challenging task. Such kind of investigations may be very valuable for us to understand how the largest solution clusters of the random K-SAT problem evolve with constraint density α. Appendix A. The 1RSB cavity equations for the random K-XORSAT problem Here we list the technical details of the 1RSB mean-field cavity theory for the random K-XORSAT problem as discussed in section 2.4. For a comprehensive description of the 1RSB mean-field theory, the reader is referred to [1,2,3,35].
When there are an exponential number of Gibbs states, the probability p + i→a as defined in section 2.3 will be different for different Gibbs states. The distribution of p + i→a among all the Gibbs states is denoted by P i→a (p + i→a ). We define p + i→a as the average value of p + i→a , i.e., p + i→a = dp + P i→a (p + )p + . We also define two auxiliary distributions as [2] Q + i→a (p + i→a |p + i→a ) = .

(A.2)
Q + i→a (p i→a |p + i→a ) is the conditional probability of p + i→a given that σ i = +1. Similarly, Q − i→a (p i→a |p + i→a ) is the conditional probability of p + i→a given that σ i = −1. At the special case of m = 1, it can be shown that the iteration equation for p + i→a has the same form as (13), with the only difference that all the p + j→b values are replaced by their corresponding mean p + j→b , i.e., p + i→a =p({p + j→b }) [3]. The iteration equations for the conditional probabilities Q + i→a (p + i→a |p + i→a ) and Q − i→a (p + i→a |p + i→a ) at m = 1 have the following expression [2,3] where σ ∂b\i denotes a spin configuration for the vertex set ∂b\i of constraint b, and is the probability of a satisfying spin assignment σ ∂b\i for constraint b given the spin value σ i of vertex i (the probability p − j→b ≡ 1 − p + j→b is the mean probability of vertex j taking spin value σ j = −1 in the absence of constraint b).
The cavity iterative equations for p + i→a and Q σ i i→a (p + i→a |p + i→a ) can be solved by the population dynamics technique [1,21]. We are interested in the solution space property at a fixed value of the overlap value q, so in the population dynamics simulation the magnitude of the coupling field x is adjusted from time to time to ensure that the mean overlap value as expressed by (14) with p + j→a replaced by p + j→a is equal to the specified value q (see, e.g., [39,21]).
The grand free-energy density g of the model system, as defined by g = (1/N) ln Z 1RSB , has the following simplified expression at m = 1: The mean free energy density f of a Gibbs state is expressed as where ∆f i and ∆f a are, respectively, the free energy increase caused by vertex i and constraint a. These two free energy increases are expressed by the following expressions at m = 1: and ∆f a = σ ∂a   w a (σ ∂a ) i∈∂a Q σ i i→a (p + i→a |p + i→a )   ln 1 + i∈∂a (2p + i→a − 1) 2 .
(A. 8) In (A.7)), p + i is expressed as The equality s total = s cm + Σ holds at m = 1. It appears that at some overlap values q, more than one mean-field 1RSB solutions can be produced by the population dynamics simulation at m = 1. These different mean-field solutions have the same value of s total but different values of Σ and s cm . In such a case we choose the mean-field solution with the largest value of Σ (similar situations of multiple mean-field solutions were observed earlier in the random K-SAT problem at m < 1 [40,35]).