text:
# Renormalization group and nonlinear susceptibilities of cubic ferromagnets at criticality

## I Introduction

Three decades ago Wilson and Fisher discovered that in the vicinity of the critical point the order parameter fluctuations may change the effective anisotropy of the system drastically . Having studied the anisotropic XY model within the $`ϵ`$ expansion they found that approaching $`T_c`$ this model either becomes effectively isotropic or developes its anisotropy further  until the fluctuation-induced first-order phase transition takes place . What scenario the system chooses depends on the value and sign of the anisotropic coupling constant in the initial Landau-Wilson Hamiltonian. In the course of studying of the generalized cubic model, the crucial role of the order parameter dimensionality $`n`$ was revealed: for $`n<n_c`$ the model undergoing the continuous phase transition demonstrates isotropic critical behavior while for $`n>n_c`$ it remains anisotropic at criticality . The numerical value of the marginal spin dimensionality $`n_c`$ separating these two regimes is of prime physical importance since it fixes the true mode of the critical behavior of real cubic ferromagnets and of some other systems of interest.

First estimates of $`n_c`$ for the cubic model were obtained within the lower-order renormalization-group (RG) calculations  and turned out to be in favor of the inequality $`n_c>3`$ indicating that cubic ferromagnets should belong to the class of universality of the three-dimensional (3D) Heisenberg model. The classical technique of the high-temperature expansions seemed to support this conclusion . However, the resummation of the three-loop RG expansions  and subsequent higher-order RG calculations performed both in three  and $`(4ϵ)`$ (Ref.) dimensions shifted the numerical estimate for $`n_c`$ downwards, fixing it below 3 . Remarkable consensus between different field-theoretical approaches and resummation techniques was achieved in the course of this study: 3D RG calculations , resummed $`ϵ`$ expansion , biased $`ϵ`$ expansion technique , and pseudo-$`ϵ`$ expansion analysis  lead to $`n_c=2.89`$, $`n_c=2.855`$, $`n_c=2.87`$, and $`n_c=2.862`$, respectively. Moreover, the predicted numerical value of $`n_c`$ was found to be stable with respect to the order of perturbation theory: starting from the four-loop approximation, the 3D RG estimate stays practically unchanged varying around 2.9 within the interval $`\pm 0.02`$ . Hence, according to the up to date, most accurate theoretical data, cubic ferromagnets should demonstrate the specific – cubic – critical behavior with a special set of critical exponents.

Since $`n_c`$ is very close to the physical value $`n=3`$, the fixed point governing the cubic critical behavior lies very near the Heisenberg fixed point at the RG flow diagram and the critical exponents for both fixed points almost coincide. For instance, the susceptibility exponent $`\gamma `$ is equal to 1.3895(50) for the 3D Heisenberg model  while at the cubic fixed point $`\gamma =1.390(6)`$ . It is clear that measuring the critical exponents one cannot distinguish between the cubic and Heisenberg critical behaviors.

In such a situation, some alternative physical quantities should be addressed to clear up how the systems with the cubic symmetry behave approaching $`T_c`$. In this paper, the nonlinear susceptibilities of the cubic model are studied in the framework of the field-theoretical RG approach in three dimensions, with particular attention paid to their anisotropy at criticality. It will be shown that the anisotropy of nonlinear susceptibilities is sensitive to the type of critical asymptotics and its measurement can be used, in principle, for detection of the cubic (non-Heisenberg) critical behavior, provided the anisotropy is as strong as the modern higher-order RG calculations predict.

## II Nonlinear susceptibilities and effective coupling constants

In the critical region, the expansion of the free energy of the cubic model in powers of the magnetization components $`M_\alpha `$ may be written in the form:

$`F(M_\alpha ,m)`$ $`=`$ $`F(0,m)+{\displaystyle \frac{1}{2}}m^{2\eta }M_\alpha ^2+m^{12\eta }(u_4+v_4\delta _{\alpha \beta })M_\alpha ^2M_\beta ^2`$ (2)
$`+m^{3\eta }(u_6+q_6\delta _{\alpha \beta }+v_6\delta _{\alpha \beta }\delta _{\alpha \gamma })M_\alpha ^2M_\beta ^2M_\gamma ^2+\mathrm{}`$
where $`\eta `$ is a Fisher exponent, $`m`$ being an inverse correlation length, and $`u_4`$, $`v_4`$, $`u_6`$, $`q_6`$, and $`v_6`$ are dimensionless effective coupling constants acquiring, under $`TT_c`$, certain universal values. These coupling constants are related to the nonlinear susceptibilities defined in a conventional way:

$$\chi _{\alpha \beta \gamma \delta }^{(4)}=\frac{^3M_\alpha }{H_\beta H_\gamma H_\delta }|_{H=0},\chi _{\alpha \beta \gamma \delta \mu \nu }^{(6)}=\frac{^5M_\alpha }{H_\beta H_\gamma H_\delta H_\mu H_\nu }|_{H=0}.$$
(3)
Of prime importance are the nonlinear susceptibilities in two particular cases, when (i) an external magnetic field is parallel to a cubic axis ($`\chi _c^{(i)}`$) and (ii) it is oriented along a cell space diagonal ($`\chi _d^{(i)}`$). For these two highly symmetric directions the differences between the corresponding values of the nonlinear susceptibilities are known to be maximal: i. e., the anisotropy is most pronounced. It is easy to show that

$$\chi _c^{(4)}=24\frac{\chi ^2}{m^3}(u_4+v_4),\chi _d^{(4)}=24\frac{\chi ^2}{m^3}\left(u_4+\frac{v_4}{3}\right).$$
(4)
$`\chi _c^{(6)}=720{\displaystyle \frac{\chi ^3}{m^6}}\left[8(u_4+v_4)^2u_6q_6v_6\right],`$ (5)
$`\chi _d^{(6)}=720{\displaystyle \frac{\chi ^3}{m^6}}\left[8\left(u_4+{\displaystyle \frac{v_4}{3}}\right)^2u_6{\displaystyle \frac{q_6}{3}}{\displaystyle \frac{v_6}{9}}\right].`$ (6)
where $`\chi `$ is a linear susceptibility. To characterize the anisotropy strength, we define the reduced parameters

$$\delta ^{(4)}=\frac{|\chi _c^{(4)}\chi _d^{(4)}|}{\chi _c^{(4)}},\delta ^{(6)}=\frac{|\chi _c^{(6)}\chi _d^{(6)}|}{\chi _c^{(6)}}.$$
(7)
Below, they will be estimated at criticality.

Let us start from the lower-order nonlinear susceptibility $`\chi ^{(4)}`$. The coordinates of the cubic fixed point $`u_4^{}`$ and $`v_4^{}`$ in three dimensions are known from the higher-order RG calculations, with resummed four-, five-, and six-loop RG expansions yielding very close numerical results . Considering recent six-loop RG estimates as the most reliable ones, we accept $`u_4^{}=0.755\pm 0.010`$, $`v_4^{}=0.067\pm 0.014`$ . The substitution of these numbers into Eqs. (3) and (5) gives

$$\delta ^{(4)}=0.054\pm 0.012.$$
(8)
In the next section, the universal critical values of the sextic coupling constants will be evaluated and an estimate for $`\delta ^{(6)}`$ will be obtained.

## III RG expansions for sextic couplings and anisotropy of sixth-order susceptibility

Let us proceed to determination of $`\chi _c^{(6)}`$, $`\chi _d^{(6)}`$, and $`\delta ^{(6)}`$. As seen from Eq. (4), we are in a position to calculate the effective coupling constants $`u_6`$, $`q_6`$, and $`v_6`$. They will be found perturbatively, using the field-theoretical RG approach in three dimensions. Our analysis is based on the well-known Landau-Wilson Hamiltonian of the 3D $`n`$-vector cubic model:

$$H=\frac{1}{2}d^3x\left[m_0^2\phi _\alpha ^2+(\phi _\alpha )^2+\frac{u_0}{12}\phi _\alpha ^2\phi _\beta ^2+\frac{v_0}{12}\phi _\alpha ^4\right],$$
(9)
where $`m_o^2`$ is the reduced deviation from the mean-field transition temperature. All the RG calculations are carried out within a massive theory under zero-momentum normalizing conditions for the renormalized Green function $`G_R(p,m)`$ and the four-point vertices $`U_R(p_i,m,u,v,)`$, $`V_R(p_i,m,u,v,)`$:

$`G_R^1(0,m)=m^2,{\displaystyle \frac{G_R^1(p,m)}{p^2}}|_{p^2=0}=1,`$ (10)
$`U_R(0,m,u,v)=mu,V_R(0,m,u,v)=mv.`$ (11)
Here, the value of the one-loop vertex graph including the factor $`(n+8)`$ is absorbed in $`u`$ and $`v`$ in order to make the coefficient for the $`u^2`$ term in the $`\beta _u`$ function equal to unity. The quartic effective couplings $`u`$ and $`v`$ thus defined are related to $`u_4`$ and $`v_4`$ entering Eqs. (1), (3), and (4) in an obvious way:

$$u=\frac{n+8}{2\pi }u_4,v=\frac{n+8}{2\pi }v_4.$$
(12)
We limit ourselves to the four-loop RG approximation that proved to lead to quite a good numerical estimate for the universal value of the sextic effective coupling for the 3D O(n)-symmetric model . Since the symmetry factors and integrals for all the relevant Feynman diagrams have been found earlier , what we have to do here is the calculation of the tensor (field) factors for the six-point vertex graphs and mass insertions generated by the O(n)-symmetric and cubic interactions. Performing this calculation and then renormalizing the ”bare” perturbative expansions for the six-point vertices precisely in the same way as has been done for the isotropic n-vector model , we arrive at the following four-loop RG expansions:

$`u_6`$ $`=`$ $`{\displaystyle \frac{8\pi ^2}{3}}{\displaystyle \frac{u^2}{(n+8)^3}}[(n+26)u+9v{\displaystyle \frac{1}{(n+8)}}({\displaystyle \frac{34n+452}{3}}u^2+124uv+18v^2)`$ (18)
$`+{\displaystyle \frac{1}{(n+8)^2}}[(1.065025n^2+157.42454n+1323.0960)u^3`$
$`+(19.382741n+1873.9825)u^2v+(2.3101121n+679.33934)uv^2+100.90652v^3]`$
$`{\displaystyle \frac{1}{(n+8)^3}}[(0.0638n^3+52.451n^2+2314.99n+14387.6)u^4`$
$`+(2.3434n^2+1006.56n+29601.0)u^3v+(100.156n+18985.4)u^2v^2`$
$`+(11.286n+5872.61)uv^3+782.865v^4]].`$
$`q_6`$ $`=`$ $`{\displaystyle \frac{8\pi ^2}{3}}{\displaystyle \frac{uv}{(n+8)^3}}[72u+27v{\displaystyle \frac{1}{(n+8)}}((4n+520)u^2+468uv+108v^2)`$ (24)
$`+{\displaystyle \frac{1}{(n+8)^2}}[(1.1357879n^2+149.54004n+5366.1581)u^3`$
$`+(32.454138n+8267.9231)u^2v+4107.2194uv^2+775.79040v^3]`$
$`{\displaystyle \frac{1}{(n+8)^3}}[(0.3986n^319.2521n^2+3351.86n+66591.7)u^4`$
$`+(11.6195n^2+115.223n+148119)u^3v+(419.561n+119322)u^2v^2`$
$`+46560.3uv^3+7520.89v^4]].`$
$`v_6`$ $`=`$ $`{\displaystyle \frac{8\pi ^2}{3}}{\displaystyle \frac{v^2}{(n+8)^3}}[54u+27v{\displaystyle \frac{1}{(n+8)}}[(6n+492)u^2+540uv+162v^2]`$ (30)
$`+{\displaystyle \frac{1}{(n+8)^2}}[(1.3828190n^250.171743n+5947.5257)u^3`$
$`+(88.612482n+10696.342)u^2v+6632.1372uv^2+1481.5855v^3]`$
$`{\displaystyle \frac{1}{(n+8)^3}}[(0.4873n^3+6.2132n^2+114.760n+83872.7)u^4`$
$`+(25.5088n^22440.17n+212729)u^3v+(1281.62n+205263)u^2v^2`$
$`+93009.1uv^3+16755.0v^4]].`$
These RG expansions should obey some exact relations appropriate to systems with cubic anisotropy. Indeed, for $`n=2`$ the model (7) is known to possess a specific symmetry property. If the field $`\phi _\alpha `$ undergoes the transformation

$$\phi _1\frac{\phi _1+\phi _2}{\sqrt{2}},\phi _2\frac{\phi _1\phi _2}{\sqrt{2}},$$
(31)
the quartic coupling constants are also transformed,

$$uu+\frac{3}{2}v,vv,$$
(32)
but the structure of the Hamiltonian itself remains unchanged . Since the RG functions of the problem are completely determined by the structure of the Hamiltonian, the RG equations should be invariant with respect to any transformation conserving this structure . This should be also true for all the expressions relating various universal quantities to each other. If, under $`n=2`$, we apply the transformation (13) to the magnetization components in Eq. (1), the free-energy expansion remains the same, provided $`u_4`$, $`v_4`$ are replaced according to Eq. (14) while the sextic coupling constants are transformed as prescribed below:

$$u_6+\frac{q_6}{2}+\frac{v_6}{4}u_6+q_6+v_6,\frac{2}{3}q_6+v_6\frac{2}{3}q_6v_6.$$
(33)
It means that the following relations between the effective coupling constants should hold:

$`u_6(u,v)+{\displaystyle \frac{q_6(u,v)}{2}}+{\displaystyle \frac{v_6(u,v)}{4}}`$ $`=`$ $`u_6(u+{\displaystyle \frac{3}{2}}v,v)+q_6(u+{\displaystyle \frac{3}{2}}v,v)+v_6(u+{\displaystyle \frac{3}{2}}v,v),`$ (34)
$`{\displaystyle \frac{2}{3}}q_6(u,v)+v_6(u,v)`$ $`=`$ $`{\displaystyle \frac{2}{3}}q_6(u+{\displaystyle \frac{3}{2}}v,v)v_6(u+{\displaystyle \frac{3}{2}}v,v).`$ (35)
As one can check, the expansions (10), (11), and (12) do really satisfy these relations. Moreover, consideration of all the remaining limiting cases (see, e. g. Refs.) shows that these RG series are in accord with their counterparts obtained earlier for the Ising  and O(n)-symmetric  models.

Let us use the expansions just found for estimation of the universal critical values of the sextic effective coupling constants for the cubic ferromagnets, i. e., in the physical case $`n=3`$. To obtain the numbers of interest from the divergent (asymptotic) RG series, a proper resummation procedure should be applied. Here we use the Padé-Borel-Leroy resummation technique which demonstrates high numerical effectiveness both for simple \[Ising and O(n)-symmetric\] models  and for complicated anisotropic systems preserving their internal symmetries (see, e. g. Ref. for detail). Since the expansions of quantities depending on two variables $`u`$ and $`v`$ are dealt with, the Borel-Leroy transformation is taken in a generalized form:

$$f(u,v)=\underset{ij}{}c_{ij}u^iv^j=\underset{0}{\overset{\mathrm{}}{}}e^tt^bF(ut,vt)𝑑t,F(x,y)=\underset{ij}{}\frac{c_{ij}x^iy^j}{(i+j+b)!}.$$
(36)
To perform an analytical continuation, the resolvent series

$$\stackrel{~}{F}(x,y,\lambda )=\underset{n=0}{\overset{\mathrm{}}{}}\lambda ^n\underset{l=0}{\overset{n}{}}\underset{m=0}{\overset{nl}{}}\frac{c_{l,nl}x^ly^{nl}}{n!}$$
(37)
is constructed which is a series in powers of $`\lambda `$ with the coefficients being uniform polynomials in $`u`$, $`v`$ and then Padé approximants $`[L/M]`$ in $`\lambda `$ at $`\lambda =1`$ are used.

With the four-loop RG expansions in hand, we can construct, in principle, three different Padé approximants: \[2/1\], \[1/2\], and \[0/3\]. To obtain proper approximation schemes, however, only diagonal \[L/L\] and near-diagonal Padé approximants should be employed . That is why further, when estimating $`u_6^{}`$, $`q_6^{}`$, and $`v_6^{}`$, we limit ourselves to the approximants \[2/1\] and \[1/2\]. Moreover, the diagonal Padé approximant \[1/1\] will be also dealt with although this corresponds, in fact, to the use of the lower-order, three-loop RG approximation.

The algorithm for estimating of the universal critical values of sextic effective couplings is as follows. Since, in fact, we have the power series for the ratios $`R_u=u_6/u^2`$, $`R_q=q_6/(uv)`$, and $`R_v=v_6/v^2`$ rather than for $`u_6`$, $`q_6`$, and $`v_6`$ themselves, we work with the RG series for these ratios. They are resummed in three different ways based on the Borel-Leroy transformation and the Padé approximants just mentioned. The Borel-Leroy integral is evaluated as a function of the shift parameter $`b`$ under $`u=u^{}`$, $`v=v^{}`$. For the cubic fixed-point coordinates the values given by the six-loop RG analysis at $`n=3`$ are adopted : $`u^{}=1.321`$, $`v^{}=0.117`$. The optimal value of $`b`$ providing the fastest convergence of the iteration scheme is then determined. It is deduced from the condition that the Padé approximants employed should give, for $`b=b_{opt}`$, the values of $`R_u`$ ($`R_q`$, $`R_v`$) which are as close as possible to each other. Finally, the average over three estimates for $`R_u`$ ($`R_q`$, $`R_v`$) is found and claimed to be a numerical value of this universal ratio.

To demonstrate how such a procedure works, in Table 1 the results of the corresponding calculations are presented. The empty cells in this table reflect the fact that for some values of the shift parameter $`b`$ Padé approximant \[1/2\] turns out to be spoiled by the ”dangerous” poles, i. e., by the poles at positive or small negative $`t`$. As one can see, for $`u_6^{}`$ (an integer) $`b_{opt}`$, providing a maximal closeness of the estimates given by three working Padé approximants, is equal to 2, while for $`q_6^{}`$ and $`v_6^{}`$, $`b_{opt}=3`$. So the results of our four-loop RG analysis are as follows:

$$u_6^{}=0.842,q_6^{}=0.175,v_6^{}=0.0108.$$
(38)
How close to the exact universal values of sextic couplings may these numbers be? Obviously, the accuracy of our estimates is limited both by the accuracy achieved in the course of evaluating the cubic fixed-point coordinates $`u^{}`$, $`v^{}`$ themselves and by the speed of convergence of the iteration procedure employed. The influence of the latter factor may be characterized by the sensitivity of the numerical results with respect to variation of $`b`$. As is seen from Table 1, the values of $`u_6^{}`$, $`q_6^{}`$, and $`v_6^{}`$, averaged over the working Padé approximants, vary by no more than $`\pm 0.005`$, $`\pm 0.002`$, and $`\pm 0.0001`$, respectively, when $`b`$ runs from 0 to 20. On the other hand, these error bars are too small to be referred to as realistic ones. Their conservative counterparts may be deduced, accepting that the true universal values of renormalized sextic couplings should lie within the intervals of the variation of the estimates given by the leading Padé approximant, namely, by the higher-order pole-free approximant \[2/1\]. These intervals are easily extracted from Table 1 and are as follows: $`\pm 0.03(u_6^{})`$, $`\pm 0.01(q_6^{})`$, $`\pm 0.001(v_6^{})`$. Contrary to the aforementioned tiny ranges, the error bars thus determined look quite realistic. Moreover, they are found to be large enough to cover the inaccuracy produced by the uncertainty of the cubic fixed-point location. Hence, we adopt these error bars as the final ones.

With the universal critical values of the quartic and sextic couplings in hand, we are able to estimate the reduced anisotropy of the sixth-order nonlinear susceptibility. Combining Eqs. (4), (5), and (19), we obtain

$$\delta ^{(6)}=0.102\pm 0.02.$$
(39)
## IV Eighth-order nonlinear susceptibility and discussion

Let us discuss the role played by the sextic coupling constants in forming the magnitudes of $`\chi _c^{(6)}`$, $`\chi _d^{(6)}`$, and $`\delta ^{(6)}`$. It may be shown that, in fact, their contrubutions to the sixth-order susceptibility and the anisotropy parameter are very small. Indeed, neglecting $`u_6`$, $`q_6`$, and $`v_6`$ in Eq. (4) changes the critical value of $`\delta ^{(6)}`$ by about 4%. In other words, the anisotropy of $`\chi ^{(6)}`$ near $`T_c`$ is fixed practically by the quartic coupling constants only. A similar situation should take place for the higher-order (eighth-order, tenth-order, etc.) susceptibilities. To demonstrate this, we consider the eighth-order susceptibility $`\chi ^{(8)}`$ for the isotropic (Heisenberg) model. As $`\chi ^{(4)}`$ and $`\chi ^{(6)}`$, it may expressed via the effective coupling constants

$$\chi ^{(8)}=40320\frac{\chi ^4}{m^9}\left(96u_4^324u_4u_6+u_8\right).$$
(40)
Up to now, the universal critical value of the octic coupling constant was determined using the field-theoretical RG approach in three dimensions , the biased $`ϵ`$ expansion technique , and the exact renormalization group machinery . The methods employed lead to the estimates $`u_8^{}=0.168`$, $`u_8^{}=0.36`$, and $`u_8^{}=0.145`$, respectively. Although these numbers are considerably scattered, all of them are several times smaller than $`u_4^{}=0.794`$ (Ref.) and $`u_6^{}=0.951`$ (Ref.). It means that, because of the big numerical coefficients of the first two terms in brackets of Eq. (21), the contribution of $`u_8`$ to $`\chi ^{(8)}`$ is negligible at criticality. Since the cubic fixed point is located very near the Heisenberg one at the RG flow diagram, the same conclusion is valid for the 3D cubic model. This enables us to estimate the critical anisotropy of $`\chi ^{(8)}`$ without the calculation of the universal values of the octic coupling constants. Making in Eq. (21) the substitutions $`u_4u_4+v_4`$, $`u_6u_6+q_6+v_6`$ and $`u_4u_4+v_4/3`$, $`u_6u_6+q_6/3+v_6/9`$ with $`u_8`$ omitted, we obtain expressions for $`\chi _c^{(8)}`$ and $`\chi _d^{(8)}`$, respectively. Using then the known coordinates of the cubic fixed point and the critical values of the sextic coupling constants (19), we obtain the reduced anisotropy of the eighth-order susceptibility:

$$\delta ^{(8)}=\frac{|\chi _c^{(8)}\chi _d^{(8)}|}{\chi _c^{(8)}}=0.144\pm 0.04.$$
(41)
The estimates (6), (20), (22) for $`\delta ^{(4)}`$, $`\delta ^{(6)}`$, and $`\delta ^{(8)}`$ at criticality do not look too small to prevent the detection of the universal cubic anisotropy in physical and computer experiments or by means of a thorough analysis of the high-temperature expansions . It is well known, however, that the approach of the universal critical regime is controlled by the tiny exponent $`\omega 0.01`$ , making both an experimental study and simulation of the true mode of critical behavior rather difficult . We believe that the results obtained in this paper will help those searching for the cubic asymptotic regimes in real systems and lattice models to choose a proper strategy and the quantities to be measured.

## V Conclusion

To summarize, we have evaluated the nonlinear susceptibilities of the forth, sixth, and eighth orders for the 3D cubic ferromagnet and estimated the parameters $`\delta ^{(4)}`$, $`\delta ^{(6)}`$, and $`\delta ^{(8)}`$ characterizing their anisotropy at the critical point. In the course of this study, we have calculated the four-loop RG expansions for the effective coupling constants $`u_6`$, $`q_6`$, and $`v_6`$ entering the small-magnetization free-energy expansion and found, resumming the RG series by the Padé-Borel-Leroy technique, the universal critical values of the renormalized sextic couplings at the cubic fixed point. The values of the anisotropy parameters $`\delta ^{(i)}`$ have turned out to be appreciable at criticality, indicating that the cubic (non-Heisenberg) critical behavior should be, in principle, detectable in physical and computer experiments.

###### Acknowledgements.

This work was supported by the Russian Foundation for Basic Research (Grant No. 01-02-17048), by the Ministry of Education of Russian Federation via the Grant Center of Natural Sciences (Grant No. E00-3.2-132), and by the Federal Program ”Integratsiya” (Project No. A 0150). One of the authors (A. I. S.) gratefully acknowledges also the support of the International Science Foundation via Grant p00–1493.


format: text

text:
# Phase diagram of a random-anisotropy mixed-spin Ising model

## I Introduction

Apart from their relevance to the description of ferrimagnetic materials, mixed-spin models are interesting from a purely theoretical point of view, being among the simplest models to exhibit tricritical behavior. So they are especially convenient for studying the effects of inhomogeneities on the phase diagrams and the multicritical behavior of magnetic systems. From a few exact  and several approximate  calculations, we now have a good picture of the phase diagrams of mixed spin-$`\frac{1}{2}`$–spin-$`1`$ Ising models in the presence of a crystal field. Our aim in this work is to use this model Hamiltonian to investigate the effects of disorder on the location of the transition lines and the tricritical point.

The mixed-spin Ising model is described as a two-sublattice system, with spin variables $`\sigma =\pm 1`$ and $`S=0,\pm 1`$, on the sites of sublattices $`A`$ and $`B`$, respectively. Restricting the interactions to nearest-neighbors (belonging to different sublattices) and single-ion terms, the most general spin Hamiltonian in even spin space is written as

$$H=J\underset{iA,jB}{}\sigma _iS_j+D\underset{jB}{}S_j^2,$$
(1)
where the first sum is over nearest-neighbor pairs, the second sum is over the sites of sublattice $`B`$, and the parameter $`J`$ is assumed to be positive (ferromagnetic exchange). For $`D>0`$, the crystal field favors the $`S_j=0`$ states; the competition between anisotropy and ferromagnetic exchange terms leads to the appearance of a tricritical point. It should be pointed out that one needs three parameters to describe the even space of the better known spin-$`1`$ Blume-Emery-Griffiths (BEG) model. In the present case of a mixed-spin Ising system, from the point of view of the calculations, the reduction of the analysis to a two-parameter space is a particularly attractive feature.

There are exact calculations for the thermodynamic functions associated with the model Hamiltonian given by Eq. (1) on a simple chain and on some three-coordinated two-dimensional structures. On a honeycomb lattice, the problem can be mapped onto a spin-$`\frac{1}{2}`$ Ising model on a triangular lattice, which does not display a tricritical point. This mixed-spin model can also be exactly solved on a Bethe lattice  (the deep interior of Cayley tree), leading to the same results of a recent cluster-variational calculation. The results on a Bethe lattice with coordination $`q`$ indicate the absence of a tricritical point for $`q<5`$, as confirmed by Migdal–Kadanoff renormalization-group calculations. In the infinite-coordination limit of the Bethe lattice, one regains the well-known results for the tricritical point displayed by the Curie–Weiss (mean-field) version of the model. An earlier approximate effective-field calculation  predicted a tricritical point for $`q4`$, but this result has been recently challenged.

In order to analyze the effects of disorder, we consider the Hamiltonian

$$H=J\underset{iA,jB}{}\sigma _iS_j+\underset{jB}{}D_jS_j^2,$$
(2)
where $`\left\{D_j\right\}`$ is a set of independent, identically distributed random variables associated with the binary probability distribution

$$\mathrm{}(D_j)=p\delta (D_j)+(1p)\delta (D_jD).$$
(3)
With this choice of disorder, and for $`D>qJ`$, the ground state can be mapped onto a percolation problem in which the dilution affects the sites belonging to only one of the sublattices (corresponding to spin $`S=1`$). This association is easy to see if we note that a uniform crystal field $`D>qJ`$ leads to $`S_j=0`$ for all $`j`$, which breaks the connectivity between the spin-$`1/2`$ variables. The presence of a randomly located distribution of $`D=0`$ crystal fields recovers that connectivity and, for sufficiently high values of $`p`$, leads to the formation of a percolating cluster. In the rather artificial case of annealed disorder, on the honeycomb lattice, there is also an exact solution  for the thermodynamic properties of the mixed-spin model described by Eqs. (2) and (3). It is interesting to remark that this solution in the annealed case reproduces the critical concentration of the percolation problem associated with the ground state of the model with quenched (frozen) disorder, which is equivalent to the usual percolation problem on the triangular lattice. For the physically more relevant case of quenched disorder, there are approximate calculations using an effective-field theory with correlations, which point to the (expected) weakening of the tricritical behavior due to the presence of disorder.

In the present work, we first analyze the temperature–anisotropy ($`T\times D`$) phase diagram of the Curie–Weiss version (mean-field limit) of the spin Hamiltonian given by Eqs. (2) and (3). Depending on the concentration $`p`$, there appear novel transition lines and multicritical points. To include the effects of thermal fluctuations, we then resort to a standard self-consistent Bethe–Peierls approximation (which is analogous, in the case of the corresponding uniform model, to performing an exact Bethe-lattice calculation). In section II, we present the Curie–Weiss results. The Bethe–Peierls approximation is discussed in section III. In the appendices, we report an exact solution in one-dimension as well as some low-temperature expansions to supplement the numerical results from the Bethe–Peierls approximation.

## II Curie–Weiss version

The Curie–Weiss version of the mixed-spin Ising model is given by the Hamiltonian

$$H=\frac{2J}{N}\underset{iA}{}\sigma _i\underset{jB}{}S_j+\underset{jB}{}D_jS_j^2,$$
(4)
where the sums are over all sites belonging to each one of the sublattices.

For a given disorder configuration $`\{D_j\}`$, we calculate the partition function by performing a partial trace over the set of spin variables $`\{S_j\}`$. In the thermodynamic limit, we use the saddle-point method and average over disorder to obtain the free energy functional

$$\mathrm{\Psi }(\sigma )=\frac{1}{2\beta }\left[\mathrm{ln}2\frac{1}{2}(1+\sigma )\mathrm{ln}(1+\sigma )\frac{1}{2}(1\sigma )\mathrm{ln}(1\sigma )\right]\frac{1}{2\beta }\mathrm{}(D_B)\mathrm{ln}\left[1+2e^{\beta D_B}\mathrm{cosh}(\beta J\sigma )\right]𝑑D_B.$$
(5)
From the minimization of $`\mathrm{\Psi }(\sigma )`$ with respect to $`\sigma `$, we have the $`A`$-sublattice magnetization,

$$\sigma =\mathrm{tanh}\left[\beta J\mathrm{}(D_B)\frac{2\mathrm{sinh}(\beta J\sigma )}{e^{\beta D_B}+2\mathrm{cosh}(\beta J\sigma )}𝑑D_B\right],$$
(6)
where the random variable $`D_B`$ satisfies the probability distribution in Eq. (3). We can now calculate various expectation values. For example, we have

$`Q`$ $`=`$ $`{\displaystyle \mathrm{}(D_B)S_B^2𝑑D_B}`$ (7)
$`=`$ $`{\displaystyle \mathrm{}(D_B)\frac{2\mathrm{cosh}(\beta J\sigma )}{e^{\beta D_B}+2\mathrm{cosh}(\beta J\sigma )}𝑑D_B}.`$ (8)
The critical line comes from the condition

$$\frac{^2\mathrm{\Psi }}{\sigma ^2}|_{\sigma =0}=0e^\mathrm{\Delta }=2\frac{(K1)\frac{1}{3}pK^2}{1\frac{2}{3}pK^2},$$
(9)
where $`\mathrm{\Delta }=\beta D`$ and $`K=\beta J`$. The thermodynamic stability of the critical line depends on the sign of the fourth derivative of $`\mathrm{\Psi }(\sigma )`$ at $`\sigma =0`$. There is then the possibility of a tricritical point, given by the additional condition

$$\frac{^4\mathrm{\Psi }}{\sigma ^4}|_{\sigma =0}=0K^2=\frac{3+9p+\sqrt{9186p+177p^2}}{8p}.$$
(10)
The stability of the tricritical point is determined by

$$\frac{^6\mathrm{\Psi }}{\sigma ^6}|_{\sigma =0}0pp_m=0.04485\mathrm{}.,$$
(11)
which means that the tricritical behavior is suppressed by disorder concentrations larger than approximately $`4.5\%`$.

In Fig. 1, we plot some $`D\times T`$ phase diagrams for a set of typical values of the concentration $`p`$. In the uniform case ($`p=0`$), there is just a tricritical point, $`P_t`$. For $`0<pp_m=0.04485\mathrm{}`$, a tricritical point still exists (see Fig. 1, for $`p=0.04`$). However, at low temperatures and sufficiently large values of $`D`$, there appears a low-density ($`Qp`$ as $`T0`$) ferromagnetic phase, which we call ferro-II phase; at fixed values of $`D`$, an increase of temperature induces a second-order transition from the ferro-II to the paramagnetic phase. This transition is represented by a critical line that meets the first-order line at a critical end point, $`P_{ce}`$. This critical end point separates the first-order line into two distinct regions: (i) at higher temperatures, there are transitions between the usual, high-density ($`Q1`$ as $`T0`$) ferromagnetic phase (ferro-I) and the paramagnetic phase; (ii) at lower temperatures, the transitions are between the ferro-I and ferro-II phases, the first-order boundary ending at a finite-temperature critical point, $`P_{cs}`$.

For $`p_m=0.04485\mathrm{}<p<3/59=0.05084\mathrm{}`$, the tricritical point is replaced by a critical end point and a simple critical point, separated by a first-order transition line between the ferromagnetic phases (see inset in Fig. 1, for $`p=0.05`$).

For $`p3/59`$, the critical line is fully stable (see Fig. 1, for $`p=0.08`$). However, for $`p0.1`$, there still exists a small finite-temperature region where there are (first-order) transitions between the ferromagnetic phases.

## III Bethe–Peierls approximation

To give an estimate of the effects of thermal fluctuations, which are not accounted for by the Curie–Weiss calculations, we now resort to a standard self-consistent Bethe–Peierls approximation. As the model is defined on a bipartite lattice, we have to consider two distinct clusters of coordination $`q`$ (see Fig. 2). In one of them, which we call $`A`$, the central site is occupied by a spin $`\sigma =\frac{1}{2}`$, connected to $`q`$ spins of the $`S=1`$ type. In the other cluster, called $`B`$, there is a central $`S=1`$ spin surrounded by $`q`$ spin-$`\frac{1}{2}`$ variables. According to the standard prescription of the Bethe–Peierls approximation, we assume that the boundary spins of cluster $`A`$ are under the action of an effective magnetic field $`\stackrel{~}{h}_B`$ and an effective crystal field $`\stackrel{~}{D}`$, while the boundary spins of cluster $`B`$ are in an effective magnetic field $`\stackrel{~}{h}_A`$. The crystal field acting on the central site of a cluster $`B`$ is a random variable $`D_B`$. We also consider external magnetic fields, $`h_A`$ and $`h_B`$, acting on the central sites of clusters $`A`$ and $`B`$, respectively.

For the sake of simplicity, we assume that the same effective crystal field $`\stackrel{~}{D}`$ acts on all boundary spins of cluster $`A`$ (and impose self-consistency between both thermal and disorder averages associated with the two clusters). In a more refined approach, we might introduce different effective crystal fields to mimic the extended disorder of real materials. However, we expect that, in the absence of external magnetic fields, the assumption of a single $`\stackrel{~}{D}`$ is reasonable, at least in the paramagnetic phase, where there is no long-range order; in particular, we expect that the paramagnetic phase boundaries obtained by the two approaches are equivalent within the Bethe–Peierls approximation. We will see that this is indeed supported by the calculations.

The partition functions associated with the two clusters are given by

$$Z_A=e^{\gamma _A}\left[1+2e^{\stackrel{~}{\mathrm{\Delta }}}\mathrm{cosh}(\stackrel{~}{\gamma }_B+K)\right]^q+e^{\gamma _A}\left[1+2e^{\stackrel{~}{\mathrm{\Delta }}}\mathrm{cosh}(\stackrel{~}{\gamma }_BK)\right]^q$$
(12)
and

$$Z_B=\left[2\mathrm{cosh}(\stackrel{~}{\gamma }_A)\right]^q+e^{\mathrm{\Delta }_B}\left\{e^{\gamma _B}\left[2\mathrm{cosh}(\stackrel{~}{\gamma }_A+K)\right]^q+e^{\gamma _B}\left[2\mathrm{cosh}(\stackrel{~}{\gamma }_AK)\right]^q\right\},$$
(13)
where $`\gamma =\beta h`$, $`\mathrm{\Delta }=\beta D`$ and $`K=\beta J`$. The effective fields $`\stackrel{~}{\gamma }_A`$, $`\stackrel{~}{\gamma }_B`$ and $`\stackrel{~}{\mathrm{\Delta }}`$ are determined by the consistency equations

$$\sigma =\left[\sigma _j\right]_{\text{av}}=\frac{\mathrm{ln}Z_A}{\gamma _A}=\frac{1}{q}\mathrm{}(D_B)\frac{\mathrm{ln}Z_B}{\stackrel{~}{\gamma }_A}𝑑D_B,$$
(14)
$$S=\left[S_j\right]_{\text{av}}=\frac{1}{q}\frac{\mathrm{ln}Z_A}{\stackrel{~}{\gamma }_B}=\mathrm{}(D_B)\frac{\mathrm{ln}Z_B}{\gamma _B}𝑑D_B,$$
(15)
and

$$Q=\left[S_j^2\right]_{\text{av}}=\frac{1}{q}\frac{\mathrm{ln}Z_A}{\stackrel{~}{\mathrm{\Delta }}}=\mathrm{}(D_B)\frac{\mathrm{ln}Z_B}{\mathrm{\Delta }_B}𝑑D_B,$$
(16)
where $`\mathrm{}`$ and $`\left[\mathrm{}\right]_{\text{av}}`$ indicate thermal and disorder averages, respectively. We point out that the introduction of the effective crystal field $`\stackrel{~}{D}`$ is essential to achieve consistency between the equations for the two clusters.

In order to analyze the critical behavior, it is convenient to choose the magnetization $`\sigma `$, the temperature $`T`$, and the external fields $`h_B`$ and $`D_B`$, as the independent thermodynamic variables. Thus, the external field $`h_A`$ is written as a function of those variables, and the second-order transitions in zero external field ($`h_A=h_B=0`$) are given by

$$\frac{\gamma _A}{\sigma }|_{\sigma =0}=0Q_0=\frac{1}{(q1)^2\mathrm{tanh}^2K},$$
(17)
where the derivative is taken for fixed values of the remaining independent variables, and

$$Q_0Q|_{\sigma =0}=\mathrm{}(D_B)\frac{2\mathrm{cosh}^qK}{e^{\mathrm{\Delta }_B}+2\mathrm{cosh}^qK}𝑑D_B.$$
(18)
To calculate the derivative in Eq. (17), we take the (implicit) derivative of the consistency equations with respect to $`\sigma `$, imposing the condition $`\sigma =0`$ and eliminating the derivatives involving $`S`$, $`Q`$ and the effective fields. Also, for $`\sigma =0`$, we have $`S=\stackrel{~}{\gamma }_A=\stackrel{~}{\gamma }_B=0`$, since those variables are odd functions of $`\sigma `$ for $`h_A=h_B=0`$. Thus, the consistency equation for $`Q`$ leads to

$$\stackrel{~}{\mathrm{\Delta }}|_{\sigma =0}=\mathrm{ln}\left(2\frac{1Q_0}{Q_0}\mathrm{cosh}K\right),$$
(19)
and the final result is

$$\frac{\stackrel{~}{\gamma }_A}{\sigma }|_{\sigma =0}=\frac{1+\left[2(q1)q^2\right]V_0+(q1)^2V_0^2}{1+(q2)V_0+(q1)^2V_0^2},$$
(20)
where $`V_0=Q_0\mathrm{tanh}^2K`$. For $`q=2`$, Eq. (20) reproduces the exact one-dimensional expression for the $`A`$-sublattice susceptibility (see Appendix A). In fact, for $`q=2`$, it is not difficult to check that we regain all the exact one-dimensional results.

It is easy to see that, in the uniform case, corresponding to $`\mathrm{}(D_B)=\delta (D_BD)`$, the critical line is given by

$$\mathrm{\Delta }=\mathrm{ln}\left\{2\left(\mathrm{cosh}K\right)^{q2}\left[q(q2)\mathrm{cosh}^2K(q1)^2\right]\right\},$$
(21)
in agreement with the results from the Bethe-lattice  and the cluster-variational  calculations.

Using the binary distribution in Eq. (3), we obtain

$$Q_0=p\frac{2\mathrm{cosh}^qK}{1+2\mathrm{cosh}^qK}+(1p)\frac{2\mathrm{cosh}^qK}{e^\mathrm{\Delta }+2\mathrm{cosh}^qK}.$$
(22)
Therefore, the critical line is given by

$$e^\mathrm{\Delta }=2\frac{(1p)\varphi (K)}{\varphi (K)}\mathrm{cosh}^qK,$$
(23)
where

$$\varphi (K)=\frac{1}{(q1)^2}\frac{\mathrm{cosh}^2K}{\mathrm{cosh}^2K1}p\frac{2\mathrm{cosh}^qK}{1+2\mathrm{cosh}^qK}.$$
(24)
In the $`T0`$ ($`K\mathrm{}`$) limit, we have

$$e^\mathrm{\Delta }\frac{e^{qK}}{2^{q1}}\frac{(q1)^21}{1p(q1)^2},$$
(25)
which has a real solution for $`\mathrm{\Delta }`$ if

$$1p(q1)^2>0p<p_{cr}=\frac{1}{(q1)^2}.$$
(26)
This last result should be anticipated for a Bethe lattice, as we can see from the following arguments. Consider a Cayley tree where the sites belonging to every other shell, say those on odd-numbered shells, are occupied with probability $`p`$, while the remaining sites are always occupied. If $`q`$ is the coordination of the tree, the average number of paths from the root (shell $`0`$) to the first shell is given by $`p(q1)`$, while we have $`p(q1)^2`$ paths from there to the second shell. According to this reasoning, the average number of paths from the root to the $`(2n)`$th shell is given by $`p^n(q1)^{2n}`$. In order to have at least one path to the surface of the tree ($`n\mathrm{}`$), it is required that $`p(q1)^21`$, which is just the condition in Eq. (26). This result, together with the reproduction of the exact one-dimensional solution, would suggest that the present treatment also gives exact results on the Bethe lattice even in the presence of disorder. However, as remarked in previous similar treatments, this works for the paramagnetic phase only, because only then it is correct to assume that all boundary sites are under the action of the same (zero) effective field. The existence of a percolating cluster, which we do not take into account in this treatment, prevents this approximation from still giving correct results for the ordered phases.

We now consider Eq. (23), in the infinite coordination limit ($`q\mathrm{}`$, $`K0`$, $`qK=\stackrel{~}{K}`$). We then have

$$e^\mathrm{\Delta }=2\frac{\left(\stackrel{~}{K}1\right)\frac{1}{3}p\stackrel{~}{K}^2}{1\frac{2}{3}p\stackrel{~}{K}^2},$$
(27)
which agrees with Eq. (9) for the Curie–Weiss version of the model.

The tricritical points are determined by Eq. (17) supplemented by the condition

$`{\displaystyle \frac{^3\gamma _A}{\sigma ^3}}|_{\sigma =0}=0,`$

which is equivalent to

$$\frac{2q^210q+6}{(q1)^5\mathrm{tanh}^2K}+3qW_0\mathrm{tanh}^2K=\frac{(q2)(q3)}{(q1)^3},$$
(28)
where $`W_0`$ is given by

$$W_0=\mathrm{}(D_B)\left(\frac{2\mathrm{cosh}^qK}{e^{\mathrm{\Delta }_B}+2\mathrm{cosh}^qK}\right)^2𝑑D_B.$$
(29)
The tricritical points are stable if

$`{\displaystyle \frac{^5\gamma _A}{\sigma ^5}}|_{\sigma =0}>0.`$

To calculate this derivative, we again take implicit derivatives of the consistency equations (up to fifth order) with respect to $`\sigma `$, at $`\sigma =0`$, and eliminate all derivatives involving $`S`$, $`Q`$ and the effective fields. Unlike the previous analysis, we have not been able to obtain closed-form expressions for the stability condition of the tricritical point, but it is not difficult to perform a number of numerical calculations.

For the uniform model, we have $`W_0=Q_0^2`$. Therefore, Eq. (28) takes the form

$$\mathrm{tanh}K=\frac{1}{q1}\sqrt{\frac{5q3}{q3}},$$
(30)
which is again identical to the result obtained from the Bethe-lattice  and the cluster-variational  calculations. Notice that this equation has real solutions only if $`q>4.561553\mathrm{}`$. Thus, the Bethe–Peierls approximation does not predict a tricritical point for the square lattice ($`q=4`$).

For the binary distribution in Eq. (3), we have

$$W_0=Q_0^2\left[1+\frac{p}{1p}\left(1\frac{1}{Q_0}\frac{2\mathrm{cosh}^qK}{1+2\mathrm{cosh}^qK}\right)^2\right].$$
(31)
In the infinite-coordination limit we can write

$$W_0=\frac{1}{\stackrel{~}{K}^4}\left[1+\frac{p}{1p}\left(1\frac{2}{3}\stackrel{~}{K}^2\right)^2\right],$$
(32)
which leads to the equation

$$\stackrel{~}{K}^23\left[1+\frac{p}{1p}\left(1\frac{4}{3}\stackrel{~}{K}^2+\frac{4}{9}\stackrel{~}{K}^4\right)\right]2=0,$$
(33)
at the tricritical point. Indeed, one of the solutions of this equation corresponds to Eq. (10), for the Curie–Weiss version of the model, while the other solution represents a thermodynamically unstable situation.

In Table I, for various values of the coordination number $`q`$, and using the binary distribution given by Eq. (3), we give the corresponding values of the concentration $`p_m`$, at which the tricritical point becomes unstable, and the critical percolation concentration $`p_{cr}`$. We see that, for $`q10`$, the tricritical behavior is suppressed for $`p_m<p_{cr}`$, while, for $`q11`$, that suppression occurs for $`p_m>p_{cr}`$. As shown in Table I, $`p_m`$ increases with $`q`$, which indicates that disorder is more effective for small coordination numbers.

As the effects of binary disorder strongly depend on the coordination, we now discuss the phase diagrams for the typical cases.

For $`q=3`$ and $`4`$, there are no tricritical points. The $`D\times T`$ phase diagram displays just a fully stable critical line. The main effect of disorder is to make the paramagnetic phase unstable at $`T=0`$, regardless of the value of $`D`$, for $`p`$ larger than the critical percolation concentration $`p_{cr}`$. The phase diagrams in Fig. 3, for $`q=3`$, are in qualitative agreement with the exact results for the honeycomb lattice (which is also three-coordinated) under annealed disorder. At $`T=0`$, there is even quantitative agreement with the value of the critical crystal field at $`p_{cr}`$, given by $`D_{cr}=5J/3`$, although of course this agreement does not extend to the value of $`p_{cr}`$ itself. Our results for $`q=3`$ and $`q=4`$ are also in qualitative agreement with those obtained by a real-space renormalization-group approach for the two-dimensional Blume-Emery-Griffiths model in a random crystal field.

For $`5q10`$, the concentration $`p_m`$ above which the tricritical point becomes unstable is lower than $`p_{cr}`$. For $`p<p_m`$, disorder depressses the tricritical temperature, and shortens the first-order transition line. For $`p_m<p<p_{cr}`$, the tricritical point is replaced by a critical end point, $`P_{ce}`$, and a simple critical point, $`P_{cs}`$, as in the Curie–Weiss version of the model. However, the paramagnetic phase is stable at $`T=0`$ if $`D>qJ`$, and the first-order line reaches $`D=qJ`$ at $`T=0`$. As $`p`$ increases, first the critical end point $`P_{ce}`$ and then the simple critical point $`P_{cs}`$ reach the $`T=0`$ axis, at values of $`p`$ which can be determined by a low-temperature expansion of the consistency equations (see Appendix B). In Fig. 4, we plot the $`D\times T`$ phase diagram for $`q=6`$ and $`p=0.011`$. To determine the first-order lines shown in that figure, we numerically solve the consistency equations to obtain the conditions $`h_A\left(\sigma _1\right)=h_A\left(\sigma _2\right)=0`$ and

$$\underset{\sigma _1}{\overset{\sigma _2}{}}h_A\left(\sigma \right)𝑑\sigma =0,$$
(34)
which correspond to a Maxwell construction.

For $`q11`$, we have $`p_m>p_{cr}`$, so the behavior of the system is quite similar to the predictions of the Curie–Weiss version of the model.

## IV Conclusions

We performed detailed calculations for the phase diagram of a random-anisotropy mixed-spin Ising model both in the mean-field limit (Curie–Weiss version of the model), in which thermal fluctuations are neglected, and according to a standard self-consistent Bethe–Peierls approximation (which turns out to be exact in one dimension). For a binary distribution of crystal fields, we obtained closed-form expressions for the critical lines and the location of the tricritical points. Depending on the concentration $`p`$, the mean-field results for the $`D\times T`$ phase diagrams predict novel first-order lines and multicritical points (besides a ferromagnetic region, for all values of the crystal field, extending down to the lowest temperatures). The Bethe–Peierls approximation shows that this additional ferromagnetic region is suppressed for concentrations below a certain percolation threshold. Also, the Bethe–Peierls results point out to the absence of a tricritical behavior for lattices with coordination $`q4`$. All results reported in this paper are in agreement with general predictions for the effects of disorder on first-order transitions and multicritical points (for a recent review, see a paper by Cardy ).

###### Acknowledgements.

We thank T. A. S. Haddad for useful discussions. This work was partially financed by the Brazilian agencies FAPESP and CNPq.

## A Exact solution in one-dimension

For an open chain with $`N+1`$ sites ($`N`$ even), and in zero external field, the Hamiltonian of the mixed-spin Ising model can be written as

$$H=J\underset{j=1}{\overset{N/2}{}}\left(\sigma _jS_j+S_j\sigma _{j+1}\right)+\underset{j=1}{\overset{N/2}{}}D_jS_j^2.$$
(A1)
Given a disorder configuration $`\{D\}=\{D_1,\mathrm{},D_{N/2}\}`$, we perform a partial trace over the spin variables $`\{S_j\}`$ to write

$`Z\{D\}`$ $`=`$ $`{\displaystyle \underset{\{\sigma \}}{}}{\displaystyle \underset{\{S\}}{}}e^{\beta H}`$ (A2)
$`=`$ $`{\displaystyle \underset{\{\sigma \}}{}}{\displaystyle \underset{j=1}{\overset{N/2}{}}}\left\{1+2e^{\mathrm{\Delta }_j}\mathrm{cosh}[K(\sigma _j+\sigma _{j+1})]\right\},`$ (A3)
where $`K=\beta J`$ and $`\mathrm{\Delta }_j=\beta D_j`$. Introducing a prefactor $`A_j`$,

$$A_j^2=\left(1+2e^{\mathrm{\Delta }_j}\right)\left[1+2e^{\mathrm{\Delta }_j}\mathrm{cosh}(2K)\right]$$
(A4)
and an effective interaction $`\stackrel{~}{K}_j`$, such that

$$e^{2\stackrel{~}{K}_j}=\frac{1+2e^{\mathrm{\Delta }_j}\mathrm{cosh}(2K)}{1+2e^{\mathrm{\Delta }_j}},$$
(A5)
the partition function can be written as the factorized form

$`Z\{D\}`$ $`=`$ $`{\displaystyle \underset{\{\sigma \}}{}}{\displaystyle \underset{j=1}{\overset{N/2}{}}}A_je^{\stackrel{~}{K}_j\sigma _j\sigma _{j+1}}`$ (A6)
$`=`$ $`{\displaystyle \underset{j=1}{\overset{N/2}{}}}2\left[1+2e^{\mathrm{\Delta }_j}\mathrm{cosh}^2K\right].`$ (A7)
From Eq. (A7), we obtain the thermal average

$$S_j^2_{\{D\}}=\frac{\mathrm{ln}Z}{\mathrm{\Delta }_j}=\frac{2e^{\mathrm{\Delta }_j}\mathrm{cosh}^2K}{1+2e^{\mathrm{\Delta }_j}\mathrm{cosh}^2K},$$
(A8)
which depends on the value of the crystal field on the $`j`$th site only. Since we are considering a nearest-neighbor one-dimensional model in zero field, the thermal averages $`S_j`$ and $`\sigma _j`$ are zero. Performing the disorder average, we obtain the expectation value

$$Q=S_j^2\underset{i=1}{\overset{N/2}{}}\mathrm{}(D_i)dD_i=\mathrm{}(D_j)S_j^2_{\{D\}}𝑑D_j.$$
(A9)
For a given disorder configuration, the magnetic susceptibilities of the $`\sigma `$ and $`S`$ sublattices are given by

$$\chi _\sigma \{D\}=\frac{1}{T}\underset{N\mathrm{}}{lim}\frac{2}{N+2}\underset{j=1}{\overset{\frac{N}{2}+1}{}}\underset{k=1}{\overset{\frac{N}{2}+1}{}}\sigma _j\sigma _k_{\{D\}}$$
(A10)
and

$$\chi _s\{D\}=\frac{1}{T}\underset{N\mathrm{}}{lim}\frac{2}{N}\underset{j=1}{\overset{N/2}{}}\underset{k=1}{\overset{N/2}{}}S_jS_k_{\{D\}}.$$
(A11)
The two-spin correlation functions,

$$\sigma _j\sigma _k_{\{D\}}=\frac{1}{Z\{D\}}\underset{\{\sigma \}}{}\underset{\{S\}}{}\sigma _j\sigma _ke^{\beta H}$$
(A12)
and

$$S_jS_k_{\{D\}}=\frac{1}{Z\{D\}}\underset{\{\sigma \}}{}\underset{\{S\}}{}S_jS_ke^{\beta H},$$
(A13)
can be calculated if we introduce the transformation

$$\tau _j=\sigma _j\sigma _{j+1}\text{with}\tau _0=\sigma _1.$$
(A14)
After some algebraic manipulations, for $`j<k`$, we have

$$\sigma _j\sigma _k_{\{D\}}=\underset{i=j}{\overset{k1}{}}\frac{2\mathrm{sinh}^2K}{e^{\mathrm{\Delta }_i}+2\mathrm{cosh}^2K}$$
(A15)
and

$`S_jS_k_{\{D\}}`$ $`=`$ $`\left({\displaystyle \underset{i=j,k}{}}{\displaystyle \frac{\mathrm{sinh}2K}{e^{\mathrm{\Delta }_i}+2\mathrm{cosh}^2K}}\right)`$ (A17)
$`\times {\displaystyle \underset{i=j+1}{\overset{k1}{}}}{\displaystyle \frac{2\mathrm{sinh}^2K}{e^{\mathrm{\Delta }_i}+2\mathrm{cosh}^2K}},`$
from which we obtain the expectation values

$`g_\sigma (\left|kj\right|)`$ $`=`$ $`{\displaystyle \sigma _j\sigma _k_{\{D\}}\underset{i=1}{\overset{N/2}{}}\mathrm{}(D_i)dD_i}`$ (A18)
$`=`$ $`\left(Q\mathrm{tanh}^2K\right)^{\left|kj\right|}`$ (A19)
and

$`g_s(\left|kj\right|)`$ $`=`$ $`{\displaystyle S_jS_k_{\{D\}}\underset{i=1}{\overset{N/2}{}}\mathrm{}(D_i)dD_i}`$ (A20)
$`=`$ $`Q\left(Q\mathrm{tanh}^2K\right)^{\left|kj\right|},`$ (A21)
which depend on the distance between sites $`j`$ and $`k`$. The expectation values of the susceptibilities are given by

$$\left[\chi _\sigma \right]_{\text{av}}=\frac{1}{T}\left[1+2\underset{r=1}{\overset{\mathrm{}}{}}g_\sigma \left(r\right)\right]=\frac{1}{T}\frac{1+Q\mathrm{tanh}^2K}{1Q\mathrm{tanh}^2K}$$
(A22)
and

$$\left[\chi _s\right]_{\text{av}}=\frac{1}{T}\left[Q+2\underset{r=1}{\overset{\mathrm{}}{}}g_s\left(r\right)\right]=\frac{Q}{T}\frac{1+Q\mathrm{tanh}^2K}{1Q\mathrm{tanh}^2K},$$
(A23)
where $`Q`$ is determined by Eq. (A9).

## B Low-temperature expansion

For the binary distribution, in the low-temperature limit ($`K=\beta J1`$), if we neglect terms of order $`\mathrm{exp}\left(2K\right)`$ and higher, the consistency equations (14)-(16) for cluster $`A`$ lead to the expressions

$$\gamma _A=\frac{1}{2}\mathrm{ln}\frac{1+\sigma }{1\sigma }+\frac{q}{2}\mathrm{ln}\frac{1C_+}{1C_{}},$$
(B1)
$$S=\frac{1+\sigma }{2}C_+\frac{1\sigma }{2}C_{},$$
(B2)
and

$$Q=\frac{1+\sigma }{2}C_++\frac{1\sigma }{2}C_{},$$
(B3)
where

$$C_\pm =\frac{e^{\pm \stackrel{~}{\gamma }_B}}{e^{\stackrel{~}{\mathrm{\Delta }}K}+e^{\pm \stackrel{~}{\gamma }_B}}.$$
(B4)
For cluster $`B`$, we have

$`\sigma `$ $`=`$ $`p\mathrm{tanh}\left(q\stackrel{~}{\gamma }_A\right)`$ (B6)
$`+(1p){\displaystyle \frac{\tau (\stackrel{~}{\gamma }_A)\mathrm{tanh}(\stackrel{~}{\gamma }_A)+\delta \mathrm{tanh}(q\stackrel{~}{\gamma }_A)}{\tau (\stackrel{~}{\gamma }_A)+\delta }},`$
$$S=p\mathrm{tanh}\left(q\stackrel{~}{\gamma }_A\right)+(1p)\frac{\delta }{\tau (\stackrel{~}{\gamma }_A)+\delta }\mathrm{tanh}\left(\stackrel{~}{\gamma }_A\right),$$
(B7)
$$Q=p+(1p)\frac{\delta }{\tau (\stackrel{~}{\gamma }_A)+\delta },$$
(B8)
where

$$\delta =\mathrm{exp}(qK\mathrm{\Delta }),$$
(B9)
and

$$\tau (x)=\frac{2^q}{\left(1+\mathrm{tanh}x\right)^q+\left(1\mathrm{tanh}x\right)^q}.$$
(B10)
Solving Eqs. (B2) and (B3) for $`C_\pm `$ in terms of $`\sigma `$, $`S`$ and $`Q`$, and using Eqs. (B6)-(B8), we can write Eq. (B1) in the form

$$\gamma _A(\sigma )=\frac{1q}{2}\mathrm{ln}\frac{1+\sigma }{1\sigma }+q\stackrel{~}{\gamma }_A(\sigma ),$$
(B11)
where $`\stackrel{~}{\gamma }_A(\sigma )`$ is determined from the solution of Eq. (B6). Notice that, according to Eqs. (B11) and (B6), $`\gamma _A(\sigma )`$ and $`\stackrel{~}{\gamma }_A(\sigma )`$ depend on the temperature through the parameter $`\delta `$ only. As $`T0`$, this parameter goes to zero (if $`D>qJ`$), or infinity (if $`D<qJ`$), except in the vicinity of the point $`P_0`$ with coordinates $`D=qJ`$, and $`T=0`$, where $`\delta `$ can assume any value.

Since the equation of state (B11) becomes asymptotically exact as $`T0`$, it can be used to determine the values of $`p`$ at which the critical end point and the simple critical point reach $`P_0`$, and thus disappear. To do that calculation, we impose the conditions

$$\gamma _A(\sigma _e)=\frac{\gamma _A}{\sigma }|_{\sigma =\sigma _e}=\frac{^2\gamma _A}{\sigma ^2}|_{\sigma =\sigma _e}=0,$$
(B12)
from which we obtain the values of $`\sigma _e`$, $`\delta _e`$ and $`p_e`$ at which the critical end point reaches $`P_0`$, and the conditions

$$\gamma _A(\sigma _s)=\frac{\gamma _A}{\sigma }|_{\sigma =\sigma _s}=_0^{\sigma _s}\gamma _A(\sigma )𝑑\sigma =0,$$
(B13)
which give the corresponding values $`\sigma _s`$, $`\delta _s`$ and $`p_s`$ for the simple critical point.


text:
# Information losses in continuous variable quantum teleportation

## Abstract

It is shown that the information losses due to the limited fidelity of continuous variable quantum teleportation are equivalent to the losses induced by a beam splitter of appropriate reflectivity.

Quantum teleportation allows the transmission of an unknown quantum state by combining the nonlocal quantum coherence of entangled states with the transmission of classical information obtained in a joint measurement of the unknown input state and one part of the entangled pair . Ideally, the classical information transmitted contains no information about the input state and the output state is exactly identical with the input state. However, this ideal form of quantum teleportation requires maximal entanglement. In continuous variable quantum teleportation , only non-maximal entanglement is available. As a consequence, the output state is not perfectly identical to the input state while the statistics of the joint measurement depend on the properties of the input state . The classical information channel then carries information on the input state which may be extracted, e.g. in order to eavesdrop on a quantum communication channel . The relationship between the measurement information extracted and the change of the quantum state can be described in terms of a measurement dependent transfer operator . In the following, the transfer operator describing the continuous variable teleportation process is derived and the equivalence with a feedback compensated beam splitter is established. The information obtained in quantum teleportation can then be identified with the reflected amplitude at the beam splitter while the output state of the teleportation corresponds to the transmitted beam, displaced by the feedback. The loss of quantum information due to the limited fidelity of the teleportation process is thus shown to be equivalent to the loss of quanta at a beam splitter.

As illustrated in figure 1 (a), quantum teleportation transfers an unknown quantum state in mode A using an entangled state of a reference mode R and an output mode B. For continuous variables, this is achieved by measuring the difference $`\widehat{x}_{}=\widehat{x}_A\widehat{x}_R`$ and the sum $`\widehat{y}_+=\widehat{y}_A+\widehat{y}_R`$ of the orthogonal quadrature components of the input mode A and the reference mode R. The quantum state $`\psi _\text{B}(\beta )`$ of the output mode B is then conditioned by the input state $`\psi _\text{A}`$ in mode A and the measurement result $`\beta =x_{}+iy_+`$ which ideally defines an eigenstate $`\beta (A,R)`$ of modes A and R. Reallistically, the finite resolution due to limited detector efficiencies and other technical problems of the measurement will result in some additional noise, which can be simulated by a classical Gaussian error. In the following, however, we will focus on the ideal quantum limit of the information transfer in order to determine the quantum state distortions originating from non-maximal entanglement. Using $`q(R,B)`$ to denote the initial entangled state of modes R and B, this conditional state in B can be written as

$$\psi _\text{B}(\beta )=\beta (A,R)\psi _Aq(R,B).$$
(1)
Note that the output state $`\psi _\text{B}(\beta )`$ is not normalized, since the probability $`P(\beta )`$ of the measurement outcome $`\beta `$ is given by

$$P(\beta )=\psi _\text{B}(\beta )\psi _\text{B}(\beta ).$$
(2)
Making use of the displacement operator $`\widehat{D}(\beta )`$ and the photon number expansion of entanglement, the eigenstates $`\beta (A,R)`$ and the entangled state $`q(R,B)`$ can be expressed as

$`\beta (A,R)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\pi }}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\widehat{D}_A(\beta )n;n_{A,R}`$ (3)
$`q(R,B)`$ $`=`$ $`\sqrt{1q^2}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}q^nn;n_{R,B},`$ (4)
where the entanglement coefficient $`q`$ provides a quantitative measure characterizing the degree of entanglement obtained by parametric amplification. The quantum state $`\psi _\text{B}(\beta )`$ of the output mode B conditioned by the measurement of modes A and R is then given by

$$\psi _B(\beta )=\sqrt{\frac{1q^2}{\pi }}\underset{n=0}{\overset{\mathrm{}}{}}q^nnn\widehat{D}_A(\beta )\psi _A.$$
(5)
For $`q=1`$ this state is a copy of the input state displaced by a field difference of $`\beta `$. Therefore the final step of quantum teleportation is the reversal of this displacement by the addition of a coherent field amplitude $`g\beta `$ to obtain the final output state, $`\psi _{\text{out}}(\beta )=\widehat{D}(g\beta )\psi _B(\beta )`$. The gain factor $`g`$ allows an adjustment of the teleported amplitude . $`g=1`$ reproduces the average input amplitude in the output and thus optimizes the fidelity for the teleportation of high amplitude coherent states. As discussed in a previous paper , the conditional output of quantum teleportation can be described using the transfer operator $`\widehat{T}_g(\beta )`$,

$`\psi _{\text{out}}(\beta )`$ $`=`$ $`\widehat{T}_g(\beta )\psi _A`$ (6)
$`P(\beta )`$ $`=`$ $`\psi _A\widehat{T}_g^{}(\beta )\widehat{T}_g(\beta )\psi _A`$ (7)
with $`\widehat{T}_g(\beta )=\sqrt{{\displaystyle \frac{1q^2}{\pi }}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}q^n\widehat{D}(g\beta )nn\widehat{D}(\beta ).`$ (8)
This transfer operator characterizes the teleportation process of an arbitrary quantum state by correlating the extracted information $`\beta `$ with the quantum information in the output state $`\psi _{\text{out}}(\beta )`$. When applied to a coherent state $`\alpha `$, the result reads

$`\widehat{T}_g(\beta )\alpha `$ $`=`$ $`\sqrt{{\displaystyle \frac{1q^2}{\pi }}}\mathrm{exp}\left((1q^2){\displaystyle \frac{|\alpha \beta |^2}{2}}\right)`$ (10)
$`\times \mathrm{exp}\left((1gq){\displaystyle \frac{\alpha \beta ^{}\beta \alpha ^{}}{2}}\right)q\alpha +(gq)\beta .`$
For coherent states, the output is also a coherent state with an amplitude give by the sum of the attenuated input amplitude $`q\alpha `$ and a measurement dependent displacement of $`(gq)\beta `$. Since the quantum state in the output depends on the randomly varying measurement result $`\beta `$, the teleportation output is generally a statistical mixture of different coherent states. Only in the special case of $`g=q`$, the amplitude $`q\alpha `$ of the coherent output state does not depend on $`\beta `$ and the output is a well defined pure state even if the output is averaged over all measurement results $`\beta `$.

The attenuation of the signal amplitude described by equation (10) corresponds to the losses induced by a beam splitter with a reflectivity of $`1q^2`$. For the special case of $`g=q`$, this property of teleportation has been pointed out previously by Polkinghorne and Ralph , based on an analysis of the quantum fluctuations in the teleportation. In the following, we will generalize this analogy by deriving the proper transformation of quantum states in a beam splitter measurement and considering the possibility of compensating the beam splitter losses by feedback. This formalism includes all the details necessary for an evaluation of the information obtained on the system and the minimal back-action on the signal field. If the reflected amplitude of $`\sqrt{1q^2}\alpha `$ is measured by eight port homodyne detection, the corresponding positive operator valued measure is given by projections onto the non-normalized, non-orthogonal coherent states

$`P(\beta )`$ $`=`$ $`\sqrt{{\displaystyle \frac{1q^2}{\pi }}}\sqrt{1q^2}\beta `$ (11)
with $`{\displaystyle d^2\beta P(\beta )P(\beta )}=\widehat{1},`$ (12)
and the transmitted state reads

$`\psi _{\text{trans}}(\beta )`$ $`=`$ $`\sqrt{{\displaystyle \frac{1q^2}{\pi }}}\sqrt{1q^2}\beta \sqrt{1q^2}\alpha q\alpha `$ (13)
$`=`$ $`\sqrt{{\displaystyle \frac{1q^2}{\pi }}}\mathrm{exp}\left((1q^2){\displaystyle \frac{|\alpha \beta |^2}{2}}\right)`$ (15)
$`\times \mathrm{exp}\left((1q^2){\displaystyle \frac{\alpha \beta ^{}\beta \alpha ^{}}{2}}\right)q\alpha ,`$
where only the probability amplitude of $`\psi _{\text{trans}}(\beta )`$ depends on the outcome $`\beta `$. This result corresponds to equation (10) if the gain $`g`$ is equal to the entanglement coefficient $`q`$. The effects of quantum teleportation at a gain of $`g=q`$ are therefore identical to the effects of a field measurement by eight port homodyne detection performed on the reflected part of the signal field using a beam splitter of reflectivity $`R=1q^2`$. At other gain coefficients, quantum teleportation is equivalent to a feedback compensated beam splitter measurement in which the losses induced in the transmitted beam are compensated using a linear feedback based on the measurement result $`\beta `$ obtained from the reflected light . At a feedback amplitude of $`f\beta `$, the output state of the feedback compensated beam splitter reads

$$\widehat{D}\left(f\beta \right)\psi _{\text{trans}}(\beta )=\widehat{T}_{g=f+q}(\beta )\alpha .$$
(16)
In particular, a gain of $`g=1`$ corresponds to a beam splitter feedback amplitude of $`(1q)\beta `$. In this case, the measurement operator of the beam splitter setup is hermitian, minimizing the back-action of the measurement to the minimal noise required by the Heisenberg principle .

Since all quantum states may be expanded in terms of coherent states, equation (16) proofs the equivalence of continuous variable quantum teleportation and feedback compensated beam splitting with respect to both the changes in the transmitted state and the information obtained in the measurement of $`\beta `$. In the special case of $`g=q`$, no additional photons are created in the teleportation process, indicating that all photons emitted into the output field B by the parametric amplifier are reabsorbed in the displacement transformation. This effect allows a teleportation of the vacuum with a fidelity of one, making a more reliable distinction of signal pulses from a vacuum background possible. The transmission probability for photons teleported at $`g=q`$ is equal to $`q^2`$. The loss of quantum information in continuous variable teleportation can thus be expressed in terms of photon losses. In an experimental realization of quantum teleportation, the case of $`g=q`$ can be used to characterize the performance of the setup. Specifically, the point at which $`g=q`$ can be found by minimalizing the output intensity at a vacuum input. The remaining intensity at that point arises from the finite resolution of the measurement, imperfect phase matching, and similar technical problems in the optical setup. It is then possible to separate quantum noise effects from the classical noise contributions in the teleportation setup.

At $`g>q`$, the loss of quantum information is compensated by the classical information obtained from the measurement of $`\beta `$. However, the original quantum state cannot be restored by this purely classical manipulation, limiting the achievable fidelity to a value well below one. At $`g<q`$, the displacement actually reduces the quantum information in the output further, until at $`g=0`$, there is no correlation between the input state and the output density matrix formed by integrating over all measurement results $`\beta `$. In general, the effect of quantum teleportation on the one photon state $`1`$ is given by

$$\widehat{T}_g(\beta )1=\sqrt{\frac{1q^2}{\pi }}\mathrm{exp}\left((1q^2)\frac{|\beta |^2}{2}\right)\widehat{D}\left((gq)\beta \right)\left((1q^2)\beta ^{}0+q1\right).$$
(17)
This representation of the teleported state has only two components, corresponding to a displaced vacuum and a displaced photon number state, respectively. At $`g=q`$, the displacement is zero and the two components correspond to the actual loss or transmission of the photon. At other gain factors, the coherent displacement can generate photon numbers $`n>1`$ in the output. Figure 1 shows a schematic comparison of the experimental setups for quantum teleportation setup and for the compensated beam splitter. Both methods employ linear transformations on the input field mode and two vacuum modes, extracting information on the unknown input field from the homodyne detection measurements on two of the output modes. However, in quantum teleportation, the only physical connection between the input field and the output field is given by the measurement dependent displacement. While the beam splitter transmits the attenuated input field by a direct physical interaction, quantum teleportation achieves the same result by combining the entanglement with the classical information $`\beta `$. The entanglement coefficient $`q`$ is the measure of non-maximal entanglement which corresponds to the attenuation of the transmitted amplitude at an equivalent beam splitter .

In conclusion, the analysis of the transfer operator $`\widehat{T}_g(\beta )`$ shows that the information transfer in quantum teleportation is essentially equivalent to a feedback compensated beam splitter. This result clarifies the nature of information losses in quantum teleportation and allows an assessment of the information extracted with respect to applications such as continuous variable eavesdropping . Moreover, the analogy provides a quantification of the information transfer properties of non-maximal entanglement in terms of the photon transmission probability $`q^2`$ and simplifies the derivation of quantum coherent transfer properties for few photon inputs, e.g. for the entanglement swapping scheme discussed in .

One of us (HFH) would like to acknowledge support from the Japanese Society for the Promotion of Science, JSPS.