text:
# Critical parameters and universal amplitude ratios of two-dimensional spin-𝑆 Ising models using high- and low-temperature expansions

## I Introduction

The properties of the spin $`S=1/2`$ two-dimensional Ising model with nearest-neighbor interactions in zero magnetic field, have been extensively explored in the last six decades. Much more modest efforts have been devoted to the study of the simplest generalizations of the model to spin $`S>1/2`$. The main reason is probably that these models are not known to be solvable or, at least, to have any simple duality property which can help to extend the small body of information coming from numerical methods of limited accuracy such as stochastic simulations, series expansions or transfer-matrix calculations.

The first important result from a comparative study of Ising models for different values of the spin came from pioneering work by Domb and Sykes dosy . They analysed the high-temperature(HT) expansion of the susceptibility $`\chi (K;S)`$ through $`O(K^6)`$ in the three-dimensional case and conjectured that the value of the critical exponent $`\gamma (S)`$ is independent of the spin magnitude. This was the first step towards the modern formulation of the critical-universality hypothesis. Similar analyses were soon repeated by other authors using both HTmoore ; camp  and low-temperature (LT) expansionsfoxgu  for two-dimensional systems. Unfortunately, the series derived in those years were rather short and therefore the results of the analyses could not reach a sufficient accuracy or were inconclusive. It was only in 1980 that Nickelnick21 ; nr90  finally extended through $`O(K^{21})`$ the HT series in two dimensions on the square (sq) lattice and in three dimensions on the body-centered cubic lattice. The expansions of $`\chi (K;S)`$ and of the second moment of the spin-spin correlation function $`\mu _2(K;S)`$ were then published only for $`S=1/2,1,2,\mathrm{}`$. More recently also the LT expansions on the sq lattice for $`S=1/2,1,3/2,2,5/2,3`$ were considerably extended by Enting, Jensen and Guttmannejg . We have summarized in Table 1 and Table 2 the state of HTcamp ; bakin  and LTejg  expansions before our work.

In spite of the very large number of LT expansion coefficients now available, the analysis of the series remains arduous due to occurrenceejg  of numerous unphysical singularities in the complex temperature planefisherzeroes  which are closer to the origin than the physical singularity and whose structure becomes increasingly complicated with $`S`$. As a consequence, the LT study of Ref.ejg  has been an alarming lesson on the subtleties in the analysis of slowly convergent series more than a source of accurate estimates of the critical parameters of the models.

Many intriguing indications and conjectures about the structure of these unphysical singularities also came in the same period from work by Matveev and Shrockms  who examined the spin $`S`$ models on various two-dimensional lattices using transfer-matrix methods.

Here we discuss some results of an analysis of HT series for the sq lattice recently extendedbcesse25  by linked-cluster expansion techniques.

For the nearest-neighbor correlation function $`G(K;S)`$, for $`\chi (K;S)`$, and $`\mu _2(K;S)`$ our series reach order $`K^{25}`$, while for the second field derivative of the susceptibility $`\chi _4(K;S)`$ they extend through $`O(K^{23})`$. In order to make alternative analyses possible, our vast collection of series data both for two- and three-dimensional lattices was made easily accessiblebcol  on-line for $`S=1/2,1,3/2,2,5/2,3,7/2,4,5,\mathrm{}`$. It should be noted that HT and LT expansions as extensive as those obtained by Nickel et al. in Ref.nick21  and more recently in Ref.orri  (only for $`\chi `$) in the very special case of the (partially solvable) two-dimensional $`S=1/2`$ Ising model seem presently beyond reach for $`S>1/2`$.

The HT series show somewhat simpler and faster convergence properties than the LT series, because the behavior of the coefficients is dominated by the physical singularity. Although, even in this case, these favorable properties slightly deteriorate for $`S>2`$, we can hope to determine basic HT critical parameters with a reasonable accuracy for various values of $`S`$. Moreover it is also worthwhile to reconsider the LT expansions of Ref.ejg  for the sq lattice, because by relying on the results of our HT analysis, we can improve some estimates of the LT critical parameters and thus obtain new determinations of universal combinations aha  of LT and HT amplitudes. No theoretical surprises are expected from this analysis, however we believe it is still useful to improve the rather modest numerical precision presently available even for basic critical parameters like the critical temperatures, to determine various critical amplitudes for which no estimates are yet known and to use our results to test with higher accuracy the validity of hyperscaling and of universality with respect to the magnitude of the spin.

Almost all the computational effort in extending series for the two-dimensional Ising model for $`S>1/2`$ has been devoted to square-lattice series. However by making use of the theory of lattice-lattice scaling, as developed by Betts et al.bgj  and extended by Gaunt and Guttmanngg , using our estimates of the critical amplitudes on the square lattice, we are able to calculate the corresponding amplitudes on other two-dimensional lattices to precisely the same precision as they are known for the square lattice.

## II  The spin-$`S`$ Ising models

The spin-$`S`$ Ising models with nearest-neighbor interaction are defined by the Hamiltonian:

$$H\{s\}=\frac{J}{2}\underset{(\stackrel{}{x},\stackrel{}{x}^{})}{}s(\stackrel{}{x})s(\stackrel{}{x}^{})h\underset{\stackrel{}{x}}{}s(\stackrel{}{x})$$
(1)
where $`J`$ is the exchange coupling, and $`s(\stackrel{}{x})=s^z(\stackrel{}{x})/S`$ with $`s^z(\stackrel{}{x})`$ a classical spin variable at the lattice site $`\stackrel{}{x}`$, taking the $`2S+1`$ values $`S,S+1,\mathrm{},S1,S`$. The sum runs over all nearest-neighbor pairs of sites. We shall restrict ourselves to the square lattice and consider expansions either in the usual HT variable $`K=J/k_BT`$ and in the natural LT variable $`u(S)=exp(K/S^2)`$. Here $`T`$ is the temperature, $`k_B`$ the Boltzmann constant, and $`K`$ will be called “inverse temperature” for brevity. In the critical region we shall also refer to the standard reduced-temperature variable $`t(S)=T/T_c(S)1=K_c(S)/K1`$.

In the HT phase, the basic observables are the connected $`2n`$-spin correlation functions. Our seriesbcol  cover quantities related to the two-spin correlation functions $`s(\stackrel{}{x})s(\stackrel{}{y})_c`$ and to the four-spin correlation functions $`s(\stackrel{}{x})s(\stackrel{}{y})s(\stackrel{}{z})s(\stackrel{}{t})_c`$.

In the LT phase the symmetry is broken and the $`n`$-spin correlations are non trivial also for odd $`n`$. In particular, we shall reconsider the LT expansions of the magnetization, the susceptibility and the specific heat derived for $`S=1,3/2,2,5/2,3`$ in Ref.ejg .

The spontaneous magnetization is defined by

$$M(T;S)=\underset{h0+}{lim}s(\stackrel{}{0}).$$
(2)
The internal energy per spin is given in terms of the nearest-neighbor correlation function by

$$U(T;S)=\frac{qJ}{2}s(\stackrel{}{0})s(\stackrel{}{\delta })$$
(3)
where $`\stackrel{}{\delta }`$ is a nearest-neighbor lattice vector and $`q`$ is the lattice coordination number.

The specific heat is the temperature-derivative of the internal energy at fixed zero external field

$$C_H(T;S)=\frac{dU(T;S)}{dT}.$$
(4)
In terms of $`\chi (T;S)`$, the zero-field reduced susceptibility,

$$\chi (T;S)=\underset{\stackrel{}{x}}{}s(\stackrel{}{0})s(\stackrel{}{x})_c$$
(5)
and of $`\mu _2(T;S)`$, the second moment of the correlation function,

$$\mu _2(T;S)=\underset{\stackrel{}{x}}{}\stackrel{}{x}^2s(\stackrel{}{0})s(\stackrel{}{x})_c$$
(6)
the “second-moment correlation length” $`\xi (T;S)`$ is defined by

$$\xi ^2(T;S)=\frac{\mu _2(T;S)}{4\chi (T;S)}.$$
(7)
The second field-derivative of the susceptibility $`\chi _4(T;S)`$ is defined by

$$\chi _4(T;S)=\underset{\stackrel{}{x},\stackrel{}{y},\stackrel{}{z}}{}s(\stackrel{}{0})s(\stackrel{}{x})s(\stackrel{}{y})s(\stackrel{}{z})_c.$$
(8)
## III  Definitions of critical parameters

In terms of the asymptotic behavior of these observables, we can now define the critical parameters, amplitudes and exponents that we are going to estimate using HT and LT series.

The spontaneous magnetization has the asymptotic behavior

$$M^{()}(T;S)B^{()}(S)|t(S)|^{\beta (S)}(1+a_M^{()}(S)|t(S)|^{\theta (S)}+..)$$
(9)
as $`t(S)0`$.

The asymptotic behaviour of the susceptibility as $`t(S)0\pm `$, is expected to be

$$\chi ^{(\pm )}(T;S)C^{(\pm )}(S)|t(S)|^{\gamma (S)}\left(1+a_\chi ^{(\pm )}(S)|t(S)|^{\theta (S)}+\mathrm{}+b_\chi ^{(\pm )}(S)t(S)+\mathrm{}\right)$$
(10)
The correlation length

$$\xi ^{(\pm )}(T;S)f^{(\pm )}(S)|t(S)|^{\nu (S)}\left(1+a_\xi ^{(\pm )}(S)|t(S)|^{\theta (S)}+\mathrm{}+b_\xi ^{(\pm )}(S)t(S)+\mathrm{}\right)$$
(11)
the specific heat

$$C_H^{(\pm )}(T;S)/k_BA^{(\pm )}(S)\mathrm{ln}|t(S)|\left(1+a_C^{(\pm )}(S)|t(S)|^{\theta (S)}+\mathrm{}+b_C^{(\pm )}(S)t(S)+\mathrm{}\right)$$
(12)
and the second field-derivative of the susceptibility $`\chi _4(K;S)`$

$$\chi _4^{(\pm )}(T;S)C_4^{(\pm )}(S)|t(S)|^{\gamma _4(S)}\left(1+a_4^{(\pm )}(S)|t(S)|^{\theta (S)}+\mathrm{}+b_4^{(\pm )}(S)t(S)+\mathrm{}\right)$$
(13)
have analogous asymptotic behaviours.

Different (universal) critical exponents $`\beta (S),\gamma (S),\nu (S),\gamma _4(S)`$ and different (non-universal) critical amplitudes $`B^{()}(S)`$, $`C^{(\pm )}(S)`$, $`f^{(\pm )}(S)\mathrm{}`$, $`a_\chi ^{(\pm )}(S)`$, $`a_\xi ^{(\pm )}(S)`$, etc. are associated with the various observables. We have reported in such detail our definitions of the critical amplitudes, because they differ significantly from those of other authors and it is necessary to use a consistent normalisation convention when comparing models expected to belong to the same universality class. Let us notice in particular that the estimates reported in the tables of Ref.ejg  for the critical amplitudes of the susceptibility $`\chi ^{()}(u;S)`$ are related to ours by the factor $`S^2`$ (-ln$`u_c(S))^\gamma /u_c(S)^{4S}`$. A similar remark applies to the specific heat amplitudes $`C^{()}(u;S)`$ for which the conversion factor is $`1/u_c(S)^{4S}`$(ln$`u_c(S))^2`$. Finally, the magnetization amplitudes of Ref.ejg  should be multiplied by the factor $`(`$-ln$`u_c(S))^{1/8}/S`$ to agree with ours. Of course, the amplitudes of the conformal field theory considered in the study of Ref.delf  are not comparable to our series quantities.

As indicated in eqs.(9) - (13), for a given spin $`S`$, all asymptotic forms are moreover expectedweg  to contain leading non-analytic confluent corrections characterized by the same exponent $`\theta (S)`$. Higher order corrections are also expected to contain logarithmicweg  factors. If universality holds, all exponents have to be $`S`$-independent.

The presence and the value of the confluent exponent has been discussedahafi ; case ; blonij ; barfish ; adler  several times. From RG calculationszinn , both in the $`ϵ`$-expansion and in the fixed-dimension approach, it was conjectured that $`\theta 4/3`$ for the universality class of the two-dimensional Ising model. Aharony and Fisher and later Blöte and den Nijs arguedahafi ; blonij  that $`a_\chi ^{(\pm )}(S)=0`$ for $`S=1/2`$ and indeed no such correction was revealed by the later very accurate studywu ; orri  of the critical asymptotic expansion for $`\chi (K;1/2)`$. However, in the absence of more general results, the reliable assessment of the subleading asymptotic critical behaviour remained an open problem when $`S>1/2`$.

## IV  Estimates of universal amplitude combinations

In terms of $`\chi ^{(+)}(K;S)`$, $`\xi ^{(+)}(K;S)`$ and $`\chi _4^{(+)}(K;S)`$, a “hyper-universal” combination of critical amplitudes denoted by $`g_r^{(+)}(S)`$ and usually called the “dimensionless renormalized coupling constant”, can be defined by

$$g_r^{(+)}(S)=\frac{3vC_4^{(+)}(S)}{8\pi (af^{(+)}(S))^2C^{(+)}(S)^2}$$
(14)
Here the normalisation factor $`\frac{3}{8\pi }`$ is chosen in order to match the usual field theoretic definitionzinn  of $`g_r^{(+)}(S)`$ and $`v`$ denotes the volume per lattice site, measured in units of the square of a lattice constant. For all lattices one has $`v=\sigma a^2`$, with $`a`$ the lattice constant. For the triangular lattice, $`\sigma =\sqrt{3}/2,`$ for the honeycomb lattice $`\sigma =3\sqrt{3}/4,`$ and for the kagomé lattice $`\sigma =2/\sqrt{3}.`$

We have also studied the hyper-universal combination usually denoted as

$$R_\xi ^{(+)}(S)=(A^{(+)}(S)/v)^{1/2}(af^{(+)}(S))$$
(15)
and the Watson combinationwats

$$R_C(S)=A^{(+)}(S)C^{(+)}(S)/B^{()}(S)^2$$
(16)
The other frequently considered universal combination

$$R_4(S)=C_4^{(+)}(S)B^{()}(S)^2/C^{(+)}(S)^3,$$
(17)
is not independent of the previous ones, since $`R_4(S)=\frac{8}{3\pi }g_r^{(+)}(S)R_\xi ^{(+)}(S)^2/R_C(S)`$.

All of these quantities are accurately known in the $`S=1/2`$ case. As indicated in Ref.wu ; aha , it is known that $`A^{(+)}(1/2)=\frac{2}{\pi }`$ln$`[\mathrm{tan}(\pi /8)]^20.4945385895`$, $`C^{(+)}(1/2)0.962581732`$, and $`B^{()}(1/2)=2^{5/16}`$ln$`[1+\sqrt{2}]^{1/8}1.22240995`$. In Refs.saso ; case , we find the very accurate estimates $`C_4^{(+)}(1/2)=4.379095(8)`$ and $`f^{(+)}(1/2)0.5670683`$.

Therefore we can conclude $`g_r^{(+)}(1/2)=1.754364(2)`$, $`R_C(1/2)0.31856939`$, $`R_\xi ^{(+)}(1/2)0.39878194`$ and $`R_4(1/2)=7.336744(10)`$.

We have also considered the ratio $`C^{(+)}(S)/C^{()}(S)`$. In Ref.wu , for $`S=1/2`$, this ratio was computed with arbitrary precision to be $`C^{(+)}(1/2)/C^{()}(1/2)37.693652`$.

Finally, we have estimated the ratio $`A^{(+)}(S)/A^{()}(S)`$ for various values of $`S.`$ This ratio equals unity for $`S=1/2`$ by self-duality. This was arguedande  in greater generality for the $`q`$-state ($`0q4`$) Potts model on the square lattice, which, for $`q=2`$, reduces to the $`S=1/2`$ Ising model.

In what follows, we determine the values of these universal amplitude combinations for $`S>1/2.`$ The preliminary part of our series analysis is aimed at estimating the critical temperatures using the expansions of $`\chi ^{(+)}(K;S)`$ for $`S1`$. We employed a variety of methods: Zinn-Justin improved-ratio formulazinn81 , Padé approximants (PA) and inhomogeneous differential approximants (DA)gutda . The best results with DA’s were obtained from approximants such that the polynomial coefficient of the highest derivative is even. (As a consequence, the approximants always contain an additional anti-ferromagnetic singularity at $`K_c(S)`$, beside the one at $`K_c(S)`$). Similarly to the LT analysis, but to a much smaller extent, the accuracy of our results tends to deteriorate with increasing $`S`$. In spite of this, our final HT estimates of the critical points, reported in Table 3, show significant improvement in apparent accuracy and sizable discrepancies from the previous LT determinationsejg .

For general values of $`S`$, less accurate estimates of $`K_c(S)`$ have been obtained in Ref.burk  from the ten term susceptibility series of Ref.camp  and from a renormalization group method. More recently other estimatesmonr  were obtained by a generalized cluster method. To our knowledge other accurate determinations of the critical points are available only for $`S=1`$. They have been obtained either by analysingadler  the 21 term HT series of Ref.nr90  for the susceptibility or by various transfer-matrix methodsbln ; yuri ; lipo . Some of these results have also been cited in Table 3.

We have then turned to the critical exponents $`\gamma (S)`$, $`\nu (S)`$ and $`\gamma _4(S)`$ and have evaluated them from the log-derivatives of the appropriate HT expansions by first order DA’s biased with our HT estimates of the critical temperatures. This computation shows that the relative variation of the exponents is smaller than $`10^3`$, in the worst case, for $`S`$ varying between $`1/2`$ and $`\mathrm{}`$. We report these results in Table 4 without further details and simply conclude that universality and hyperscaling appear to be well supported for the leading critical exponents.

It is perhaps also worth noticing that assuming the universality of $`\gamma `$ we can bias and therefore refine the determination of $`K_c(S)`$. This procedure does not change the central values of the critical points with respect to the unbiased one, but reduces the error bars.

On the other hand, the estimate of the exponent $`\theta (S)`$ of the leading singular confluent corrections to scaling in the various observables remains quite elusive. Performing either a Baker-Hunterhb  or a Zinn-Justinzinn81  analysis, we can conclude that, at the level of accuracy made possible by the present extension of the HT series, the amplitudes of these corrections are very small, (or perhaps vanishing) for all values of $`S`$. We should mention that a similar conclusion was suggested for $`S=1`$ in Ref.blonij , while the opposite conclusion was advocated in Ref.adler .

Once we have estimated the critical temperatures and verified the universality of the leading exponents, we can proceed with the analysis simply assuming that, for all values of $`S`$, these exponents take exactly the values expected for the universality class of the two-dimensional spin-1/2 Ising model and using them along with our estimated critical temperatures to bias the evaluation of the HT and LT critical amplitudes defined by eqs.(9) - (13).

Our estimates of the critical amplitudes are reported in Table 5.

We have employed quasi-diagonal non-defective PA’s or DA’s for extrapolating to $`K_c(S)`$ the effective amplitudes of the susceptibility and of the derivative of the specific heat, from the HT and the LT side of the critical points. We have similarly studied the effective amplitudes of the correlation-length (available only in the HT region) and of the magnetization. For proper comparison, in the same Table we have also cited the LT estimates of the critical amplitudes for the spontaneous magnetization, the specific heat and the susceptibility previously obtained in Ref.ejg . These quantities have been multiplied by the above indicated conversion factors to agree with our normalisation conventions. The uncertainties we have reported, which allow for the observed spreads in the approximant values, provide a subjective assessment of residual trends in the sequence of estimates and for the (unbiased) uncertainties of the critical points. The HT amplitudes can be determined with a relative accuracy ranging from $`10^3`$ in the case of the susceptibility, to $`10^2`$ in the case of the specific heat. The LT amplitudes are subject to larger relative uncertainties, increasing with $`S`$, and reaching up to $`50\%`$ for $`S>2`$. In some cases, in order to improve the accuracy of the estimates of the LT amplitudes for $`S2`$, we have based our extrapolations only on the data for $`|t(S)|0.020.04`$. This unconventional but reasonable procedure reduces the sensitivity of the approximants to the unphysical nearby singularities. Unfortunately, even this prescription fails to work satisfactorily for $`S>2`$.

Using only the HT series, we can evaluate $`g_r^{(+)}(S)`$, either directly, in terms of the amplitudes reported in Table 5, or by extrapolating to the critical points via DA’s the HT expansion of the inverse effective coupling $`1/g_r^{(+)}(K;S)=\frac{8\pi }{3}\xi (K;S)^2\chi (K;S)^2/\chi _4(K;S)`$. A third approach consists in studying the residua at $`x=1`$ of the series with coefficients $`a_n(S)=c_n(S)/d_n(S)`$, where $`c_n(S)`$ are the HT coefficients of $`\xi ^2(K;S)`$ and $`d_n(S)`$ are the coefficients of the quantity $`\chi _4(K;S)/\chi ^2(K;S)`$. In Table 5 we have reported the results of the latter procedure since it yields estimates with smaller spreads.

Several other estimatesbb ; lai ; bc96 ; kim ; case ; balo  obtained by a variety of methods are also available in the literature.

Using also the LT series, we have evaluated, directly in terms of the single amplitudes, the other mentioned universal combinations, for a range of values of $`S`$. We have reported in Table 6, our series estimates of all these quantities for $`S>1/2`$. In conclusion, whenever only HT amplitudes are involved, our estimates, within a precision up to $`0.1\%`$, are independent of $`S`$, in full agreement with universality. On the other hand, our reanalysis of the LT series has been only partially successful: whenever LT amplitudes also enter into the combinations, universality appears to be fairly well respected for $`S<5/2`$, but the uncertainties grow notably larger for larger values of the spin.

In the next section we describe the theory of lattice-lattice scaling, and show how it can be used to extend our estimates of the critical amplitudes from the square lattice to other two-dimensional lattices.

## V Lattice-lattice scaling

The theory of lattice-lattice scaling was developed by Betts, Guttmann and Joycebgj  in the early ’70s. It explains how amplitudes change within a given universality class, as one moves from one lattice to another. It can also be viewed as a generalisation of the law of corresponding states. In this section we give a terse development of the theory, and apply it to the problem at hand.

In order to review the general ideas let us first consider the Weiss theory or mean field theory of a magnetic system. The equation of state is well known to be

$$h=\frac{1}{3}m^3(1+3t/m^2).$$

Here $`t=T/T_c1,`$ $`h=\mu H/kT,`$ and $`m=M(T)/M(0)`$ are the reduced temperature, magnetic field and magnetization, respectively.

Then the law of corresponding states says that the equation of state is the same for all lattices. That is,

$$m_\mathrm{X}(t,h)=m_\mathrm{Y}(t,h),$$

where $`\mathrm{X}`$ and $`\mathrm{Y}`$ denote two lattices. That is to say, the lattice dependence is entirely contained in the critical temperature $`T_c.`$

A more complex model is the three-dimensional spherical model, for which the critical equation of state is:

$$h=D_\mathrm{X}m^5(1+t/m^2)^2.$$

Here both $`T_c`$ and the amplitude $`D`$ are lattice dependent. Thus

$$m_\mathrm{X}(t,h_\mathrm{X})=m_\mathrm{Y}(t,h_\mathrm{Y}).$$

We see that we must scale the field variable, so that $`\frac{h_\mathrm{X}}{D_\mathrm{X}}=\frac{h_\mathrm{Y}}{D_\mathrm{Y}},`$ but that there is no need to scale the reduced temperature.

Let us now consider the case of the (zero-field) spin $`S=1/2`$ Ising model on the triangular ($`\mathrm{T}`$) and hexagonal ($`\mathrm{H}`$) lattices. The star-triangle relationons ; fis  allows us to relate the free-energy, susceptibility and spontaneous magnetization between the lattices:

$`2f_\mathrm{H}(K_\mathrm{H})=f_\mathrm{T}(K_\mathrm{T}),`$
$`M_{\mathrm{T},0}(t_\mathrm{T})=M_{\mathrm{H},0}(t_\mathrm{H}),`$
$`2\chi _\mathrm{T}(v_\mathrm{T})=\chi _\mathrm{H}(v_\mathrm{H})+\chi _\mathrm{H}(v_\mathrm{H}),`$
where $`f=1/N\mathrm{ln}Z,`$ $`K_\mathrm{X}=J_\mathrm{X}/kT,`$ and $`v_\mathrm{X}=\mathrm{tanh}(J_\mathrm{X}/kT).`$ Here we see that the reduced temperature needs to be re-scaled for the free-energy to be universal. This is not restricted to the triangular-honeycomb pair, but in that case it is easy to be explicit.

All these examples can be encapsulated in the following expression for the singular part of the free-energy:

$$n_\mathrm{X}f_\mathrm{X}(t_\mathrm{X},h_\mathrm{X})=n_\mathrm{Y}f_\mathrm{Y}(t_\mathrm{Y},h_\mathrm{Y})=f(t,h),$$

where the reduced temperature and field are scaled by

$$n_\mathrm{X}h_\mathrm{X}=n_\mathrm{Y}h_\mathrm{Y}=h$$

and

$$g_\mathrm{X}t_\mathrm{X}=g_\mathrm{Y}t_\mathrm{Y}=t.$$

The singular part of the free-energy, $`f(t,h),`$ is then a universal (lattice independent) function for a given model.

Equivalently, by differentiation we obtain

$$m_\mathrm{X}(t_\mathrm{X},h_\mathrm{X})=m_\mathrm{Y}(t_\mathrm{Y},h_\mathrm{Y})=m(t,h),$$

and

$$\chi _\mathrm{X}(t_\mathrm{X},h_\mathrm{X})/n_\mathrm{X}=\chi _\mathrm{Y}(t_\mathrm{Y},h_\mathrm{Y})/n_\mathrm{Y}=\chi (t,h),$$

where $`m(t,h)`$ and $`\chi (t,h)`$ are universal functions for a given model.

Writing $`m_\mathrm{X}=B_\mathrm{X}^{()}(t_\mathrm{X})^\beta ,`$ it follows that

$$\frac{g_\mathrm{X}}{g_\mathrm{Y}}=\left(\frac{B_\mathrm{X}^{()}}{B_\mathrm{Y}^{()}}\right)^{1/\beta }.$$

Using this result, and the exact scaling parameters $`g_\mathrm{X}`$ and $`g_\mathrm{Y}`$ given below, it is a trivial matter to calculate the magnetization amplitudes for the other lattices we consider (triangular, hexagonal and kagomé ($`\mathrm{K}`$)), taking as input the square lattice amplitudes given in Table V. These amplitude estimates are given in Tables VII, VIII and IX.

Writing $`\chi _\mathrm{X}=C_\mathrm{X}^{(\pm )}t_\mathrm{X}^\gamma `$ it follows that

$$C_\mathrm{X}^{(+)}/C_\mathrm{Y}^{(+)}=\frac{n_\mathrm{X}}{n_\mathrm{Y}}\left(\frac{g_\mathrm{X}}{g_\mathrm{Y}}\right)^\gamma =C_\mathrm{X}^{()}/C_\mathrm{Y}^{()}.$$

Similarly, it is a trivial matter to calculate the susceptibility amplitudes for the other lattices we consider, taking as input the $`S=1/2`$ square lattice amplitudes given in Table V.

Further differentiation gives the corresponding relationship for higher field derivatives, and we readily obtain

$$C_{2l,\mathrm{X}}^{(+)}/C_{2l,\mathrm{Y}}^{(+)}=\left(\frac{n_\mathrm{X}}{n_\mathrm{Y}}\right)^{2l1}\left(\frac{g_\mathrm{X}}{g_\mathrm{Y}}\right)^{\gamma _{2l}},$$

for the high-temperature field derivatives, (where only the even-order derivatives are non-zero). The corresponding result for low-temperature field derivatives is

$$C_{l,\mathrm{X}}^{()}/C_{l,\mathrm{Y}}^{()}=\left(\frac{n_\mathrm{X}}{n_\mathrm{Y}}\right)^{l1}\left(\frac{g_\mathrm{X}}{g_\mathrm{Y}}\right)^{\gamma _l^{}}.$$

Taking temperature derivatives, one readily establishes that the specific heat amplitudes satisfy

$$A_\mathrm{X}^{(+)}/A_\mathrm{Y}^{(+)}=\frac{n_\mathrm{Y}}{n_\mathrm{X}}\left(\frac{g_\mathrm{X}}{g_\mathrm{Y}}\right)^{2\alpha }=A_\mathrm{X}^{()}/A_\mathrm{Y}^{()}.$$

As $`\alpha =0`$ for the two-dimensional Ising model, this simplifies to

$$A_\mathrm{X}^{(+)}/A_\mathrm{Y}^{(+)}=\frac{n_\mathrm{Y}}{n_\mathrm{X}}\left(\frac{g_\mathrm{X}}{g_\mathrm{Y}}\right)^2=A_\mathrm{X}^{()}/A_\mathrm{Y}^{()}.$$

We have similarly calculated the specific-heat amplitudes for the other lattices we consider, taking as input the square lattice amplitudes given in Table VI. These are also given in Tables VII, VIII, and IX.

Finally, the correlation function amplitudes were calculated exactly from various star-triangle transformations by Thompson and Guttmann TG75  for the true correlation length. It follows from universality that these same transformations should hold also for the second-moment correlation length amplitudes. We show below that this is equivalent to the universality of $`R_\xi ^{(+)}(S),`$ which is a conclusion of this work.

In TG75  it was shown by explicit calculation that

$$f_\mathrm{T}^{(+)}(1/2)K_{c,\mathrm{T}}^{(+)}(1/2)/f_\mathrm{H}^{(+)}(1/2)K_{c,\mathrm{H}}^{(+)}(1/2)=1/3.$$

$$f_\mathrm{S}^{(+)}(1/2)K_{c,\mathrm{S}}^{(+)}(1/2)/f_\mathrm{T}^{(+)}(1/2)K_{c,\mathrm{T}}^{(+)}(1/2)=\sqrt{3}.$$

$$f_\mathrm{K}^{(+)}(1/2)K_{c,\mathrm{K}}^{(+)}(1/2)/f_\mathrm{H}^{(+)}(1/2)K_{c,\mathrm{H}}^{(+)}(1/2)=2/3.$$

The universality of $`R_\xi ^{(+)}(S)=\left(\frac{A^+(S)}{v}\right)^{1/2}af^+(S),`$ taken together with the above lattice-lattice scaling relation for the specific heat amplitude $`A^+(S)`$ implies the lattice-lattice scaling relation

$$f_\mathrm{X}^{(+)}/f_\mathrm{Y}^{(+)}=\frac{g_\mathrm{Y}}{g_\mathrm{X}}\left(\frac{\sigma _\mathrm{X}n_\mathrm{X}}{\sigma _\mathrm{Y}n_\mathrm{Y}}\right)^{1/2},$$

where $`\sigma `$ is the area per site, and $`\sigma _\mathrm{X}=1,\sqrt{3}/2,3\sqrt{3}/4`$ and $`2/\sqrt{3}`$ for the square, triangular, honeycomb and kagomé lattices respectively. This is equivalent to the explicit amplitude scaling reported in TG75 , and confirms the expectation that the true correlation length amplitude and the second-moment correlation length amplitude scale similarly.

For the 2d Ising model we can calculate the scaling parameters $`n_\mathrm{X}`$ and $`g_\mathrm{X}`$ exactly from known spin-1/2 spontaneous magnetization and specific heat amplitudes. The critical points are also exactly known. These are given in the table below:

| $`\mathrm{X}`$ | $`K_c`$ | $`n_\mathrm{X}`$ | $`g_\mathrm{X}K_{c,\mathrm{T}}/K_{c,\mathrm{X}}`$ |
| --- | --- | --- | --- |
| Tr. | $`\mathrm{ln}(3)/4`$ | 1 | 1 |
| Sq. | $`\mathrm{ln}(\sqrt{2}1)/2`$ | $`3\sqrt{3}/4`$ | $`1/\sqrt{2}`$ |
| Ho. | $`\mathrm{ln}(2\sqrt{3})/2`$ | 2 | $`1/\sqrt{3}`$ |
| Ka. | $`\mathrm{ln}(2/\sqrt{3}1)/4`$ | $`3^5(\sqrt{3}1)^{16}`$ | $`9(\sqrt{3}1)^8`$ |

In addition to the above results, one can also derive scaling relations for the amplitudes of sub-dominant singularities. For example, consider the susceptibility of the 2d Ising model on lattice X. Writing

$$\chi _\mathrm{X}C_0^\mathrm{X}t^\gamma +C_1^\mathrm{X}t^{1\gamma },$$

from lattice-lattice scaling we can derive the following amplitude relations:

$$\frac{C_0^\mathrm{X}}{C_0^\mathrm{T}}=\frac{n_\mathrm{X}}{n_\mathrm{Y}}\left(\frac{g_\mathrm{Y}}{g_\mathrm{X}}\right)^\gamma ,$$

$$\frac{C_1^\mathrm{X}}{C_1^\mathrm{T}}=\frac{n_\mathrm{X}}{n_\mathrm{Y}}\left(\frac{g_\mathrm{Y}}{g_\mathrm{X}}\right)^{\gamma 1}.$$

The second expression is false for the kagomé lattice. It is corrected by the theory of extended lattice-lattice scaling developed by Gaunt and Guttmanngg  in 1978. In that theory, a third scaling parameter needs to be introduced.

It is widely accepted, and in complete agreement with the results of the first part of this paper, that the spin-$`S`$ Ising model is in the same universality class as the spin-$`1/2`$ Ising model. This is the only assumption we require in order to apply the theory of lattice-lattice scaling to the square lattice data given in the previous section. The only subtlety is whether “universality” really extends to lattice-lattice universality. While this assumption seems natural, we did attempt to verify it by estimating amplitudes for other lattices from the rather short data available in camp . The longest effective series is the triangular lattice series. We found the high-temperature susceptibility amplitude, as estimated by Padé approximants, to be $`C^{(+)}(1)=0.529`$ from this series, in complete agreement with the more precise value $`C^{(+)}(1)=0.5294(2)`$ found by lattice-lattice scaling, and reported in Table VII.

Note too that there is no loss of accuracy, as the conversions from lattice to lattice are exact. For example, based on the recent estimate of the leading susceptibility amplitude of the spin $`1/2`$ square lattice Ising model given inorri , application of lattice-lattice scaling gives the corresponding amplitude on the triangular lattice to the same precision, viz:

$$C_+^T=0.9242069582451643296971575778559317176696261520028389.$$

Similarly accurate results for other lattices can be readily written down, as can equally accurate sub-dominant amplitudes.

###### Acknowledgements.

This work was completed before the untimely death of Marco Comi, our lifelong friend and coauthor. To his dear memory we dedicate this paper. This work has been partially supported by the Ministry of Education, University and Research (PB, MC), and the Australian Research Council (AJG).


format: text

text:
# Dynamics of Scale Free Random Threshold Network

## Abstract

We study the dynamics of Random Threshold Network (RTN) on scale free networks, with asymmetric links, some interaction rules where propagation of local perturbations depends on in-degree $`k`$ of the nodes. We find that there is no phase transition with respect to average connectivty independently of network topology for the case temperature $`T=0`$, threshold $`h=0`$ and the probability distribution of indegree $`P(k)`$ satisfies $`P(0)=D=0`$. We have investigated the emergence of phase transition involving three parameters, i.e. $`T,h`$ and $`D`$. RTN can be continuously connected to Random Boolean Network (RBN) in $`T\mathrm{}`$, and we find moderate thermal noise extends the regime of ordered dynamics, compared to RTN in $`T=0`$ regime and RBN. Furthermore, we discuss the dynamic properties from another point of view, dynamical mean field reaction rate equation.

Recently there has been attracted considerable interest in the formation of random networks, in the connectivity and especially in the dynamics of these networks WaSt ; AlBa ; BaAlRev . Many naturally occuring networks such as the Internet AlBar , WWW CaMa , gene regulatory networkAg , protein networks Je ; MaSn  exhibit a power-law; i.e. scale free degree distribution for the links. A scale free network have a distinctive property that a small fraction of the nodes which is often called ’hub’ are highly connected whereas the majority of the nodes have low connectivities. The absence of a typical scale of the connectivity is often related to the organization of the network as a hierarchy. Although their topological properties have been studied in detail, dynamical behavior provided with asymmetric links and some interaction rules are not sufficiently explored. Especially in neural networks and gene regulatory networks, interaction rules between nodes plays an important role to organize the dynamical properties of the network. Such a dynamical network model, e.g. Random Boolean Network (RBN) is first introduced by Kauffman Ka ; Ka2  to model the gene regulatory netoworks in the biological system. In this models the state $`\sigma _i\{0,1\}`$ of a network node $`i`$ is a logical function of the states of in-degree $`k_i`$ other nodes chosen at random. The logical functions are chosen at random, with a suitably biased probability. RBN model shows a phase transition with respect to the number of inputs per node $`K`$ at a critical average connectivity $`K_c=[2p(1p)]^1`$ from ordered to chaotic state, where $`p`$ is the probability that the randomly chosen output of node $`i`$ are unity. In chaotic regime ($`K>K_c`$), a small perturbation in the initial state propagates across the entire system, whereas all perturbation in the initial state dies out in ordered regime ($`K<K_c`$). The dynamics of RBN on scale free is recently investigated Aldana . However, its dynamical property is the same as that of classical RBN model because the output state of the nodes in RBN is independent of their in-degree distribution. In order to study the dynamical effects of the in-degree (scale free) distribution, we focus on Random Threshold Network (RTN), first investigated as diluted and asymmetric spinglass models Derrida , and asymmetric neural networks DeGa ; KrZi . It is a subset of Boolean Network which is thought to show dynamical behavior similar to RBN BoRo .

We consider a network of $`N`$ randomly interconnected binary nodes with state $`\sigma _i=\pm 1`$. In this model, probability distribution of outputs strongly depends on their in-degree parameter. At time $`t`$, the fields $`f_i(t)`$ of node $`i`$ are computed:

$`f_i(t)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}c_{ij}\sigma _j(t)+h`$ (1)
where $`h(0)`$ is the threshold parameter. Threshold characteristics can be often shown in neural networks and gene regulatory networks, i.e. threshold for neuron firing and that for gene expression. In this brief report, The interaction weight $`c_{ij}`$ take discrete values $`c_{ij}=\pm 1`$, with equal posiibility. If node $`i`$ does not take a signal from $`j`$, one has $`c_{ij}=0`$. We assume $`c_{ii}`$ is always zero to avoid $`self`$ $`connection`$. The in-degree $`k_i`$ of node $`i`$ is defined as $`k_i=_j^N|c_{ij}|`$. The interaction weights $`c_{ij}`$ is chosen randomly from all nodes but follows the in-degree distribution $`P(k)`$. For each node $`i`$, its state at time $`t+1`$ are functions of the inputs it receives from other nodes at time $`t`$.

$`\sigma _i(t+1)`$ $`=`$ $`\{\begin{array}{cc}\hfill 1,& \text{with probability }G(f_i(t))\hfill \\ \hfill 1,& \text{with probability }G(f_i(t))\hfill \end{array}`$ (4)
where $`G(x)=[1+\text{exp}(2x/T)]^1`$. The connection between RBN and RTN is continuously made by introducing thermal noise to probability distribution of outputs function $`G(x)`$. The parameter $`T`$ defines the temperature of the system. The $`N`$ network nodes are idealized by synchronized updating. Generalizing so-called Annealing Approximation method introduced by Derrida and Pomeau DePo , we approach the system analytically in the limit of a large system size $`N`$. One finds that normalized overlap function $`x`$ obey the equation

$`x(t+1)`$ $`=`$ $`F(x(t))`$ (5)
where the mapping function $`F(x)`$ is given by

$`F(x)`$ $``$ $`1{\displaystyle \underset{k}{\overset{\stackrel{~}{K}}{}}}p_s(k,h,T)P(k)(1x^k)`$ (6)
Normalized overlap function $`x(t)=1d(\stackrel{}{\sigma (t)},\stackrel{}{\sigma ^{}(t)})/N`$ is a subtraction of the normalized Hamming distance from unity where Hamming distance $`d(\stackrel{}{\sigma (t)},\stackrel{}{\sigma ^{}(t)})`$ denotes overlap between time independent two distinctive configurations. In the $`N\mathrm{}`$ limit, $`x`$ represent the probability of two arbitary configurations to be equal. The possibility $`p_s`$ is output state reversal by changing the state of a single input $`j`$. Using combinatorial methods RoBo , the stochastic distribution of $`p_s`$ is derived.

$`p_s(k,h,T)=`$ $`{\displaystyle \frac{1}{k2^k}}{\displaystyle \underset{i=0}{\overset{k}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{i}}\right)(G(k+h2i)(k(ki)G(k+h22i)iG(k+h+22i))`$ (7)
$`+(1G(k+h2i))((ki)G(k+h22i)+iG(k+h+22i)))`$
Stationary value of overlap $`x^{}`$ is obtained by means of fixed point equation $`x=F(x)`$. The existence of a phase transition depends uniquely on the nature of the fixed point $`x=1`$. If it is attractive, two initial configurations differing by an infinitesimal fraction of nodes will become almost identical, whereas if it is repulsive, they will produce diverging trajectories. A possible critical point is then determined by the equation

$`{\displaystyle \frac{dF}{dx}}|_{x=1}=1`$ (8)
Let us recall briefly sufficient conditions for attractive fixed point $`x^{}`$ , that is written as, Eq.(6) satisfies (a)$`dF/dx0`$ for all $`0x1,`$ (b)$`F(0)>0`$ and (c)$`dF/dx|_{x=1}>1`$. For $`T=0`$, $`p_s`$ is given by,

$`p_s(k,h)=\{\begin{array}{cc}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{(k+h)/2}}\right)/2^k,\hfill & \text{for }k+h\text{ is even}\hfill \\ p_s(k1,h),\hfill & \text{for }k+h\text{ is odd}\hfill \end{array}`$ (11)
Moreover, $`h`$ is also set to zero. Item (a), (b) are manifest, and (c) is surprisingly satisfied for any $`P(k)`$ if $`P(0)=0`$. It is straight-forward to prove (c), that we show $`kp_s(k)10,\text{for }k1`$, using asymptotically $`p_s(k)k^{1/2}`$ for large $`k`$. In this case, fixed point $`x=1`$ is always repulsive and another fixed point $`x^{}<1`$ is stable. That is, independent of the in-degree distribution $`P(k)`$, there is no phase transition and a small perturbation/damage could propagate over the entire system.

We now introduce the scale free in-degree distribution $`k_i`$ given by,

$`P(k)`$ $`=`$ $`\{\begin{array}{cc}D,\hfill & k=0\hfill \\ (1D)\eta (\gamma ,\stackrel{~}{K})^1k^\gamma ,\hfill & 1k\stackrel{~}{K}\hfill \end{array}`$ (14)
where $`\eta (\gamma ,\stackrel{~}{K})=_1^{\stackrel{~}{K}}k^\gamma `$ with cutoff $`\stackrel{~}{K}`$, $`0D1`$ and $`\gamma `$ is usually called scale free exponent. Substituting Eq.(14) to Eq.(8), we obtain the critical condition for $`D`$.

$`D_c={\displaystyle \frac{A(\gamma )1}{A(\gamma )}},A(\gamma )=\eta (\gamma ,\stackrel{~}{K})^1{\displaystyle \underset{k=1}{\overset{\stackrel{~}{K}}{}}}p_s(k)k^{1\gamma }`$ (15)
Since $`1p_s(k,0)0`$ for $`k1`$, critical average connectivity $`K_c=\eta (\gamma 1,\stackrel{~}{K})/(A(\gamma )\eta (\gamma ,\stackrel{~}{K}))`$ has a minimum value in the $`\gamma \mathrm{}`$, i.e. minimum critical average $`K_c=1`$ is realized when $`D_c=0`$.

In Fig.1, we plot the stationary overlap function $`x^{}`$ for various threshold parameters, provided that $`T=D=0`$. The overlap function for the case $`h=0`$ is shown to asymptotically approach $`x^{}1`$ with increasing $`\gamma `$, i.e. average connectivity is approaching $`K=1`$. It is obviously different from RBN in which exists well known chaotic to ordered phase transition. For finite thresholds $`h`$, we observed equilibrium overlap is decreased as threshold incereases, which suggests they eliminates the memory of initial condition and two different initial condition become identical. For the comparison, The curve of equilibrium overlap for RBN is denoted. One observes a pertuabation is rather suppressed in RTN with zero-threshold than in RBN for low $`\gamma `$ regime. The opposite is observed for high $`\gamma `$ regime. The critical point where the strength of perturbation reverse is given by $`\gamma _s1.754`$ in this model. This result supports the numerical simulation observed in Fig.8 of RoBo .

Fig.2 shows the curve of critical scale free exponent $`\gamma _c`$ as a function of temperature for $`h=D=0`$, which distinguish ordered state from chaotic state. For finite temperature $`T`$, $`dF/dx|_{x=1}`$ in large $`\gamma `$ limit is given by,

$`{\displaystyle \frac{dF}{dx}}|_{x=1_{}}=`$ $`{\displaystyle \frac{1}{(1+\mathrm{exp}(2/T))^2}}`$
$`+`$ $`{\displaystyle \frac{1}{(1+\mathrm{exp}(2/T))^2}}+𝒪(\text{1}/𝒩)`$ (16)
The existence of infinitesimal $`T`$ which satisfies $`dF/dx|_{x=1}<1`$ results in critical temperature $`T_c0`$. The critical line would diverge in the zero temperature limit. An interesting feature of the $`\gamma _c`$ curve is that it has a minimum value $`\gamma _c=2.1397`$ , i.e. maximum $`K_c=`$2.611, at $`T=1.97\pm .003`$. The same phenomenon can be observed in other connectivity models, for example, Poisson distribution model where the in-degree distribution follows $`P(k)=e^{K_p}K_p^k/k!`$. Parameter $`K_p`$ is the average connectivity of the network. It has the critical average connectivity $`K_c=1.8494`$ at $`T=0`$ RoBo . In the low $`T`$ regime, the higher the temperature $`T`$, the higher the critical average connectivity $`K_c`$, and has a maximum value $`K_c=2.1994`$ around $`T=1.683\pm .001`$ (subsequently goes to $`K_c=2`$). Both results show that moderate thermal noise can extend the regime of ordered dynamics in RTN to higher average network connectivities, compared to the zero temperature and the infinite temperature regime. The infinite temprature $`T\mathrm{}`$ asymptotically leads to $`p_s=1/2`$, which is identical to RBN with $`p=1/2`$. Substituting it to Eq.(8), the critical condition in RBN is given by,

$`{\displaystyle \frac{\eta (1+\gamma ,\stackrel{~}{K})}{2\eta (\gamma ,\stackrel{~}{K})}}=1`$ (17)
Since average connectivity of the network is given by $`K=\eta (1+\gamma ,\stackrel{~}{K})/\eta (\gamma ,\stackrel{~}{K})`$, one finds $`K_c=2`$ from Eq.(17). Notice this value is independent of $`\stackrel{~}{K}`$ and $`\gamma `$. Indeed, it is independent of any topological structures. In the $`\stackrel{~}{K}\mathrm{}`$ limit, Eq.(17) is represented by Riemann $`\zeta `$ function, which gives the critical value $`\gamma _c2.4788`$.

In the following, we approach the system from another point of view, by describing the single-site equation governing the time evolution of the node state. We show a simple example of agent’s strategy game. We first take the node state $`\sigma _i(t)`$ as agent’s strategy at time $`t`$, i.e. $`\sigma _i(t)=1`$ as strategy $`A`$, and $`\sigma _i(t)=1`$ as strategy $`B`$. One can define $`\rho _k(t)`$ as a density of agents which have connectivity $`k`$ and take strategy $`B`$ at time $`t`$. We assume the interaction weight $`c_{ij}`$ takes only $`1`$ and $`0`$ and threshold $`h`$ to zero to simplify the dynamics. This alteration does not change the dynamical overlap equation Eq.(7). This model can be also interpreted as a ’locally majority selection game’ which means each agent choose the strategy at time $`t+1`$ the same as the majority of their inputs at time $`t`$. If the number of inputs strategies $`A`$ and $`B`$ is equal, the agent choose either one with equal probability. The dynamical mean field reaction rate equation for $`\rho _k(t)`$ is written as MaDi ; SatVes

$`_t\rho _k(t)=\delta _k\rho _k(t)+\nu _k(1\rho _k(t))`$ (18)
Agents with connectivity $`k`$ and strategy $`A`$ would change its strategy to $`B`$ with rate $`\nu _k`$, whereas agents with strategy $`B`$ would change to $`A`$ with rate $`\delta _k`$. By using $`\rho =_k^{\stackrel{~}{K}}P(k)\rho _k`$ which gives the average density of agents with strategy $`B`$, these rate parameters are written as

$`\delta _k`$ $`=`$ $`{\displaystyle \underset{i=0}{\overset{k}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{i}}\right)\rho ^i(1\rho )^{ki}G(k2i)`$
$`\nu _k`$ $`=`$ $`1\delta _k`$ (19)
Imposing the stationary condition on Eq.(18), we find stationary densities $`\rho _k^{}=1\delta _k^{}`$ since $`\rho _k`$ is on its turn function of $`\delta _k`$. A self consistency equation that allows us to find $`\rho ^{}`$ is given by:

$`\rho =1{\displaystyle \underset{k=1}{\overset{\stackrel{~}{K}}{}}}{\displaystyle \underset{i=0}{\overset{k}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{i}}\right)\rho ^i(1\rho )^{ki}G(k2i)P(k)`$ (20)
From Eq.(20), An asymptotic $`\gamma T`$ phase diagram similar to Fig.2 is obtained. In $`\gamma <\gamma _c(T)`$ regime, the solution of Eq.(20) has 3 fixed points $`\rho =\rho _1(<1/2),1/2,\rho _2(>1/2)`$. However, the fixed point at $`\rho =1/2`$ is an unstable solution because the derivative of right hand in Eq.(20) has more than unity. Whether the density $`\rho `$ goes to $`\rho _1`$ or $`\rho _2`$ in the stationary state depends on its initial condition. In $`\gamma >\gamma _c(T)`$ regime, the equation has only one fixed and stable point at $`\rho =1/2`$, which leads that any initial configurations asymptotically reach $`\rho =1/2`$ in the stationary state.

In Fig.3, we show the density $`\rho `$ in the upper plane ($`\rho >1/2`$) where the density $`\rho _2`$ reaches $`\rho =1/2`$ as $`\gamma `$ increases. On the other hand, in the $`T0`$ limit, stationary density $`\rho `$ is given by,

$`\underset{T0}{lim}\rho _1=0,\underset{T0}{lim}\rho _2=1,\text{for finite }\gamma `$ (21)
The time interval required to reach the equilibrium state $`t^{}`$ is estimated as $`t^{}e^\gamma `$. Since it takes long time to reach stationary state for large $`\gamma `$, the network would be practically fixed to its initial state for large $`\gamma `$ limit. In a small $`\gamma `$ area, the cutoff $`\stackrel{~}{K}`$ become crucial, which leads to $`lim_T\mathrm{}\rho _2(\gamma =0)=\stackrel{~}{K}/(\stackrel{~}{K}+1),lim_{T0}\rho _2(\gamma =0)=1`$. For finite $`h`$ case, $`\rho =0`$ is the only stable point for $`\gamma >\gamma _c`$, whereas there is another stable point for $`\gamma <\gamma _c`$.

In conclusion, we have studied the properties of Random Threshold Network model on scale free networks. In $`h=0,T=0,D=0`$ case, we found there is no phase transition with respect to average connectivity in any network topology. , i.e. propagation of perturbation/damage independent of $`\gamma `$. Regarding to $`D`$, it is quite interesting that minimum critical average connectivity is obtained when the rate of $`k=0`$ in-degree distribution is zero. We also clarifies the relationship between RTN and RBN, our findings that moderate thermal noise extends the regime of ordered dynamics to higher average connectivities, i.e. suppresses the propagation of perturbation/damage, suggest a new aspect to e.g. the theory of neural networks. Tuning thermal noise parameter in networks with asymmetric connections (as usually found in neural networks) where the basic units (neuron) follows threshold dynamics, thermal noise, i.e. non-zero error rate, may in fact improve information processing, given hierachical structure of the network topology.


text:
# Dynamical Nuclear Polarization by Electrical Spin Injection in Ferromagnet-Semiconductor Heterostructures

## Abstract

Electrical spin injection from Fe into AlxGa1-xAs quantum well heterostructures is demonstrated in small ($`<500`$ Oe) in-plane magnetic fields. The measurement is sensitive only to the component of the spin that precesses about the internal magnetic field in the semiconductor. This field is much larger than the applied field and depends strongly on the injection current density. Details of the observed hysteresis in the spin injection signal are reproduced in a model that incorporates the magnetocrystalline anisotropy of the epitaxial Fe film, spin relaxation in the semiconductor, and the dynamical polarization of nuclei by the injected spins.

The injection of spin from a conventional ferromagnetic metal into a semiconductor is a prerequisite for several proposed magneto-electronic devices Datta:Das:SpinTransistor . Although spin transport across the ferromagnet-semiconductor (FM-S) interface has recently been demonstrated Hammar:Johnson:1st-FM-SC\_spinInjection ; Ploog:Zhu:InGaAsPRL ; Jonker:Hanbicki:AlGaAsAPL ; IMEC:Motsnyi:OHE , most injection experiments on metallic FM-S systems have required relatively large magnetic fields, in excess of several kilogauss, to produce a spin component perpendicular to the FM-S interface. The most useful properties of typical ferromagnetic thin films, however, such as low-field switching and hysteresis, can be exploited only by coupling to the in-plane component of the magnetization OhnoY:Awschalom:FirstGaMaAs\_spinLED . In the case of metallic FM-S structures, in-plane coupling has been observed only as a small change in transport properties Hammar:Johnson:1st-FM-SC\_spinInjection  or using optically pumped carriers Kawakami:Awschalom:FM\_imprinting\_in\_GaAs ; Epstein:Awschalom:FPP\_in\_n-GaAs .

In this Letter we report a demonstration of electrical spin injection in FM-S heterostructures using small ($`<500`$ Oe) in-plane magnetic fields. We measure only the component of the spin that precesses after injection into the semiconductor using electroluminescence polarization (ELP) as a detection technique Fiederling:Molenkamp:FirstBeMnZnSe\_spinLED ; OhnoY:Awschalom:FirstGaMaAs\_spinLED . The effective magnetic field inducing the precession depends strongly on the electrical bias conditions and is dramatically enhanced at the highest injection current densities. The origin of the hysteresis in the spin polarization signal is magnetization reversal in the ferromagnet, but the magnitude and shape of the observed loops depend on the effective field in the semiconductor. Modeling based on the results of optical pumping experiments demonstrates that the origin of the large effective field is dynamical nuclear polarization due to the spin-polarized current injected from the ferromagnet Johnson:ElectricalDNP-proposal . This approach to dynamical nuclear polarization in semiconductors is a simple alternative to the use of optical pumping or high magnetic fields as sources of spin-polarized electrons Smet:GateControl ; Barrett:QWNMR .

We report results from two heterostructures with different quantum well (QW) spin detectors. The samples are grown by molecular beam epitaxy on p+ GaAs (100) substrates and consist of p-AlxGa1-xAs/QW/n-AlxGa1-xAs/Fe/Al Isakovic:SpinInjection . Intrinsic setback layers are placed on both sides of the QW, and the resulting p-i-n structure forms a light-emitting diode (LED). The 50 Å thick iron and the 25 Å Al layers are grown in situ at $`0^{}`$ C. A Si $`\delta `$-doped layer ($`\delta `$ = $`2\times 10^{13}`$ atoms/cm2 for sample A and $`3\times 10^{13}`$ atoms/cm2 for sample B) is grown 25 Å below the Fe-semiconductor interface in order to form a tunneling contact. For sample A, $`x=0.2`$ and the QW is 100 Å Al0.1Ga0.9As with p-doped barriers. For sample B, $`x=0.1`$ and the QW is 100 Å GaAs with intrinsic barriers. The important difference between the samples is the effective $`g`$-factor $`g^{}`$ for electrons in the QW. For sample B, $`g^{}0.4`$, while sample A was designed to have an electron $`g`$-factor close to zero Chadi:AlGaAs\_g-factor . The magnetocrystalline anisotropy of the epitaxial Fe film results in easy and hard axes along  and \[01$`\overline{1}`$\], as shown in Fig. 1.

The samples are processed into 200 $`\mu `$m wide mesas. The ferromagnetic contact is 80 $`\mu `$m wide, and a 40 $`\mu `$m wide gold contact is evaporated over the center of the bar, which is cleaved into pieces 1–1.5 mm long. The devices are operated with the Schottky contact under reverse bias and the p-i-n LED under forward bias. Electrons tunnel from the Fe into the n-layer and recombine in the QW with holes supplied by the substrate. Electroluminescence (EL) is collected along the sample normal $`\widehat{𝐳}`$, and is dominated by heavy-hole recombination. Under these conditions the circular polarization of the EL, which we will refer to as the ELP signal, is equal to the electron spin polarization along $`\widehat{𝐳}`$ at the time of recombination OO:DyakonovPerel .

In a typical ELP measurement in the Faraday geometry, $`\widehat{𝐤}\widehat{𝐳}𝐁_a`$, where $`\widehat{𝐤}`$ is the direction of light propagation, and $`𝐁_a`$ is the applied magnetic field Ploog:Zhu:InGaAsPRL ; Jonker:Hanbicki:AlGaAsAPL . The inset of Fig. 1 shows the ELP signal measured in the Faraday geometry for sample A. The observed signal is proportional to the magnetization of the Fe film and is significantly larger than the background determined by either the photoluminescence (PL) polarization measured at a bias just below the EL threshold or the magnetoabsorption (MA) of the Fe film. This result indicates a steady-state spin polarization in the QW of at least 4% for applied fields above 2.1 T, the saturation field for Fe.

In the transverse field configuration ($`\widehat{𝐤}\widehat{𝐳}𝐁_a`$) shown in the left inset of Fig. 1, the magnetocrystalline anisotropy of the Fe film allows for a significant angle between the electron spin $`𝐒`$ and $`𝐁_a`$. The magnetization remains in the sample plane and $`𝐒`$ precesses out of the plane after injection into the semiconductor. This approach is related philosophically to the oblique Hanle effect OO:DyakonovPerel , in which optically pumped spins are created at an angle $`0<\theta <\pi /2`$ with respect to the applied field, an approach followed in the recent electrical spin injection experiment of Motsnyi et al. IMEC:Motsnyi:OHE . An important advantage of our approach is that we detect only the precessing component of the spin, making the measurement essentially immune to background effects.

Figures 2(a) and (b) show the ELP signal as a function of the transverse magnetic field, which is applied in the sample plane at a small angle (4$`{}_{}{}^{}\pm 2^{}`$) from the \[01$`\overline{1}`$\] axis. Data are shown at several different current densities for samples A and B. At the lowest current densities, just above the threshold for light emission, the ELP signal is nearly independent of the applied magnetic field and shows no hysteresis, although there is a polarization offset of 0 - 2.5%. As the current density is increased, a double loop structure develops. In the case of sample A \[Fig. 2(a)\], the magnitude of the polarization signal increases with increasing bias, reaching a maximum value of approximately 4%, determined from the peak-to-peak height of the loop. In the case of sample B \[Fig. 2(b)\], both the magnitude and the shape of the observed loops depend on the current density. A maximum in the signal is observed at a current density of approximately 3.5 A/cm2, above which the signal decreases and the central dip in the loop becomes much narrower. Photoluminescence measurements on the same structures show no polarization signal in a transverse magnetic field.

The origin of the ELP signal observed in Fig. 2 is precession of the spin injected from the ferromagnet. The observed polarization depends on the angle between $`𝐁_a`$ and the magnetization $`𝐌`$ of the Fe film as well as the precession frequency and spin relaxation time in the QW. Under the assumption that the magnetization reversal in the Fe film occurs by coherent rotation StonerWohlfarth , the angle between $`𝐌`$ and $`𝐁_a`$, and hence the orientation of the injected spin $`𝐒_0`$ can be determined from the magnetization data. The spin dynamics in the semiconductor are modeled using the formalism developed for the optical Hanle effect OO:DyakonovPerel :

$$\frac{d𝐒}{dt}=\frac{𝐒_0}{\tau }\frac{𝐒}{\tau _s}\frac{𝐒}{\tau }+𝛀\times 𝐒,$$
(1)
where $`\tau _s^1`$ is the spin relaxation rate, $`\tau ^1`$ is the recombination rate, and $`𝛀=g^{}\mu _B𝐁/\mathrm{}`$, where $`\mu _B`$ is the Bohr magneton.

The steady-state spin $`𝐒`$ that is measured by ELP can be found from Eq. 1 by setting $`d𝐒/dt=0`$. We will be interested in an upper bound for $`𝐒`$ and will hence assume that $`𝐒_0=ϵ\widehat{𝐦}`$, where $`ϵ=0.42`$ is the spin polarization of Fe Meservey:Tedrow:spin\_pol\_tunneling  and $`\widehat{𝐦}`$ is a unit vector parallel to $`𝐌`$. The parameters in Eq. 1 are thus reduced to two: the ratio $`\alpha =\tau _s/\tau `$, and a characteristic field scale $`B_{1/2}=\mathrm{}/(\mu _Bg^{}T_s)`$, where $`T_s^1=\tau _s^1+\tau ^1`$. We determine $`\alpha `$ and $`B_{1/2}`$ by fitting the results of Hanle effect measurements, in which electron spins are optically injected along $`\widehat{𝐳}`$ using circularly polarized light tuned near the band-edge of the QW barriers. The circular polarization of the PL as a function of the transverse magnetic field is shown for both samples in Fig. 2(f). Fits to the expected Lorentzian shape OO:DyakonovPerel  are shown as the solid curves. The dramatic difference in the widths of the two Hanle curves is due primarily to the different electron $`g`$-factors in the two quantum wells.

Once the magnetization of the ferromagnet and the Hanle curve are known for each sample, it is then possible to model the ELP signal without any free parameters. The result for sample A is shown as the solid curve in Fig. 2(c) along with data obtained at a current density of 3.6 A/cm2. (A field-independent offset has been removed from the experimental data.) The ELP signal vanishes at large applied fields because the injected spin is parallel to the field, and the signal vanishes at zero field because there is no precession. At intermediate fields, the magnetization is between the easy and hard axes, and the non-zero angle between $`𝐒_0`$ and $`𝐁_a`$ results in an ELP signal.

In spite of the general qualitative agreement in Fig. 2(c), the magnitude of the signal determined from our model is smaller than the experimental result by a factor of $`20`$, even though we have already assumed a spin injection efficiency of unity. More significantly, the model of Eq. 1 cannot replicate the current-driven change in the shape of the ELP loops observed for sample B without requiring either the spin lifetimes or $`g^{}`$ to increase by over an order of magnitude with increasing bias.

The most probable explanation for the discrepancy between the simple dynamical model and the experiment is that the effective magnetic field in the semiconductor is larger than the applied field and depends on the current density. This hypothesis is motivated by the observation of large internal magnetic fields in optical pumping experiments due to dynamical nuclear polarization (DNP) by electron spins Paget:GaAsNuclearSpinCoupling ; Flinn:Kerr:1stOptNMR\_single\_GaAsQW . In our experiment, the spins are injected electrically, and the magnitude of the nuclear polarization is determined by the current density Johnson:ElectricalDNP-proposal . The total magnetic field for electron spins is then the sum of the applied field $`𝐁_a`$ and a nuclear field Paget:GaAsNuclearSpinCoupling

$$𝐁_N=b_N\frac{(\widehat{𝐬}𝐁_a)𝐁_a}{B_a^2+B_0^2},$$
(2)
where $`\widehat{𝐬}`$ is a unit vector parallel to the injected spin, $`B_0`$ is the nuclear dipolar field and $`b_N`$ is the maximum nuclear field. We set $`B_0`$ to 2 G, a typical value for GaAs Paget:GaAsNuclearSpinCoupling , and allow $`b_N`$ to vary with current density. The total field $`𝐁=𝐁_a+𝐁_N`$ is then used in Eq. 1 to determine $`𝐒`$ and the expected ELP signal. The results of modeling with a nuclear field, adjusting $`b_N`$ as the only free parameter, are shown in Figs. 2(d) and (e). Agreement with the magnitude of the experimental results is obtained using values of $`b_N`$ between 0 and 9 kG, which are indicated next to each curve. These are comparable to the internal fields obtained in optical pumping experiments and are well below the maximum possible internal field of 5.3 T Paget:GaAsNuclearSpinCoupling ; Salis:OptNMR-PRB .

Including the effects of DNP also accounts for the major differences between the ELP curves observed for samples A and B. Although the shape of the ELP loops is nearly independent of current density for sample A, the position of the maximum signal observed for sample B shifts toward zero field as the current density increases, and the jump observed at larger applied fields disappears. The same trend is observed in the theoretical curves. The difference in behavior reflects the fact that the magnetization reversal in the two samples is essentially the same, while the widths of the Hanle curves differ by a factor of six. In the case of sample A, the total internal field is less than the half-width of the Hanle curve. In this case, the average spin precesses only a fraction of a cycle before recombination, and the shape of the ELP loop is determined only by the rotation of the magnetization in the ferromagnet. In contrast, the maximum nuclear field in sample B exceeds the half-width of the Hanle curve, and hence spin precession in the semiconductor changes the shape of the ELP signal as well as determining its magnitude.

Given that only one parameter was varied to obtain the theoretical curves shown in Fig. 2, the overall agreement with the experiment provides strong support for DNP by electrical spin injection. The most important discrepancy between the data and the model occurs near zero field, where the suppression of the experimental signal occurs over a much wider field range than in the model. The narrow region observed in the model is based on the nuclear dipole field of 2 G for bulk GaAs. However, the suppression of DNP over a much wider field range has been observed in optical pumping studies of ferromagnet-semiconductor interfaces Kawakami:Awschalom:FM\_imprinting\_in\_GaAs  and quantum dots Gammon:QD\_hyperfine\_interaction . Furthermore, the behavior near zero field is very sensitive to the details of the magnetic reversal process which may not be described completely by the simple coherent rotation model.

Additional evidence for the existence of a nuclear field comes from data obtained for fields applied near the $`[011]`$ axis. As can be seen from the square hysteresis loop shown in Fig. 1, the magnetization along this direction is essentially parallel to $`𝐁_a`$ for any field, and in this case the torque term in Eq. 1 is zero. As a result, no precession is expected. This is consistent with the behavior observed in the middle curve of Fig. 3(a), which was obtained for a magnetic field applied in the (100) plane at an angle of $`5^{}`$ with respect to .

The top and bottom curves in Fig. 3(a) were obtained for rotations of the field axis $`\pm 5^{}`$ out of the (100) plane. The curves in Fig. 3(b) are the results of the model of Eq. 1 using only the applied field (left) and including the nuclear field from Eq. 2 (right). The effective field $`b_N`$ of 8 kG is the same as that used to fit the hard axis data at the same current density, and the raw magnetization data were used to determine the direction of $`𝐒_0`$. As can be seen in Figs. 3(a) and (b), there is good agreement with the experiment only if the nuclear field is included. The extreme sensitivity of the modeling results to the presence of the nuclear field in the easy axis case follows from the dot product $`\widehat{𝐬}𝐁_a`$ in Eq. 2. This compensates for the fact that the torque term in Eq. 1 nearly vanishes for small angles between $`𝐒`$ and $`𝐁_a`$. The dip in the signal at zero field shown in the inset of Fig. 3(a) is due to dipolar relaxation.

The efficiency of the DNP process should decrease with increasing temperature NucMag:Abragam . A strong suppression of the ELP signal with increasing temperature is seen in both samples and is shown for sample B in Fig. 3(c). The change in the shape of the ELP hysteresis loop with temperature also mirrors the dependence on current density seen in Fig. 2.

The approach to spin injection outlined in this Letter realizes the possibility of using conventional ferromagnets, small magnetic fields, and DC electrical currents to create and manipulate spin polarized carriers in a semiconductor device. There remain several open questions, including a microscopic description of the DNP mechanism and its relation to the injection process. This will require a more detailed study of the transport properties of these systems, with the goal of achieving a purely electronic means of spin detection.

We acknowledge L. J. Sham and C. Ciuti for helpful discussions. This work was supported by ONR, the DARPA/ONR SPINS program, the University of Minnesota MRSEC (NSF DMR-0212032), and the Institute for Rock Magnetism. A.F.I. and B.D.S thank 3M for fellowship support.


text:
# Optimal strategy for a single-qubit gate and trade-off between opposite types of decoherence

## Abstract

We study reliable quantum information processing (QIP) under two different types of environment. First type is Markovian exponential decay, and the appropriate elementary strategy of protection of qubit is to apply fast gates. The second one is strongly non-Markovian and occurs solely during operations on the qubit. The best strategy is then to work with slow gates. If the two types are both present, one has to optimize the speed of gate. We show that such a trade-off is present for a single-qubit operation in a semiconductor quantum dot implementation of QIP, where recombination of exciton (qubit) is Markovian, while phonon dressing gives rise to the non-Markovian contribution.

Already at the early stage of quantum information theory, the implementation of quantum computing was found to be extremely challenging due to decoherence processes Unruh (1995). It was claimed, nevertheless, that a quantum computer could reliably perform quantum algorithms once the error per operation is below some threshold level which is approximately $`\delta 10^5`$ Aharonov and Ben-Or  (for some caveats coming from memory of environment see Alicki et al. (2002)). In any case, minimizing the error per gate is a reasonable strategy. Usually, one considers the decoherence time of the system $`\tau _\mathrm{d}`$ and the time of performing operation (quantum gate) $`\tau _\mathrm{g}`$, so that the error is given by the ratio $`\stackrel{~}{\delta }=\tau _\mathrm{g}/\tau _\mathrm{d}`$ DiVincenzo (1995). This suggests to search for systems with $`\stackrel{~}{\delta }`$ as small as possible, and the obvious strategy to diminish the error is to apply fast gates. One then tacitly assumes that the process of decoherence is independent of running the gate, which formally means that it is a Markovian process, where the error, to first order, grows linearly in time. Gate speed-up may be achieved by selecting materials to provide favorable spectrum characteristics De Rinaldis et al. (2002) or by applying techniques reducing unwanted transitions within the register space Pochung Chen et al. (2001) as well as outside this space (leakage) Tian and Lloyd (2000).

The assumption of Markovian character of noise is by no means obvious and, indeed, it has recently been found Alicki et al. (2002) that due to non-Markovian effects it may prove to be completely invalid. The notion of minimal decoherence model was introduced where the error occurs solely during gate operation, and it grows for fast gates as $`\delta 1/\tau _\mathrm{g}^2`$ (for the spectral density of the reservoir $`\omega ^3`$). An example of such a model is a degenerate system interacting with a bosonic field at vacuum via dipole interaction. This kind of error favors the counterintuitive strategy of slowing down the gate operation.

In real systems, this type of decoherence competes with the Markovian damping, although this effect may be covered by other errors (e.g. leakage). In this paper we point out that such a competitive effect is essential e.g. for the solid-state qubit implementation using excitonic (charge) states in quantum dots (QDs) Biolatti et al. (2000), with computational states defined by the absence ($`|0`$) or presence ($`|1`$) of one exciton in the ground state of the dot, operated by resonant coupling to laser light. In such a system, Markovian decoherence (exponential damping) is related to exciton recombination on 1 ns timescale, while strongly non-Markovian effects result from lattice inertia. The interplay of these two decoherence mechanisms, favoring opposite strategies (fast vs. slow gates), leads to a kind of trade-off, resulting in optimal speed of gate for most reliable operation on the qubit.

To investigate the effect, we analyze decoherence in a general spin-boson system van Hove (1952); Alicki , and find decoherence to be strongly non-Markovian, fitting into “minimal decoherence model”. Assuming additional, Markovian damping, we obtain the trade-off formula for the error caused by decoherence, averaged over input states of the qubit

$$\delta =\frac{\gamma _{\mathrm{nM}}}{\tau _\mathrm{g}^2}+\gamma _\mathrm{M}\tau _\mathrm{g}$$
(1)
(actually, in the solid state example, singularity of the first term is lifted by the presence of upper cut-off). The constants $`\gamma _{\mathrm{M},\mathrm{nM}}`$ express the strength of the Markovian and non-Markovian decoherence, respectively. We determine their values for the exciton in a QD with typical parameters and show that the gate duration leading to minimal overall error is of the order of $`1`$ ps. Since level spacing in such a system allows even 100 fs gating Pochung Chen et al. (2001) the restrictions imposed by non-Markovian phonon mechanisms are decisive for the gate speed. Our result has both practical as well as general implications for QIP: (i) it suggests the proper direction in research towards semiconductor implementation of QIP (ii) it provides a non-trivial dynamical strategy to minimize decoherence in quantum systems. for a class of reservoirs.

Minimal decoherence in the spin-boson model. We consider a qubit described by means of the spin-boson model. The Hamiltonian $`H^{\prime \prime }`$ of the system plus reservoir is

$$H^{\prime \prime }=H_S^0+H_S(t)+H_{SR}+H_R,$$
(2)
where $`S`$ is our qubit system, $`R`$ is (bosonic) reservoir. $`H_S^0=(1/2)\mathrm{\Omega }\sigma _z`$ is the self Hamiltonian of the qubit, $`H_S(t)=(1/2)ϵ(t)(e^{i\mathrm{\Omega }t}\sigma _++e^{i\mathrm{\Omega }t}\sigma _{})`$ is the gate Hamiltonian (in rotating wave approximation), $`ϵ(t)`$ being the shape of the pulse. $`H_R=\omega _𝒌a_𝒌a_𝒌^{}`$ is the self Hamiltonian of the reservoir, where $`\omega _𝒌`$ is the boson energy. Finally, $`H_{SR}`$ is the interaction with reservoir

$$H_{SR}=\sigma _z[f_𝒌^{}a_𝒌+\mathrm{H}.\mathrm{c}.],$$
(3)
where $`f_𝒌`$ is the coupling constant for the mode $`𝒌`$.

In the interaction picture with respect to $`H_S^0`$ we get $`H^{}=ϵ(t)\sigma _x+H_{SR}+H_R.`$ Let us now represent $`H^{}`$ in the basis of eigenstates of the total system (dressed states). This is achieved by the unitary operation $`U=|00|W|11|W^{}`$, where $`W=\mathrm{exp}[(f_𝒌^{}/\omega _𝒌)a_𝒌\mathrm{H}.\mathrm{c}]`$ and $`|0,|1`$ are the eigenvectors of $`\sigma _z`$. In this basis the Hamiltonian $`H=U^{}H^{}U`$ up to a constant has the form

$`H`$ $`=`$ $`{\displaystyle \frac{1}{2}}ϵ(t)\sigma _x\mathrm{cos}(2i{\displaystyle \underset{𝒌}{}}{\displaystyle \frac{f_𝒌}{\omega _𝒌}}a_𝒌^{}+\mathrm{H}.\mathrm{c}.)`$
$`+{\displaystyle \frac{1}{2}}ϵ(t)\sigma _y\mathrm{sin}(2i{\displaystyle \underset{𝒌}{}}{\displaystyle \frac{f_𝒌}{\omega _𝒌}}a_𝒌^{}+\mathrm{H}.\mathrm{c}.)+H_R.`$
To the first order in $`f_𝒌/\omega _𝒌`$ we obtain

$$Hϵ(t)\sigma _x+ϵ(t)\sigma _y(i\underset{𝒌}{}\frac{f_𝒌}{\omega _𝒌}a^{}+\mathrm{H}.\mathrm{c}.)+H_R.$$
(5)
We define the qubit in terms of the dressed states Alicki . With such a choice, the reservoir is decoupled and the interaction term is present only during gate operation. Thus our system may be reduced to the minimal decoherence model: decoherence is not present at all, if the qubit is not active. Now we want to calculate the resulting error and relate it to the speed of gate. Since the regime is strongly non-Markovian we solve the Master equation in the Born approximation and compute the fidelity $`F=1\delta `$ of a single-qubit operation (see Alicki et al. (2002) for details).

The error, averaged over the initial qubit states can be represented as the overlap of two functions

$$\delta =\frac{\mathrm{d}\omega }{\omega ^2}R(\omega )S(\omega ).$$
(6)
Here $`R(\omega )`$ is the spectral density of the reservoir

$$R(\omega )=[n_\mathrm{B}(\omega )+1]\underset{𝐤}{}\left[\delta (\omega _𝐤\omega )+\delta (\omega _𝐤+\omega )\right]|f_𝒌|^2,$$
(7)
where $`n_\mathrm{B}(\omega )`$ is the Bose-Einstein distribution. The function $`S(\omega )`$ fully represents the spectral characteristics of the system and is given by

$`S(\omega )`$ $`=`$ $`\left[\psi |YY^{}|\psi \psi |Y^{}|\psi \psi |Y|\psi \right]_{\mathrm{av}}`$
$``$ $`{\displaystyle \frac{1}{3}}(|F_{}(\omega )|^2+|F_+(\omega )|^2),`$
where $`[]_{\mathrm{av}}`$ denotes averaging over the states $`\psi `$ and

$`Y=i(F_+|+|+F_{}|+|),`$ $`|\pm ={\displaystyle \frac{|0\pm |1}{\sqrt{2}}};`$
$`F_\pm =\pm {\displaystyle _{\mathrm{}}^+\mathrm{}}due^{\pm i\varphi (u)}ϵ(u)e^{i\omega u};`$ $`\varphi (t)={\displaystyle _{\mathrm{}}^t}duϵ(u).`$
A complete minimization of $`\delta `$ would require full optimization of the pulse shape. However, in order to demonstrate the idea of the trade-off in simple terms, let us restrict the discussion to qubit rotations performend by Gaussian pulses, $`ϵ(t)=\alpha /(\sqrt{2\pi }\tau _\mathrm{g})e^{\frac{1}{2}\left(t/\tau _\mathrm{g}\right)^2}`$. Here $`\tau _\mathrm{g}`$ is the gate duration, while $`\alpha `$ is the angle determining the gate, e.g. $`\alpha =\frac{\pi }{2}`$ is $`\sqrt{\mathrm{NOT}}`$, while $`\alpha =\pi `$ is $`\sigma _x`$ (bit flip).

The functions $`|F_\pm (\omega )|^2`$ may be written as

$$|F_\pm (\omega )|^2\alpha ^2e^{\tau _g^2\left(\omega \pm \frac{\alpha }{\sqrt{2\pi }\tau _g}\right)^2}.$$
(9)
Although the original minimal decoherence model with its $`\omega ^3`$ dependence appears in many physical situations, other characteristics are obviously possible. In general, as may be seen from (6), (7) and (9), for a spectral density $`R(\omega )\omega ^n`$ the error scales with the gate duration as $`\tau _\mathrm{g}^{n+1}`$ and $`\tau _\mathrm{g}^{n+2}`$ at low and high temperatures, respectively. Therefore, for $`n>2`$ (typical e.g. for various types of phonon reservoirs Krummheuer et al. (2002)) the error grows for faster gates. Assuming the spectral density of the form $`R(\omega )=R_0\omega ^3`$ (see the QD example below), we obtain

$$\delta _{\mathrm{nM}}=\frac{1}{3}\alpha ^2R_0\tau _\mathrm{g}^2,\text{at}T=0$$
(10)
This leading order formula holds for $`\delta <<1`$. Also, if we introduce the upper cut-off, the error will be finite even for an infinitely fast gate Jacak et al. (2003a) (see Fig. 1).

Trade-off between two types of decoherence. As we have shown, in our model the error grows as the speed of gate increases. This could result in obtaining arbitrarily low error by choosing suitably low speed of gates. However, if the system is also subject to other types of noise this becomes impossible. Indeed, assuming an additional contribution growing with rate $`\gamma _\mathrm{M}`$, the total error per gate is

$$\delta =\frac{\gamma _{\mathrm{nM}}}{\tau _\mathrm{g}^2}+\gamma _\mathrm{M}\tau _\mathrm{g},\gamma _{\mathrm{nM}}=\frac{1}{3}\alpha ^2R_0,\gamma _\mathrm{M}=\frac{1}{\tau _\mathrm{r}},$$
(11)
where $`\tau _\mathrm{r}`$ is the characteristic time of Markovian decoherence (recombination time in the excitonic case). As a result, the overall error is unavoidable and optimization is needed. The above formulas lead to the optimal values of the form (for $`T=0`$)

$$\delta _{\mathrm{min}}=\frac{3}{2}\left(\frac{2\alpha ^2R_0}{3\tau _\mathrm{r}^2}\right)^{1/3},\text{for}\tau _\mathrm{g}=\left(\frac{2}{3}\alpha ^2R_0\tau _\mathrm{r}\right)^{1/3}.$$
(12)
We will see that the trade-off is present also for finite temperatures.

Example: exciton in a quantum dot (QD). Let us now estimate the error magnitude for the specific semiconductor QD qubit implementation Biolatti et al. (2000). The reservoir is then constituted by phonons, the most important branch for our present purpose being the longitudinal acoustical (LA) one, characterized by linear dispersion, $`\omega _k=ck`$, coupled via deformation potential (DP) Mahan (2000). In fact, the relatively simple model (2) is very accurate for a description of the QD system for $`1`$ ps timescales relevant here, as confirmed by the excellent agreement between the theoretical calculations Vagov et al. (2002) and experimental results Borri et al. (2001). This is due to the fact that neither the high-frequency optical phonons nor direct or phonon-induced leakage to higher states contribute considerably to the decoherence. Also effects related to piezoelectric coupling to LA and TA phonons may be neglected in weakly piezoelectric systems (e.g. GaAs). Moreover, the qubit control is all-optical, eliminating the need for additional noise-inducing device structures. On the other hand, details of the QD structure (shape, stress, composition) may lead only to quantitative modifications of secondary importance.

The non-Markovian error has a dynamical origin and is related to a which path trace left by exciting the phonon modes rather than to an influence of noise. The spectral density is uniquely defined by the lattice response characteristics: the mode-dependent coupling strength and the density of states. Due to fundamental restrictions (global charge neutrality, translational invariance), the frequency dependence is always strongly super-ohmic Krummheuer et al. (2002). On the other hand, the approximate momentum conservation holding for weakly confined carrier states leads to exponential cut-off of carrier-phonon interaction at high frequencies Jacak et al. (2002). In the specific case of DP coupling, the coupling constants are Mahan (2000)

$$f_𝒌=\frac{1}{2}\sqrt{\frac{k}{2\varrho vc}}\left(\sigma _\mathrm{e}_𝒌^{(\mathrm{e})}\sigma _\mathrm{h}_𝒌^{(\mathrm{h})}\right),$$
(13)
where $`v`$ is the volume of the unit cell, $`\rho `$ is the crystal density, $`\sigma _{\mathrm{e},\mathrm{h}}`$ are deformation potential constants for the electron and the hole (material parameters are taken as in Krummheuer et al. (2002)) and $`_𝒌^{(\mathrm{e},\mathrm{h})}`$ are formfactors for the corresponding wavefunctions (approximated by Gaussians, on the grounds of numerical diagonalization in parabolic confinement Jacak et al. (2003a)),

$$_𝒌^{(\mathrm{e},\mathrm{h})}=e^{\frac{1}{4}(k_{}^2l_{\mathrm{e},\mathrm{h}}^2+k_z^2l_z^2)},$$
(14)
where $`l_{\mathrm{e},\mathrm{h}}`$ and $`l_z`$ are the wavefunction widths in the dot plane, for electron and hole, and in the growth direction, respectively, and $`k_{}`$ and $`k_z`$ are the corresponding components of the phonon wavevector.

Note that the most effectively coupled phonon modes correspond to wavelengths around the dot size. For gapless bosons with linear dispersion, their frequency determines the reservoir memory times and sets up the timescale of non-Markovian effects which, in the case of the QD system, may be interpreted as “dressing” the localized carriers with coherent lattice deformation field (cf. Alicki ). This kind of dynamics has been described in the limit of very rapid gating pulses Jacak et al. (2003a); Krummheuer et al. (2002) and has been shown to lower the degree of coherent control over the system and to destroy the coherence of polarization oscillations Krummheuer et al. (2002); Vagov et al. (2002). Due to relatively large dot sizes and low sound speed the cut-off frequency ($`1`$ meV) is lower than the transition energy between different carrier states in a small, self-assembled dot and the corresponding 1 ps times are experimentally accessible.

For large enough gate duration only the low-frequency sector contributes and the coupling (13) may be approximated by

$$f_𝒌\frac{1}{2}(\sigma _\mathrm{e}\sigma _\mathrm{h})\sqrt{\frac{k}{2\varrho vc}},$$
(15)
leading, according to (7), to

$$R(\omega )R_0\omega ^3,R_0=\frac{(\sigma _\mathrm{e}\sigma _\mathrm{h})^2}{16\pi ^2\varrho c^5}.$$
(16)
Hence, Eqs. (12) are applicable and one finds for the specific material parameters of GaAs,

$$\tau _\mathrm{g}=\alpha ^{2/3}1.47\mathrm{ps},\delta _{\mathrm{min}}=\alpha ^{2/3}0.0035.$$
(17)
The full solution within the proposed model, taking into account the phonon cut-off for an anisotropic shape ($`l_z<l_{e,h}`$), according to Eqs. (14), (13), and allowing finite temperatures by numerically calculating the spectral density (7), is shown in Fig. 1. The size-dependent cut-off is reflected by a shift of the optimal parameters for the two dot sizes: larger dots admit slightly faster gates and lead to lower error. Interestingly, the trade-off becomes more apparent at nonzero temperature.

The presence of the upper cut-off could suggest applying the dynamical decoupling (DD) technique Viola and Lloyd (1998) to diminish decoherence. However, for high frequencies many other mechanisms of decoherence become relevant. On the other hand, combining the bounded-control version of DD Viola and Knill (2003) with the optimization proposed here might lead to some reduction of the resulting error. Such techniques might also be useful for eliminating other sources of noise, not included in our discussion.

Another possibility to reduce the decoherence effect in the QD system might be to encode qubits in decoherence-free subspaces Zanardi and Rossi (1998). However, the phonon wavelength corresponding to the optimal range of gating times is comparable to the single dot size, precluding the necessary collective interaction with the whole QD system. Thus, for decoherence effects which do not involve real transitions and taking into account the feasible system geometry and actual nature of phonon coupling, only a small decrease of the minimal decoherence may be expected.

In conclusion, we have exhibited the trade-off between two opposite types of decoherence: usual Markovian damping and dynamically induced non-Markovian decoherence for a realistic super-ohmic reservoir. To protect the qubit, opposite elementary strategies are needed: fast and slow gates, respectively. The minimization of the overall error leads to optimal speed of gate. We have shown that the trade-off is present in a semiconductor implementation of quantum information processing. The Markovian error is caused by recombination, while the non-Markovian one occurs if the gate operation is not adiabatic with respect to lattice modes. We have evaluated the minimal error in this case for a single qubit gate, showing that the two processes indeed compete. The optimal gating time ($`1`$ ps) sets up the limit beyond which any further gate speed-up is unfavorable. Even at this optimal point the trade-off gives rise to a significant error. For two qubit gates involving single-qubit rotations (as proposed e.g. in Biolatti et al. (2000)), the present result gives a rough lower bound for the error, which is of 1-2 orders of magnitude higher than that admitted by fault-tolerant schemes Aharonov and Ben-Or  known so far ($`10^5`$). However, possible improvements of the latter schemes cannot be excluded. It follows also that diminishing the responsible constants $`\gamma _\mathrm{M}`$ (e.g. by elimination of radiative losses Troiani et al. (2003)) and $`\gamma _{\mathrm{nM}}`$ (by optimizing the system parameters) is the most important task towards semiconductor implementation of a quantum computer. It is also important to explore whether the same effect can occur in other implementations, as well as to what extent the error avoiding techniques may be helpful here.

We thank Peter Knight and Martin Plenio for useful feedback to the first version of this paper. L. J. and P. M. are grateful to T. Kuhn and V. M. Axt for discussions. This work was supported by EC grants EQUIP (IST-1999-11053), RESQ (IST-2002-37559), SQID (IST-1999-11311), QUPRODIS (IST-2001-38877), and by the Polish government grant PBZ-MIN-008/P03/2003. P.M. is grateful to the Humboldt Foundation for support.