Systems with Symmetry Breaking and Restoration

Abstract

Statistical systems, in which spontaneous symmetry breaking can be accompanied by spontaneous local symmetry restoration, are considered. A general approach to describing such systems is formulated, based on the notion of weighted Hilbert spaces and configuration averaging. The approach is illustrated by the example of a ferroelectric with mesoscopic fluctuations of paraelectric phase. The influence of the local symmetry restoration on the system characteristics, such as sound velocity and Debye-Waller factor, is discussed.

1. Introduction

It is generally accepted that phase transitions are related to symmetry changes, so that the transformation from a disordered to an ordered phase is accompanied by spontaneous symmetry breaking [1,2,3]. There also exist more complicated cases, when the symmetry, being broken in the major part of the sample, at the same time, is restored in other parts of the same system, which is called symmetry restoration [4]. Reciprocally, in the generally disordered phase, there can appear local regions with broken symmetry. In dynamical modeling, the arising regions of symmetry that differs from the symmetry of the surrounding matrix are related to the generation of solitonic clusters and local coherent structures [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]. It turns out that the systems with such local structures are more stable, as compared to homogeneous systems. Generally, the systems that become more stable by spontaneously changing their properties pertain to the class of self-optimizing systems [20]. There are numerous examples of condensed matter, where the sample is not homogeneous, but consists of regions of different symmetry. One tells that such systems display mesoscopic phase separation, and they are termed heterophase. This, for instance, concerns many high-temperature superconductors, in which superconducting regions coexist with normal regions [21,22,23,24,25], and some low-temperature superconductors [26,27]. Such a coexistence of superconducting and normal phases can occur even in atomic nuclei [28]. In many magnetic materials, magnetically ordered phase includes paramagnetic clusters [29,30,31,32,33]. Around the points of structural phase transitions, there exist regions, where phases with different symmetry coexist [34,35,36,37,38,39,40]. A number of ferroelectric materials displays the coexistence of ferroelectric and paraelectric phases [41,42,43,44,45,46,47,48,49,50,51,52,53,54], which also concerns such novel materials as relaxor ferroelectrics [55,56,57,58]. Much more examples of heterophase matter can be found in the review articles [59,60].

The basic difficulty in the description of heterophase systems is the necessity of dealing with two or more different phases, possessing principally different symmetry properties, but coexisting inside the same volume. The problem is aggravated by the fact that the location of the germs of different phases in space is chaotic. The situation is drastically different from the case of a sample consisting of several domains with well defined spatial locations and structure [61]. A typical heterophase system is a mixture of regions, randomly located in space and having various and often very irregular shapes. Moreover, in many cases, the heterophase regions are not static, but can move in space, vary their shapes, and even appear and disappear.

To describe such a complicated matter, it has been necessary to develop an approach allowing for the treatment of these inhomogeneous and nonequilibrium systems. More precisely, such systems are to be locally equilibrium, since the notion of phase requires the existence of at least local equilibrium. A general approach for treating heterophase systems has been advanced [62,63,64,65,66] and summarized in reviews [59,60].

One of the main problems in treating the systems with coexisting regions of different symmetry is how to separate the states, corresponding to different symmetry properties, in the space of microscopic quantum states. The standard situation is when the system as a whole is characterized by a given Hilbert space, with a prescribed symmetry. Then how would it be possible to describe subsystems of different symmetry in the frame of the same Hilbert space? This can be done, for instance, by imposing different conditions on the equations characterizing different phases [62,67] or by invoking the method of restricted trace [68]. Probably, the most powerful, convenient, and rigorous method of separating different phases is the approach based on the notion of weighted Hilbert spaces, advanced in Refs. [59,60]. In the latter publications, however, this method was not formulated in the most mathematically general way. The aim of the present paper is to further develop the method of weighted Hilbert spaces, formulating it in the most general form. As an illustration, the method will be applied to deriving a general model of a heterophase ferroelectric. The choice of ferroelectrics for illustration is caused by several reasons. First, there are numerous examples of the heterophase systems of this type [41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60], hence their correct description is of great importance. Second, the general model related to these materials is generic for many other substances, hence it serves as a good example for extending the approach to other heterophase systems.

2. Single-Phase Systems

Before developing an approach for treating multiphase systems, it is necessary to briefly recall the main points of describing the standard case of a single-phase system, which will be used in what follows. For generality, quantum systems are considered throughout the paper.

Let us have a set $\{{\phi }_n\}$ of states forming a basis for the closed linear envelope

$$\overline{\mathcal{E}}\equiv {\mathrm{Span}}_n\{{\phi }_n\}.$$

(1)

Here the index n implies a multi-index that can pertain to either discrete or continuous set. The discrete set can be infinite. Let the scalar product $<f|h>$ be defined for each pair $f,h\in \mathcal{E}$. The norm, generated by the scalar product, is

$$||f||\equiv \sqrt{\langle f|f\rangle }(f\in \overline{\mathcal{E}}).$$

(2)

Completing the linear envelope (1) by the norm yields the complete normed space, that is, the Hilbert space

$$\mathcal{H}\equiv \{\overline{\mathcal{E}},||f||\}.$$

(3)

The so-defined Hilbert space can be separable or not, depending on the physics of a concrete problem and, respectively, on the nature of the multi-index n enumerating the basis $\{{\phi }_n\}$. The existence of a Hilbert space, associated with the considered physical system, is the necessary prerequisite for characterizing the system.

The basis $\{{\phi }_n\}$ can be taken to be orthonormalized, such that

$$\langle {\phi }_m|{\phi }_n\rangle ={\delta }_{mn}.$$

(4)

And let it be complete. Then for each $f,h\in \overline{\mathcal{E}}$, one can write their expansions over the basis,

$$f=\sum _n{f}_n{\phi }_n,h=\sum _n{h}_n{\phi }_n,$$

(5)

with the expansion coefficients

$$f_n=\langle {\phi }_n|f\rangle ,h_n=\langle {\phi }_n|h\rangle .$$

Therefore, the scalar product can be represented as

$$\langle f|h\rangle =\sum _n{f}_n^{\ast }{h}_n.$$

(6)

Suppose the algebra $\mathcal{A}$ of local observables, represented by Hermitian operators $\widehat{A}\in \mathcal{A}$, be given on $\mathcal{H}$. The matrix elements of $\widehat{A}$ over the basis $\{{\phi }_n\}$ are

$$A_{mn}\equiv \langle {\phi }_m|\widehat{A}|{\phi }_n\rangle .$$

(7)

The system statistics are characterized by a statistical operator $\widehat{\rho }$ that is a trace-one operator, acting on $\mathcal{H}$. The operator averages are defined as

$$\langle \widehat{A}\rangle \equiv {\mathrm{Tr}}_{\mathcal{H}}\widehat{\rho }\widehat{A}=\sum _n\langle {\phi }_n|\widehat{\rho }\widehat{A}|{\phi }_n\rangle .$$

(8)

The set of all available averages is termed the statistical state:

$$\langle \mathcal{A}\rangle =\{\langle \widehat{A}\rangle \}.$$

(9)

Characterizing thermodynamics phases, a special role is played by the order operator $\widehat{\eta }\in \mathcal{A}$, whose average defines the system order parameter

$$\eta =\langle \widehat{\eta }\rangle .$$

(10)

In general, this can be a scalar, vector, or matrix quantity, so that one talks of an order parameter just for short. The main point is that this order parameter is specific for each thermodynamic phase, such that different phases, possessing different symmetries, enjoy different order parameters.

3. Symmetry Breaking

In the language of symmetries, the existence of phase transitions assumes the following. If there is a statistical state that is invariant with respect to a symmetry group, then it can be decomposed into a sum of several terms describing different pure phases [1,69,70]. To select a state with a particular type of symmetry, one employs one of the variants of the quasiaveraging techniques [1,67,71]. When the state with a broken symmetry is thermodynamically more stable than the invariant state, one calls this the spontaneous breaking of symmetry.

The methods of quasiaverages allow one to select a particular state with the desired symmetry only in the case of thermodynamically equilibrium systems, when all the system is characterized by one and the same symmetry. But our aim is to describe the situation, when inside the considered system there appear regions with different types of symmetry. How could we proceed in such a principally different case? For this purpose, it is necessary to develop a more general method of symmetry breaking, which could be used not only for equilibrium systems, but also for quasiequilibrium, metastable, or even for arbitrary nonequilibrium systems. Such a general method is developed below.

Let us consider a statistical system that, under different thermodynamic conditions, could be in different thermodynamic phases, enumerated by the index $\nu =1,2,\dots$. We keep in mind that the space of microscopic states, related to the considered system, is the Hilbert space (3), with a basis $\{{\phi }_n\}$, as described in Sec. 2. Let us put into correspondence to a vector ${\phi }_n$ a probability

$$p_n^{\nu }\equiv p^{\nu }({\phi }_n).$$

(11)

The set $\{p_n^{\nu }\}$ of these probabilities forms a probability measure with the standard normalization property

$$\sum _n{p}_n^{\nu }=1(0\le p_n^{\nu }\le 1).$$

(12)

Now, let us introduce the weighting operator

$${\widehat{W}}_{\nu }\equiv \sum _n{p}_n^{\nu }|{\phi }_n\rangle \langle {\phi }_n|.$$

(13)

Using this, we define the weighted scalar product

$${\langle f|h\rangle }_{\nu }\equiv \langle f|{\widehat{W}}_{\nu }|h\rangle .$$

(14)

With the weighting operator (13), we have

$${\langle f|h\rangle }_{\nu }=\sum _n{p}_n^{\nu }\langle f|{\phi }_n\rangle \langle {\phi }_n|h\rangle .$$

(15)

Under expansions (5), the latter reads as

$${\langle f|h\rangle }_{\nu }=\sum _n{p}_n^{\nu }{f}_n^{\ast }{h}_n.$$

(16)

The scalar product (15) generates the norm

$${||f||}_{\nu }\equiv \sqrt{{\langle f|f\rangle }_{\nu }}(f\in \overline{\mathcal{E}}),$$

(17)

which, in view of form (16), can be written as

$${||f||}_{\nu }\equiv \sqrt{{\sum }_n{p}_n^{\nu }{|f_n|}^2}.$$

(18)

The closed linear envelope (1), equipped with norm (17), is the weighted Hilbert space

$${\mathcal{H}}_{\nu }\equiv \{\overline{\mathcal{E}}{,||f||}_{\nu }\}.$$

(19)

The representation ${\mathcal{A}}_{\nu }$ of the algebra of local observables $\mathcal{A}$, acting on the weighted Hilbert space (19), consists of the operators ${\widehat{A}}_{\nu }$ defined through their matrix elements

$$\langle {\phi }_m|{\widehat{A}}_{\nu }|{\phi }_n\rangle \equiv \langle {\phi }_m\left|\frac{1}{2}{\left[\widehat{A},{\widehat{W}}_{\nu }\right]}_{+}\right|{\phi }_n\rangle .$$

(20)

Taking into account the weighting operator (13) gives

$$\langle {\phi }_m|{\widehat{A}}_{\nu }|{\phi }_n\rangle =\frac{1}{2}\left(p_m^{\nu }+p_n^{\nu }\right)\langle {\phi }_m|\widehat{A}|{\phi }_n\rangle .$$

(21)

The operator averages are

$$\langle {\widehat{A}}_{\nu }\rangle \equiv {\mathrm{Tr}}_{{\mathcal{H}}_{\nu }}{\left(\widehat{\rho }\widehat{A}\right)}_{\nu }=\sum _n\langle {\phi }_n|{\left(\widehat{\rho }\widehat{A}\right)}_{\nu }|{\phi }_n\rangle .$$

(22)

With the matrix elements (21), we get

$$\langle {\widehat{A}}_{\nu }\rangle =\sum _n{p}_n^{\nu }\langle {\phi }_n|\widehat{\rho }\widehat{A}|{\phi }_n\rangle .$$

(23)

Similarly, the average of the order operator ${\widehat{\eta }}_{\nu }\in {\mathcal{A}}_{\nu }$ becomes

$$\langle {\widehat{\eta }}_{\nu }\rangle =\sum _n{p}_n^{\nu }\langle {\phi }_n|\widehat{\eta }|{\phi }_n\rangle .$$

(24)

Respectively, the order parameter is

$${\eta }_{\nu }\equiv \langle {\widehat{\eta }}_{\nu }\rangle .$$

(25)

The probability measure $\{p_n^{\nu }\}$ is to be defined so that to guarantee the order parameter, specifying the chosen thermodynamics phase. This means that the probabilities (11) should select the states typical of the considered phase [72,73], which does not need to be equilibrium. By this construction, it is clear that the scheme, based on the weighted Hilbert spaces, includes as a particular case the selection of phases by means of the quasiaveraging method, since the latter also chooses the states typical of the desired phase, but provided this phase corresponds to a stable equilibrium system.

4. Multiphase Systems

First of all, let us stress the difference of the case we try to describe, as compared to the Gibbs phase mixture [74]. In the latter case, the system is spatially separated into several macroscopic regions filled by different thermodynamics phases. But in the case we are interested in, the system is a heterophase mixture, comprising the mesoscopic germs of several thermodynamic phases, which are randomly intermixed and coexist in a region of thermodynamic parameters in the same volume [59,60]. The space of states for such a heterophase system is the fiber space

$$\tilde{\mathcal{F}}\equiv \underset{\nu }{⨂}{\mathcal{H}}_{\nu }.$$

(26)

The basis of the latter, $\{{\tilde{\phi }}_{\mathbf{n}}\}$, is made of the tensor products

$${\tilde{\phi }}_{\mathbf{n}}\equiv \underset{\nu }{⨂}{\phi }_{n_{\nu }},$$

(27)

in which the notation for the multi-index

$$\mathbf{n}\equiv \{n_1,n_2,n_3,\dots \}$$

is used. Any state $\tilde{f}\in \tilde{\mathcal{F}}$ can be decomposed over this basis:

$$\tilde{f}=\sum _{\mathbf{n}}{f}_{\mathbf{n}}{\tilde{\phi }}_{\mathbf{n}}\left(f_{\mathbf{n}}=\langle {\tilde{\phi }}_{\mathbf{n}}|\tilde{f}\rangle \right).$$

Then the scalar product of any pair $\tilde{f},\tilde{h}\in \tilde{\mathcal{F}}$ is given by

$$\langle \tilde{f}|\tilde{h}\rangle =\sum _{\mathbf{n}}{f}_{\mathbf{n}}^{\ast }{h}_{\mathbf{n}}.$$

(28)

The scalar product generates the norm

$$||\tilde{f}||\equiv \sqrt{\langle \tilde{f}|\tilde{f}\rangle }.$$

(29)

The fiber space (26) can be equivalently represented as the closed linear envelope

$$\tilde{\mathcal{F}}=\left\{{\mathrm{Span}}_{\mathbf{n}}\{{\tilde{\phi }}_{\mathbf{n}}\},||\tilde{f}||\right\},$$

(30)

equipped with norm (29).

The operators of observables $\tilde{A}$, acting on $\tilde{\mathcal{F}}$, are defined as the direct sums

$$\tilde{A}={\displaystyle \underset{\nu }{⨁}}{\widehat{A}}_{\nu }.$$

(31)

The related matrix elements are given by the expressions

$$\langle {\tilde{\phi }}_{\mathbf{m}}|\tilde{A}|{\tilde{\phi }}_{\mathbf{n}}\rangle =\sum _{\nu }\langle {\phi }_{m_{\nu }}|{\widehat{A}}_{\nu }|{\phi }_{n_{\nu }}\rangle \prod _{\mu (\ne \nu )}{\delta }_{m_{\mu }{n}_{\mu }}.$$

(32)

The operator averages are defined as

$$\langle \tilde{A}\rangle \equiv {\mathrm{Tr}}_{\tilde{\mathcal{F}}}\tilde{\rho }\tilde{A}=\sum _{\mathbf{n}}\langle {\tilde{\phi }}_{\mathbf{n}}|\tilde{\rho }\tilde{A}|{\tilde{\phi }}_{\mathbf{n}}\rangle .$$

(33)

This, in view of form (27), yields

$$\langle \tilde{A}\rangle =\sum _{\nu }\sum _n\langle {\phi }_n|{(\widehat{\rho }\widehat{A})}_{\nu }|{\phi }_n\rangle .$$

(34)

Comparing this with (22), we obtain

$$\langle \tilde{A}\rangle =\sum _{\nu }\langle {\tilde{A}}_{\nu }\rangle .$$

(35)

The set $\{<\tilde{A}>\}$ of all observable averages is the statistical state of the heterophase system.

5. Phase Configurations

When the system is inhomogeneous, being composed of many mesoscopic regions of different phases, we need, first of all, to describe the spatial distribution of these regions inside the system. For this purpose, the space $\mathbb{V}$, occupied by the system, can be decomposed into the subregions, whose set $\{{\mathbb{V}}_{\nu }\}$ forms an orthogonal covering:

$$\mathbb{V}=\bigcup _{\nu }{\mathbb{V}}_{\nu },V=\sum _{\nu }{V}_{\nu },$$

(36)

such that

$${\mathbb{V}}_{\mu }\bigcap {\mathbb{V}}_{\nu }={\delta }_{\mu \nu }{\mathbb{V}}_{\nu },$$

(37)

where

$$V\equiv {\int }_{\mathbb{V}}d\mathbf{r},V_{\nu }\equiv {\int }_{{\mathbb{V}}_{\nu }}d\mathbf{r}.$$

(38)

The set of the regions, occupied by a $\nu$-phase is described by the manifold characteristic function [75], or manifold indicator

$$\begin{array}{c}\hfill {\xi }_{\nu }(\mathbf{r})\equiv \left\{\begin{array}{cc}1,\hfill & \mathbf{r}\in {\mathbb{V}}_{\nu }\hfill \\ 0,\hfill & \mathbf{r}\notin {\mathbb{V}}_{\nu }\hfill \end{array}\right..\end{array}$$

(39)

These functions satisfy the properties

$$\sum _{\nu }{\xi }_{\nu }(\mathbf{r})=1(\mathbf{r}\in \mathbb{V})$$

(40)

and

$${\int }_{\mathbb{V}}{\xi }_{\nu }(\mathbf{r})d\mathbf{r}=V_{\nu }.$$

(41)

The collection of all manifold indicators defines the phase configuration

$$\xi \equiv \{{\xi }_{\nu }(\mathbf{r}):\nu =1,2,\dots ;\mathbf{r}\in \mathbb{V}\}.$$

(42)

In turn, the space, occupied by a $\nu$-phase, can be decomposed into subregions ${\mathbb{V}}_{\nu }$, forming an orthogonal subcovering $\{{\mathbb{V}}_{\nu }\}$, such that

$${\mathbb{V}}_{\nu }=\bigcup _{i=1}^{z_{\nu }}{\mathbb{V}}_{\nu i},$$

(43)

and

$${\mathbb{V}}_{\mu i}\bigcap {\mathbb{V}}_{\nu j}={\delta }_{\mu \nu }{\delta }_{ij}{\mathbb{V}}_{\nu i}.$$

(44)

Each subregion ${\mathbb{V}}_{\nu i}$ is described by its manifold indicator

$$\begin{array}{c}\hfill {\xi }_{\nu i}(\mathbf{r}-{\mathbf{a}}_{\nu i})\equiv \left\{\begin{array}{cc}1,\hfill & \mathbf{r}\in {\mathbb{V}}_{\nu i}\hfill \\ 0,\hfill & \mathbf{r}\notin {\mathbb{V}}_{\nu i}\hfill \end{array}\right.,\end{array}$$

(45)

in which ${\mathbf{a}}_{\nu \mathbf{i}}$ is a fixed vector pertaining to ${\mathbb{V}}_{\nu i}$. Then the manifold indicator (39) writes as

$${\xi }_{\nu }(\mathbf{r})=\sum _{i=1}^{z_{\nu }}{\xi }_{\nu i}(\mathbf{r}-{\mathbf{a}}_{\nu i}).$$

(46)

Under a given configuration, the relative volume, occupied by a $\nu$-phase, is characterized by its geometrical fraction

$$x_{\nu }\equiv \frac{1}{V}{\int }_{\mathbb{V}}{\xi }_{\nu }(\mathbf{r})d\mathbf{r}=\frac{V_{\nu }}{V}.$$

(47)

As is clear, the latter satisfy the normalization condition

$$\sum _{\nu }{x}_{\nu }=1(0\le x_{\nu }\le 1).$$

(48)

This shows that the set $\{x_{\nu }\}$ of all admissible geometrical fractions, with the given normalization condition, can be regarded as a probability measure.

Since the the phase regions are randomly distributed inside the considered system, the configuration (42) is to be treated as a random variable. The related differential measure can be defined as

$$\mathcal{D}\xi =\underset{\{z_{\nu }\to \infty \}}{\lim }\delta \left(\sum _{\nu }{x}_{\nu }-1\right)\prod _{\nu }d{x}_{\nu }\prod _{\nu }\prod _{i=1}^{z_{\nu }}\frac{d{\mathbf{a}}_{\nu i}}{V}.$$

(49)

All possible configurations constitute a topological space

$$\mathcal{X}\equiv \{\xi ,\mathcal{D}\xi \}.$$

(50)

The above constructions describe the situation, when the system is separated into several regions filled by different phases, and these regions are randomly distributed in space.

6. Configuration Averaging

An inhomogeneous heterophase system is not necessarily in complete equilibrium. But it must be at least quasi-equilibrium in order that it would be admissible to talk about the germs of phases. As is known, the notion of phases is not strictly defined for finite systems. The mathematically rigorous definition of phases assumes the introduction of thermodynamic limit [2,3,72,74,76]. However, in practice, it is possible to speak about germs of phases already when each of such germs consists of a large number of particles $N\gg 1$. Monte Carlo simulations show that thermodynamics phases can be well defined already for $10-100$ particles in a finite cluster [77]. Respectively, though the symmetry, related to a phase, may be not strictly defined for a finite cluster, but it is possible to talk about an asymptotic symmetry that is approximately defined for a large number of particles $N\gg 1$, keeping in mind that the symmetry becomes exact in the thermodynamic limit.

A quasi-equilibrium system is described by a quasi-Hamiltonian [59] depending on a given phase configuration and having the operator structure as in (31):

$$\widehat{Q}(\xi )={\displaystyle \underset{\nu }{⨁}}{\widehat{Q}}_{\nu }({\xi }_{\nu }).$$

(51)

The partition function

$$Z\equiv {\mathrm{Tr}}_{\tilde{\mathcal{F}}}\int \exp \left\{-\widehat{Q}(\xi )\right\}\mathcal{D}\xi$$

(52)

includes the quantum averaging over the given quantum variables and the averaging over phase configurations. This function defines the quasi-equilibrium thermodynamic potential

$$y\equiv -\frac{1}{N}\ln Z.$$

(53)

Considering asymptotically large systems, we, as usual, keep in mind the thermodynamic limit

$$N\to \infty ,V\to \infty ,\frac{N}{V}\to const.$$

(54)

Similarly to (31), the operators of observables, under the given phase configuration, have the form

$$\widehat{A}(\xi )={\displaystyle \underset{\nu }{⨁}}{\widehat{A}}_{\nu }({\xi }_{\nu }).$$

(55)

And for what follows, it is convenient to introduce the notation

$${\widehat{A}}_{\nu }(x_{\nu })\equiv \underset{{\xi }_{\nu }\to x_{\nu }}{\lim }{\widehat{A}}_{\nu }({\xi }_{\nu }).$$

(56)

The statistical operator of a multi-phase system, with a fixed phase configuration, is

$$\widehat{\rho }(\xi )=\frac{1}{Z}\exp \left\{-\widehat{Q}(\xi )\right\}.$$

(57)

The observable quantities, related to the operators (55), are given by the averages

$$\langle \tilde{A}\rangle \equiv {\mathrm{Tr}}_{\tilde{\mathcal{F}}}\int \widehat{\rho }(\xi )\widehat{A}(\xi )\mathcal{D}\xi .$$

(58)

The following theorem is valid [59,78,79].

Theorem 1.

If ${\widehat{Q}}_{\nu }({\xi }_{\nu })$ can be represented as an expansion in powers of ${\xi }_{\nu }$, then the thermodynamic potential (53), in the thermodynamic limit (54), asymptotically equals

$$y=\mathrm{abs}\underset{\mathbf{w}}{\min }y(\mathbf{w}),$$

(59)

with the set

$$\mathbf{w}\equiv \{w_1,w_2,\dots \}$$

(60)

forming the probability measure enjoying the standard properties

$$\sum _{\nu }{w}_{\nu }=1(0\le w_{\nu }\le 1),$$

(61)

and where

$$y(\mathbf{w})=\sum _{\nu }{y}_{\nu }(w_{\nu })$$

(62)

is the sum of the terms

$$y_{\nu }(w_{\nu })=-\frac{1}{N}\ln \mathrm{Tr}{Z}_{\nu },$$

(63)

in which

$$Z_{\nu }={\mathrm{Tr}}_{{\mathcal{H}}_{\nu }}\exp \left\{-{\widehat{Q}}_{\nu }(w_{\nu })\right\}.$$

(64)

The quantities $w_{\nu }$ are the phase geometric probabilities, showing the corresponding weights of the coexisting thermodynamic phases.

Corollary 1

From this theorem, it follows that the thermodynamic potential (59) can also be represented in the form

$$y=-\frac{1}{N}\ln \mathrm{Tr}\tilde{Z},$$

(65)

where

$$\tilde{Z}\equiv \prod _{\nu }{Z}_{\nu }.$$

The above theorem defines the thermodynamic potential of a heterophase system. The observable quantities, corresponding to the averages of the operators from the algebra of local observables, are given by the following theorem [59,78,79].

Theorem 2. 

Assume that ${\widehat{Q}}_{\nu }({\xi }_{\nu })$ and ${\widehat{A}}_{\nu }({\xi }_{\nu })$ can be expanded in powers of ${\xi }_{\nu }$, then the averages (58), for asymptotically large N, take the form of the sum

$$\langle \tilde{A}\rangle =\sum _{\nu }\langle {\widehat{A}}_{\nu }\rangle ,$$

(66)

with the terms

$$\langle {\widehat{A}}_{\nu }\rangle ={\mathrm{Tr}}_{{\mathcal{H}}_{\nu }}{\widehat{\rho }}_{\nu }{\widehat{A}}_{\nu },$$

(67)

in which

$${\widehat{\rho }}_{\nu }=\frac{1}{Z_{\nu }}\exp \left\{-{\widehat{Q}}_{\nu }(w_{\nu })\right\},$$

(68)

and where the notation ${\widehat{A}}_{\nu }\equiv {\widehat{A}}_{\nu }(w_{\nu })$ is used.

Corollary 2. 

Equivalently, the averages (66) can be represented in another way by introducing the operators

$$\tilde{A}\equiv {\displaystyle \underset{\nu }{⨁}}{\widehat{A}}_{\nu },$$

(69)

for which

$$\langle \tilde{A}\rangle ={\mathrm{Tr}}_{\tilde{\mathcal{F}}}\tilde{\rho }\tilde{A},$$

(70)

where

$$\tilde{\rho }=\underset{\nu }{⨂}{\widehat{\rho }}_{\nu }.$$

The second theorem defines the method of calculating the averages for a heterophase system.

7. Effective Hamiltonians

The quasi-Hamiltonian ${\widehat{Q}}_{\nu }({\xi }_{\nu })$ is connected with the local Hamiltonian ${\widehat{H}}_{\nu }(\mathbf{r},{\xi }_{\nu })$ through the relation

$${\widehat{Q}}_{\nu }({\xi }_{\nu })=\int {\beta }_{\nu }(\mathbf{r},{\xi }_{\nu }){\widehat{H}}_{\nu }(\mathbf{r},{\xi }_{\nu })d\mathbf{r},$$

(71)

in which the inverse local temperatures ${\beta }_{\nu }(\mathbf{r},{\xi }_{\nu })$ play the role of the Lagrange multipliers.

If the considered system is such that the phase fluctuations can appear randomly in any part of the system, the latter is called phase-uniform on average [59]. For such a system, the configuration averaging of the inverse temperatures gives

$${\beta }_{\nu }\equiv \int {\beta }_{\nu }(\mathbf{r},{\xi }_{\nu })\mathcal{D}\xi .$$

(72)

The mixture of different phases is in thermal equilibrium, when the temperatures of these phases coincide:

$${\beta }_{\nu }=\beta =\frac{1}{T}(\forall \nu ).$$

(73)

Then the quasi-Hamiltonian ${\widehat{Q}}_{\nu }(w_{\nu })$, entering the partition function (64), becomes

$${\widehat{Q}}_{\nu }(w_{\nu })=\beta {\widehat{H}}_{\nu },$$

(74)

with the renormalized Hamiltonian

$${\widehat{H}}_{\nu }\equiv \int {\widehat{H}}_{\nu }(\mathbf{r},w_{\nu })d\mathbf{r}.$$

(75)

The partial statistical operators (68) acquire the form

$${\widehat{\rho }}_{\nu }=\frac{1}{Z_{\nu }}\exp \left(-\beta {\widehat{H}}_{\nu }\right),$$

(76)

with the partition functions

$$Z_{\nu }={\mathrm{Tr}}_{{\mathcal{H}}_{\nu }}\exp \left(-\beta {\widehat{H}}_{\nu }\right).$$

(77)

It is convenient to introduce the effective Hamiltonian

$$\tilde{H}\equiv {\displaystyle \underset{\nu }{⨁}}{\widehat{H}}_{\nu },$$

(78)

using which, the thermodynamic potential (65) can be represented as

$$y=-\frac{1}{N}\ln \mathrm{Tr}{e}^{-\beta \tilde{H}}.$$

(79)

The terms of sum (78) are called the phase-replica Hamiltonians, since they have a similar mathematical structure, but are associated with different phases.

In order to connect the thermodynamic potential (79) with the free energy, let us recall that the latter is defined as

$$F\equiv -T\ln \mathrm{Tr}{e}^{-\beta \tilde{H}}.$$

(80)

Therefore, the potential y is directly connected to the free energy by means of the relations

$$F=NTy=\sum _{\nu }{F}_{\nu }$$

and

$$F_{\nu }=-T\ln \mathrm{Tr}{Z}_{\nu }=NT{y}_{\nu }.$$

To specify the consideration, let us take the local Hamiltonians in the usual form

$${\widehat{H}}_{\nu }(\mathbf{r},{\xi }_{\nu })={\xi }_{\nu }(\mathbf{r}){\psi }_{\nu }^{†}(\mathbf{r})\left(-\frac{{\nabla }^2}{2m}+U\right){\psi }_{\nu }(\mathbf{r})+$$

$$+\frac{1}{2}\int {\xi }_{\nu }(\mathbf{r}){\xi }_{\nu }({\mathbf{r}}^{\prime }){\psi }_{\nu }^{†}(\mathbf{r}){\psi }_{\nu }^{†}({\mathbf{r}}^{\prime })\mathsf{\Phi }(\mathbf{r}-{\mathbf{r}}^{\prime }){\psi }_{\nu }({\mathbf{r}}^{\prime }){\psi }_{\nu }(\mathbf{r})d\mathbf{r}.$$

(81)

Here $U=U(\mathbf{r})$ is an external potential and ${\psi }_{\nu }(\mathbf{r})$ are the field operators of the particles forming the system.

Then the phase-replica Hamiltonians are

$${\widehat{H}}_{\nu }=w_{\nu }\int {\psi }_{\nu }^{†}(\mathbf{r})\left(-\frac{{\nabla }^2}{2m}+U\right){\psi }_{\nu }(\mathbf{r})d\mathbf{r}+$$

$$+\frac{w_{\nu }^2}{2}\int {\psi }_{\nu }^{†}(\mathbf{r}){\psi }_{\nu }^{†}({\mathbf{r}}^{\prime })\mathsf{\Phi }(\mathbf{r}-{\mathbf{r}}^{\prime }){\psi }_{\nu }({\mathbf{r}}^{\prime }){\psi }_{\nu }(\mathbf{r})d\mathbf{r}d{\mathbf{r}}^{\prime }.$$

(82)

Introducing the notation for the kinetic-energy operator

$${\widehat{K}}_{\nu }\equiv \int {\psi }_{\nu }^{†}(\mathbf{r})\left(-\frac{{\nabla }^2}{2m}+U\right){\psi }_{\nu }(\mathbf{r})d\mathbf{r}$$

(83)

and for the operator

$${\widehat{\mathsf{\Phi }}}_{\nu }\equiv \int {\psi }_{\nu }^{†}(\mathbf{r}){\psi }_{\nu }^{†}({\mathbf{r}}^{\prime })\mathsf{\Phi }(\mathbf{r}-{\mathbf{r}}^{\prime }){\psi }_{\nu }({\mathbf{r}}^{\prime }){\psi }_{\nu }(\mathbf{r})d\mathbf{r}d{\mathbf{r}}^{\prime },$$

(84)

related to the potential energy part of the Hamiltonian (82), we have for the latter

$${\widehat{H}}_{\nu }=w_{\nu }{\widehat{K}}_{\nu }+\frac{w_{\nu }^2}{2}{\widehat{\mathsf{\Phi }}}_{\nu }.$$

(85)

In this way, we have derived all basic equations for treating heterophase systems. The derivation has been based on the following three major points making it possible to separate different thermodynamic phases. First, to distinguish the phases in the space of microscopic states, the notion of weighted Hilbert spaces is introduced. Second, to separate the phases in real space, the manifold indicators were employed. And, finally, the procedure of averaging over phase configurations is accomplished, leading to the set of equations for equilibrium on average phase replicas. The idea of the averaging procedure reminds the method of averaging [80] and the scale separation approach [81,82,83], used for nonlinear equations in dynamical theory. The main difference from the latter is that here we have averaged out slow heterophase fluctuations, slow with respect to the fast microscopic motion of particles, while in dynamical theory one usually averages out fast fluctuations, leaving at the end the slow motion of guiding centers.

8. Stability Conditions

The developed theory should be complimented by an important addition discussing the stability of heterophase systems. Thermodynamic stability is characterized by the minimization of a thermodynamic potential. For concreteness, let us consider the case of a heterophase system, where two phases coexist, so that $\nu =1,2$. And let us denote

$$w_1\equiv w,w_2=1-w.$$

(86)

Minimizing the thermodynamic potential (79) with respect to w implies

$$\frac{\partial y}{\partial w}=0,\frac{{\partial }^{2}y}{\partial w^2}>0.$$

(87)

The first of these equations gives

$$\langle \frac{\partial \tilde{H}}{\partial w}\rangle =0.$$

(88)

And the second condition leads to the inequality

$$\langle \frac{{\partial }^2\tilde{H}}{\partial w^2}\rangle >\beta \langle {\left(\frac{\partial \tilde{H}}{\partial w}\right)}^2\rangle .$$

(89)

Since the right-hand side in the above inequality is nonnegative, the necessary condition of heterophase stability is

$$\langle \frac{{\partial }^2\tilde{H}}{\partial w^2}\rangle >0.$$

(90)

To specify these conditions, let us take the effective Hamiltonian (78), with the replica Hamiltonians (85). And let us use the notations for the averages

$$K_{\nu }\equiv \langle {\widehat{K}}_{\nu }\rangle ,{\mathsf{\Phi }}_{\nu }\equiv \langle {\widehat{\mathsf{\Phi }}}_{\nu }\rangle .$$

(91)

Then we have

$$\frac{1}{N}\langle \frac{\partial \tilde{H}}{\partial w}\rangle =K_1+w{\mathsf{\Phi }}_1-K_2-(1-w){\mathsf{\Phi }}_2$$

and

$$\frac{1}{N}\langle \frac{{\partial }^2\tilde{H}}{\partial w^2}\rangle ={\mathsf{\Phi }}_1+{\mathsf{\Phi }}_2.$$

This yields the equation for the phase probability

$$w=\frac{{\mathsf{\Phi }}_2+K_2-K_1}{{\mathsf{\Phi }}_1+{\mathsf{\Phi }}_2}.$$

(92)

The stability condition (89) gives

$${\mathsf{\Phi }}_1+{\mathsf{\Phi }}_2>\frac{\beta }{N}\langle {\left(\frac{\partial \tilde{H}}{\partial w}\right)}^2\rangle ,$$

(93)

and from the stability condition (90), we get

$${\mathsf{\Phi }}_1+{\mathsf{\Phi }}_2>0.$$

(94)

One more condition, follows from the definition of $w_{\nu }$ as of phase probabilities, according to which

$$0\le w\le 1.$$

This, with the use of (92), results in the inequalities

$$-{\mathsf{\Phi }}_1\le K_1-K_2\le {\mathsf{\Phi }}_2.$$

(95)

The stability conditions must be valid in order that the considered heterophase system could really exist. Note that the stability condition (90) is analogous to the condition of diffusion stability [84].

9. Heterophase Ferroelectrics

To illustrate how the developed theory works, let us consider ferroelectric materials. Ferroelectrics are known to be a good example of matter demonstrating heterophase properties. In many ferroelectrics, above the transition point from the disordered into the ordered state, there exist small polarized clusters [37]. These fluctuational embryos of the ordered phase inside the disordered phase where termed by Cook [36] antiphase fluctuations, emphasizing that they where a particular case of heterophase fluctuations, whose existence was predicted by Frenkel [85]. Such embryos of one phase inside another phase can exist in a whole region around the phase-transition point $T_c$. At temperatures below $T_c$, the embryos of the paraelectric phase inside the ferroelectric phase arise at a temperature $T_n$ called the lower nucleation point [86]. And above $T_c$, there is another temperature $T_n^{\ast }$, called the upper nucleation point [86], where the embryos of the ferroelectric phase appear inside the paraelectric phase. Thus, around the phase transition temperature, there can exist a region of temperatures

$$T_n<T_c<T_n^{\ast },$$

where the embryos of two phases coexist. The phase transition, generically, can be either of first or of second order, but the appearance of heterophase fluctuations, usually smears it into a continuous crossover [41,42,43,44].

Heterophase fluctuations were observed, for instance, in such well known ferroelectrics as HCl, DCl, mixed crystals HCl${}_{1-x}$DCl${}_x$, and RbCaF${}_3$, where they were intensively studied by Brookeman and Rigamonti using nuclear quadrupole resonance [42,43] and nuclear magnetic resonance [44]. Heterophase fluctuations in these ferroelectrics arise in a finite region around the phase transition point. The appearance of these fluctuations occurs even without external defects, although the presence of defects intensifies their nucleation [41].

Such heterophase fluctuations have also been observed in many other ferroelectrics, e.g., in C${}_4$O${}_4$H${}_2$ [41], in KH${}_2$As${}_4$ [46], in Rb${}_x$(ND${}_4$)${}_{1-x}$D${}_2$PO${}_4$ [47], in Na${}_x$Bi${}_{1-x}$TiO${}_3$ [49], and others. They also arise in such novel materials as relaxor ferroelectrics, for example, in (PbZn${}_x$Nb${}_{1-x}$O${}_3$) [55], in Ba${}_2$NdTi${}_2$Nb${}_3$O${}_{15}$ and Ba${}_2$La${}_{0.5}$Nd${}_{0.5}$Ti${}_2$Nb${}_3$O${}_{15}$ [56,58], and, with a high probability, in many other relaxors, such as PbMg${}_{1/3}$Nb${}_{2/3}$O${}_3$, (Ba${}_x$Pb${}_{1-x}$)(Zn${}_{1/3}$Nb${}_{2/3}$)O${}_3$, (Sr${}_x$Pb${}_{1-x}$)(Zn${}_{1/3}$Nb${}_{2/3}$)O${}_3$, (Ba${}_x$Pb${}_{1-x}$)(Yb${}_{0.5}$Nb${}_{0.5}$)O${}_3$, Pb${}_{1-x}$Ba${}_x$(Yb${}_{1/2}$Ta${}_{1/2}$)O${}_3$, (Fe${}_{1/2}$Nb${}_{1/2}$)O${}_3$, Pb(Fe${}_{1/2}$Ta${}_{1/2}$)O${}_3$, Pb(Yb${}_{1/2}$Ta${}_{1/2}$)O${}_3$-Pb(Fe${}_{1/2}$Ta${}_{1/2}$)O${}_3$, and Pb(Mg${}_{1/3}$Nb${}_{2/3}$)O${}_3$-PbTiO${}_3$, [58].

A simple model of a heterophase ferroelectric has been considered in Refs. [50,87,88]. In this model, phonon degrees of freedom were not taken into account. However, the latter are important because of two reasons. First, the appearance of heterophase fluctuations is frequently accompanied by heterostructural fluctuations [89,90], which are connected with phonon excitations. Second, the occurrence of the paraelectric-ferroelectric phase transition is intimately related with phonon characteristics that can be directly measured. In the present paper, we derive a generalized model of a heterophase ferroelectric, taking into account the phonon degrees of freedom. This makes it possible to find out the influence of the heterophase fluctuations on such observable quantities as the Debye-Waller factor and sound velocity. This influence is especially pronounced in the phase-transition region.

10. Basic Hamiltonian

Let us consider a ferroelectric, in which there can arise the embryos of the competing phase. So that the sample can house two coexisting phases, ferroelectric and paraelectric, being randomly intermixed with each other. The ferroelectric phase will be indexed by $\nu =1$ and the paraelectric phase will be labeled by the index $\nu =2$. For concreteness, we shall study the model of a ferroelectric, in which the order is characterized by the pseudospin operator $S_j^z$ describing the shift of a charged particle into one of the wells of a double-well potential at the cite j of the crystalline lattice [91,92]. Then the ordered ferroelectric phase corresponds to the nonzero order parameter

$$\langle S_{j1}^z\rangle \ne 0.$$

(96)

On the contrary, the disordered paraelectric phase is characterized by the zero order parameter

$$\langle S_{j2}^z\rangle \equiv 0.$$

(97)

Starting with the ferroelectric Hamiltonian, having the mathematical structure characterized by the pseudospin variables [91,92], we follow the general scheme described above, and after averaging over the random phase configurations, we come to the effective Hamiltonian (78) with the replica Hamiltonians

$${\widehat{H}}_{\nu }=w_{\nu }\sum _j\frac{{\mathbf{p}}_j^2}{2m}+\frac{w_{\nu }^2}{2}\sum _{i\ne j}A({\mathbf{r}}_{ij})-w_{\nu }\mathsf{\Omega }\sum _j{S}_{j\nu }^x+$$

$$+w_{\nu }^2\sum _{i\ne j}B({\mathbf{r}}_{ij})S_{i\nu }^x{S}_{j\nu }^x-w_{\nu }^2\sum _{i\ne j}I({\mathbf{r}}_{ij})S_{i\nu }^z{S}_{j\nu }^z.$$

(98)

Here, the first term represents kinetic energy, $A(\mathbf{r}),\mathbf{B}(\mathbf{r})$, and $I(\mathbf{r})$ are particle interactions, $\mathsf{\Omega }$ is the tunneling frequency, and the abbreviated notation

$${\mathbf{r}}_{ij}\equiv {\mathbf{r}}_i-{\mathbf{r}}_j.$$

(99)

is employed.

The phonon variables can be introduced in the standard way by defining the deviation ${\mathbf{u}}_{\mathbf{j}}$ from the lattice site with the lattice vector ${\mathbf{a}}_{\mathbf{j}}$ as

$${\mathbf{r}}_j={\mathbf{a}}_j+{\mathbf{u}}_j.$$

(100)

The lattice vectors are assumed to form an equilibrium lattice, being defined as the averages

$${\mathbf{a}}_j\equiv \langle {\mathbf{r}}_j\rangle .$$

(101)

Hence the average deviation, by definition, is zero,

$$\langle {\mathbf{u}}_j\rangle =0.$$

(102)

The interactions, as usual, are supposed to be symmetric with respect to the spatial coordinate inversion:

$$A(-{\mathbf{r}}_{ij})=A({\mathbf{r}}_{ij}),B(-{\mathbf{r}}_{ij})=B({\mathbf{r}}_{ij}),I(-{\mathbf{r}}_{ij})=I({\mathbf{r}}_{ij}).$$

(103)

In what follows, we shall also use the short-hand notation for the vector differences

$${\mathbf{a}}_{ij}\equiv {\mathbf{a}}_i-{\mathbf{a}}_j.$$

(104)

and

$${\mathbf{u}}_{ij}\equiv {\mathbf{u}}_i-{\mathbf{u}}_j.$$

(105)

The interactions are expanded in powers of the deviations as

$$A({\mathbf{r}}_{ij})\cong A_{ij}+\sum _{\alpha }{A}_{ij}^{\alpha }{u}_{ij}^{\alpha }+\frac{1}{2}\sum _{\alpha \beta }{A}_{ij}^{\alpha \beta }{u}_{ij}^{\alpha }{u}_{ij}^{\beta },$$

(106)

where

$$A_{ij}\equiv A({\mathbf{a}}_{ij}),A_{ij}^{\alpha }\equiv \frac{\partial A_{ij}}{\partial a_i^{\alpha }},A_{ij}^{\alpha \beta }\equiv \frac{{\partial }^2{A}_{ij}}{\partial a_i^{\alpha }\partial a_i^{\beta }}.$$

(107)

The same expansions are made for $B(\mathbf{r})$ and $I(\mathbf{r})$. Owing to the symmetry properties (103), we have

$$A_{ij}\equiv A_{ji},A_{ij}^{\alpha }=-A_{ji}^{\alpha },A_{ij}^{\alpha \beta }=A_{ij}^{\beta \alpha }=A_{ji}^{\beta \alpha }=A_{ji}^{\alpha \beta }.$$

(108)

In the Hamiltonian (98), the double summation over the lattice excludes the self-action terms with $i=j$. In order to simplify the notation, we can sum over all lattice sites, setting the diagonal elements

$$A_{ii}=B_{ii}=I_{ii}\equiv 0.$$

(109)

The lattice is treated as ideal, because of which

$$A\equiv \sum _j{A}_{ij}=const$$

(110)

does not depend on the index i. Using the ideality of the lattice, we get

$$\sum _j{A}_{ij}^{\alpha }=\frac{\partial A}{\partial a_i^{\alpha }}=0.$$

(111)

and

$$\sum _j{A}_{ij}^{\alpha \beta }=\frac{{\partial }^{2}A}{\partial a_i^{\alpha }\partial a_i^{\beta }}=0.$$

(112)

Consequently,

$$\sum _{ij}{A}_{ij}^{\alpha }{u}_{ij}^{\alpha }=0.$$

(113)

The same type of expansions is accomplished for both the ferroelectric and paraelectric phases. Therefore, in what follows, we shall consider the resulting transformations only for the ferroelectric phase, keeping in mind that the same is done for the paraelectric phase. And to simplify the notation, we shall not write explicitly the index $\nu =1$. Then, substituting the above expansions into Hamiltonian (98), with $\nu =1$, and invoking the notation

$$S_{ij}^{\alpha }\equiv S_i^{\alpha }{S}_j^{\alpha },$$

(114)

we obtain the Hamiltonian

$${\widehat{H}}_1=w\sum _j\frac{{\mathbf{p}}_j^2}{2m}+\frac{w^2}{2}NA+\frac{w^2}{4}\sum _{ij}\sum _{\alpha \beta }{A}_{ij}^{\alpha \beta }{u}_{ij}^{\alpha }{u}_{ij}^{\beta }-w\mathsf{\Omega }\sum _j{S}_j^x+$$

$$+w^2\sum _{ij}\left(B_{ij}+\sum _{\alpha }{B}_{ij}^{\alpha }{u}_{ij}^{\alpha }+\frac{1}{2}\sum _{\alpha \beta }{B}_{ij}^{\alpha \beta }{u}_{ij}^{\alpha }{u}_{ij}^{\beta }\right)S_{ij}^x-$$

$$-w^2\sum _{ij}\left(I_{ij}+\sum _{\alpha }{I}_{ij}^{\alpha }{u}_{ij}^{\alpha }+\frac{1}{2}\sum _{\alpha \beta }{I}_{ij}^{\alpha \beta }{u}_{ij}^{\alpha }{u}_{ij}^{\beta }\right)S_{ij}^z,$$

(115)

in which $w\equiv w_1$.

11. Pseudospin-Phonon Decoupling

The obtained Hamiltonian is yet too complicated to be treated, and some approximation is required. For any approximation, we have to keep in mind that

$$\langle u_{ij}^{\alpha }\rangle =\langle u_i^{\alpha }\rangle -\langle u_j^{\alpha }\rangle =0.$$

(116)

Because the phonon and pseudospin operators are of different nature, it is reasonable to decouple them in the second-order, with respect to the deviations, terms as

$$u_{ij}^{\alpha }{u}_{ij}^{\beta }{S}_{ij}^{\gamma }=\langle u_{ij}^{\alpha }{u}_{ij}^{\beta }\rangle S_{ij}^{\gamma }+u_{ij}^{\alpha }{u}_{ij}^{\beta }\langle S_{ij}^{\gamma }\rangle -\langle u_{ij}^{\alpha }{u}_{ij}^{\beta }\rangle \langle S_{ij}^{\gamma }\rangle .$$

(117)

At the same time, the terms linear in the deviations can be left for a while, since later they can be dealt with by using canonical transformations.

Using again the lattice ideality, we see that

$$\sum _{ij}{A}_{ij}^{\alpha \beta }{u}_{ij}^{\alpha }{u}_{ij}^{\beta }=2\sum _{ij}{A}_{ij}^{\alpha \beta }{u}_i^{\alpha }{u}_j^{\beta }.$$

(118)

The average $<S_{ij}>$ depends on the difference ${\mathbf{a}}_{\mathbf{ij}}$. Thence

$$\sum _j{B}_{ij}^{\alpha \beta }\langle S_{ij}^x\rangle =\frac{{\partial }^2}{\partial a_i^{\alpha }\partial a_i^{\beta }}\sum _j{B}_{ij}\langle S_{ij}^x\rangle =0,$$

$$\sum _j{I}_{ij}^{\alpha \beta }\langle S_{ij}^z\rangle =\frac{{\partial }^2}{\partial a_i^{\alpha }\partial a_i^{\beta }}\sum _j{I}_{ij}\langle S_{ij}^z\rangle =0.$$

(119)

Similarly,

$$\sum _{ij}{B}_{ij}^{\alpha \beta }{u}_{ij}^{\alpha }{u}_{ij}^{\beta }\langle S_{ij}^x\rangle =2\sum _{ij}{B}_{ij}^{\alpha \beta }\langle S_{ij}^x\rangle u_i^{\alpha }{u}_j^{\beta },$$

$$\sum _{ij}{I}_{ij}^{\alpha \beta }{u}_{ij}^{\alpha }{u}_{ij}^{\beta }\langle S_{ij}^z\rangle =2\sum _{ij}{I}_{ij}^{\alpha \beta }\langle S_{ij}^z\rangle u_i^{\alpha }{u}_j^{\beta }.$$

(120)

Let us introduce the notation for the dynamical matrix

$${\mathsf{\Phi }}_{ij}^{\alpha \beta }\equiv A_{ij}^{\alpha \beta }+2B_{ij}^{\alpha \beta }\langle S_{ij}^x\rangle -2I_{ij}^{\alpha \beta }\langle S_{ij}^z\rangle$$

(121)

and for the renormalized interactions

$${\tilde{B}}_{ij}\equiv B_{ij}+\sum _{\alpha \beta }{B}_{ij}^{\alpha \beta }\langle u_i^{\alpha }{u}_j^{\beta }\rangle .$$

(122)

and

$${\tilde{I}}_{ij}\equiv I_{ij}+\sum _{\alpha \beta }{I}_{ij}^{\alpha \beta }\langle u_i^{\alpha }{u}_j^{\beta }\rangle .$$

(123)

Separating the non-operator energy part

$$E_1\equiv \frac{w^2}{2}NA+w^2\sum _{ij}\sum _{\alpha \beta }\left(I_{ij}^{\alpha \beta }\langle S_{ij}^z\rangle -B_{ij}^{\alpha \beta }\langle S_{ij}^x\rangle \right)\langle u_i^{\alpha }{u}_j^{\beta }\rangle ,$$

(124)

we reduce Hamiltonian (115) to the sum

$${\widehat{H}}_1=E_1+{\widehat{H}}_{ph}+{\widehat{H}}_{ps}+{\widehat{H}}_{lin}.$$

(125)

The second term here is the effective phonon Hamiltonian

$${\widehat{H}}_{ph}=w\sum _j\frac{{\mathbf{p}}_j^2}{2m}+\frac{w^2}{2}\sum _{ij}\sum _{\alpha \beta }{\mathsf{\Phi }}_{ij}^{\alpha \beta }{u}_i^{\alpha }{u}_j^{\beta }.$$

(126)

The third term is the effective pseudospin Hamiltonian

$${\widehat{H}}_{ps}=-w\mathsf{\Omega }\sum _j{S}_j^x+w^2\sum _{ij}\left({\tilde{B}}_{ij}{S}_{ij}^x-{\tilde{I}}_{ij}{S}_{ij}^z\right).$$

(127)

And the last term is the linear pseudospin-phonon interaction Hamiltonian. The latter is obtained by invoking the properties

$$\sum _{ij}{B}_{ij}^{\alpha }{u}_{ij}^{\alpha }{S}_{ij}^x=2\sum _{ij}{B}_{ij}^{\alpha }{S}_{ij}^x{u}_i^{\alpha }=-2\sum _{ij}{B}_{ij}^{\alpha }{S}_{ij}^x{u}_j^{\alpha },$$

$$\sum _{ij}{I}_{ij}^{\alpha }{u}_{ij}^{\alpha }{S}_{ij}^z=2\sum _{ij}{I}_{ij}^{\alpha }{S}_{ij}^z{u}_i^{\alpha }=-2\sum _{ij}{I}_{ij}^{\alpha }{S}_{ij}^z{u}_j^{\alpha }$$

(128)

and denoting

$$K_{ij}^{\alpha }\equiv B_{ij}^{\alpha }{S}_{ij}^x-I_{ij}^{\alpha }{S}_{ij}^z,$$

(129)

which yields

$${\widehat{H}}_{lin}=-2w^2\sum _{ij}\sum _{\alpha }{K}_{ij}^{\alpha }{u}_j^{\alpha }.$$

(130)

Thus, the Hamiltonian parameters are renormalized due to the interactions between pseudospins and phonons.

12. Dressed Phonons

The quantization of the phonon degrees of freedom can be done in the way, similar to how it is done for pure crystalline phases. However, here we have to be careful, taking into account the presence of the factors $w=w_1$, characterizing the weight of the related phase. Thus, the eigenproblem equation for the phonon frequencies and polarization vectors takes the form

$$\frac{w}{m}\sum _{\alpha \beta }{\mathsf{\Phi }}_{ij}^{\alpha \beta }\exp \left(i\mathbf{k}·{\mathbf{a}}_{ij}\right)e_{ks}^{\beta }={\omega }_{ks}^2{e}_{ks}^{\alpha }.$$

(131)

The phonon frequencies and polarization vectors can be chosen to be symmetric with respect to the momentum inversion,

$${\omega }_{-ks}={\omega }_{ks},{\mathbf{e}}_{-ks}={\mathbf{e}}_{ks}.$$

The polarization vectors are orthonormalized, such that

$${\mathbf{e}}_{ks}·{\mathbf{e}}_{k{s}^{\prime }}={\delta }_{s{s}^{\prime }},\sum _s{e}_{ks}^{\alpha }{e}_{ks}^{\beta }={\delta }_{\alpha \beta },$$

with the momentum summation over the Brillouin zone. The eigenproblem (131) can be rewritten as

$$\frac{w}{m}\sum _{\beta }{\mathsf{\Phi }}_k^{\alpha \beta }{e}_{ks}^{\beta }={\omega }_{ks}^2{e}_{ks}^{\alpha },$$

(132)

where

$${\mathsf{\Phi }}_k^{\alpha \beta }\equiv \sum _j{\mathsf{\Phi }}_{ij}^{\alpha \beta }{e}^{i\mathbf{k}·{\mathbf{a}}_{ij}}.$$

The phonon quantization, in the presence of the linear pseudospin-phonon interactions, differs from the standard case by the necessity to involve a nonuniform canonical transformation

$${\mathbf{p}}_j=-\frac{i}{\sqrt{N}}\sum _{ks}\sqrt{\frac{m}{2}{\omega }_{ks}}{\mathbf{e}}_{ks}\left(b_{ks}-b_{-ks}^{†}\right)e^{i\mathbf{k}·{\mathbf{a}}_j},$$

$${\mathbf{u}}_j={\mathbf{v}}_j+\frac{1}{\sqrt{N}}\sum _{ks}\frac{{\mathbf{e}}_{ks}}{\sqrt{2m{\omega }_{ks}}}\left(b_{ks}+b_{-ks}^{†}\right)e^{i\mathbf{k}·{\mathbf{a}}_j}.$$

(133)

The nonuniformity comes through an additional term in the expression for $u_j$.

Hamiltonian (125) transforms to

$${\widehat{H}}_1=E_1+{\widehat{H}}_{ph}^{\prime }+{\widehat{H}}_{ps}+{\widehat{H}}_{lin}^{\prime }.$$

(134)

The effective phonon Hamiltonian is

$${\widehat{H}}_{ph}^{\prime }=w\sum _{ks}{\omega }_{ks}\left(b_{ks}^{†}{b}_{ks}+\frac{1}{2}\right),$$

(135)

with the phonon frequency given by the equation

$${\omega }_{ks}^2=\frac{w}{m}\sum _j\sum _{\alpha \beta }{\mathsf{\Phi }}_{ij}^{\alpha \beta }{e}_{ks}^{\alpha }{e}_{ks}^{\beta }{e}^{i\mathbf{k}·{\mathbf{a}}_{ij}}.$$

(136)

The latter equation can also be written as

$${\omega }_{ks}^2=\frac{w}{m}\sum _{\alpha \beta }{\mathsf{\Phi }}_k^{\alpha \beta }{e}_{ks}^{\alpha }{e}_{ks}^{\beta }.$$

The momentum $\mathbf{k}$ pertains to the Brillouin zone.

After transformation (133), the term, remaining from the renormalized linear pseudospin-phonon interaction, reads as

$${\widehat{H}}_{lin}^{\prime }=-2w^2\sum _{ij}\sum _{\alpha }{K}_{ij}^{\alpha }{v}_j^{\alpha },$$

(137)

with

$$v_f^{\alpha }=\frac{1}{2N}\sum _{ij}\sum _{\beta }{\gamma }_{jf}^{\alpha \beta }{K}_{ij}^{\beta },$$

(138)

where the notation

$${\gamma }_{ij}^{\alpha \beta }\equiv 4w\sum _{ks}\frac{e_{ks}^{\alpha }{e}_{ks}^{\beta }}{m{\omega }_{ks}^2}\exp \left(i\mathbf{k}·{\mathbf{a}}_{ij}\right)$$

(139)

is used. The latter quantity possesses the properties

$${\gamma }_{ij}^{\alpha \beta }={\gamma }_{ji}^{\alpha \beta }={\gamma }_{ji}^{\beta \alpha }={\gamma }_{ij}^{\beta \alpha }.$$

Combining (137) and (138) gives

$${\widehat{H}}_{lin}^{\prime }=-\frac{w^2}{N}\sum _{ij}\sum _{fg}\sum _{\alpha \beta }{K}_{ij}^{\alpha }{\gamma }_{jf}^{\alpha \beta }{K}_{fg}^{\beta },$$

(140)

which shows that this is an effective four-pseudospin interaction.

In the summation over momenta in (139), the main contribution comes from the term with $k=0$ because of the fast oscillations of the exponential. Therefore, expression (139) can be well approximated as

$${\gamma }_{ij}^{\alpha \beta }\cong {\gamma }^{\alpha \beta }\equiv 4w\sum _{ks}\frac{e_{ks}^{\alpha }{e}_{ks}^{\beta }}{m{\omega }_{ks}^2}.$$

(141)

Then (138) becomes

$$v_f^{\alpha }\cong \frac{1}{2N}\sum _{ij}\sum _{\beta }{\gamma }^{\alpha \beta }{K}_{ij}^{\beta }.$$

(142)

Owing to the property

$$\sum _{ij}{K}_{ij}^{\beta }=-\sum _{ij}{K}_{ij}^{\beta }=0,$$

we have

$${\widehat{H}}_{lin}^{\prime }=0,v_f^{\alpha }=0.$$

(143)

Consequently, the term $H_{lin}^{\prime }$ in Hamiltonian (134) can be omitted.

This is in agreement with the following. The existence of the linear, in the deviations, terms in Hamiltonian (130), formally, leads to the fact that, according to (133), the average $<u_j>$ could be nonzero, which, however would be in contradiction with condition (102). Therefore, the linear in deviations terms should be zero, as in Eqs. (143). Also, the appearance of the linear terms in Hamiltonian (130) breaks the symmetry of the initial Hamiltonian with respect to the inversion $u_j\to -u_j$. Generally, such linear terms should either be zero, due to symmetry properties, or are to be canceled by counterterms preserving the equilibrium condition (102).

The same conclusion could be obtained, if, when decoupling the pseudospin and phonon degrees of freedom, we would employ the decoupling

$$u_{ij}^{\alpha }{S}_{ij}^{\beta }=\langle u_{ij}^{\alpha }\rangle S_{ij}^{\beta }+u_{ij}^{\alpha }\langle S_{ij}^{\beta }\rangle -\langle u_{ij}^{\alpha }\rangle \langle S_{ij}^{\beta }\rangle ,$$

(144)

which, in view of (116), yields

$$u_{ij}^{\alpha }{S}_{ij}^{\beta }=u_{ij}^{\alpha }\langle S_{ij}^{\beta }\rangle .$$

(145)

Then, employing the properties

$$\sum _j{B}_{ij}^{\alpha }\langle S_{ij}^x\rangle =\frac{\partial }{\partial a_i^{\alpha }}\sum _j{B}_{ij}\langle S_{ij}^x\rangle =0,$$

$$\sum _j{I}_{ij}^{\alpha }\langle S_{ij}^z\rangle =\frac{\partial }{\partial a_i^{\alpha }}\sum _j{I}_{ij}\langle S_{ij}^z\rangle =0$$

(146)

and

$$\sum _{ij}{B}_{ij}^{\alpha }{u}_{ij}^{\alpha }\langle S_{ij}^x\rangle =0,\sum _{ij}{I}_{ij}^{\alpha }{u}_{ij}^{\alpha }\langle S_{ij}^z\rangle =0,$$

(147)

following from the ideality of the lattice, results in the Hamiltonian

$${\widehat{H}}_1=E_1+{\widehat{H}}_{ph}+{\widehat{H}}_{ps}.$$

(148)

The first term here is the nonoperator part (124), the second term is the effective phonon Hamiltonian

$${\widehat{H}}_{ph}=w\sum _{ks}{\omega }_{ks}\left(b_{ks}^{†}{b}_{ks}+\frac{1}{2}\right).$$

(149)

And the effective pseudospin Hamiltonian is

$${\widehat{H}}_{ps}=-w\mathsf{\Omega }\sum _j{S}_j^x+w^2\sum _{ij}{\tilde{B}}_{ij}{S}_i^x{S}_j^x-w^2\sum _{ij}{\tilde{I}}_{ij}{S}_i^z{S}_j^z.$$

(150)

The pseudospin interactions are renormalized by the existence of the phonon vibrations. And the phonon characteristics are renormalized due to the phonon interactions with pseudospins. In addition, all quantities are renormalized by the presence of the heterophase fluctuations. It is, therefore, possible to call the resulting effective phonons as dressed phonons.

13. Heterophase Fluctuations

With the derived effective Hamiltonian (148), we can explicitly calculate all phonon characteristics. For instance, the phonon distribution is

$$n_{ks}\equiv \langle b_{ks}^{†}{b}_{ks}\rangle ={\left[\exp \left(\frac{w{\omega }_{ks}}{T}\right)-1\right]}^{-1}.$$

(151)

The deviation-deviation correlation function read as

$$\langle u_j^{\alpha }{u}_j^{\beta }\rangle =\frac{1}{2N}\sum _{ks}\frac{e_{ks}^{\alpha }{e}_{ks}^{\beta }}{m{\omega }_{ks}}\coth \left(\frac{w{\omega }_{ks}}{2T}\right).$$

(152)

The mean kinetic energy becomes

$$\langle \frac{{\mathbf{p}}^2}{2m}\rangle =\frac{1}{4N}\sum _{ks}{\omega }_{ks}\coth \left(\frac{w{\omega }_{ks}}{2T}\right).$$

(153)

In order to treat the pseudospin variables, we can use the mean-field approximation for the products $S_i^{\alpha }{S}_j^{\beta }$. This is possible because the critical region in ferroelectrics is known [91,92] to be narrow due to the smallness of the Ginzburg number. Then, in the mean-field approximation for the pseudospins, we find the averages for the x-component

$$\langle S_j^x\rangle =w\frac{{\mathsf{\Omega }}_j}{2H_j}\tanh \left(\frac{H_j}{2T}\right),$$

(154)

for the y-component

$$\langle S_j^y\rangle =0,$$

(155)

and for the z-component

$$\langle S_j^z\rangle =w^2\langle S_j^z\rangle \frac{\tilde{I}}{H_j}\tanh \left(\frac{H_j}{2T}\right).$$

(156)

Here, the notations are used for the effective tunneling frequency

$${\mathsf{\Omega }}_j=\mathsf{\Omega }-2w\tilde{B}\langle S_j^x\rangle$$

(157)

and the effective field

$$H_j=w\sqrt{{\mathsf{\Omega }}_j^2+4w^2{\tilde{I}}^2\langle S_j^z\rangle },$$

(158)

in which

$$\tilde{B}\equiv \sum _j{\tilde{B}}_{ij},\tilde{I}\equiv \sum _j{\tilde{I}}_{ij}.$$

(159)

To deal further with considering the properties of the heterophase ferroelectric, we need to restore the phase indices $\nu =1,2$. The effective Hamiltonian of the heterophase system is

$$\tilde{H}={\widehat{H}}_1⨁{\widehat{H}}_2.$$

(160)

The conditions, distinguishing the ordered and disordered phases are

$$\langle S_{j1}^z\rangle \ne 0,\langle S_{j2}^z\rangle =0.$$

All necessary equations for the ferroelectric phase are written above. The equations for the paraelectric phase can be obtained from the above expressions by setting there the order parameter $<S_j^z>=0$. The equations for the phase probabilities are given in Sec. 8.

14. Phase Transition

The total system of equations, defining the heterophase ferroelectric, can be solved numerically. Below, we present some of the most interesting conclusions describing the influence of the heterophase fluctuations on the system properties. The strongest influence of these fluctuations occurs in the vicinity of the phase transition point, that is represented by the temperature

$${\tilde{T}}_c=\frac{(1-\tilde{b})\tilde{\omega }}{4\mathrm{artanh}(2\tilde{\omega })},$$

(161)

where

$$\tilde{b}\equiv \frac{\tilde{B}}{\tilde{I}+\tilde{B}},\tilde{\omega }\equiv \frac{\tilde{\mathsf{\Omega }}}{\tilde{I}+\tilde{B}}.$$

The existence of the heterophase fluctuations can be noticed and their influence measured by studying, e.g., the Debye-Waller factor $f_{DW}$ and sound velocity s. It is convenient to characterize the influence of the heterophase fluctuations by comparing their values in the presence of the latter with the related values without these fluctuations. For example, comparing the sound velocity s in the heterophase system with the sound velocity $s_0$ in a pure system without such fluctuations, it is useful to introduce the relative sound-velocity decrease

$$\delta s\equiv \frac{s-s_0}{s_0}.$$

(162)

Another useful characteristic is the Debye-Waller factor $f_{DW}$ that can be measured by x-ray scattering, coherent neutron scattering, and by Mössbauer spectroscopy. All details and definitions can be found in the books [93,94,95]. We need to compare the Debye-Waller factor ${\tilde{f}}_{DW}$ for the heterophase system and its value $f_{DW}$ for a pure system. Again, it is useful to employ the relative Debye-Waller factor decrease

$$\delta f_{DW}\equiv \frac{{\tilde{f}}_{DW}-f_{DW}}{f_{DW}}.$$

(163)

It is interesting that, as follows from numerical calculations, these quantities, the relative sound velocity decrease and the relative Debye-Waller factor decrease, are universal, weakly depending on the considered materials. For the relative sound velocity decrease at the critical transition point, we get

$$\delta s\approx -0.3(T=T_c)$$

(164)

and for the relative Debye-Waller factor decrease at this point,

$$\delta f_{DW}\approx -0.3.$$

(165)

This decrease of the sound velocity and of the Debye-Waller factor is due to the scattering caused by heterophase fluctuations.

15. Conclusion

The systems are addressed, exhibiting phase transitions between thermodynamic phases with different symmetry, in which spontaneous symmetry breaking can be accompanied by local spontaneous symmetry restoration, caused by the appearance of heterophase fluctuations. In such systems, above the phase transition point, where a disordered phase dominates, there appear the germs of the ordered phase. And below the transition point, where an ordered phase prevails, there arise the embryos of the disordered phase. Thus, around the phase transition point there is a region, where the phases with higher and lower symmetries coexist. The embryonic regions of a competing phase are distributed randomly in space. Their characteristic sizes are mesoscopic, so that the typical size $r_f$ of such an embryo is much larger than the mean interparticle distance a, but much smaller than the characteristic size of the whole system L, that is,

$$a\ll r_f\ll L.$$

The typical size $r_f$ should be understood as an average size, since the shapes of the embryonic heterophase fluctuations are not necessarily regular, but can be rather ramified. Therefore, terming these germs mesoscopic could be better applied to the average number of particles $N_f$ composing each of such embryos, as compared to the total number N of particles in the system, so that

$$1\ll N_f\ll N.$$

This type of heterophase fluctuations is common for a number of substances.

The general approach for treating these heterophase systems is developed, being based on the notion of weighted Hilbert spaces. The real-space distribution of the phases is described by means of manifold indicators. The averaging over random phase configurations reduces the problem to the consideration of an effective renormalized Hamiltonian composed of the phase-replica Hamiltonians representing the phases of different symmetry. Stability conditions define the geometric phase probabilities in a self-consistent way.

The method is illustrated by applying it to heterophase ferroelectrics that are the typical materials exhibiting the appearance of such heterophase fluctuations around their phase transition points between the paraelectric and ferroelectric phases. The influence of the heterophase fluctuations is the strongest in the vicinity of the phase transition point. Numerical calculations show that the occurrence of such fluctuations leads to the noticeable decrease of the sound velocity and Debye-Walle factor at the transition point. The relative values of this decrease turn out to be universal, only weakly depending on the material parameters.

In conclusion, it is worth mentioning that thermodynamic phases and phase transitions between them can be conveniently characterized not only by order parameters but also by order indices [96]. Another important characteristic is the measure of entanglement in the considered physical system [97]. These three characteristics are interrelated with each other [98,99,100]. The usual situation is when the increasing order is accompanied by the decreasing entanglement [100]. Since a mesoscopic mixture is a system that is between an absolutely disordered and a completely ordered phases, its entanglement should be between these two limiting cases. An additional entanglement arises in mesoscopic mixtures because a heterophase system consists of several spatially separated regions with different symmetry, and these mesoscopic regions are mutually entangled. This increases the system entanglement, as compared to the completely ordered case. The problem of studying the entanglement level of mesoscopic mixtures is extremely interesting. Since the mixture properties can be regulated by means of thermodynamic parameters, this gives an additional possibility of governing the level of entanglement in such mesoscopic mixtures, which does not exist in the case of pure systems. This novel possibility of regulating entanglement is really exciting, but it is a separate problem that goes out of the scope of the present paper. It is clear that the mesoscopic mixtures provide a nontrivial possibility for regulating the system entanglement, which can be used for quantum information processing.