text:
# Deterministic Secure Direct Communication Using Entanglement

## Abstract

A novel secure communication protocol is presented, based on an entangled pair of qubits and allowing asymptotically secure key distribution and quasi-secure direct communication. Since the information is transferred in a deterministic manner, no qubits have to be discarded. The transmission of information is instantaneous, i.e. the information can be decoded during the transmission. The security against arbitrary eavesdropping attacks is provided. In case of eavesdropping attacks with full information gain, the detection rate is 50% per control transmission. The experimental realization of the protocol is feasible with relatively small effort, which also makes commercial applications conceivable.

Introduction.— Cryptographic schemes based on quantum mechanics are usually *non-deterministic* BB84 ; Ekert ; Bruss . Alice, the sender, can encode a classical bit into a quantum state, which is then sent to Bob, but she cannot determine the bit value that Bob will finally *decode*. Inspite of that, such non-deterministic communication can be used to establish a *shared secret key* between Alice and Bob, consisting of a sequence of random bits. This secret key can then be used to encrypt a message which is sent through a classical public channel. Recently, a novel quantum communication protocol has been presented Almut1  that allows *secure direct communication*, where the message is deterministically send through the quantum channel, but can only be decoded after a final transmission of classical information. We present a communication scheme, the “ping-pong protocol”, that also allows for deterministic communication. This protocol can be used for the transmission of either a secret key or a plaintext message. In the latter case, the protocol is *quasi-secure*, i.e. an eavesdropper is able to gain a small amount of message information before being detected. In case of a *key transmission* the protocol is asymptotically secure. In contrast to other quantum cryptographic schemes, the presented scheme is *instantaneous*, i.e. the information can be decoded during the transmission and no final transmission of additional information is needed. The basic idea of the protocol, encoding information by local operations on an EPR pair, has already been raised by Bennett and Wiesner BW92 . In our protocol, we follow this idea, but abandon the *dense* coding feature in favour of a *secure* transmission.

The ping-pong protocol.— When two photons are maximally *entangled* in their polarization degree of freedom, then each single photon is not polarized at all. Denote the horizontal and vertical polarization state by $`|0`$ and $`|1`$, respectively, then the *Bell states* $`|\psi ^\pm =\frac{1}{\sqrt{2}}(|01\pm |10)`$ are maximally entangled states in the two-particle Hilbert space $`=_A_B`$. A measurement of the polarization of one photon, say $`A`$, leads to a completely random result. This is reflected by the fact that the corresponding *reduced density matrices*, $`\rho _A^\pm :=\mathrm{Tr}_B\{|\psi ^\pm \psi ^\pm |\}`$ are both equal to the complete mixture, $`\rho _A^\pm =\frac{1}{2}\mathrm{𝟙}_A`$. Hence, no experiment performed on only one photon can distinguish these states from each other. However, since the states $`|\psi ^\pm `$ are mutually orthogonal, a measurement on *both* photons can perfectly distinguish the states from each other. In other words: One bit of information can be encoded in the states $`|\psi ^\pm `$, which is completely unavailable to anyone who has only access to one of the photons. As one can easily verify, the unitary operator $`\widehat{\sigma }_z^A(\widehat{\sigma }_z\mathrm{𝟙})=(|00||11|)\mathrm{𝟙}`$ flips between the two states $`|\psi ^\pm `$,

$$\widehat{\sigma }_z^A|\psi ^\pm =|\psi ^{}.$$
(1)
Altough $`\widehat{\sigma }_z^A`$ acts *locally*, i.e. on one photon only, it has a *non-local* effect. Someone who has access to one single photon only, can *encode* one bit of information, but he cannot *decode* it, since he has no access to the other photon. This is a situation perfectly suited for a cryptographic scenario. Bob prepares two photons in the state $`|\psi ^+`$. He keeps one photon, the “*home qubit*”, and sends the other one, the “*travel qubit*”, to Alice (“ping!”). Alice decides either to perform the operation $`\widehat{\sigma }_z`$ on the travel qubit or to do nothing, i.e. to perform the operation $`\mathrm{𝟙}`$. Then she sends the travel qubit back to Bob (“pong!”). Bob, who has now both qubits again, performs a Bell measurement resulting in either $`|\psi ^+`$ or $`|\psi ^{}`$, depending on what Alice did. Thus, he has received one bit of information from Alice. One qubit travels forth and back (“ping-pong!”) and one bit of information flows from Alice to Bob. Let us introduce two communication modes, “*message mode*” and “*control mode*” (see Figs. 1,2). By default, Alice and Bob are in message mode and communicate the way described above. With probability $`c`$, Alice switches to control mode and instead of performing her operation on the travel qubit, she performs a measurement in the basis $`_z=\{|0,|1\}`$. Using the public channel, she sends the result to Bob, who then also switches to control mode and performs a measurement in the same basis $`_z`$. Bob compares his own result with Alice’s result. If both results coincide, Bob knows that Eve is in the line and stops the communication.

\[t\] Let us give an explicit algorithm for the protocol.

1. Protocol is initialized. $`n=0`$. The message to be transmitted is a sequence $`x^N=(x_1,\mathrm{},x_N)`$, where $`x_n\{0,1\}`$.
2. $`n=n+1`$. Alice and Bob are set to message mode. Bob prepares two qubits in the Bell state $`|\psi ^+=\frac{1}{\sqrt{2}}(|01+|10)`$.
3. He stores one qubit, the *home qubit*, and sends the other one, the *travel qubit*, to Alice through the quantum channel.
4. Alice receives the travel qubit. With probability $`c`$ she switches to control mode and proceeds with c.1, else she proceeds with m.1.

1. Alice measures the travel qubit in the basis $`_z`$ and obtains the result $`i\{0,1\}`$ with equal probability.
    2. She sends $`i`$ through the public channel to Bob.
    3. Bob receives $`i`$ from the public channel, switches to control mode and measures the home qubit in the basis $`_z`$ resulting in the value $`j`$.
    4. ($`i=j`$): Eve is detected. Abort transmission. ($`ij`$): Set $`n=n1`$ and Goto p.1.
    1. Define $`\widehat{C}_0:=\mathrm{𝟙}`$ and $`\widehat{C}_1:=\widehat{\sigma }_z`$. For $`x_n\{0,1\}`$, Alice performs the coding operation $`\widehat{C}_{x_n}`$ on the travel qubit and sends it back to Bob.
    2. Bob receives the travel qubit and performs a Bell measurement on both qubits resulting in the final state $`|\psi ^{}\{|\psi ^+,|\psi ^{}\}`$. He decodes the message as

    $$|\psi ^{}=\{\begin{array}{cc}|\psi ^+x_n=0\hfill & \\ |\psi ^{}x_n=1\hfill & \end{array}.$$
    (2)
    3. ($`n<N`$): Goto p.1. ($`n=N`$): Goto p.4.
5. Message $`x^N`$ is transmitted from Alice to Bob. Communication successfully terminated.

Security proof.— Eve is an evil quantum physicist able to build all devices that are allowed by the laws of quantum mechanics. Her aim is to find out which operation Alice performs. Eve has no access to Bob’s home qubit, so all her operations are restricted to the travel qubit, whose state is (to Eve) indistinguishable from the complete mixture $`\rho _A=\mathrm{Tr}_B\{|\psi ^+\psi ^+|\}=\frac{1}{2}\mathrm{𝟙}_A`$. The most general quantum operation is a *completely positive map* $`:𝒮(_A)𝒮(_A)`$ on the state space $`𝒮(_A)`$. Due to the *Stinespring dilation theorem* Stinespring , any completely positive map can be realized by a unitary operation on a larger Hilbert space. For $`_A`$ and $``$ given, there is an *ancilla space* $`_E`$ of dimension $`dim_E(dim_A)^2`$, an ancilla state $`|\chi _E`$, and a unitary operation $`\widehat{E}`$ on $`_A_E`$, such that for all states $`\rho _A𝒮(_A)`$, we have

$$(\rho _A)=\mathrm{Tr}_E\{\widehat{E}(\rho _A|\chi \chi |)\widehat{E}^{}\}.$$
(3)
In order to gain information about Alice’s operation, Eve should first perform the unitary *attack operation* $`\widehat{E}`$ on the composed system, then let Alice perform her coding operation $`\widehat{C}`$ on the travel qubit, and finally perform a measurement on the composed system (see Fig. 3). Since a probable control measurement by Alice takes place *before* Eve’s final measurement, the latter has no influence on the detection probability for Eve’s attack. All that can be detected is the attack operation $`\widehat{E}`$. Let us analyze the detection probability $`d`$, given an attack operation $`\widehat{E}`$. Since for Eve the state of the travel qubit is indistinguishable from the complete mixture, we can replace the state of the travel qubit by the *a priori* mixture $`\rho _A=\frac{1}{2}|00|+\frac{1}{2}|11|`$, which corresponds to the situation where Bob sends the travel qubit in either of the states $`|0`$ or $`|1`$, with equal probability $`p=1/2`$. Let us at first consider the case where Bob sends $`|0`$. Alice adds an ancilla in the state $`|\chi `$ and performs the unitary operation $`\widehat{E}`$ on both systems, resulting in

$`|\psi ^{}`$ $`=`$ $`\widehat{E}|0,\chi =\alpha |0,\chi _0+\beta |1,\chi _1,`$ (4)
where $`|\chi _0,|\chi _1`$ are pure ancilla states uniquely determined by $`\widehat{E}`$, and $`|\alpha |^2+|\beta |^2=1`$. In a subsequent control measurement, Alice measures the travel qubit in the basis $`_z=\{|0,|1\}`$ and sends the result to Bob. Without Eve, the result will always read “0”, hence the detection probability for Eve’s attack in a control run reads

$$d=|\beta |^2=1|\alpha |^2.$$
(5)
Now let us analize how much information Eve can maximally gain when there is no control run. After Eve’s attack operation, the state of the system reads

$`\rho ^{}`$ $`=`$ $`|\psi ^{}\psi ^{}|=|\alpha |^2|0,\chi _00,\chi _0|+|\beta |^2|1,\chi _11,\chi _1|+\alpha \beta ^{}|0,\chi _01,\chi _1|+\alpha ^{}\beta |1,\chi _10,\chi _0|,`$ (6)
which can be rewritten in the orthogonal basis $`\{|0,\chi _0,|1,\chi _1\}`$ as

$$\rho ^{}=\left(\begin{array}{cc}|\alpha |^2& \alpha \beta ^{}\\ \alpha ^{}\beta & |\beta |^2\end{array}\right).$$
(7)
Alice encodes her bit by applying the operation $`\widehat{C}_0=\mathrm{𝟙}`$ or $`\widehat{C}_1=\widehat{\sigma }_z`$ to the travel qubit, with probability $`p_0`$ and $`p_1`$, respectively. The state of the travel qubit after Eve’s attack operation and after Alice’s encoding operation reads

$$\rho ^{\prime \prime }=\left(\begin{array}{cc}|\alpha |^2& \alpha \beta ^{}(p_0p_1)\\ \alpha ^{}\beta (p_0p_1)& |\beta |^2\end{array}\right).$$
(8)
The maximal amount $`I_0`$ of classical information that can be extracted from this state is given by the *von-Neumannn entropy*, $`I_0=S(\rho ^{\prime \prime })\mathrm{Tr}\{\rho ^{\prime \prime }\mathrm{log}_2\rho ^{\prime \prime }\}`$. In order to calculate the von-Neumann entropy we need the eigenvalues $`\lambda `$ of $`\rho ^{\prime \prime }`$, which are the roots of the characteristic polynomial $`det(\rho ^{\prime \prime }\lambda \mathrm{𝟙})`$, yielding the two eigenvalues

$$\lambda _{1,2}=\frac{1}{2}\left(1\pm \sqrt{14|\alpha \beta |^2[1(p_0p_1)^2]}\right),$$
(9)
so we have

$$I_0=\lambda _1\mathrm{log}_2\lambda _1\lambda _2\mathrm{log}_2\lambda _2.$$
(10)
The maximal information gain $`I_0`$ can be expressed as a function of the detection probability $`d`$. Using (5), we have $`|\alpha \beta |^2=(1|\beta |^2)|\beta |^2=(dd^2)`$, and therefore

$$\lambda _{1,2}=\frac{1}{2}\pm \frac{1}{2}\sqrt{1(4d4d^2)[1(p_0p_1)^2]}.$$
(11)
Now assume that Bob sends $`|1`$ rather than $`|0`$. The above calculations can be done in full analogy, resulting in the same crucial relations (10,11). Eve’s task is, of course, to minimize $`d`$. Though if she chooses an eavesdropping action $`\widehat{E}`$ that provides $`d=0`$, then $`\lambda _1=1,\lambda _2=0`$, which implies $`I_0=0`$, therefore Eve can gain no information at all. Thus we have shown:

Any effective eavesdropping attack can be detected.— In the case $`p_0=p_1=1/2`$, where Alice encodes exactly 1 bit, expression (11) simplifies to $`\lambda _{1,2}=\frac{1}{2}\pm |\frac{1}{2}d|`$, or $`\lambda _1=d,\lambda _2=1d`$. Interestingly, the maximal information gain is equal to the Shannon entropy of a binary channel,

$$I_0(d)=d\mathrm{log}_2d(1d)\mathrm{log}_2(1d).$$
(12)
The function $`I_0(d)`$ has a maximum at $`d=1/2`$, and can be inversed on the interval $`[0,1/2]`$, giving a monotonous function $`0d(I_0)1/2,I_0[0,1]`$. By choosing a desired information gain $`I_0>0`$ per attack, Eve has to face a detection probability $`d(I_0)>0`$. If she wants to gain the *full* information ($`I_0=1`$), the detection probability is $`d(I_0=1)=1/2`$.

Direct communication versus key distribution.— In contrast to quantum key distribution protocols like BB84 BB84 , the ping-pong protocol provides a *deterministic* transmission of bits, hence it is possible to communicate the message *directly* from Alice to Bob. Assuming that Eve wants to gain full information in each attack, the ping-pong protocol provides a detection probability of $`d=1/2`$, which is significantly higher than the detection probability of the BB84 protocol, where we have $`d=\frac{1}{2}\times \frac{1}{2}=\frac{1}{4}`$ for the same situation. Furthermore, the BB84 protocol has a probability of $`1/2`$ that a transmitted bit has to be discarded due to the wrong choice of basis on both sides.

Taking into account the probability $`c`$ of a control run, the *effective transmission rate*, i.e. the number of message bits per protocol run, reads $`r=1c`$, which is equal to the probability for a message transfer. Say, Eve wants to eavesdrop one message transfer without being detected. The probability for this event reads

$`s(c,d)`$ $`=`$ $`(1c)+c(1d)(1c)+`$ (13)
$`c^2(1d)^2(1c)+\mathrm{}`$
$`=`$ $`{\displaystyle \frac{1c}{1c(1d)}},`$ (14)
where the terms in the (geometric) series correspond to Eve having to survive $`0,1,2,\mathrm{}`$ control runs before she gets to eavesdrop on a message run, finally yielding the desired information of $`I_0(d)`$ bits. After $`n`$ successful attacks Eve gains $`nI_0(d)`$ bits of information and survives with probability $`s^n`$, thus the probability to successfully eavesdrop $`I=nI_0(d)`$ bits reads $`s(I,c,d)=s(c,d)^{I/I_0(d)}`$, so

$$s(I,c,d)=\left(\frac{1c}{1c(1d)}\right)^{I/I_0(d)},$$
(15)
where $`I_0(d)`$ is given by (12). For $`c>0`$, $`d>0`$, this value decreases exponentially but is nonzero. In the limit $`I\mathrm{}`$ (a message or key of infinite length) we have $`s0`$, so the protocol is *asymptotically secure*, just like the BB84 protocol. Let us give an example. A convenient choice of the control parameter is $`c=0.5`$, where on the average every second bit is a control bit. Say, Eve wants to gain full information in each attack, thus $`I_0=1`$ and $`d=1/2`$. The probability that Eve successfully eavesdrops 1 character (8 bits) is already as low as $`s0.039`$. In Fig. 4 we have plotted the eavesdropping success probability as a function of the information gain $`I`$, for $`c=0.5`$ and for different detection probabilities $`d`$ that Eve can choose. (Note that for $`d<1/2`$ Eve only gets *part* of the message right and does not even know *which* part.) If desired, the security can arbitrarily be improved by increasing the control parameter $`c`$ at the cost of decreasing the transmission rate. Let us call such communication “*quasi-secure*”. If we want a *perfectly* secure communication (which is, strictly speaking, also not really perfect), we must abandon the direct transfer in favour of a *key* transfer. In this case, Alice does not transmit the message directly to Bob but rather takes a random sequence of $`N`$ bits from a secret random number generator. After a succesful transmission, the random sequence is used as a *shared secret key* between Alice and Bob. Eve has virtually no advantage in eavesdropping only a few bits, because one can choose classical privacy amplification protocols that make it *very* hard to decode parts of the message with only *some* of the key bits given. The one-time-pad scheme, by the way, is not quite a good choice, because here each eavesdropped key bit directly yields one decoded message bit. Anyway, as soon as Eve is detected, the transfer stops and she has learned nothing but a sequence of nonsense random bits.

Experimental feasibility.— The Bell state $`|\psi ^+`$ can be created by *parametric down-conversion*. Bob’s Bell measurement must only distinguish between the states $`|\psi ^\pm `$, which can be accomplished, too. The storage of one photon is necessary only for a duration corresponding to twice the distance between Alice and Bob. The encoding procedure corresponds to a controlled $`\widehat{\sigma }_z`$-operation, which can be realized by triggered optical elements. The correlation test involves a simple measurement of the linear polarization in a fixed basis. Altogether, the experimental realization of the ping-pong protocol should be feasible using nowaday’s technology. Even a commercial application could be envisaged.

We had fruitful discussions with Almut Beige, Luke Rallan, Jens Eisert, Martin Plenio, Sougato Bose, and others. This work is supported by the Deutsche Forschungsgemeinschaft (DFG) and by the European Union (EQUIP).


format: text

text:
# Thermodynamics of fiber bundles

## 1 INTRODUCTION

The fiber-bundle model with global load sharing is a simple model for failure in tension introduced almost 80 years ago \[1–3\] and having received considerable attention and extensions over the past 15 years \[4–19\] . Although this model may have little pertinence to real fibrous systems such as a rope breaking in tension, it is of interest because it possesses exact analytical properties.

In recent work \[20–22\] , we have developed a general statistical theory for determining the properties of a disordered solid that is accumulating irreversible damage due to cracking under stress. The ensembles in this theory are created by considering different realizations of the quenched disorder in the local breaking strength of the material. Such realizations are made either for the system as a whole, or of more pertinence to real systems, by dividing a given system into smaller “mesovolumes” and letting each mesovolume correspond to a different realization of the quenched disorder. The ensembles so obtained have nothing to do with thermal fluctuations (molecular dynamics). Given a system with such quenched disorder, our theory determines the probability of emergent crack states by maximizing Shannon’s measure of disorder subject to constraints coming from the energetics of the fracture process.

In the present paper, we apply this theory to the specific problem of fiber bundles with global load sharing and demonstrate that for any quenched-disorder distribution, it yields the known exact results. Furthermore, a previously unnoticed phase transition is demonstrated to exist where entropy goes through a maximum. This phase transition is distinct from the well-known stress maximum transition and was noticed in the present theory because of the prominent role played by entropy.

However, the importance of our theory is not that it yields a new result in this old model, but that it applies and yields analytical results about phase transitions for any quasi-static damage model; albeit, for models involving crack interactions, approximate treatment of the functional integrations may be required to obtain analytical results (such as renormalization or mean-field approximations). Using our theory, we recently treated the problem of how the mechanical properties of rocks change due to cracks arriving and interacting in compressive shear \[20–22\]. We analytically demonstrated that the localization transition observed in experiments is a second-order phase transition when it occurs.

## 2 THERMODYNAMICS

Our theory was originally developed for a disordered solid under the influence of a stress tensor. For fiber bundles, a much simpler scalar theory applies and so in this section, the formalism is rederived in this simplified context.

A fiber bundle is depicted in Fig. 1. A collection of $`N`$ fibers are stretched between two rigid supports. One support is held fixed, while the other is free to displace. A load $`F_N`$ is applied to the bundle through the free support so that the fibers are in a state of tension. In this paper, the load $`F_N`$ will always be normalized by $`N`$ to define the overall tension $`\tau =F_N/N`$. Each non-broken fiber in a bundle has the same length $`L`$. If, when $`\tau =0`$, this length is $`L_o`$, then the measure of strain is $`\epsilon =L/L_o1`$. Experiments may be performed on the bundle either by controlling $`\tau `$ or $`\epsilon `$.

Each of the fibers has the same Young’s modulus which is taken to be unity so that the axial strain $`\epsilon `$ of each fiber is identical to the tension in the fiber. The $`N`$ fibers have strengths $`\epsilon _1`$, $`\epsilon _2`$,…$`\epsilon _N`$ which are independent random variables sampled from a distribution $`p(\epsilon )`$, whose cumulative distribution is defined $`P(\epsilon )=_0^\epsilon p(x)𝑑x`$. As the strain $`\epsilon `$ of the bundle is increased, the individual fibers will break once their tension (strain) gets to their fixed strength threshold. The assumption of “global load sharing” is that when one of the fibers breaks at fixed load, all the other fibers will extend by the same amount thus increasing the tension in each of the surviving fibers so that the load as a whole is always supported entirely by the surviving fibers. All of this defines the fiber bundle model with global load sharing. Of interest are the mechanical properties of such bundles as averaged over all possible realizations of the fiber strengths.

We first need to know the probability $`p_j`$ of observing one of the realizations to be in a particular state $`j`$ of damage when the ensemble as a whole is at an applied strain $`\epsilon `$. In the fiber bundle model, a damage state $`j`$ is defined by which of the $`N`$ fibers are broken. One could define $`j`$ using a local order parameter that is 1 if a fiber is intact and 0 if the fiber is broken.

Our theory postulates that the fraction $`p_j`$ of all realizations observed to be in state $`j`$ is obtained by maximizing Shannon’s measure of disorder

$$S=\underset{j}{}p_j\mathrm{ln}p_j$$
(1)
subject to constraints. Such constraints must involve the independent variables of $`S`$. To identify the independent variables, we consider how both $`S`$ and the average energy in the ensemble of bundles change as the strain is increased.

When $`\epsilon `$ increases to $`\epsilon +d\epsilon `$, there is both a work carried out in reversibly stretching the fibers and an additional work carried out due to irreversible fiber breaks. Due to breaking, some of the members of the ensemble (individual realizations of the disorder) will be led out of their current damage state and into state $`j`$, while others that were in state $`j`$ will transfer to still different states. If there is a difference in the number of members entering and leaving state $`j`$, there will be a change $`dp_j`$ in the occupation probability of state $`j`$ and such changes are what cause Shannon’s disorder measure $`S`$ to change.

The average energy density (average energy normalized by $`N`$) in the ensemble is given by $`U=_jp_jE_j`$. Here, $`E_j`$ is the energy density required to create state $`j`$ at imposed strain $`\epsilon `$ and as averaged over all members that have been led to state $`j`$. Depending on the breaking strengths of a given realization, the work performed in arriving at state $`j`$ can be different. It is through $`E_j`$ that all dependence on the quenched-disorder distribution enters the problem. The change that occurs when $`\epsilon `$ increases to $`\epsilon +d\epsilon `$ is

$$dU=\underset{j}{}E_jdp_j+\underset{j}{}p_jdE_j.$$
(2)
The first term is the energy expended in changing the disorder over the collection of realizations. It is thus proportional to the disorder change and can be written

$$TdS=\underset{j}{}E_jdp_j.$$
(3)
The second term is written

$$fd\epsilon =\underset{j}{}p_jdE_j.$$
(4)
and represents both the reversible stretching energy in those members that did not experience breaks during the deformation increment, as well as the irreversible energy changes due to all the breaks that did not result in a net change in the occupation numbers of each state.

This decomposition of the energy increment can be thought of as follows. When breaks occur throughout the ensemble of realizations as $`\epsilon `$ increases to $`\epsilon +d\epsilon `$, there is a flow of members between the states. This flow involves energy changes due to fibers breaking in the interval $`d\epsilon `$. It may be resolved into a uniform “incompressible” part whose associated energy is included within $`fd\epsilon `$ as well as a non-uniform “compressible” part associated with more members arriving in a given state than are leaving that state. The energy associated with the non-uniform flow between states is entirely contained in $`TdS`$.

From these expressions it may be concluded that if $`U`$ is to be treated as a fundamental function, then $`U=U(S,\epsilon )`$, or equivalently if $`S`$ is to be treated as the fundamental function then $`S=S(U,\epsilon )`$. In other words, the independent variables that must be involved in the constraints on the maximization of $`S`$ are $`U`$ and $`\epsilon `$. The proportionality constants $`T`$ and $`f`$ are defined

$$T=\left(\frac{U}{S}\right)_\epsilon \text{and}f=\left(\frac{U}{\epsilon }\right)_S.$$
(5)
The state function $`f`$ is something different than the overall tension $`\tau `$ since we also have that $`\tau d\epsilon =dU`$. Thus, in general, $`(\tau f)d\epsilon =TdS`$ so that $`f\tau `$ due to fibers breaking in a positive increment $`d\epsilon `$. If strain were to be decreased, fibers do not break and so $`dS=0`$ and the state function $`f`$ would be defined using only the purely elastic part of the energy changes $`dE_j`$. Changes in $`f`$ in this case are equivalent to changes in $`\tau `$. Last, since we have taken $`S`$ to be extensive (proportional to $`N`$) while $`U`$ is a density independent of $`N`$, $`T`$ comes out being proportional to $`1/N`$. This choice is made so that factors of $`N`$ do not clutter the equations that follow.

The constraint involving $`\epsilon `$ is that each non-broken fiber throughout the entire ensemble has the same length which implies $`\epsilon _j=\epsilon `$. The constraint involving $`U`$ is that $`U=_jp_jE_j`$. Carrying out the maximization of $`S`$ subject to these constraints using Lagrangian multipliers gives the probability distribution as

$$p_j(\beta ,\epsilon )=\frac{e^{\beta E_j(\epsilon )}}{Z(\beta ,\epsilon )}$$
(6)
where $`\beta =1/T`$ and where the partition function $`Z`$ is defined

$$Z(\beta ,\epsilon )=\underset{j}{}e^{\beta E_j(\epsilon )}.$$
(7)
The free energy $`F=F(\beta ,\epsilon )`$ in this ensemble (for this set of constraints) is the Legendre transform of $`U`$ with respect to $`S`$; *i.e.,* $`F=UTS`$. Differentiation then gives

$$dF=\beta ^2Sd\beta +fd\epsilon .$$
(8)
Upon introducing $`U=_jp_jE_j`$ and $`S=_jp_j\mathrm{ln}p_j`$ into the Legendre transform $`F=US/\beta `$ and using that $`\mathrm{ln}p_j=\beta E_j\mathrm{ln}Z`$ we obtain the standard result

$$\beta F=\mathrm{ln}Z.$$
(9)
Thus, the standard canonical ensemble emerges in this problem where quenched disorder alone (not molecular fluctuations) is responsible for the existence of ensembles.

An important question in such an approach is whether anything precise can be said about the temperature $`T=1/\beta `$. Indeed, $`\beta `$ can formally be found as the solution of a differential equation. This differential equation is obtained from the previously stated but unused fact that $`dU=\tau d\epsilon `$ which can be written

$$\tau =\frac{dU}{d\epsilon }=\underset{j}{}\frac{dp_j}{d\epsilon }E_j+\underset{j}{}p_j\frac{dE_j}{d\epsilon }.$$
(10)
Now from Eq. (6), we obtain that

$$\frac{dp_j}{d\epsilon }=\left[\frac{dE_j}{d\epsilon }\beta E_j\frac{d\beta }{d\epsilon }\frac{d\mathrm{ln}Z}{d\epsilon }\right]p_j,$$
(11)
while from Eqs. (8) and (9)

$$\frac{d\mathrm{ln}Z}{d\epsilon }=F\frac{d\beta }{d\epsilon }\beta \frac{dF}{d\epsilon }=U\frac{d\beta }{d\epsilon }\beta f.$$
(12)
Since each member is at the same $`\epsilon `$, each member has its own $`\tau _j`$, and so $`\tau =_jp_j\tau _j`$. We then obtain the differential equation for $`\beta `$ in the form

$$a\frac{d\beta }{d\epsilon }+b\beta +c=0$$
(13)
with coefficients given by

$$a=\underset{j}{}p_jE_j(UE_j);b=\underset{j}{}p_jE_j\left(f\frac{dE_j}{d\epsilon }\right);c=\underset{j}{}p_j\left(\frac{dE_j}{d\epsilon }\tau _j\right).$$
(14)
Since $`p_j=p_j(\beta ,\epsilon )`$, this equation is non-linear and thus difficult to solve. In the present work, we demonstrate that a proposed function $`\beta =\beta (\epsilon )`$ exactly satisfies this equation and thus is the true “fiber-bundle” $`\beta `$. To make progress, we next obtain $`E_j`$ and $`\tau _j`$ for the specific problem of fiber bundles with global-load sharing.

## 3 FIBER BUNDLE MODEL

Figure 2 details the history of how the overall tension $`\tau `$ might evolve when controlled variations in $`\epsilon `$ are applied to a bundle having fiber strengths $`\epsilon _1,\epsilon _2`$$`\epsilon _N`$. In this particular example, $`n=3`$ fibers have broken when the strain is at $`\epsilon `$.

The load on the bundle is equally shared by the $`Nn`$ surviving fibers so that

$$\tau _n=\epsilon \left(1\frac{n}{N}\right)$$
(15)
which is a relation independent of the history; i.e. it depends only on the actual state of the bundle through the number of broken fibers $`n`$, and not on the breaking thresholds $`\epsilon _1,\epsilon _2`$$`\epsilon _n`$.

The total work density $`E_j^p`$ for states $`j`$ consisting of $`n`$ broken fibers is the area under the charging curve in a particular realization (the area under the solid line in Fig. 2)

$$E_j^p=_0^\epsilon \tau (x)𝑑x=\underset{m=0}{\overset{n}{}}_{\epsilon _m}^{\epsilon _{m+1}}\tau _m(x)𝑑x=\underset{m=0}{\overset{n}{}}\left(1\frac{m}{N}\right)\left(\frac{\epsilon _{m+1}^2}{2}\frac{\epsilon _m^2}{2}\right)$$

where from Eq. (15) we have used $`\tau _m(x)=x(1m/N)`$ and where by convention $`\epsilon _{n+1}=\epsilon `$ is the final applied strain. A direct recursion gives exactly

$$E_j^p=\left(1\frac{n}{N}\right)\frac{\epsilon ^2}{2}+\underset{m=1}{\overset{n}{}}\frac{\epsilon _m^2}{2N}.$$
(16)
The first term here is the elastic energy that can be reversibly recovered upon decreasing the strain while the second term represents the energy irreversibly consumed in the breaking process. Both contributions can be directly vizualized in Fig. 2.

Equation (16) is next averaged over the quenched disorder to obtain the average energy $`E_j`$ needed to create state $`j`$. Each breaking threshold $`\epsilon _m`$ is an independant variable, randomly distributed according to $`p(\epsilon _m)`$ under the condition that $`0\epsilon _m\epsilon `$. They are therefore distributed according to the probability density $`p(\epsilon _m)/P(\epsilon )`$ where the normalization factor accounts for the fact that the upper limit $`\epsilon `$ is independent of the threshold values $`\epsilon _m`$. Thus, averaging over the quenched disorder gives

$$h(\epsilon )=\frac{\epsilon _m^2}{2}_{\mathrm{q}.\mathrm{d}.}=\frac{1}{P(\epsilon )}_0^\epsilon \frac{x^2}{2}p(x)𝑑x$$
(17)
where $`h(\epsilon )`$ designates the average energy that is lost when each fiber breaks.

The average work density (Hamiltonian) required to create the state $`j`$, averaged over all realizations of the quenched disorder is then in general

$$E_j=\left(1\frac{n_j}{N}\right)\frac{\epsilon ^2}{2}+\frac{n_j}{N}h(\epsilon )=\frac{\epsilon ^2}{2}\frac{n_j}{N}\left[\frac{\epsilon ^2}{2}h(\epsilon )\right].$$
(18)
For example, under the special assumption that the breaking strengths are randomly sampled from a uniform distribution on $`0\epsilon _m1`$, we have that $`p(\epsilon )=1`$, $`P(\epsilon )=\epsilon `$, $`h(\epsilon )=\epsilon ^2/6`$, and $`E_j=\epsilon ^2/2(n_j/N)\epsilon ^2/3`$.

## 4 Average Properties

We now apply the above theory and determine the thermodynamic functions as a function of applied strain $`\epsilon `$. The following analysis is valid for any properly normalized quenched-disorder distribution $`p(\epsilon )`$.

The cumulative distribution $`P(\epsilon )=_0^\epsilon p(x)𝑑x`$ is the probability that any one fiber has broken when the strain is at $`\epsilon `$. Thus, the fraction of all possible realizations having $`n_j`$ broken fibers and $`Nn_j`$ unbroken fibers is exactly

$$p_j^{\text{exact}}=(1P)^{Nn_j}P^{n_j}.$$
(19)
This distribution may be written

$$p_j^{\text{exact}}=p_0\mathrm{exp}\left[n_j\mathrm{ln}\left(\frac{1P}{P}\right)\right]$$
(20)
where $`p_0=(1P)^N`$ is the probability of the unbroken state $`j=0`$. Upon using the Hamiltonian of Eq. (18), the postulate of entropy maximization predicts this same distribution to be given by

$$p_j=p_0\mathrm{exp}\left[n_j\frac{\beta }{N}\left(\frac{\epsilon ^2}{2}h\right)\right]$$
(21)
where $`p_0=\mathrm{exp}(\beta \epsilon ^2/2)/Z`$ is the probability of the unbroken state.

These two distributions are both Boltzmannians in the number $`n_j`$ of broken fibers and are identical if the temperature $`T=1/\beta `$ is given by

$$\beta (\epsilon )=\frac{N}{\epsilon ^2/2h(\epsilon )}\mathrm{ln}\left(\frac{1P(\epsilon )}{P(\epsilon )}\right).$$
(22)
If it can be shown that this $`\beta `$ satisfies the differential equation of Eq. (13), then our theory is exact when applied to fiber bundles with global load sharing.

From the definition $`P(\epsilon )h(\epsilon )=_0^\epsilon p(x)x^2/2𝑑x`$, one has $`h(\epsilon )<\epsilon ^2/2`$ for any distribution $`p(\epsilon )`$. Thus, $`\beta (\epsilon )`$ is a negative increasing function up to the strain point $`\epsilon =\epsilon _\beta `$ where it smoothly goes to zero. The inflection point $`\epsilon _\beta `$ is obtained from the condition that $`P(\epsilon _\beta )=1/2`$ and defines a previously unnoticed phase transition that will be shown to be distinct from the transition at peak stress. When $`\epsilon >\epsilon _\beta `$, $`\beta `$ becomes a positive increasing function of $`\epsilon `$.

There are two key average properties upon which all the thermodynamic functions depend; namely, the average fraction of broken fibers in each bundle $`n_j/N`$ and the average of this fraction squared $`(n_j/N)^2`$. Using the exact probabilities of Eq. (20) \[which is equivalent to using the $`\beta `$ of Eq. (22) in our probability law\], one obtains

$$\frac{n_j}{N}=\underset{j}{}p_j\frac{n_j}{N}=\underset{n=0}{\overset{N}{}}c_n^N\frac{n}{N}P^n(1P)^{Nn}$$
(23)
where $`c_n^N=N!/[n!(Nn)!]`$ defines the number of ways of selecting $`n`$ objects from a group of $`N`$ distinguishable items. The binomial theorem states that

$$(x+y)^N=\underset{n=0}{\overset{N}{}}c_n^Nx^ny^{Nn}.$$
(24)
Upon differentiating this equation with respect to $`x`$ and then multiplying by $`x/N`$, gives that when $`x=P`$ and $`y=1P`$

$$\frac{n_j}{N}=P$$
(25)
which is a known result consistent with the meaning of $`P`$. Differentiating a second time yields

$$\left(\frac{n_j}{N}\right)^2=P^2+\frac{P(1P)}{N}.$$
(26)
Using these two results, the other averages defining the thermodynamic variables are easily read off.

The average stress $`\tau (\epsilon )`$ is thus obtained to be

$$\tau =\underset{j}{}p_j\tau _j=(1P)\epsilon $$
(27)
which is initially an increasing function of $`\epsilon `$ having the slope

$$\frac{d\tau }{d\epsilon }=1P\epsilon p.$$
(28)
This slope goes to zero and the stress has a maximum at any strain point $`\epsilon _\tau `$ that is a solution of $`1P(\epsilon _\tau )\epsilon _\tau p(\epsilon _\tau )=0`$; i.e., at the point(s) where

$$\epsilon _\tau p(\epsilon _\tau )=1_0^{\epsilon _\tau }p(x)𝑑x$$
(29)
admits a solution. This is a known exact result . In general, we can expect that $`\epsilon _\tau \epsilon _\beta `$. The condition required for equality of these strain points is that simultaneously $`\epsilon p(\epsilon )=1/2`$ and $`_0^\epsilon p(x)𝑑x=1/2`$ which for monotonic distribution functions $`p(\epsilon )`$ can only occur with the uniform distribution $`p(\epsilon )=1`$ in which case $`\epsilon _\tau =\epsilon _\beta =1/2`$ and the two transitions coincide. If $`p(\epsilon )`$ is a monotonic increasing function of $`\epsilon `$, then $`\epsilon _\tau <\epsilon _\beta `$ while if it is a decreasing function of $`\epsilon `$, then $`\epsilon _\tau >\epsilon _\beta `$. For non-monotonic distributions, there can be an arbitrary number of stress maximas and Eq. (29) can have either no solutions or multiple solutions  . The nature of the phase transitions at the distinct strain points $`\epsilon _\beta `$ and $`\epsilon _\tau `$ is discussed in the following section.

The average energy in the ensemble is

$$U=\underset{j}{}p_jE_j=(1P)\frac{\epsilon ^2}{2}+hP.$$
(30)
Again recalling the definition of $`h`$ from Eq. (17) gives

$$\frac{dU}{d\epsilon }=(1P)\epsilon p\frac{\epsilon ^2}{2}+\frac{d(hP)}{d\epsilon }=\tau $$
(31)
which is also the equation that gives the differential equation for temperature. This is sufficient for demonstrating that the $`\beta (\epsilon )`$ of Eq. (22) satisfies its differential equation. Nonetheless, as a consistency test, the coefficients $`a`$, $`b`$, and $`c`$ defined in Eq. (14) will be derived and the differential equation will explicitly be shown to be satisfied.

To obtain the state function $`f`$, we need first

$$\frac{dE_j}{d\epsilon }=\epsilon \frac{n_j}{N}\left[\epsilon \frac{p}{P}\left(\frac{\epsilon ^2}{2}h\right)\right].$$
(32)
From Eq. (4) and the lemma of Eq. (25) we then have

$$f=\underset{j}{}p_j\frac{dE_j}{d\epsilon }=(1P)\epsilon +p\left(\frac{\epsilon ^2}{2}h\right).$$
(33)
The variation of the entropy is obtained from the energy balance as

$$\frac{dS}{d\epsilon }=\beta \left(\frac{dU}{d\epsilon }f\right)=Np\mathrm{ln}\left(\frac{1P}{P}\right).$$
(34)
Together with the initial condition $`S(0)=0`$, this is readily integrated to give

$$S=N\left[P\mathrm{ln}P+(1P)\mathrm{ln}(1P)\right].$$
(35)
This expression is the classical Shannon result for a set of $`N`$ random variables in two possible states having probabilities $`P`$ and $`1P`$, which is precisely the case of the fiber bundle with global sharing. This is another consistency check. This entropy $`S`$ increases from zero \[total certainty that every member is intact\] and goes smoothly through a maximum \[total uncertainty as to what state a member may be in\] at the same strain point $`P(\epsilon _\beta )=1/2`$ where $`\beta `$ goes to zero. After the smooth maximum, $`S`$ decreases to reach zero if $`P(\epsilon )`$ reaches $`1`$ in which case there is total certainty that each member is entirely broken. We note that Eq. (35) can also be directly obtained from the Shannon formula upon applying the binomial theorem.

Finally, the free energy $`F`$ is again obtained from its Legendre transform definition $`F=US/\beta `$ to be

$$F=\frac{\epsilon ^2}{2}\left(\frac{\epsilon ^2}{2}h\right)\frac{\mathrm{ln}(1P)}{\mathrm{ln}(1P)\mathrm{ln}P}.$$
(36)
At the point $`\epsilon _\beta `$ defined by $`P(\epsilon _\beta )=1/2`$, the free energy diverges due to the fact that $`\beta (\epsilon _\beta )=0`$ while $`S(\epsilon _\beta )`$ remains finite. So long as $`\epsilon _\beta \epsilon _\tau `$, the free energy does not diverge when (if) $`\epsilon =\epsilon _\tau `$. In passing, we also note that

$$Z=\left[(1P)^hP^{\epsilon ^2/2}\right]^{N/(\epsilon ^2/2h)}$$
(37)
is the exact expression of the partition function.

With the above results established, we now obtain the coefficients $`a,b,c`$ of the differential equation for the temperature as

$`a`$ $`=`$ $`{\displaystyle \frac{1P}{N}}\left(Ph^2{\displaystyle \frac{\epsilon ^2}{2}}\right)^2;b={\displaystyle \frac{1P}{N}}(Ph^2{\displaystyle \frac{\epsilon ^2}{2}})\left(P\epsilon p{\displaystyle \frac{\epsilon ^2}{2}}+ph\right);`$
$`c`$ $`=`$ $`p\left(h{\displaystyle \frac{\epsilon ^2}{2}}\right).`$ (38)
Using these together with Eq. (22) for $`\beta `$ and its derivative

$$\frac{d\beta }{d\epsilon }=\frac{Np}{P(1P)(h\epsilon ^2/2)}+N\mathrm{ln}(\frac{1P}{P})\left[\frac{\epsilon +(h\epsilon ^2/2)p/P}{(h\epsilon ^2/2)^2}\right]$$
(39)
shows that the differential equation $`ad\beta /d\epsilon +b\beta +c=0`$ is exactly satisfied.

## 5 Phase Transitions

### 5.1 The entropy maximum transition

The strain point $`\epsilon _\beta `$ defined by $`P(\epsilon _\beta )=1/2`$ is where simultaneously $`\beta =0`$, the entropy is a maximum, and the free energy diverges. It is distinct from the stress-maximum transition(s) at $`\epsilon _\tau `$. The interpretation of $`\epsilon =\epsilon _\beta `$ as a phase transition is natural, since the most probable configuration of a bundle suddenly jumps from being entirely intact to entirely broken. What are the measurable manifestations of this transition at $`\epsilon _\beta `$?

Since the entropy is a maximum at this point, the fluctuations between realizations should also be a maximum as we now demonstrate. Define the fraction of broken fibers in state $`j`$ to be $`\rho _j=n_j/N`$, and define the average of this fraction to be $`\rho =\rho _j=P`$ where $`P=P(\epsilon )`$ is again the probability of a fiber being broken. The quantity of interest here is the root-mean-square fluctuation $`\mathrm{\Delta }\rho `$ in the fraction of broken fibers given by

$$\mathrm{\Delta }\rho =\sqrt{\mathrm{\Delta }\rho _j^2}=\sqrt{\rho _j^2\rho ^2}=\frac{\sqrt{P(1P)}}{\sqrt{N}}$$
(40)
which indeed goes through a maximum at $`P(\epsilon _\beta )=1/2`$ as expected. This maximum is something that can be directly measured in numerical experiments on fiber bundles but has never before been commented on. The reason it has been discovered in the present theory is because entropy and temperature are explicitly present. We encourage someone to numerically measure $`\mathrm{\Delta }\rho `$ and to verify that it is a maximum at the transition point $`\epsilon _\beta `$. Recall that for monotonic quenched-disorder distributions, if $`p(\epsilon )`$ is a decreasing function (more weak fibers than strong fibers), then $`\epsilon _\beta <\epsilon _\tau `$. So the numerical observation in this case should be that $`\mathrm{\Delta }\rho `$ goes through a maximum prior to peak stress. Equivalent comments hold when $`p(\epsilon )`$ is an increasing function and $`\epsilon _\beta >\epsilon _\tau `$.

The other fluctuations that are potentially of interest include the root-mean-square stress fluctuation

$$\mathrm{\Delta }\tau =\sqrt{\mathrm{\Delta }\tau _j^2}=\sqrt{\tau _j^2\tau ^2}=\epsilon \mathrm{\Delta }\rho $$
(41)
and the root-mean-square energy fluctuation

$$\mathrm{\Delta }U=\sqrt{\mathrm{\Delta }E_j^2}=\sqrt{E_j^2U^2}=\left(\frac{\epsilon ^2}{2}h\right)\mathrm{\Delta }\rho ,$$
(42)
but since these are simply proportional to $`\mathrm{\Delta }\rho `$ it seems that the interesting signature of this phase transition is the maximum in $`\mathrm{\Delta }\rho `$.

The “order” of this transition is not classically defined in the Ehrenfest scheme; however, it seems inappropriate to call it a continuous transition because the free energy is singular at $`\epsilon _\beta `$ as shown above. Nonetheless, in the limit as $`\epsilon \epsilon _\beta `$ the singular part of $`F`$ diverges according to the scaling law

$$F_s=\frac{\mathrm{ln}2}{8}\frac{[\epsilon _\beta ^2h(\epsilon _\beta )/2]}{p(\epsilon _\beta )}(\epsilon \epsilon _\beta )^1$$
(43)
which has a (trivial) universal exponent. Since $`F`$ and its derivatives are difficult to numerically measure, the principle manifestation of this phase transition remains $`\mathrm{\Delta }\rho `$ going through a maximum at a point $`\epsilon _\beta \epsilon _\tau `$.

### 5.2 The stress maximum transition

The phase transition at peak stress $`\epsilon =\epsilon _\tau `$ is the one that researchers up to now have focused on. Since $`S`$ and its derivatives remain continuous and finite there, an attempt to classify this transition as being first-order or second-order is meaningless. The label of “continuous” transition seems the most appropriate.

There is universality at this transition due to the fact that for any monotonic analytic quenched-disorder distribution, $`\tau (\epsilon )`$ is a smooth analytic function so that upon developing this function in the neighborhood of its maximum one is guaranteed $`|\tau \tau (\epsilon _\tau )||\epsilon \epsilon _\tau |^2`$ with the exponent of course being independent of the distribution $`p(\epsilon )`$ or other model parameters. Such development of an analytic function about a critical point is the way any mean-field theory acquires a universal scaling law. However, Kloster, Hansen and Hemmer  do demonstrate that a non-analytic quenched-disorder distribution can lead to highly non-analytic stress-strain behavior that thus falls outside the quadratic universality class. Also, in other damage models in which elastic interaction between cracks (damage points) is important, the simple observation of a model exhibiting an averaged stress maximum does not by itself guarantee an exactly quadratic stress-strain relation in the neighborhood of the maximum. Such a relation may be non-analytic due to a diverging correlation length in the correlation between cracks.

The quadratic stress maximum means that $`d\epsilon /d\tau `$ diverges as $`|\tau \tau (\epsilon _\tau )|^{1/2}`$  . Accordingly, the average rate at which the fraction of broken fibers increases with stress $`d\rho /d\tau =pd\epsilon /d\tau `$ also diverges as $`|\tau \tau (\epsilon _\tau )|^{1/2}`$  . Using additional considerations not developed in this paper, one can futher show that the average size of avalanches also diverges as $`|\tau \tau (\epsilon _\tau )|^{1/2}`$ , and that the distribution $`D(\mathrm{\Delta })`$ of the size $`\mathrm{\Delta }`$ of the avalanches scales as $`D(\mathrm{\Delta })\mathrm{\Delta }^{5/2}`$  at the stress maximum. These are the principal observable characteristics of the stress-maximum transition.

## 6 Summary and Conclusions

Two principal results have been obtained in this paper. First, it was demonstrated that for any quenched-disorder distribution used in the fiber-bundle model with global load sharing, the probability of the emergent damage states may be exactly calculated by maximizing Shannon’s entropy subject to constraints. It is through the constraints that the nature of the quenched-disorder distribution enters the macroscopic thermodynamics.

Second, a previously unnoticed phase transition occurs in the fiber-bundle model with global load sharing when the entropy goes through a maximum. This phase transition is distinct from the stress-maximum transition; although, for the uniform quenched-disorder distribution, the two transitions coincide. The principle manifestation of this transition is that the root-mean-square fluctuation in the number of broken fibers will go through a maximum, which is a quantity that can be directly measured in numerical experiments.

To conclude, we postulate that the probability of emergent states can always be calculated through entropy maximization for any damage model of interest including models in which there is elastic interaction between the local damage (order) parameters. Damage models with order-parameter interactions are not normally considered amenable to analytical treament; however, using our approach a rather standard canonical-ensemble partition function emerges and the various functional integration procedures available for studying the partition function in the neighborhood of critical points may be employed.

Such generality is the principal utility of our approach. A concern is the ability to produce an expression for the model temperature. The temperature can be found in principle as the solution of a well-posed initial-value problem. Unfortunately, the differential problem is highly non-linear and thus difficult to solve. Because stress maxima and entropy maxima do not normally coincide, a stress-maximum transition may be studied by simply assuming the temperature to be well behaved in the neighborhood of the stress maxima. But if an explicit expression for the temperature is desired, the following approach can be employed.

In models involving crack interactions, there are always a subset of dilute states in which the interactions are negligible. The exact probabilities of such states can usually be determined from the quenched-disorder distribution alone. A comparison between the probabilities calculated from entropy maximization and such exact probabilities then determines the model temperature to be used in the Boltzmannian. With the temperature so defined, and with a proper model Hamiltonian in hand, the partition function can be analyzed.


text:
# RANDOM FIELD MODELS FOR RELAXOR FERROELECTRIC BEHAVIOR

## I INTRODUCTION

Shortly after the seminal work of Imry and MaIM75  on the effects of random fields in ferromagnets, it was pointed out by Halperin and VarmaHV76  that similar ideas could be used to understand the effects of atomic disorder in perovskite ferroelectrics. Halperin and Varma showed that atomic disorder which coupled linearly to the ferroelectric order parameter behaved in an equivalent fashion to the random field in a ferromagnet, and thus could cause the ferroelectric phase transition to become smeared. They seem to have been unaware, however, that ferroelectrics displaying such smeared transitions had been discovered long before.SI54  These materials, now called relaxor ferroelectrics, were originally referred to as ferroelectrics with diffuse phase transitions. The random-field ideas of Halperin and Varma did not immediately become popular in the field. In the 1984 review of IsupovSmo84  the diffuse phase transitions were modeled using only random-bond disorder. However, by the early 1990’s it was generally recognized that random-bond disorder was not enough to explain all of the observed effects, and models including both random bonds and random fields became widely used.WKG92 ; Kle93 ; PTB87 ; PB99  Until now it has been assumedPTB87 ; PB99  that in the presence of a strong cubic crystal-field anisotropy one could treat the different Ising components independently. In this work we will show that for a certain type of random field new phenomena occur which cannot be understood without including coupling between the components.

The identification of relaxor ferroelectrics as a broad class of materials in which electric dipoles behave in a manner analogous in many respects to the magnetic dipoles in spin-glasses was made forcefully by Burns and Dacol,BD83  who generalized and extended ideas of Courtens.Cou82  It is essential, however, to note that, as shown by Halperin and Varma,HV76  in the case of electric dipoles it is almost inevitable that the disorder which is an essential feature of these materials will produce random fields. The presence of these random fields limits the development of long-range spin-glass order in three dimensions.HFA79  Therefore, it is not expected that the glassy behavior which is seen in the disordered electric dipole materials represents a true phase transition into a state with long-range spin-glass-type order.

We would like to study the development of spatial correlations among the electric dipoles. In order to do this, we will work with a model which does not retain all of the atomic details. In the model studied here, we crudely represent each perovskite unit cell by only one dynamical variable. We will refer to this variable as a spin, although it actually represents atomic displacements in the unit cell. The model includes a cubic crystal anisotropy, a nearest-neighbor two-spin interaction, and a random field.

As discussed by Pirc and Blinc,PB99  if we assume that the cubic crystal anisotropy is so strong that each spin points along one of the  directions and the random field is represented by a simple vector coupling linearly to the spin, then the different components of the spin act independently.Suz67  Thus the model reduces to a set of decoupled random-field Ising models. We may also include randomness in the bonds. (To perform this decoupling on Eq. (1), we use a coordinate system rotated by 45.) This model has been studied in detail over a number of years. It is believed to be generally well-understood, although there remains some controversy about the values of the critical exponents.Kle02 ; Ye02

Similar models have been studied in which the spin components do not act independently. Then the behavior does not reduce to that of an Ising model. For example, there are the three-state Potts model with a random field,EB96  and the four-state Potts model with a random field,QB96  which are still qualitatively similar to the random-field Ising model. When the random field becomes strong enough, the ground state of the system breaks up into Imry-Ma domains, and the long-range order is destroyed.

There is also the cubic model with random anisotropy,Fis93  which, when the randomness is strong, shows a new type of domain phase with long-range order. This phase may be thought of as one in which each of the spin components orders independently, but different parts of the lattice are dominated by different components. Thus, in this phase, there is a network of 90 domain walls embedded in the ordered phase. A somewhat similar domain state was found by Kartha, Castan, Krumhansl and SethnaKCKS91  in a model of shape-memory alloys.

The existence of such a domain state with long-range order has been proposedBDKB95 ; BKB96  as the explanation of relaxor ferroelectric behavior. Further, EgamiEg00  has emphasized that the network of 90 domain walls is probably an essential element leading to the large piezoelectric response of relaxor ferroelectrics. Therefore, we would like to study a random-field model which has such a phase.

## II THE MODEL

The simplest model which has the desired properties is a four-state vector Potts model (four-state clock model) with a random Potts field, on a simple cubic lattice. The Hamiltonian of this model is

$$H=J\underset{ij}{}𝐒_i𝐒_j+h_r\underset{i}{}\delta _{𝐒_i,𝐧_i}.$$
(1)
Each spin $`𝐒_i`$ is a dynamical variable which has four allowed states: (1,0), (0,1), (-1,0) and (0,-1). Each $`𝐧_i`$ is an independent quenched random variable which assumes one of these four allowed states with equal probability. The $`ij`$ indicates a sum over nearest neighbors of a simple cubic lattice. The extension to three-component spins is conceptually straightforward, but will not be discussed in detail here.

The sign of $`h_r`$ is chosen so that for large positive values the probability of spin $`𝐒_i`$ in state $`𝐧_i`$ is strongly suppressed. Density-functional calculations for a typical relaxor ferroelectric alloyCGMR02  have shown that the most significant interactions created by the alloy disorder are repulsive and short-ranged. Thus we believe that this is a reasonable first approximation to the interactions produced by the disorder.

When $`h_r`$ is negative, the random Potts field favors one state on each site, just as a vector field would. In that case there is no qualitative difference between the random Potts field and the random vector field.GH96 ; Fis97  Thus we will not present any calculations for negative $`h_r`$ here.

When $`h_r`$ is large and positive, we can approximate the random field term as a projection operator which forbids $`𝐒_i`$ from occupying the state $`𝐧_i`$. Then within mean-field theory it becomes simple to calculate the approximate ground states. There are eight of these mean-field ground states. For example, in one such ground state 75% of the spins are in state (1,0), and the other 25% of the spins, which are prevented from being in this spin state by their local $`𝐧_i`$ are in state (0,-1). The energy of this state for a simple cubic lattice is easily calculated to be $`1.875J`$ per spin. Additional details of the mean-field theory are discussed in the Appendix.

As we will see, the energies of the ground states found by computer calculation are less than $`2J`$, and they differ qualitatively from the mean-field ground states. The problem with the mean-field ground states is easy to understand. If we consider the above example, on a simple cubic lattice the 25% of spins in the state (0,-1) are broken up into finite clusters. If we move all of the spins in such a finite cluster into the state (0,1), the total energy does not change. One might expect, therefore, that for large $`h_r`$, the long-range order in the ground state would still be along one of the four  directions, since the 25% of the spins in the finite clusters should not be able to exhibit long-range order.

However, a true ground state on the simple cubic lattice actually has a very different structure. It turns out to be possible to find states with, for example, 50% of the spins in state (1,0) and 50% of the spins in state (0,-1) which have a much lower energy than the mean-field ground states. As we will show, to find states which look like the true ground states for large $`h_r`$ within mean-field theory, we need to add another term to Eq. (1).

The author finds it helpful to compare the model studied here to the results for the three dimensional Ashkin-Teller model,DBGK80  whose Hamiltonian is

$$H_{AT}=J\underset{ij}{}(𝝈_i𝝈_j+𝝉_i𝝉_j)2J_4\underset{ij}{}(𝝈_i𝝈_j)(𝝉_i𝝉_j),$$
(2)
where $`𝝈=\pm \frac{1}{2}(1,1)`$ and $`𝝉=\pm \frac{1}{2}(1,1)`$ are Ising variables. (Note that our notation and normalization differ from those of Ditzian et al.DBGK80 )

Each of the three linearly independent components of a random field of a four-state model can now be identified with one of the terms of Eq. (2), and thus one of the mean-field order parameters. Expressing the random-field term using the $`𝝈`$ and $`𝝉`$ variables, and adding it to Eq. (2), we get

$$H_{RFAT}=H_{AT}+(h_r/4)\underset{i}{}[1+2(𝝈_i𝐧_i)+2(𝝉_i𝐧_i)+4\eta (𝝈_i𝐧_i)(𝝉_i𝐧_i)].$$
(3)
If we set $`\eta =1`$ and $`J_4=0`$, then Eq. (3) becomes just Eq. (1) expressed in the coordinate system defined by $`𝝈`$ and $`𝝉`$. If we then set $`\eta =0`$, the random Potts field term is replaced by the usual random vector field term.

In a real alloy we would not expect that exactly one of the local states was blocked on each site. Instead, we would expect to find a random distribution, with some sites having no blocked states, some with one blocked state, etc. It would be straightforward to generalize the random field term to allow this, by letting $`\eta `$ in Eq. (3) become a function of i. For example, we could use a probability distribution for $`\eta _i`$ which had weight at $`\pm 1`$. Such a generalization is likely to be useful in modeling the properties of specific alloy systems. However, a more general model has more parameters. Exploring its entire parameter space would be tedious, and seems premature at this time.

## III PHASE DIAGRAM

If $`h_r/J`$ is chosen to be an integer, then the energy of any state is an integer multiple of $`J`$. Then it becomes straightforward to design a heat bath Monte Carlo computer program to study Eq. (1) which uses integer arithmetic for the energies, and a look-up table to calculate Boltzmann factors. This substantially improves the performance of the computer program, and was used for all the calculations reported here. (If desired, one could do almost as well for half-integer values of $`h_r/J`$.) Lattices with periodic boundary conditions were used throughout. Two different linear congruential random number generators were used: one to choose the random $`𝐧_i`$ and one to choose the Boltzmann factors for the Monte Carlo spin flips.

A series of $`L\times L\times L`$ samples of various lattice sizes and values of $`h_r`$ was studied to map out the phase diagram, which is shown in Fig. 1. For $`h_r/J3`$ the behavior of the random Potts field model remains generally similar to the behavior of the random vector field model. At low temperature the system develops long-range order, with the order parameter aligned along one of the spin-coordinate axes. It is not required, however, that the critical exponents for the phase transition here will be those of the random-field Ising model, as happens for the vector random-field case.

The $`𝝈`$ phase found by Ditzian, Banavar, Grest and KadanoffDBGK80  has many of the features of the domain FE phase which we find in our model for large $`h_r`$ and low $`T`$ . The $`𝝈`$ phase of Eq. 2 only occurs for $`J_4<0`$. This indicates that one of the effects of the random Potts field in Eq. (1) is to generate a negative effective value of the $`J_4`$ coupling. This effect was also seen in the random-anisotropy cubic model.Fis93  It is the random-anisotropy part (the $`\eta `$ term of Eq. (3)) which generates the $`J_4`$ effective coupling under rescaling, and prevents the decoupling of Eq. (1) into two Ising models. If we want to study a mean-field theory for Eq. (1) we should include a $`J_4`$ term, thus arriving at Eq. (3). The effective value of $`J_4`$ will depend on $`h_r`$. In this way we can obtain a mean-field theory which reproduces the phase diagram of Fig. 1.

A necessary condition for the existence of stable long-range  domain FE order is that there be infinite connected clusters of both of the two dominant spin states. Since it becomes impossible to fulfill this condition when the amount of long-range order is small, the transition from the domain FE phase to the paraelectric phase must be first order. This, however, is specific to our model with only first neighbor interactions on the simple cubic lattice. It is possible that this phase transition could become continuous if, for example, second neighbor exchange is included in the Hamiltonian.

As discussed by Ditzian et al.DBGK80 , even in the absence of randomness the details of what happens as we move from the $`𝝈`$ phase to the Baxter phase are not clear. There may be a small area in the region of the  FE to  domain FE phase boundary in which an intermediate phase is stable. This would be similar to the intermediate phase which has recently been found experimentally.NCSGCP99  We have indicated the uncertainty about what is going on in this region of the phase diagram by the question mark in Fig. 1. A generalization of the $`𝝈`$ phase also exists for the three-component Ashkin-Teller model in three dimensions.GW81

The dotted line in Fig. 1 represents the approximate location of the onset of thermal hysteresis. Above the dotted line, the system will quickly relax to a state which is independent of initial conditions, but below this line the system retains a memory of initial conditions for a long time. This line does not represent any true thermodynamic singularity. Its precise location depends somewhat on the size of $`L`$ and the length of the Monte Carlo runs. Thus the dotted line should be identified as the glass temperature, $`T_g`$. For small values of $`L`$ the system can develop order on the scale of $`L`$ above $`T_g`$. For $`L=64`$ and $`h_r3`$ it is not possible to find the ground states in a reasonable time by cooling the system from high temperature.

## IV SPECIFIC HEAT

In Fig. 2 we display the specific heat for Eq. (1) using four different values of $`h_r`$. These curves were calculated by numerically differentiating and averaging Monte Carlo data for the energy of runs for four $`L=64`$ lattices for each value of $`h_r`$, starting at high temperature and cooling slowly. For each run, the temperature was changed in steps of $`0.05J`$. 40,960 Monte Carlo passes through the lattice were performed at each value of $`T`$, with data being collected after each 20 passes. The first 248 data points at each temperature were discarded, and the remaining 1,800 data points were averaged. A similar procedure was used for heating runs.

Note that for $`h_r=0`$ (not displayed), the specific heat for Eq. (1) reduces to a that of an ordinary Ising model on the simple cubic lattice, with $`T_c/J=2.256`$. This model has a very sharp singularity at $`T_c`$, whose height diverges to infinity with $`L`$. We see in Fig. 2 that as $`h_r`$ increases, the peak in the specific heat broadens and shifts to lower energy, with approximately half of the change occurring by $`h_r=3`$. By comparing with Fig. 1, it is seen that for $`h_r>0`$ the peak in the specific heat is centered above the region of glassy behavior. The average heights and widths of the peaks shown in Fig. 2 are essentially the same for $`L=32`$ (not shown) as for $`L=64`$. As one would expect, there is more variation from sample to sample for the smaller lattices.

In Fig. 3 we compare data from the cooling runs with $`h_r=6`$ with heating runs on the same lattices. The initial condition for each heating run was a state with long-range order. The temperature at which the long-range order collapses is identified as $`T_g`$. We see that below $`T_g`$ there is a region of $`T`$ where the specific heat is slightly higher for the cold-start sample. The crossing of the two curves near $`T/J=1.2`$ gives us a crude estimate for the ordering temperature $`T_c`$.

It is often possible to estimate $`T_c`$ for a first-order phase transition with high precision. In order to do this with a computer simulation, however, it is necessary to run the simulation close to $`T_c`$, with $`L`$ larger than the correlation length, for a time which is long enough so that the sample can spontaneously order and disorder several times. We are far from being able to satisfy that condition here, due to the glassy dynamics for $`T<T_g`$.

The integrated area between the cooling curve data and the heating curve data gives us an estimate of the value of the latent heat at $`T_c`$. We can also look directly at the difference in energy between the heating run and the cooling run at $`T=1.2J`$. Thus we find that the latent is about $`0.01J`$. The fact that the slope of the heating curve is so similar to that of the cooling curve and the latent heat at $`T_c`$ is so small indicates that there is very little difference in the local structure of this model between the ordered states of the  domain ferroelectric phase and the metastable states which are found by slow cooling. A small latent heat is a natural consequence of a large (but finite) correlation length at $`T_c`$.

## V CORRELATION FUNCTIONS

In order to obtain information about the two-spin correlations, we calculate the angle-averaged structure factor $`S(k)`$. In Fig. 4 we display a log-log plot of $`S(k)`$, obtained by averaging over the data from the cooling runs for the $`L=64`$ lattices, as a function of $`T/J`$. We see that for all these values of $`h_r3`$, the qualitative behavior of a sample cooled from high $`T`$ without any external ordering field is the same. As we lower $`T`$, the peak near $`k=0`$ first grows and then saturates. The rate of growth of this peak is maximal at the temperature where the specific heat has its maximum. As $`h_r`$ decreases, the height of the peak and the value of the correlation length increase. For $`h_r=3`$, the correlation length at low $`T`$ becomes larger than our sample size. However, these zero-field-cooled samples never show true long-range order. With true long-range order, the size of the small-$`k`$ peak would decrease as $`T`$ decreases below $`T_c`$, because the intensity in the long-range-order $`\delta `$-function will not appear on the log-log plot.

In Fig. 5 we show data for $`S(k)`$ from the heating runs, where the samples begin in an ordered state. Under these conditions, the size of the small-$`k`$ peak continues to grow as $`T`$ is increased above $`T_c`$, until $`T`$ reaches $`T_g`$ and the sample can equilibrate. In Fig. 5(c) we show data from runs which are ten times longer than the standard runs. The $`h_r=4`$ samples are able to equilibrate at $`T/J=1.7`$ during the long run (i.e. they become indistinguishable from the data from the corresponding cooling run), but the changes at $`T/J=1.5`$ are small. By $`T/J=1.4`$ increasing the length of the run by a factor of ten has a negligible effect on the observed state of the sample.

## VI GROUND STATES

Somewhat surprisingly, it turned out to be possible to find approximate ground states of the $`h_r=\mathrm{}`$ samples by a simulated annealing procedure. In Fig. 6 we show $`S(k)`$ for the approximate ground states of the $`L=64`$ lattices with $`h_r=\mathrm{}`$, averaging over each of the four  directions for each of the four samples. The data shown in the figure are fit quite well with by a Lorentzian line shape, with a correlation length of about 10 lattice units. This diffuse-scattering part of $`S(k)`$ contains about 54% of the spectral weight. The other 46% of the weight is in the $`\delta `$-function at $`k=0`$, which is due to the long-range  order.

What this means in detail is that for large $`h_r`$ a ground state contains approximately 49% each of the two majority spin states, and about 1% each of the two minority states. The energy per spin of the lowest ground states found for each sample was about $`2.010J`$, with the energies of the ”ground states” of the same sample in the other three  directions lying typically about $`0.003J`$ higher. Since $`L=64`$ is large compared to a correlation length of 10, this result for the ground state energy should represent the large $`L`$ limit, while the differences in energies between alternate ground state directions of the same sample should scale like $`1/L^{3/2}`$. Note that these actual ground state energies are substantially below the energy of the simple mean-field ground state of the $`h_r=\mathrm{}`$ model.

For $`L=64`$ and $`h_r=\mathrm{}`$ the energy of a zero-field-cooled state at $`T=0`$, which has no long-range order, is only about $`0.01J`$ higher than the energy of a true ground state. Thus for lattices smaller than about $`L=32`$ the energy differences between the ground states of different  directions become larger than the energy difference between a zero-field-cooled state and a true ground state. Under these conditions the zero-field-cooled sample is able to reach a long-range-ordered state. The zero-field-cooled $`L=64`$ samples can be decomposed into domains of different  ground states, with 90 walls between the domains.

The probabilities of having effective fields of magnitude 6, 5, … and 0 in a ground state, averaged over the four ground states for each of these $`L=64`$ $`h_r=\mathrm{}`$ samples are 0.1121, 0.2780, 0.3148, 0.1685, 0.0751, 0.0383, and 0.0131, respectively. Thus, in a ground state only about one spin in nine is surrounded by fully aligned nearest neighbors. More than 1% of the spins are in zero effective field.

Some of the spins can flip freely between two positions in a ground state with no cost in energy. Since the spins in this model are discrete and two-dimensional, being in zero effective field is not a necessary condition for a spin to be ”free” in this way. The bulk of these free spins are flipping between the two majority spin states for that particular ground state. It is easy to find that the fraction of free spins in a particular ground state is about 6.5%. However, a spin which is not free in one ground state may become free by the flipping of other free spins. The fraction of spins which can become free by the flipping of other free spins, one at a time, is 11.9%. Since 0.119 is much less than the critical concentration for uncorrelated site percolation on a simple cubic lattice, it is not surprising that a set of free spins defined in this way consists of small isolated clusters.

The probability that a free spin may occupy all three of its allowed spin states with no energy cost is negligible to this level of accuracy. Thus a lower bound on the ground state entropy per spin for large $`h_r`$ is $`0.065ln(2)`$, and $`0.119ln(2)`$ is an upper bound. Calculating the ground state entropy of this model is a complex problem, because it requires finding all of the ground states.

Thus each of the four ” ground states” of a sample is actually a cluster of ground states in the phase space of the model. The simulated annealing procedure works because, if we use one of the simple mean-field ground states, which has 75% of the spins in one spin state and 25% of the spins in a second spin state, as an initial condition, the direction of the order parameter will rapidly relax to the closest of the  energy minima, even at temperatures well below $`T_c`$. The fact that this is possible is partly due to the large number of spin rearrangements which can be made with no energy cost. This prevents the system from being easily trapped in a metastable minimum of the free energy which is close to the initial mean-field ground state.

Although a ground state energy of $`2.01J`$ means that about 33% of nearest neighbor spin pairs are oriented at 90 from one another, the fraction of nearest neighbor spin pairs which are pointing in opposite directions in one of these ground states is only 0.00011. The small number of spins in the minority states exist in compact blobs whose diameter is approximately the correlation length of 10 lattice units.

For $`h_r/J`$ = 3, every sample studied had four ground states, each one fully aligned along one of the four  spin states. There is an exponentially small probability of having an unusual local configuration of the $`𝐧_i`$ which would misalign a small cluster of spins in a ground state, but this was never observed. Thus, ignoring the statistical fluctuations in the sample average of $`𝐧_i`$, the ground state energy in the  ferroelectric phase is found to be $`E_0=3J+0.25h_r`$. The value of $`E_0`$ for any value of $`h_r`$ cannot be greater than $`2.01J`$, its value for $`h_r=\mathrm{}`$. From this we expect that the boundary between the  phase and the  domain phase should be at $`h_r/J`$ slightly less than 4, as observed.

## VII DISCUSSION

In the past, various kinds of evidence have been presented to argue that relaxor ferroelectrics represent a distinct class of materials, and are not merely the inevitable consequence of adding some alloy disorder to any ferroelectric. De Yoreo, Pohl and BurnsDPB85  studied the low-temperature properties of a variety of ferroelectric alloys, and found that they could be separated into two classes. Relaxor ferroelectrics showed glassy low temperature behavior and normal ferroelectrics did not, even when the normal ferroelectrics were random alloys. Recently, Viehland et al.DXV94 ; XKLV96 ; VPCL01  have argued that in some materials they find a well-defined transition between a normal ferroelectric phase and a relaxor ferroelectric phase. It would be remarkable if the low temperature thermal properties of such a crystal could be switched from glassy behavior to crystalline behavior by poling the sample. To the author’s knowledge, this experiment has not yet been tried.

For the two-spin exchange interaction of Eq. (1), a 180 domain wall has twice the energy of a 90 domain wall. As a result of this, 180 domain walls become rare at low temperature in this model. As pointed out by Bürgel, Kleemann and Bianchi,BKB96  ferroelastic interactions raise the energy of 90 walls, but do not have much effect on the energy of 180 walls.

To incorporate this effect into our model, we would add a biquadratic term to Eq. (1) of the form $`K(𝐒_i𝐒_j)^2`$. When $`K/J`$ becomes large, the structure of the domain state will change,BKB96  because 180 domain walls now have a lower energy than 90 domain walls. It is difficult to study this Bürgel, Kleemann, Bianchi domain state with computer simulations if the size of a ferroelastic domain is large. A sample which is the size of a single ferroelastic domain (or smaller) will behave essentially as a random-field Ising model.

## VIII CONCLUSION

In this work we have used Monte Carlo computer simulations to study a simple model of a ferroelectric alloy. In this model, the alloy disorder causes some of the positions of the polarizable ions to be strongly disfavored. We have found that this model displays a type of ferroelectric domain phase which does not exist in the model of Pirc and Blinc.PB99  This phase may be thought of as resulting from the addition of a random field to the $`𝝈`$ phaseDBGK80  of an Ashkin-Teller model. It seems likely that this model, and some natural generalizations of it, will help in the understanding of relaxor ferroelectric behavior.

*

## Appendix A MEAN FIELD THEORY

The presence of the random field greatly complicates the mean-field theory for Eq. (1). It is necessary to distinguish the four classes of spins corresponding to the four different allowed values of $`𝐧`$. Label the four directions (1,0), (0,1), (-1,0) and (0,-1) as direction 1, 2, 3 and 4, respectively. Define $`p_{\mu \nu }`$ to be the probability that a spin of class $`\mu `$ points in the $`\nu `$ direction. Then $`_\nu p_{\mu \nu }=1`$ for any $`\mu `$. But, in the general case, we have to deal with twelve independent variables in the mean-field theory. If we take the limit that $`h_r`$ becomes infinite, there are only eight independent variables, because we can then set $`p_{\mu \nu }`$ to 0 when $`\mu =\nu `$.

Let $`z`$ be the number of neighbors of any spin. Then the mean-field expression for the energy per spin is

$$E=(zJ/32)\underset{\mu ,\nu ,\mu ^{},\nu ^{}}{}p_{\mu \nu }p_{\mu ^{}\nu ^{}}\mathrm{cos}(\pi (\nu \nu ^{})/2).$$
(4)
In addition to the usual factor of 1/2 to avoid double-counting, there is factor of 1/16 arising from the assumption that each value of $`\mu `$ (i.e. each $`𝐧`$) is equally weighted. The entropy per spin is

$$S=1/4\underset{\mu ,\nu }{}p_{\mu \nu }\mathrm{ln}(p_{\mu \nu }),$$
(5)
and the free energy is, as usual, $`F=ETS`$. Then solving the mean-field theory requires finding the minimum of $`F(T)`$ in the phase space of the $`p_{\mu \nu }`$.

###### Acknowledgements.

The author is grateful to T. Egami, I-W. Chen, S. Vakhrushev, A. B. Harris, I. Grinberg, and A. M. Rappe for helpful discussions. This work was supported by the Center for Piezoelectrics by Design under ONR Grant N00014-01-1-0365.


text:
# Counting curves via lattice paths in polygons

## 1. Introduction: the numbers $`N^{\mathrm{\Delta },\delta }`$

Let $`\mathrm{\Delta }^2`$ be a convex polygon with integer vertices. It defines a finite-dimensional linear system $``$ of curves in $`(^{})^2`$, where $`^{}=\{0\}`$. These curves are the zero loci in $`(^{})^2`$ of the (Laurent) polynomials

$$f(z,w)=\underset{(j,k)\mathrm{\Delta }^2}{}a_{jk}z^jw^k,$$

$`a_{jk}`$. The polynomials $`f`$ themselves form the vector space $``$. Recall that the Newton polygon of $`f`$ is $`\mathrm{Convexhull}\{(j,k)|a_{jk}0\}`$. Thus $``$ contains polynomials whose Newton polygon is contained in $`\mathrm{\Delta }`$. Clearly $``$ is a complex projective space of dimension

$$m=\mathrm{\#}(\mathrm{\Delta }^2)1.$$

Curves with the Newton polygon $`\mathrm{\Delta }`$ form an open dense set $`𝒰`$. A generic curve in $``$ is a smooth curve of genus

$$l=\mathrm{\#}(\mathrm{Int}\mathrm{\Delta }^2).$$

Let $`C`$ be a curve. Even if $`C`$ is not irreducible we can define its genus $`g(C)`$. Consider the decomposition $`C=C_1\mathrm{}C_n`$ into the irreducible components $`C_j`$. We define $`g(C)=\underset{j=1}{\overset{n}{}}g(C_j)+1n`$. Note that this definition of genus allows for negative values (cf. ). If $`C`$ is singular then its genus is strictly smaller than $`l`$.

The curves of genus $`l\delta `$ and with the Newton polygon $`\mathrm{\Delta }`$ form a subvariety $`\mathrm{\Sigma }_\delta ^{}𝒰`$ of dimension $`m\delta `$. Let $`\mathrm{\Sigma }_\delta `$ be the projective closure of $`\mathrm{\Sigma }_\delta ^{}`$. We define $`N^{\mathrm{\Delta },\delta }`$ to be the degree of the $`(m\delta )`$-dimensional subvariety $`\mathrm{\Sigma }_\delta `$ in $``$. The degree is the intersection number with a projective subspace of codimension $`m\delta `$. Curves from $``$ passing through a point $`z(^{})^2`$ form a hyperplane.

The number $`N^{\mathrm{\Delta },\delta }`$ has the following enumerative interpretation. Let $`z_1,\mathrm{},z_{m\delta }(^{})^2`$ be generic points. The number $`N^{\mathrm{\Delta },\delta }`$ equals to the number of algebraic curves of Newton polygon $`\mathrm{\Delta }`$ and genus $`l\delta `$ passing through $`z_1,\mathrm{},z_{m\delta }`$. Note that $`N^{\mathrm{\Delta },0}=1`$ for any $`\mathrm{\Delta }`$. These numbers get more interesting when $`\delta >0`$.

###### Remark 1.

Another way to look at the same problem is to consider the compactification of the torus $`(^{})^2`$. Recall that the polygon $`\mathrm{\Delta }`$ defines a compact toric surface $`T_\mathrm{\Delta }`$, see e.g. . (Some readers may be more familiar with the definition of toric surfaces by fans, in our case the fan is formed by the dual cones at the vertices of $`\mathrm{\Delta }`$.) The surface $`T_\mathrm{\Delta }`$ may have isolated singularities that correspond to some vertices of $`\mathrm{\Delta }`$.

In addition to a complex structure (which depends only on the dual fan) the polygon $`\mathrm{\Delta }`$ defines a holomorphic linear bundle $``$ over $`T_\mathrm{\Delta }`$. We have a canonical identification $`\mathrm{\Gamma }()=`$, where $`\mathrm{\Gamma }()`$ is the space of the sections of $``$. The projective space $``$ can be also considered as the space of all holomorphic curves in $`T_\mathrm{\Delta }`$ such that their homology class is Poincaré dual to $`c_1()`$. The number $`N^{\mathrm{\Delta },\delta }`$ is the number of holomorphic curves $`\overline{C}T_\mathrm{\Delta }`$ such that $`z_1,\mathrm{},z_{m_\delta }\overline{C}`$, the homology class $`[\overline{C}]`$ is dual to $`c_1()`$, the Euler characteristic of the normalization of $`\overline{C}`$ is $`22(l\delta )`$ and no irreducible component of $`\overline{C}`$ is contained in $`T_\mathrm{\Delta }(^{})^2`$.

###### Remark 2.

Note that in the set-up of Remark 1 the number $`N^{\mathrm{\Delta },\delta }`$ appears related to the Gromov-Witten invariant of $`T_\mathrm{\Delta }`$ (see ) corresponding to $`c_1()`$. The difference is that the corresponding Gromov-Witten invariant also has a contribution from the curves of genus $`l\delta `$ which have components contained in $`T_\mathrm{\Delta }(^{})^2`$. This contribution is zero (by the dimension reasons) if $`T_\mathrm{\Delta }`$ is smooth and does not have exceptional divisors. Thus if $`\mathrm{\Delta }=\mathrm{\Delta }_d=\mathrm{ConvexHull}\{(0,0),(d,0),(0,d)\}`$ or $`\mathrm{\Delta }=[0,r]\times [0,s]`$ then the number $`N^{\mathrm{\Delta },\delta }`$ is the multicomponent Gromov-Witten invariant of genus $`(l\delta )`$ and degree $`d`$ in $`^2`$ or of bidegree $`(r,s)`$ in $`^1\times ^1`$.

###### Special Case .

Suppose $`\mathrm{\Delta }=\mathrm{\Delta }_d`$ so that $`T_\mathrm{\Delta }=^2`$. We have $`m=\frac{d(d+3)}{2}`$ and $`l=\frac{(d1)(d2)}{2}`$. The number $`N^{\mathrm{\Delta }_d,\delta }=N_{g,d}`$ is the number of genus $`g`$, degree $`d`$ (not necessarily irreducible) curves passing through $`3d1+g`$ generic points in $`^2`$, $`g=\frac{(d1)(d2)}{2}\delta `$.

The formula $`N^{\mathrm{\Delta }_d,1}=3(d1)^2`$ is well-known as the degree of the discriminant (cf. ). An elegant recursive formula for the number of irreducible rational curves (the one-component part of $`N_{0,d}`$) was found by Kontsevich . An algorithm for computing $`N_{g,d}`$ for arbitrary $`g`$ is due to Caporaso and Harris . See  for computations for some other rational surfaces, in particular, the Hirzebruch surfaces (this corresponds to the case when $`\mathrm{\Delta }`$ is a trapezoid).

## 2. Lattice paths and their multiplicities

A path $`\gamma :[0,n]^2`$, $`n`$, is called a lattice path if $`\gamma |_{[j1,j]}`$, $`j=1,\mathrm{},n`$ is an affine-linear map and $`\gamma (j)^2`$, $`j0,\mathrm{},n`$. Clearly, a lattice path is determined by its values at the integer points. Let us choose an auxiliary linear map $`\lambda :^2`$ that is irrational, i.e. such that $`\lambda |_^2`$ is injective. Let $`p,q\mathrm{\Delta }`$ be the vertices where $`\alpha |_\mathrm{\Delta }`$ reaches its minimum and maximum respectively. A lattice path is called $`\lambda `$-increasing if $`\lambda \gamma `$ is increasing.

The points $`p`$ and $`q`$ divide the boundary $`\mathrm{\Delta }`$ into two increasing lattice paths

$$\alpha ^+:[0,n_+]\mathrm{\Delta }\text{and}\alpha ^{}:[0,n_{}]\mathrm{\Delta }.$$

We have $`\alpha _+(0)=\alpha _{}(0)=p`$, $`\alpha _+(n_+)=\alpha _{}(n_{})=q`$, $`n_++n_{}=ml+3`$. To fix a convention we assume that $`\alpha _+`$ goes clockwise around $`\mathrm{\Delta }`$ wile $`\alpha _{}`$ goes counterclockwise.

Let $`\gamma :[0,n]\mathrm{\Delta }^2`$ be an increasing lattice path such that $`\gamma (0)=p`$ and $`\gamma (n)=q`$. The path $`\gamma `$ divides $`\mathrm{\Delta }`$ into two closed regions: $`\mathrm{\Delta }_+`$ enclosed by $`\gamma `$ and $`\alpha _+`$ and $`\mathrm{\Delta }_{}`$ enclosed by $`\gamma `$ and $`\alpha _{}`$. Note that the interiors of $`\mathrm{\Delta }_+`$ and $`\mathrm{\Delta }_{}`$ do not have to be connected.

We define the positive (resp. negative) multiplicity $`\mu _\pm (\gamma )`$ of the path $`\gamma `$ inductively. We set $`\mu _\pm (\alpha _\pm )=1`$. If $`\gamma \alpha _\pm `$ then we take $`1kn1`$ to be the smallest number such that $`\gamma (k)`$ is a vertex of $`\mathrm{\Delta }_\pm `$ with the angle less than $`\pi `$ (so that $`\mathrm{\Delta }_\pm `$ is locally convex at $`\gamma (k)`$).

If such $`k`$ does not exist we set $`\mu _\pm (\gamma )=0`$. If $`k`$ exist we consider two other increasing lattice paths connecting $`p`$ and $`q`$ $`\gamma ^{}:[0,n1]\mathrm{\Delta }`$ and $`\gamma ^{\prime \prime }:[0,n]^2`$. We define $`\gamma ^{}`$ by $`\gamma ^{}(j)=\gamma (j)`$ if $`j<k`$ and $`\gamma ^{}(j)=\gamma (j+1)`$ if $`jk`$. We define $`\gamma ^{\prime \prime }`$ by $`\gamma ^{\prime \prime }(j)=\gamma (j)`$ if $`jk`$ and $`\gamma ^{\prime \prime }(k)=\gamma (k1)+\gamma (k+1)\gamma (k)^2`$. We set

$$\mu _\pm (\gamma )=2\mathrm{Area}(T)\mu _\pm (\gamma ^{})+\mu _\pm (\gamma ^{\prime \prime }),$$

where $`T`$ is the triangle with the vertices $`\gamma (k1)`$, $`\gamma (k)`$ and $`\gamma (k+1)`$. The multiplicity is always integer since the area of a lattice triangle is half-integer.

Note that it may happen that $`\gamma ^{\prime \prime }(k)\mathrm{\Delta }`$. In such case we use a convention $`\mu _\pm (\gamma ^{\prime \prime })=0`$. We may assume that $`\mu _\pm (\gamma ^{})`$ and $`\mu _\pm (\gamma ^{\prime \prime })`$ is already defined since the area of $`\mathrm{\Delta }_\pm `$ is smaller for the new paths. Note that $`\mu _\pm =0`$ if $`n<n_\pm `$ as the paths $`\gamma ^{}`$ and $`\gamma ^{\prime \prime }`$ are not longer than $`\gamma `$.

We define the multiplicity of the path $`\gamma `$ as the product $`\mu _+(\gamma )\mu _{}(\gamma )`$. Note that the multiplicity of a path connecting two vertices of $`\mathrm{\Delta }`$ does not depend on $`\lambda `$. We only need $`\lambda `$ to determine whether a path is increasing.

###### Example 1.

Consider the path $`\gamma :[0,8]\mathrm{\Delta }_3`$ depicted on the extreme left of Figure 1. This path is increasing with respect to $`\lambda (x,y)=xϵy`$, where $`ϵ>0`$ is very small.

Let us compute $`\mu _+(\gamma )`$. We have $`k=2`$ as $`\gamma (2)=(0,1)`$ is a locally convex vertex of $`\mathrm{\Delta }_+`$. We have $`\gamma ^{\prime \prime }(2)=(1,3)\mathrm{\Delta }_3`$ and thus $`\mu _+(\gamma )=\mu _+(\gamma ^{})`$, since $`\mathrm{Area}(T)=\frac{1}{2}`$. Proceeding further we get $`\mu _+(\gamma )=\mu _+(\gamma ^{})=\mathrm{}=\mu _+(\alpha _+)=1`$.

Let us compute $`\mu _{}(\gamma )`$. We have $`k=3`$ as $`\gamma (3)=(1,2)`$ is a locally convex vertex of $`\mathrm{\Delta }_{}`$. We have $`\gamma ^{\prime \prime }(3)=(0,0)`$ and $`\mu _{}(\gamma ^{\prime \prime })=1`$. To compute $`\mu _{}(\gamma ^{})=1`$ we note that $`\mu _{}((\gamma ^{})^{})=0`$ and $`\mu _{}((\gamma ^{})^{\prime \prime })=1`$. Thus the full multiplicity of $`\gamma `$ is 2.

## 3. The formula

In the previous section we fixed an auxiliary linear function $`\lambda :^2`$ which determines the extremal vertices $`p,q`$ of $`\mathrm{\Delta }`$.

###### Theorem 1.

The number $`N^{\mathrm{\Delta },\delta }`$ equals to the number (counted with multiplicities) of $`\lambda `$-increasing lattice paths $`[0,m\delta ]\mathrm{\Delta }`$ connecting $`p`$ and $`q`$.

This theorem is proved in  (to appear). The proof is based on the application of the so-called tropical algebraic geometry (see e.g. Chapter 9 of ). The relation between the classical enumerative problem and the corresponding tropical problem is provided by passing to the “large complex limit” as suggested by Kontsevich (see  for these ideas in a more general setting).

Note that an immediate corollary of Theorem 1 is that the number of $`\lambda `$-increasing lattice paths of a fixed length does not depend on the choice of $`\lambda `$.

###### Example 2.

Let us compute $`N^{\mathrm{\Delta },1}=5`$ for the polygon $`\mathrm{\Delta }`$ depicted on Figure 2 in two different ways. Using $`\lambda (x,y)=x+ϵy`$ for a small $`ϵ>0`$ we get the left two paths depicted on Figure 2. Using $`\lambda (x,y)=x+ϵy`$ we get the three right paths. The corresponding multiplicities are shown under the path. All other $`\lambda `$-increasing paths have zero multiplicity.

In the next two examples we use $`\lambda (x,y)=xϵy`$ as the auxiliary linear function.

###### Example 3.

Figure 3 shows a computation of the well-known number $`N^{\mathrm{\Delta }_3,1}=N_{0,3}`$. This is the number of rational cubic curves through 8 generic points in $`^2`$.

###### Example 4.

Figure 4 shows a computation of a less well-known number $`N^{\mathrm{\Delta }_4,2}=N_{1,4}`$. This is the number of genus 1 quartic curves through 12 generic points in $`^2`$.

## 4. Real aspects of the count

Suppose that $`z_1,\mathrm{},z_{m\delta }(^{})^2(^{})^2`$ are generic real points. We may ask how many of the $`N^{\mathrm{\Delta },\delta }`$ relevant complex curves are real, i.e. defined over $``$. Note that this number depends on the configuration of real points.

Theorem 1 can be modified to give the relevant count of real curves. In order to do this we need to define the real multiplicity of a lattice path $`\gamma :[0,n]\mathrm{\Delta }`$ connecting the vertices $`p`$ and $`q`$. We introduce the sequence of the pairs of signs $`\sigma _1,\mathrm{},\sigma _n_2_2`$ (this sequence will record the quadrants of the points $`z_j(^{})^2`$). The sign $`\sigma _j`$ is prescribed to the edge $`\gamma [j1,j]`$. We make a convention that $`\sigma _j`$ and $`\sigma _j^{}`$ are equivalent if $`\sigma _j\sigma _j^{}\gamma (j)\gamma (j1)(mod2).`$

We set

(1)
$$\mu _\pm ^{}(\gamma )=a(T)\mu _\pm ^{}(\gamma ^{})+\mu _\pm ^{}(\gamma ^{\prime \prime }).$$

The definition of the new paths $`\gamma ^{}`$, $`\gamma ^{\prime \prime }`$ and the triangle $`T`$ is the same as in section 2. The sign sequence for $`\gamma ^{\prime \prime }`$ is $`\sigma _j^{\prime \prime }=\sigma _j,jk,k+1`$, $`\sigma _k^{\prime \prime }=\sigma _{k+1}`$, $`\sigma _{k+1}^{\prime \prime }=\sigma _k`$. The sign sequence for $`\gamma ^{}`$ is $`\sigma _j^{}=\sigma _j,j<k`$, $`\sigma _j^{}=\sigma _{j+1},j>k`$. We define the sign $`\sigma _k^{}`$ and the function $`a(T)`$ as follows.

* If all sides of $`T`$ are odd we set $`a(T)=1`$ and define the sign $`\sigma _k^{}`$ (up to the equivalence) by the condition that the three equivalence classes of $`\sigma _k`$, $`\sigma _{k+1}`$ and $`\sigma _k^{}`$ do not share a common element.
* If all sides of $`T`$ are even we set $`a(T)=0`$ if $`\sigma _{k1}\sigma _k`$. In this case we can ignore $`\gamma ^{}`$ (and its sequence of signs). We set $`a(T)=4`$ if $`\sigma _k=\sigma _{k+1}`$. In this case we define $`\sigma _k^{}=\sigma _k=\sigma _{k+1}`$.
* Otherwise we set $`a(T)=0`$ if the equivalence classes of $`\sigma _k`$ and $`\sigma _{k+1}`$ do not have a common element. We set $`a(T)=2`$ if they do. In the latter case we define the equivalence class of $`\sigma _k^{}`$ by the condition that $`\sigma _k`$, $`\sigma _{k+1}`$ and $`\sigma _k^{}`$ have a common element. There is one exception to this rule. If the even side is $`\gamma (k+1)\gamma (k1)`$ then there are two choices for $`\sigma _k^{}`$ satisfying the above condition. In this case we replace $`a(T)\mu _\pm ^{}(\gamma ^{})`$ in (1) by the sum of the two multiplicities of $`\gamma ^{}`$ equipped with the two allowable choices for $`\sigma _k^{}`$ (note that this agrees with $`a(T)=2`$ in this case).

Similar to section 2 we define $`\mu _\pm ^{}(\alpha _\pm )=1`$ and $`\mu ^{}(\gamma )=\mu _+^{}(\gamma )\mu _{}^{}(\gamma )`$. As before $`\lambda :^2`$ is a linear map injective on $`^2`$ and $`p`$ and $`q`$ are the extrema of $`\lambda |_\mathrm{\Delta }`$.

###### Theorem 2.

For any choice of $`\lambda `$ and $`\sigma _j,j=1,\mathrm{},m\delta `$ there exists a configuration of $`m\delta `$ of generic points in the respective quadrants such that the number of real curves among the $`N^{\mathrm{\Delta },\delta }`$ relevant complex curves is equal to the number of $`\lambda `$-increasing lattice paths $`\gamma :[0,m\delta ]\mathrm{\Delta }`$ connecting $`p`$ and $`q`$ counted with multiplicities $`\mu ^{}`$.

###### Example 5.

Here we use the choice $`\sigma _j=(+,+)`$ so all the points $`z_j`$ are in the positive quadrant $`(_{>0})^2(^{})^2`$. The first count of $`N^{\mathrm{\Delta },1}`$ from Example 2 gives a configuration of 3 real points with 5 real curves. The second count gives a configuration with 3 real curves as the real multiplicity of the last path is 1. Note also that the second path on Figure 2 changes its real multiplicity if we reverse its direction.

Example 3 gives a configuration of 9 generic points in $`^2`$ with all 12 nodal cubics through them real. Example 4 gives a configuration of 12 generic points in $`^2`$ with 217 out of the 225 quartics of genus 1 real. The path in the middle of Figure 4 has multiplicity 9 but real multiplicity 1. A similar computation shows that there exists a configuration of 11 generic points in $`^2`$ such that 564 out of the 620 irreducible quartic through them are real.

###### Remark 3.

Real nodal curves have three types of nodes: hyperbolic, elliptic and imaginary. Theorem 2 can be refined to count curves with different types of nodes separately. In accordance with  let us prescribe a sign $`(1)^e`$ to a real nodal curve, where $`e`$ is the number of its elliptic nodes. To compute the corresponding algebraic number of curves we introduce the multiplicity $`\nu ^{}`$ by replacing (1) with $`\nu _\pm ^{}(\gamma )=b(T)\nu _\pm ^{}(\gamma ^{})+\nu _\pm ^{}(\gamma ^{\prime \prime }).`$ Here we define $`b(T)=0`$ if at least one side of $`T`$ is even and $`b(T)=(1)^{\mathrm{\#}(\mathrm{Int}T^2)}`$ otherwise. It can be shown with the help of this formula and a combinatorial observation made by Itenberg, Kharlamov and Shustin (to appear) that in the case $`\mathrm{\Delta }=\mathrm{\Delta }_d`$ the algebraic number of irreducible curves counted by $`\nu ^{}`$ is positive for any genus $`0g\frac{(d1)(d2)}{2}`$ if $`\lambda (x,y)=yϵx`$.

Note that unlike $`\mu ^{}`$ the multiplicity $`\nu ^{}`$ does not depend on the quadrant choices $`\sigma _j`$. Furthermore, in  Welschinger stated that if $`g=0`$ then this number is independent of the configuration of generic real points. Corollary 1.2 of  combined with Remark 3 implies the following statement (which answers a question asked e.g. by Rokhlin and Kharlamov). For any configuration of generic $`3d1`$ points in $`^2`$ there exists a real rational curve of degree $`d`$ passing through this configuration.