4 we will assume the order parameter is uniform: Formal extension to a non-uniform g ( r ) is straightfor- ward. In the following, it will be often convenient to present the six independent entries of the symmetric tensors u$_{ij}$ and d$_{ij}$ as 6-component vectors, which is sometimes re- ferred to as the Voigt (or Voigt-Mandel) notation. 26,27 Here we choose the specific realization of the Voigt nota- tion, in which the tensor Λ is a second rank (positively de- fined) tensor of size 6 × 6. Upon a suitable linear transfor- mation, see Appendix A, we can formally rewrite Eq.(6) as where and Throughout the article, we will consistently indicate 6- component vectors with arrows: ⃗ g = ∑ 6 α =1 g$_{α}$⃗ e$_{α}$ , where ⃗ e$_{α}$⃗ e$_{β}$ = δ$_{αβ}$ ; and three component vectors with bold sym- bols: u = u$_{i}$ e$_{i}$ , where e$_{i}$e$_{j}$ = δ$_{ij}$ . We will not apply the Einstein summation convention to Greek indeces. Even though we have written ⃗ g ( r ) (or d$_{ij}$ ( r )) as a con- tinuum field, it is understood that ⃗ g ( r ) actually stands for a discrete vector located at the site of the lattice of the non-linear Einstein oscillators that is closest to the point r in space. Switching between discrete summation over lattice sites and continuous integration over space is straightforward: If the volume occupied by a single vector ⃗ g is a 3 - implying k$_{micro}$ ≃ π/a - we may interchange the summation and integration according to a 3 ∑ ↔ ∫ dV . Since vectors ⃗ g correspond to 6 degrees of freedom, the re- gion corresponding to one vector ⃗ g must contain at least two atoms. ⃗ φ ( r ) (or u$_{ij}$ ( r )) is of course, too, an ap- proximation to acoustic modes of a discrete lattice, at wave-lengths exceeding 2 π/k$_{micro}$ . Note that the number of these acoustic modes is three per vector ⃗ g ( r ), corre- sponding to the three translational degrees of freedom of the vector-containing region. For the sake of argument we will consider the simplest case of an isotropic reference state: where λ and µ are the Lam´ e coefficients. 25 µ is also called the shear modulus. In the isotropic case, the components of vectors ⃗ φ and ⃗ g have a particularly lucid meaning. One component - called “dilatational” - corresponds to uniform compression: and the rest five - called “isochoric” - to pure shear: where K ≡ [ λ + (2 / 3) µ ] is the usual bulk modulus. 25 We inquire about the possibility of (meta-stable) states with non-zero g by testing for non-concavity of F ( g ) from Eq.(3). The latter peculiarity is only possible if F$_{0}$ ( g ) has an inflection point. An analytic evaluation of F$_{0}$ ( g ) ap- pears difficult (except in the mean-field limit, see next Section). Despite the simple bilinear form in Eq.(6), the functional integration over u ( r ) in Eq.(4) at a fixed con- figuration d$_{ij}$ ( r ) could not be performed by simply chang- ing variables u = u ′ − d and performing a Gaussian inte- gration with respect to the new variable u $^{′}$: In contrast to the purely harmonic, long-wavelength motions u , the new displacement field u ′ would now contain short-wavelength motions, which, by construction, are anharmonic. It will suffice for our purposes here to extract the small and large g asymptotics of F$_{0}$ ( g ), which one may accom- plish by first averaging the Boltzmann weight e − β ∫ dV f ( r ) from Eq.(4) over the directions of the ⃗ g vectors on each site. By the aforementioned prescription for interchange between discrete summation over the oscillator sites and integration in space, we break up the integral in the exponent into a sum ∑ i a $^{3}$f ( r$_{i}$ ). Note the formula 〈 exp( ⃗ x⃗ y ) 〉 = I$_{m/}$$_{2}$$_{−}$$_{1}$ ( xy )( xy/ 2) 1 − m/ $^{2}$Γ( m/ 2), where the averaging 〈 . . . 〉 is with respect to the mutual angle be- tween two m -component vectors ⃗ x and ⃗ y (see Ref. 28 or formula 9.6.18 of Ref. $^{29}$). Here, I$_{ν}$ ( x ) is the modified Bessel function of the first kind of order ν and Γ is the gamma function. 29 Using this formula, we average e − βa 3 ⃗ φ⃗ g over the directions of ⃗ g on each site and obtain the following free energy at site r , per non-linear oscilla- tor: up to an additive constant. The subscripts φ, g at f$_{φ,g}$ ( r ) above signify that the latter function is derived from the original free energy density f ( r ) from Eq.(8), but now depends only on the absolute values of the 6- component vectors ⃗ φ and ⃗ g characterizing the long- and short-wavelength distortions respectively. We note the expression (14) represents the exact free energy of our ansatz, up to the small numerical ambiguities related to the detailed shape of the cut-off surface in the Fourier space and the lattice of the spins. f$_{φ,g}$ ( r ) = 1 2 [ φ $^{2}$( r ) + g $^{2}$] − 1 βa 3 ln I$_{2}$ [ βa $^{3}$φ ( r ) g ] [ βa $^{3}$φ ( r ) g ] 2 , (14) Λ$_{ijkl}$ = λδ$_{ij}$ δ$_{kl}$ + µ ( δ$_{ik}$δ$_{jl}$ + δ$_{il}$δ$_{jk}$ ) , (11) g 2 ≡ 6 ∑ α =1 g 2 α = d$_{ij}$ Λ$_{ijkl}$ d$_{kl}$. (10) φ 2 ≡ 6 ∑ α =1 φ 2 α = u$_{ij}$ Λ$_{ijkl}$ u$_{kl}$, (9) f ( r ) ≡ f$_{⃗}$ φ,⃗ $_{g}$( r ) = 1 2 ( ⃗ φ + ⃗ g ) $^{2}$, (8) g $^{2}$( r ) = g $^{2}$. (7) g 2 1 = Kd 2 $_{ii}$, (12) 6 ∑ α =2 g 2 α = 2 µ ( d$_{ij}$ − 1 3 δ$_{ij}$d$_{ll}$ ) $^{2}$, (13) II. GLASS AS A FROZEN-IN, MULTIPLY DEGENERATE STRESS PATTERN