format: text

text:
# 1 Introduction

## 1 Introduction

In QED quantized in covariant gauge, longitudinally polarized on–shell photons are present, but due to gauge invariance decouple, order by order in perturbation theory, in expressions for physical quantities. For the virtual photon with nonzero virtuality 111In this paper the virtuality $`\tau `$ of a particle with four–momentum $`k`$ and mass $`m`$ is defined as $`\tau m^2k^2`$. In this convention, $`P^2>0`$ in the space–like region relevant for hard collisions involving photons in the initial state. its longitudinal polarization, denoted $`\gamma _L^{}`$, does, however, give nonzero contributions to physical quantities and gauge invariance merely requires that these contributions vanish as $`P^20`$. In this paper we discuss hard collisions 222Characterized by some “hard scale” denoted generically as $`Q^2`$. In practice “$`Q^2`$” may be standard $`Q^2`$ in DIS, $`E_T^2`$ in jet studies, $`M_Q^2`$ in heavy quark production, etc. of virtual photons and, moreover, restrict our attention to kinematical region $`P^2Q^2`$ where the concept of virtual photon structure makes good sense.

This work has been motivated by the lack of the general attitude toward the role of $`\gamma _L^{}`$ in hard collisions 333The relevance of $`\gamma _L^{}`$ has recently been pointed out in  as well. and in particular by our disagreement with the statements of a recent paper , where the treatment of virtuality dependence of physical quantities is based on the following two claims 444$`f_{\gamma (P^2)/e}^L`$ in the notation of  corresponds to $`f_L^\gamma (P^2)`$ in our formula (3).

(i) effects of $`f_{\gamma (P^2)/e}^L`$ should be neglected since the corresponding longitudinal cross–sections are suppressed by powers of $`P^2/Q^2`$, and

(ii) cross–sections of partonic subprocesses involving $`\gamma (P^2)`$ should be calculated as if $`P^2=0`$ due (partly) to the $`P^2/Q^2`$ suppression of any additional terms.

In the next Section we first analyze these claims and point out why they are wrong. Then we recall the reasons for introducing the concept of virtual photon structure and recollect basic formulae concerning the structure of $`\gamma _L^{}`$. Numerical results illustrating the importance of including the contributions of $`\gamma _L^{}`$ are presented in Section 3. The feasibility of extracting the information of partonic content of $`\gamma _L^{}`$ from jet production at HERA is addressed in Section 4, followed by the summary and conclusions in Section 5.

## 2 Theoretical considerations

### 2.1 Virtuality dependence of $`\sigma (\gamma _L^{})`$

Before recalling practical usefulness of the concept of partonic structure of virtual photons, let us show why the claims made in  and mentioned in the Introduction are incorrect. The fact that in the resolved photon channel the cross–sections of $`\gamma _L^{}`$ are not suppressed by $`P^2/Q^2`$ follows directly from analysis of the formula (E.1) in  for the cross–section $`\sigma _{TL}`$ (denoted $`\sigma _{TS}`$ there), which shows that for small $`P^2m_q^2`$ ($`m_q`$ being the quark mass) its contribution to $`F_2^\gamma (x,P^2,Q^2)`$ behaves as 555In the notation of  the first and second subscripts in $`\sigma _{ij}`$ refer to polarizations of the probing and target photons respectively, with wirtualities $`Q^2`$ and $`P^2`$. Most of the terms in the expression for $`\sigma _{TL}`$ do, indeed, behave as $`P^2/Q^2`$, but there is one, proportional to $`(\mathrm{\Delta }tq_1^4/T)`$, which does not and which yields (1).

$$F_{TL}^\gamma (x,P^2,Q^2)=\frac{P^2}{m_q^2}\frac{\alpha }{\pi }4x^3(1x)^2=\frac{P^2}{m_q^2}2x\frac{\alpha }{2\pi }4x^2(1x)^2.$$
(1)
This expression coincides, apart from the factor $`N_ce_q^2`$ appropriate to a quark with electric charge $`e_q`$ and $`N_c`$ colors, with the QED expression for distribution function of quarks inside $`\gamma _L^{}`$ in our formulae (9) below 666In practical applications the factorization scale $`M^2`$ in (5) is identified with the generic hard scale $`Q^2`$. multiplied by $`2x(\alpha /2\pi )`$. For $`P^2m_q^2`$, on the other hand, the distribution function (8) is proportional to $`4x(1x)`$ with no $`P^2/Q^2`$ suppression. As a result, in the region $`P^2m_q^2`$, $`\gamma _L^{}`$ supplies finite contribution to $`F_2^\gamma (x,P^2,Q^2)`$ equal to $`(\alpha /\pi )(N_ce_q^2)4x(1x)`$ . Similarly, also the second claim (ii) is incorrect, because also part of the contribution of $`\sigma _{TT}`$ has the same $`P^2`$ behaviour as in (1). Physical explanation of this behaviour is simple: even for large values of $`Q^2`$ the virtuality $`\tau `$ of the quarks (antiquarks) from the primary splitting $`\gamma ^{}q\overline{q}`$ of the target photon comes predominantly from the region close to its minimal value $`\tau ^{\mathrm{min}}=xP^2+m_q^2/(1x)`$ and therefore the threshold behaviour is governed by the quark mass $`m_q`$ rather than $`Q^2`$.

On the other hand, virtuality dependence of the contributions of $`\gamma _T^{}`$ and $`\gamma _L^{}`$ can be safely neglected in the LO direct photon hard processes, for instance in large $`E_T`$ jet production via the photon–gluon fusion subprocess $`\gamma ^{}Gq\overline{q}`$. In these processes virtuality of the exchanged quark (or antiquark) is forced by kinematics to be proportional to jet transverse energy $`E_T`$ and therefore the virtuality dependent part is suppressed by powers of $`P^2/E_T^2`$. Of course, in realistic QCD the onset of quark distribution functions of $`\gamma _L^{}`$ is not expected to be determined directly by quark masses, but rather by some nonperturbative parameter related to confinement, but the basic features of the dependence on $`P^2`$, exemplified in (1), are likely to persist.

### 2.2 Equivalent photon approximation

Most of the present knowledge of the structure of the photon comes from experiments at the ep and e+e- colliders, where the incoming leptons act as sources of transverse and longitudinal virtual photons. To order $`\alpha `$ their respective unintegrated fluxes are given as

$`f_T^\gamma (y,P^2)`$ $`=`$ $`{\displaystyle \frac{\alpha }{2\pi }}\left({\displaystyle \frac{1+(1y)^2)}{y}}{\displaystyle \frac{1}{P^2}}{\displaystyle \frac{2m_\mathrm{e}^2y}{P^4}}\right),`$ (2)
$`f_L^\gamma (y,P^2)`$ $`=`$ $`{\displaystyle \frac{\alpha }{2\pi }}{\displaystyle \frac{2(1y)}{y}}{\displaystyle \frac{1}{P^2}}.`$ (3)
Phenomenological analyses of interactions of virtual photons have so far concentrated on its transverse polarization. The same holds for available parameterizations of parton distribution functions (PDF) of virtual photons. Neglecting longitudinal photons is in general a good approximation for $`y1`$, where the flux $`f_L^\gamma (y,P^2)0`$, as well as for very small virtualities $`P^2`$, where PDF of $`\gamma _L^{}`$ vanish by gauge invariance. But how small is “very small” in fact? For instance, should we take into account the contribution of $`\gamma _L^{}`$ to jet cross–section in the region $`E_T5`$ GeV, $`P^21`$ GeV2, where most of the data on virtual photons obtained in ep collisions at HERA come from? The rest of this paper is devoted to addressing this and related questions.

### 2.3 Who needs the concept of partonic structure of virtual photons?

Let us briefly recall the virtue of extending the concept of partonic “structure” to virtual photons. The arguments for it were discussed in detail in  and we therefore merely summarize the most important points:

* In principle, the concept of partonic structure of (sufficiently) virtual photons can be dispensed with because higher order perturbative QCD corrections to cross–sections of processes involving virtual photons in the initial state are well–defined and finite even for massless partons.
* In practice, however, this concept is extraordinarily useful as it allows us to include the resummation of higher order QCD effects that come from physically well–understood region of (almost) parallel emission of partons off the quark (or antiquark) coming from the primary $`\gamma ^{}q\overline{q}`$ splitting and subsequently participating in hard processes.

In other words, for the virtual photon, as opposed to the real one, its PDF can be regarded as “merely” describing higher order perturbative effects and not the “true” structure. Although this distinction between the content of PDF of real and virtual photons does exist, it does not affect the extraordinary phenomenological usefulness of PDF of the virtual photon. As shown in  the nontrivial part of the resolved photon contributions to NLO calculations of dijet production at HERA obtained with JETVIP  is large and affects significantly the conclusions of phenomenological analyses of existing experimental data.

### 2.4 Structure of $`\gamma _L^{}`$ in QED

The definition and evaluation of quark distribution functions of the virtual photon in QED serves as a guide to QCD improved parton model predictions of virtuality dependence of their pointlike parts. In pure QED and to order $`\alpha `$ the probability of finding inside $`\gamma _T^{}`$ or $`\gamma _L^{}`$ of virtuality $`P^2`$ a quark with mass $`m_q`$, electric charge $`e_q`$, momentum fraction $`x`$ and virtuality $`\tau M^2`$, is given, in units of $`3e_q^2\alpha /2\pi `$, as  ($`k=T,L`$)

$$q_k^{\mathrm{QED}}(x,m_q^2,P^2,M^2)=f_k(x)\mathrm{ln}\left(\frac{M^2}{\tau ^{\mathrm{min}}}\right)+\left[f_k(x)+\frac{g_k(x)m_q^2+h_k(x)P^2}{\tau ^{\mathrm{min}}}\right]\left(1\frac{\tau ^{\mathrm{min}}}{M^2}\right),$$
(4)
where $`\tau ^{\mathrm{min}}=xP^2+m_q^2/(1x)`$. The quantity defined in (4) has a clear physical interpretation: it describes the flux of quarks and antiquarks that are almost collinear with the incoming photon and “live” longer 777In fact most of these quarks live much longer than $`1/M`$. than $`1/M`$. For $`\tau ^{\mathrm{min}}M^2`$ the expression (4) simplifies to

$$q_k^{\mathrm{QED}}(x,m_q^2,P^2,M^2)=f_k(x)\mathrm{ln}\left(\frac{M^2}{xP^2+m_q^2/(1x)}\right)f_k(x)+\frac{g_k(x)m_q^2+h_k(x)P^2}{xP^2+m_q^2/(1x)},$$
(5)
which for $`x(1x)P^2m_q^2`$ reduces further to

$$q_k^{\mathrm{QED}}(x,0,P^2,M^2)=f_k(x)\mathrm{ln}\left(\frac{M^2}{xP^2}\right)f_k(x)+\frac{h_k(x)}{x}.$$
(6)
The functions $`f_k,g_k,h_k`$ are given as

$$\begin{array}{ccc}f_T(x)=x^2+(1x)^2,\hfill & g_T(x)=\frac{1}{1x},\hfill & h_T(x)=0,\hfill \\ f_L(x)=0,\hfill & g_L(x)=0,\hfill & h_L(x)=4x^2(1x).\hfill \end{array}$$
(7)
For $`M^2x(1x)P^2`$ the quark distribution function of $`\gamma _L^{}`$ has a simple form

$`q_L^{\mathrm{QED}}(x,m_q^2,P^2,M^2)={\displaystyle \frac{4x^2(1x)^2P^2}{x(1x)P^2+m_q^2}}`$ $``$ $`4x(1x);x(1x)P^2m_q^2`$ (8)
$``$ $`{\displaystyle \frac{P^2}{m_q^2}}4x^2(1x)^2;x(1x)P^2m_q^2`$ (9)
### 2.5 QCD corrections

For $`\gamma _T^{}`$ QCD corrections to QED formula (6) are well understood. Though important, in particular for large and very small $`x`$, they do not change its basic features and the main nontrivial effect comes from the emergence of gluons inside $`\gamma _T^{}`$. For $`\gamma _L^{}`$ the effects of collinear parton radiation off the quarks/antiquarks from the $`\gamma _L^{}q\overline{q}`$ splitting result in factorization scale dependence that resembles those of hadrons and will be discussed in separate paper. For the purpose of this exploratory study we use the QED formula (8) throughout this paper.

## 3 Numerical results

### 3.1 DIS on $`\gamma ^{}`$ in QED

The cleanest evidence of the importance of taking into account the contribution of $`\gamma _L^{}`$ has been provided by the L3 and OPAL measurements  of the QED structure function $`F_2^{\gamma ,\mathrm{QED}}`$ at LEP. In these measurements, based on the analysis of $`\mu ^+\mu ^{}`$ final states, the average target photon virtuality is small ($`P^2=0.033`$ GeV2 in  and $`P^2=0.05`$ GeV2 in ) but still sufficiently large with respect to $`m_\mu ^20.01`$ GeV2 to see the decrease of $`F_2^{\gamma ,\mathrm{QED}}(x,P^2,Q^2)`$ with respect to the QED prediction for the real photon. To order $`\alpha `$ these predictions were calculated exactly in  and contain contributions of both transverse and longitudinal polarizations of the target photon. In the region $`m_e^2P^2Q^2`$ experiments at LEP actually measure the following sum of $`\gamma ^{}\gamma ^{}`$ cross–sections, the first and second indices corresponding to probe and target photon respectively,

$$F_{\mathrm{eff}}^\gamma (x,P^2,Q^2)\frac{Q^2}{4\pi ^2\alpha }\left(\sigma _{TT}+\sigma _{LT}+\sigma _{TL}+\sigma _{LL}\right)=\frac{Q^2}{4\pi ^2\alpha }\sigma (P^2,Q^2,W^2),$$
(10)
where all cross–sections $`\sigma _{jk}`$ are functions of $`W^2,P^2`$ and $`Q^2`$ and $`x=Q^2/(W^2+Q^2+P^2)`$. As shown in  the data are in very good agreement with QED prediction for (10) provided the dependence on target photon virtuality $`P^2`$ is taken into account. For OPAL kinematical region the QED predictions for $`f(x,P^2,Q^2)(2\pi /\alpha )F_{\mathrm{eff}}^\gamma (x,P^2,Q^2)`$ as well as the individual contributions $`f_{ij}`$, are shown in Fig. 1a, together with the results (shown as dotted curves) corresponding to the real photon and the approximations using formulae of the preceding Section (dashed curves). The variations of $`f_{jk}(x,P^2,Q^2)`$ and $`f(x,P^2,Q^2)`$ with respect to real photon, defined as ($`i,j=T,L`$)

$$\mathrm{\Delta }f_{jk}(x,P^2,Q^2)f_{jk}(x,P^2,Q^2)f_{jk}(x,0,Q^2)),$$
(11)
are plotted in Fig. 1b. The contribution $`\mathrm{\Delta }f_{TL}=f_{TL}`$ to the variation $`\mathrm{\Delta }F_{\mathrm{eff}}^\gamma (x,P^2,Q^2)`$ coming from target $`\gamma _L^{}`$ is clearly comparable in magnitude to $`\mathrm{\Delta }f_{TT}`$ and $`\mathrm{\Delta }f_{LT}`$ coming from target $`\gamma _T^{}`$. Neglecting $`\mathrm{\Delta }f_{TL}`$ would thus lead to serious disagreement between QED predictions and data.

### 3.2 DIS on $`\gamma ^{}`$ in QCD

In LO QCD the structure function $`F_2^\gamma `$ is given in terms of quark distribution functions by the same expression as for hadrons 888In the present paper we disregard the consequences of the reformulation of QCD analysis of $`F_2^\gamma `$ proposed in  as they do not concern the main point of our discussion.

$$F_2^\gamma (x,P^2,Q^2)=\underset{i}{}2xe_i^2\left(q_i(x,P^2,Q^2)+\overline{q}_i(x,P^2,Q^2)\right).$$
(12)
In all existing phenomenological analyses of experimental data only target $`\gamma _T^{}`$ has been taken into account and to the best of our knowledge no attempt has been made to extract PDF of $`\gamma _L^{}`$ therefrom. In this exploratory study we compare the results for $`F_2^\gamma `$ obtained with Schuler–Sjöstrand (SaS) parameterization  of $`q_T(x,P^2,M^2)`$ with the QED prediction (8) for $`q_L(x,P^2,M^2)`$.

In Fig. 2 this comparison is performed for typical values of $`P^2`$ and $`Q^2`$ accessible at LEP and $`m_q^2=1,0.1`$ GeV2 and $`m_q^2=0`$. The importance of the contributions of $`\gamma _L^{}`$ with respect to those of $`\gamma _T^{}`$ depends sensitively on the value of $`m_q`$: whereas for $`m_q1`$ GeV, $`\gamma _L^{}`$ is largely irrelevant, for $`m_q0.3`$ GeV, medium values of $`x`$ and $`Q^2100`$ GeV2, its contributions in the considered region of $`P^2`$ and $`Q^2`$ are comparable to those of SaS1D parameterization of $`\gamma _T^{}`$. Only for very large $`Q^2`$ does $`\gamma _L^{}`$ become really negligible with respect to $`\gamma _T^{}`$. For fixed $`Q^2`$ the relative importance of $`\gamma _L^{}`$ with respect to $`\gamma _T^{}`$ grows with $`P^2`$, but to retain clear physical meaning of PDF we stay throughout this paper in the region where $`P^2Q^2`$.

The comparison of the contributions of $`\gamma _L^{}`$ and $`\gamma _T^{}`$ at the same values of $`P^2`$ is one measure of the relevance of $`\gamma _L^{}`$. If we are interested in virtuality dependence of $`F_2^\gamma (x,P^2,Q^2)`$, the appropriate comparison is with difference

$$\mathrm{\Delta }F_2^{\gamma _T}(x,P^2,Q^2)F_2^{\gamma _T}(x,0,Q^2)F_2^{\gamma _T}(x,P^2,Q^2),$$
(13)
of SaS results for $`\gamma _T^{}`$, denoted in Fig. 2 by thin solid curves. At small to moderate $`x`$, lower hard scales $`Q^2`$ and larger virtualities $`P^2`$, the contributions of $`\gamma _L^{}`$ appear by this measure less important than when compared to $`F_2^\gamma (x,P^2,Q^2)`$ itself. However, this is due largely to the fact that $`F_2^\gamma `$ of the real photon gets a large contribution from its VDM component, whereas the parameterization of $`q_L`$ used in this comparison corresponds to purely pointlike expression (8). Compared to the difference (13) of the pointlike parts of $`\gamma _T^{}`$ only, denoted by dashed curves in Fig. 2, the contributions of $`\gamma _L^{}`$ are, at least for $`m_q0.3`$ GeV, again quite significant throughout large part of the kinematical range considered.

### 3.3 LO calculations of dijet production in ep and e+e- collisions

The measurement of dijet production in ep and e+e- collisions provides another way of investigating interactions of the virtual photon . In general the cross–sections for dijet production are given as sums of contributions of all possible parton level subprocess. To demonstrate the importance of including the contributions of target $`\gamma _L^{}`$ it is, however, sufficient to use the approximation of the single effective subprocess  in which dijet cross–sections are expressed in terms of the so called effective parton distribution function of the target photon

$$D_{\mathrm{e}ff}(x,P^2,M^2)\underset{i=1}{\overset{n_f}{}}\left(q_i(x,P^2,M^2)+\overline{q}_i(x,P^2,M^2)\right)+\frac{9}{4}G(x,P^2,M^2).$$
(14)
In Fig. 3 we perform for this quantity the same comparisons as we did in Fig. 2 for $`F_2^\gamma `$, including the comparison with the difference $`\mathrm{\Delta }D_{\mathrm{eff}}(x,P^2,Q^2)`$, defined analogously to (13).

The fact that in QED $`\gamma _L^{}`$ contains no gluons is reflected in substantially smaller relative importance of $`\gamma _L^{}`$ for $`D_{\mathrm{eff}}`$ at small values of $`x`$. Otherwise, however, the messages of Figs. 2 and 3 are the same: in hard processes the relative importance of the contributions of target $`\gamma _L^{}`$ with respect to those of $`\gamma _T^{}`$

* depends sensitively on the value of $`m_q`$,
* peaks around $`x0.6`$ and vanish for $`x0`$ and $`x1`$,
* grows with target photon virtuality $`P^2`$ and
* decreases with factorization scale $`M^2`$.

For physically reasonable value $`m_q=0.3`$ GeV, Figs. 2 and 3 suggest that at least in part of the kinematical range accessible at HERA $`\gamma _L^{}`$ should definitely be taken into account.

### 3.4 NLO calculations of dijet production in ep collisions

In the preceding subsections we have discussed the importance of including the contributions of $`\gamma _L^{}`$ to QED or the LO QCD quantities $`F_{\mathrm{eff}},F_2^\gamma `$ and $`D_{\mathrm{eff}}`$. In this subsection we shall address the same question within the NLO QCD parton level calculations of dijet cross–sections in ep collisions, obtained with JETVIP , currently the only NLO parton level MC program that includes both direct and resolved photon contributions 999In specifying the powers of $`\alpha `$ corresponding to various diagrams we discard one common power of $`\alpha `$ coming from the vertex where the incoming electron emits the virtual photon. This vertex is also left out in diagrams of Fig. 4. JETVIP contains the full set of partonic cross–sections for the direct photon contribution up the order $`\alpha \alpha _s^2`$. Examples of such diagrams are in Fig. 4a ($`\alpha \alpha _s`$ tree diagram) and Fig. 4b ($`\alpha \alpha _s^2`$ tree diagram). To go one order of $`\alpha _s`$ higher and perform complete calculation of the direct photon contributions up to order $`\alpha \alpha _s^3`$ would require evaluating tree diagrams like that in Fig. 4e, as well as one–loop corrections to diagrams like in Fig. 4b and two–loop corrections to diagrams like in Fig. 4a. So far, such calculations are not available.

In addition to complete $`O(\alpha \alpha _s^2)`$ direct photon contribution JETVIP includes also the resolved photon one with partonic cross–sections up to the order $`\alpha _s^3`$, exemplified by diagrams in Fig. 4c,d. The justification for including in the resolved channel terms of the order $`\alpha _s^3`$ are discussed in detail in . Once the concept of virtual photon structure is introduced, part of the direct photon contribution (which for the virtual photon is actually nonsingular) is subtracted and included in the definition of PDF appearing in the resolved photon contribution. For $`\gamma _T^{}`$ the subtracted term is given as the convolution of the splitting function 101010JETVIP works with massless quarks and includes in (15) additional function of $`x`$.

$$q_T^{\mathrm{s}plit}(x,P^2,M^2)=q_T^{\mathrm{QED}}(x,P^2,Q^2)=\frac{\alpha }{2\pi }3e_q^2\left(x^2+(1x)^2\right)\mathrm{ln}\frac{M^2}{xP^2}.$$
(15)
with $`\alpha _s^2`$ partonic cross–sections. To avoid misunderstanding we shall henceforth use the term “direct unsubtracted” (DIRuns) to denote NLO direct photon contributions before this subtraction and reserve the term “direct” for the results after it. In this terminology the complete JETVIP calculations are given by the sum of direct and resolved parts and denoted DIR$`+`$RES. In JETVIP only the terms defining quark distribution function of the transverse virtual photon are subtracted from DIRuns calculations.

In  we discussed dijet cross–sections calculated by means of JETVIP in the kinematical region typical for HERA experiments

$$E_T^{(1)}E_T^c+\mathrm{\Delta },E_T^{(2)}E_T^c,E_T^c=5\mathrm{GeV},\mathrm{\Delta }=2\mathrm{GeV}$$

$$2.5\eta ^{(i)}0,i=1,2,$$

in four windows of photon virtuality

$$1.4P^22.4\mathrm{GeV}^2;2.4P^24.4\mathrm{GeV}^2;4.4P^210\mathrm{G}eV^2;10P^225\mathrm{GeV}^2$$

and for $`0.25y0.7`$. The whole analysis has been performed in $`\gamma ^{}`$p CMS. The cuts on $`E_T`$ were chosen in such a way that in all $`P^2`$ windows $`P^2E_T^2`$, thereby ensuring that the virtual photon lives long enough for its “structure” to develop before the hard scattering takes place. The asymmetric cut in $`E_T`$ is appropriate for our decision to plot the sums of $`E_T`$ and $`\eta `$ distributions of the jets with highest and second highest $`E_T`$. In JETVIP jets are defined by means of the standard cone algorithm with jet momenta defined using the $`E_T`$–weighting recombination procedure and supplemented with the $`R_{\mathrm{sep}}`$ parameter. All calculations presented below correspond to $`R_{\mathrm{sep}}=2`$ and were obtained setting the renormalization scale $`\mu `$ as well as the factorization scale $`M`$ equal to jet transverse energy. The sensitivity to these parameters as well as other ambiguities are discussed in detail in .

Beside the splitting term (15), which generates quark distribution function of $`\gamma _T^{}`$, one can subtract from NLO direct photon calculations also the integral over the term proportional to $`h_L(x)`$ and put it into the definition of quark distribution function of $`\gamma _L^{}`$. To do that properly would, however, require modifying the original code in order to take into account different $`y`$ dependence of the fluxes of $`\gamma _T^{}`$ and $`\gamma _L^{}`$ in (2-3). In this exploratory study we neglect this difference and fake the contributions of $`\gamma _L^{}`$ simply by running JETVIP in the resolved photon channel using (8) with $`m_q=0`$ as the input PDF. As in the considered region $`y0.4`$, the error incurred by this approximation does not exceed $`16`$%.

But does it make any sense to introduce the concept of PDF of $`\gamma _L^{}`$? Admittedly, for interactions of virtual photons we can stay solely within the framework of DIRuns calculations and thus dispense with the concept of PDF of virtual photons at all. On the other hand, as argued in , the effects incorporated in the transverse part of resolved photon component of JETVIP are numerically large. In particular, we have emphasized the importance of including in the resolved photon component of JETVIP the $`\alpha _s^3`$ partonic cross–sections. These are not included in exact $`\alpha \alpha _s^2`$ DIRuns calculations and in part of accessible kinematical range more than double the resolved photon contribution to dijet production at HERA compared to the contribution of the $`\alpha _s^2`$ partonic cross–sections.

The same effect can be expected for $`\gamma _L^{}`$ as the NLO DIRuns calculations contain at the order $`\alpha \alpha _s^2`$ exact matrix elements which include both transverse and longitudinal polarization of the target photon. To illustrate the importance of including the effects of $`\gamma _L^{}`$ we compare in Fig. 5 the convolutions $`q_L^{\mathrm{QED}}\sigma (\alpha _s^2)`$ and $`q_L^{\mathrm{QED}}\sigma (\alpha _s^3)`$ with the convolution $`q_T^{\mathrm{QED}}\sigma (\alpha _s^2)`$. In addition, we overlay the complete NLO DIR$`+`$RES and DIRuns results, from which, however, the LO direct photon contribution has been subtracted. Fig. 5 shows that the contributions of $`\gamma _L^{}`$, though smaller, are nevertheless comparable to those of $`\gamma _T^{}`$, in particular for $`\eta `$ close to $`\eta 0`$. Moreover, in the region $`\eta 1.75`$ the sum of the contributions $`(q_T^{\mathrm{QED}}+q_L^{\mathrm{QED}})\sigma (\alpha _s^2)`$ approximates remarkably well the exact $`\alpha \alpha _s^2`$ DIRuns calculations. The excess of the exact results over this sum in the region $`\eta 1.75`$ is primarily due to the fact that the $`\alpha \alpha _s^2`$ DIRuns calculations contain beside the tree level diagrams describing the production of three final state partons, also one loop corrections to two parton final states, which contribute predominantly at large negative $`\eta `$.

The message of Fig. 5 is quantified by plotting in Fig. 6 the ratia

$$r_k(\eta ,P^2)\frac{q_k^{\mathrm{QED}}\sigma ^{\mathrm{res}}(\alpha _s^2)}{\sigma _{\mathrm{uns}}^{\mathrm{DIR}}(\alpha \alpha _s^2)},k=T,L$$
(16)
of the contributions of $`\gamma _T^{}`$ and $`\gamma _L^{}`$, as well as their sum, to the $`\alpha _s^2`$ part of the DIRuns results. The ratio of the contributions of $`\gamma _L^{}`$ and $`\gamma _T^{}`$ is above $`1/4`$ throughout the considered $`\eta `$ range and above $`1/2`$ in the region $`\eta 0`$. Within the DIRuns calculations at the order $`\alpha \alpha _s^2`$ in the kinematical region relevant for HERA, $`\gamma _L^{}`$ is thus comparable in importance to $`\gamma _T^{}`$.

The preceding discussion illustrates the importance of the contributions of $`\gamma _L^{}`$, but as the $`\alpha \alpha _s^2`$ DIRuns calculations include them exactly, the genuine nontrivial effect of introducing the concept of PDF of $`\gamma _L^{}`$ is given in Fig. 5 by the thin dashed-dotted curves, denoting the convolutions $`q_L^{\mathrm{QED}}\sigma (\alpha _s^3)`$, which, similarly to those of $`\gamma _T^{}`$, are not included in NLO DIRuns calculations. However, as the current version of JETVIP takes into account in the resolved channel only the transverse virtual photons, they are not included even in the full DIR$`+`$RES calculations. The net nontrivial effect of introducing the concept of PDF of $`\gamma _L^{}`$ into JETVIP is then quantified by plotting in Fig. 7a the ratio

$$r_{\mathrm{NLO}}(\eta ,P^2)\frac{q_L^{\mathrm{QED}}\sigma ^{\mathrm{res}}(\alpha _s^3)}{\sigma ^{\mathrm{DIR}+\mathrm{RES}}}.$$
(17)
Also by this measure the contributions of $`\gamma _L^{}`$ are sizable. This net effect is much larger when the convolution $`q_L^{\mathrm{QED}}\sigma ^{\mathrm{res}}(\alpha _s^3)`$ is compared to the difference of DIR$`+`$RES and DIRuns JETVIP results, measuring the nontrivial aspects of the concept of PDF of $`\gamma _T^{}`$ and corresponding to the gap between the thick solid and dashed curves in Fig. 5. As shown in Fig. 7b, the corresponding ratio, denoted $`r_{\mathrm{nontriv}}(\eta ,P^2)`$, is large, particularly for $`\eta `$ close to lower edge $`\eta =2.5`$.

## 4 How to measure partonic content of $`\gamma _L^{}`$

In principle there is no obstacle to extracting partonic content of the virtual photon from experimental data by analyzing dijet production at two different values of $`y`$. This procedure is analogous to that involved in measuring the longitudinal structure function $`F_L^\mathrm{p}(x,Q^2)`$ of real hadrons, which requires performing the measurement at two different collisions energies. Although straightforward in principle, no such direct measurement of $`F_L^\mathrm{p}`$ has been performed at HERA, primarily for technical reasons related to changing the proton energy. For extraction of the partonic content of $`\gamma _L^{}`$ no such change of beam energies is necessary and it suffices to perform the analysis of dijet cross–sections at two different values of $`y`$. In practice, however, the separation of the contributions of $`\gamma _T^{}`$ and $`\gamma _L^{}`$ is not that simple, because it relies on different $`y`$ dependencies of the corresponding fluxes (2-3) at large $`y`$. This in turn requires measuring jet cross–sections in narrow bins centered at two different values $`y_1`$ and $`y_2`$ instead of integrating over the whole interval of accessible $`y`$, which at HERA spans typically $`0.05y0.9`$. Optimizing the bin width and choice of the values $`y_1,y_2`$ is crucial for the success of such extraction.

## 5 Summary and conclusions

We have demonstrated the importance of including in hard collisions the contributions of the longitudinal polarization of the target virtual photon. In QED these contributions are fully calculable and their onset is determined by the ratio $`P^2/m^2`$ of photon virtuality $`P^2`$ and fermion mass $`m^2`$. The inclusion of target $`\gamma _L^{}`$ is indispensable for good quantitative agreement of QED predictions with existing LEP data.

In QCD gluon radiation off the quarks or antiquarks coupling to $`\gamma _L^{}`$ is expected to modify simple QED formulae and, in addition, generate gluons inside $`\gamma _L^{}`$. In this exploratory study we, nevertheless, neglected these effects and used the purely QED formula for quark distribution function of $`\gamma _L^{}`$. The numerical relevance of $`\gamma _L^{}`$ has been illustrated within the framework of LO analysis of observables $`F_2^\gamma `$ and $`D_{\mathrm{eff}}`$ as well as within the NLO calculations of dijet production at HERA. Better theoretical understanding of the structure of $`\gamma _L^{}`$ is, however, needed for more reliable evaluation of these effects.

Acknowledgment: We are grateful to J. Cvach, Ch. Friberg, B. Pötter, I. Schienbein and A. Valkárová for interesting discussions concerning the structure and interactions of longitudinal virtual photons. This work was supported in part by the Grant Agency of the Academy od Sciences of the Czech Republic under grants No. A1010821 and B1010005.


text:
# Roughening and superroughening in the ordered and random two-dimensional sine-Gordon models

## I Introduction

The location and characteristics of phase transitions in disordered media constitute a long-standing and controversial question, particularly in more than one spatial dimension Plischke . The question becomes even more difficult if the disorder is not very weak; then, new, non-trivial behavior is commonly found, involving features such as aging, ergodicity breaking, extremely slow dynamics, complicated energy landscapes, etc.; major examples of this are spin glasses and structural glasses Plischke ; Young . In this context, the properties of crystal surfaces growing on disordered substrates, frequently described by a two-dimensional random-phase sine-Gordon model (RsGM), have attracted a lot of attention in the last decade Toner ; Tsai ; Batrouni ; Riegercom ; Cule ; Cule2 ; Krug ; Scheidl ; Marinari ; Lancaster ; Zeng ; Tsai2 ; Rieger1 ; Coluzzi ; us1 ; jjfises ; Rieger2 ; Raul ; Hwanew ; Rieger3 . Interestingly, the same model describes many other relevant physical problems, such as randomly pinned flux lines confined in a plane Batrouni ; Hwanew ; Fisher ; Hwa2 ; Blatter , vortex lines in planar Josephson junctions Vinokur , charge density waves Fukuyama , and commensurate-incommensurate transitions Villain .

In spite of those efforts, the phase diagram and main features of the RsGM are not clear yet. To summarize what is known, we refer for comparison to the ordered sine-Gordon model (OsGM), which is rather well understood (see, e.g., Plischke ; Weeks ; Beijeren ; Barabasi ; Krug2 ; Villain2  and references therein). The hamiltonian for the OsGM and the RsGM is

$$\begin{array}{cc}\hfill H=\underset{𝒓}{}[& \frac{1}{2}\underset{n.n.}{}[h(𝒓)𝒉(𝒓^{})]^\mathit{2}+\hfill \\ & +V_0(1\mathrm{cos}[h(𝒓)𝒉^{(\mathit{0})}(𝒓)])].\hfill \end{array}$$
(1)
where $`n.n.`$ stands for the nearest neighbors of site $`𝒓`$. The OsGM corresponds to $`h^{(0)}(𝒓)\mathit{0}`$, whereas the RsGM is defined by choosing the quenched disorder variables $`h^{(0)}(𝒓)`$ randomly from a uniform distribution in $`[0,2\pi ]`$. We will discuss our work in terms of growth on flat (OsGM) or disordered (RsGM) substrates (see Hwa2 ; Blatter ; Vinokur ; Fukuyama ; Villain  for other physical interpretations): Accordingly, $`h(𝒓)`$ is a continuous variable representing the height of the growing surface at site $`𝒓`$ of the lattice, and the cosine term is a potential energy making integer (i.e., multiples of the crystalline lattice constant) values of the height energetically favorable. We will consider two-dimensional (2D) square lattices, i.e., $`𝒓=(𝒓_\mathit{1},𝒓_\mathit{2})`$, with $`N=L\times L`$ sites. As first shown by Chui and Weeks Chui  (see also Knops ) by means of a renormalization group (RG) approach Golden , the OsGM posseses a Kosterlitz-Thouless KT  type topological transition between a low temperature, flat phase and a high temperature, rough phase, the latter being described by the Edwards-Wilkinson (EW) equation EW , i.e., the diffusion equation with additive white noise (see below). Above the so called roughening temperature ($`T_R`$), thermal fluctuations effectively suppress the effect of the cosine potential, and the surface becomes free, described only by the kinetic part of the hamiltonian (1). As will be shown below, the most important measurable consequence of this is that the width of the interface,

$$w^2=\frac{1}{N}\underset{𝒓}{}(h(𝒓)\overline{𝒉})^\mathit{2},\overline{𝒉}=\frac{\mathit{1}}{𝑵}\underset{𝒓}{}𝒉(𝒓),$$
(2)
scales (in 2D) as $`w^2\mathrm{ln}L`$ in the asymptotic regime.

In contrast with the clear picture for the OsGM, there are very few generally accepted results for the RsGM. One of them is that there must be a roughening temperature above which the potential effectively vanishes (much as in the case of the OsGM) leading to EW behavior. Apart from this, theoretical predictions about the low temperature phase largely disagree (a good summary is given in the third paper in Ref. Cule ): While RG calculations predict a superrough low temperature phase, with $`w^2\mathrm{ln}^2L`$, replica-symmetry breaking variational approaches lead to $`w^2\mathrm{ln}L`$ independently of temperature. Numerical simulations were not very conclusive either: Batrouni and Hwa Batrouni  did not find evidence for an equilibrium phase transition in Langevin dynamics, although Monte Carlo simulations by Rieger Riegercom  showed a transition from a superrough phase to a EW phase for stronger potentials \[larger $`V_0`$ in Eq. (1)\] than those used by Batrouni and Hwa. Subsequent numerical work Marinari ; Lancaster ; Coluzzi ; jjfises  presented more evidence of superrough ($`\mathrm{ln}^2L`$) behavior, albeit with large quantitative discrepancies with the predictions of RG theories. Finally, a number of works using special optimization algorithms Zeng ; Rieger1 ; Rieger2  or direct numerical simulations us1  strongly supported superrough behavior at zero temperature. Very recently Hwanew , simulations of a related model provided more evidence of $`\mathrm{ln}^2`$ behavior at finite temperatures, although this model did not allow study of the transition. In summary, most researchers believe that there is indeed a superrough low temperature phase in the RsGM, but its nature (glassy or not), the transition temperature, and its dependence on the model parameters remain unclear.

In this paper, we attempt to shed light on these issues by simultaneously studying the OsGM and the RsGM in different regions of their phase diagram. As we will show below, it turns out that the potential strength, $`V_0`$, crucially determines the model features. In addition, it is also natural to ask about the intensity of the disorder: How does the model phenomenology change if the disorder takes values in $`[0,ϵ]`$ with $`ϵ<2\pi `$? The importance of these points can be clearly seen in the Langevin equation,

$$\frac{dh(𝒓,𝒕)}{dt}=\frac{\delta H}{\delta h(𝒓,𝒕)}+\eta (𝒓,𝒕),$$
(3)
where $`\eta (𝒓,𝒕)`$ is a gaussian white noise of zero mean and correlations obeying the fluctuation-dissipation theorem,

$$\eta (𝒓^{},𝒕^{})\eta (𝒓,𝒕)=\mathit{2}𝑻\delta ^{(\mathit{2})}(𝒓𝒓^{})\delta (𝒕𝒕^{}),$$
(4)
where $`T`$ stands for the temperature of the system, $`\mathrm{}`$ indicates thermal averages (over $`\eta `$) and Boltzmann’s constant is set to $`k_B=1`$. For our hamiltonian, Eq. (1), the Langevin equation is

$$\begin{array}{cc}\hfill \frac{h(𝒓,𝒕)}{t}=& ^2h(𝒓,𝒕)\hfill \\ & V_0\mathrm{sin}[h(𝒓,𝒕)𝒉^{(\mathit{0})}(𝒓)]+\eta (𝒓,𝒕),\hfill \end{array}$$
(5)
and changing variables according to $`\stackrel{~}{h}(𝒓,𝒕)=𝒉(𝒓,𝒕)𝒉^{(\mathit{0})}(𝒓)`$ (i.e., the height referred to the substrate), we find

$$\begin{array}{cc}\hfill \frac{\stackrel{~}{h}(𝒓,𝒕)}{t}=& ϵ\stackrel{~}{F}(𝒓)+^\mathit{2}\stackrel{~}{𝒉}(𝒓,𝒕)\hfill \\ & V_0\mathrm{sin}(\stackrel{~}{h}(𝒓,𝒕))+\eta (𝒓,𝒕)\hfill \end{array}$$
(6)
where $`F(r)=^2h^{(0)}(𝒓)/ϵ`$. In this form, the disorder is a random (correlated in space) chemical potential acting on a surface growing on a flat substrate. If we think of the roughening transition for the OsGM in terms of the interplay between the temperature $`T`$ and the energy scale introduced by the potential, $`V_0`$, when $`ϵ0`$, we have another energy scale, $`ϵ`$, which can modify the universal features of the roughening transition or even give rise to novel thermodynamical transitions.

Having the above issues in mind, we discuss our work according to the following scheme: Section II presents an analytical study of the EW equation and other statistical mechanics results about the energy and roughness. Section III deals with the main part of our work, namely Langevin dynamics simulations of the OsGM and the RsGM, beginning with $`V_0=1`$ and $`ϵ=2\pi `$ and subsequently analyzing the model behavior for different $`V_0`$ and $`ϵ`$. Finally, in Sec. IV we present and discuss our conclusions and indicate future directions.

## II Analytical results

### II.1 Linear theory: The Edwards-Wilkinson equation

According to the RG approach Weeks , the high temperature phase of the OsGM obeys the EW equation EW ,

$$\frac{h(𝒓,𝒕)}{t}=\mathbf{}^2h(𝒓,𝒕)+\eta (𝒓,𝒕).$$
(7)
Equation (7) can be solved by means of the Fourier decomposition (see Krug2  for details),

$$\widehat{h}_𝒒=\frac{1}{L}\underset{𝒓}{}\mathrm{e}^{i𝒒𝒓}h(𝒓,t),$$
(8)
where $`𝒒=\frac{\mathit{2}\pi }{𝑳}𝒌`$, $`k_i=0,\mathrm{},L1`$ is the reciprocal vector. The structure factor can then be shown to be

$$S(𝒒)=\widehat{𝒉}_𝒒\widehat{𝒉}_𝒒=𝑻\frac{\mathit{1}\mathrm{e}^{\mathit{2}\omega _𝒒𝒕}}{\omega _𝒒},$$
(9)
$`\omega _𝒒`$ being the 2D EW discrete dispersion relation

$$\omega _𝒒=4\mathrm{sin}^2\left(\frac{q_1}{2}\right)+4\mathrm{sin}^2\left(\frac{q_2}{2}\right).$$
(10)
From $`S(𝒒)`$ we can obtain the relevant magnitudes, such as the total roughness

$$w^2(t)=\frac{1}{L^d}\underset{𝒓}{}[h(𝒓,𝒕)\overline{𝒉}]^\mathit{2}=\frac{1}{L^d}\underset{𝒒\mathit{0}}{}S(𝒒),$$
(11)
the correlation function,

$`C(𝒓)`$ $`=`$ $`{\displaystyle \frac{1}{L^2}}{\displaystyle \underset{𝒔}{}}[h(𝒔+𝒓)𝒉(𝒔)]^\mathit{2}`$ (12)
$`=`$ $`{\displaystyle \frac{2}{L^2}}{\displaystyle \underset{𝒒}{}}S(𝒒)[\mathit{1}\mathrm{cos}(𝒒𝒓)].`$
the total slope,

$$s^2(t)=\frac{1}{L^2}\underset{𝒒\mathit{0}}{}S(𝒒)\omega _𝒒,$$
(13)
and the slope difference correlation function,

$$\begin{array}{cc}\hfill G(𝒓)& =\frac{1}{L^2}\underset{𝒔}{}[\mathbf{}h(𝒔+𝒓)\mathbf{}𝒉(𝒔)]^\mathit{2}=\hfill \\ & =\frac{4}{L^2}\underset{𝒒}{}S(𝒒)\underset{𝒊=\mathit{1}}{\overset{\mathit{2}}{}}[\mathit{1}\mathrm{cos}(𝒒_𝒊𝒓_𝒊)][\mathit{1}\mathrm{cos}(𝒒_𝒊)].\hfill \end{array}$$
(14)
From Eq. (9) we can find an estimate for the time needed to reach saturation, $`t_\times `$, and the dynamic exponent, $`z`$: We compute the time that the structure factor needs to be within a 1% of its saturated form for the slowest Fourier mode, that with the lowest $`|𝒒|`$, $`|𝒒|=\mathit{2}\pi /𝑳`$:

$$t_\times 3\times 10^2L^2,$$
(15)
implying $`z=2`$. For the saturated roughness, we obtain

$$w^2(t\mathrm{},L)=\frac{T}{4L^2}\underset{k_1,k_2=1}{\overset{L1}{}}\left[\mathrm{sin}^2\left(\frac{q_1}{2}\right)+\mathrm{sin}^2\left(\frac{q_2}{2}\right)\right]^1,$$
(16)
which cannot be computed exactly but, for large $`L`$, can be approximated changing the sum by an integral and the sine functions by their arguments, arriving at

$`w^2(t\mathrm{},L)`$ $``$ $`{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle _{2\pi /L}^\pi }{\displaystyle _{2\pi /L}^\pi }𝑑q_x𝑑q_y{\displaystyle \frac{T}{(q_x^2+q_y^2)}}`$ (17)
$``$ $`{\displaystyle \frac{T}{2\pi }}{\displaystyle _{2\pi /L}^\pi }{\displaystyle \frac{dq}{q}}={\displaystyle \frac{T}{2\pi }}\mathrm{ln}L,`$
yielding a roughness exponent $`\alpha =0`$. As for the total slope, it tends to a value independent of $`L`$,

$$s^2T,$$
(18)
whereas for the correlation functions, we find for large $`r`$

$`C(𝐫)`$ $``$ $`{\displaystyle \frac{T}{\pi }}\mathrm{ln}r,`$ (19)
$`G(𝐫)`$ $``$ $`2s^2+{\displaystyle \frac{T}{4\pi }}\mathrm{ln}{\displaystyle \frac{r^8}{[(r^2+1)^24r_1^2][(r^2+1)^24r_2^2]}}`$ (20)
$``$ $`2T.`$
### II.2 Other results for the energy and the roughness

At equilibrium, the partition function for the OsGM is

$$Z=\left[\underset{𝒓}{}dh(𝒓)\right]𝒆^{\beta [\frac{\mathit{1}}{\mathit{2}}_𝒓(𝒉(𝒓))^\mathit{2}+𝑽_\mathit{0}(\mathit{1}\mathrm{cos}𝒉(𝒓))]},$$
(21)
where $`\beta =T^1`$. Expanding for high temperatures berez  means rewriting Eq. (21) as a series in powers of $`\beta V_0`$:

$`Z`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\beta V_0)^n}{n!}}{\displaystyle }\left[{\displaystyle \underset{𝒓}{}}dh(𝒓)\right]\mathrm{e}^{\beta _𝒓\frac{\mathit{1}}{\mathit{2}}(𝒉(𝒓))^\mathit{2}}\times `$
$`\times \left({\displaystyle \underset{𝒓}{}}\mathrm{cos}h(𝒓)\right)^𝒏=`$
$`=`$ $`Z_0{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\beta V_0)^n}{n!}}\left({\displaystyle \underset{𝒓}{}}\mathrm{cos}h(𝒓)\right)^𝒏_{𝑯_\mathit{0}},`$ (22)
where $`Z_0`$ is the partition function of $`H_0`$, the free hamiltonian \[i.e., Eq. (1) without the potential term\]. By means of this expansion we obtain

$$\begin{array}{cc}\hfill \mathrm{cos}h(𝒓^{})_𝑯=& \frac{Z_0}{Z}\underset{n=0}{\overset{\mathrm{}}{}}\frac{(\beta V_0)^n}{n!}\times \hfill \\ & \times \left(\underset{𝒓}{}\mathrm{cos}h(𝒓)\right)^𝒏\mathrm{cos}𝒉(𝒓^{})_{𝑯_\mathit{0}},\hfill \end{array}$$
(23)
where subdominant terms such as $`\beta ^n\mathrm{exp}(A/\beta )`$ have been neglected, and only terms of the order $`\beta ^n`$ have been kept. The expression above can be put in the form

$$\begin{array}{cc}\hfill \mathrm{cos}h(𝒓^{})_𝑯& =\left(\underset{n=1}{\overset{\mathrm{}}{}}\frac{2^{2n1}(n!)^2}{n(\beta V_0)^{2n1}}\right)\times \hfill \\ & \times \left(\underset{n=0}{\overset{\mathrm{}}{}}\frac{(\beta V_0)^{2n}}{2^{2n}(n!)^2}\right)^1=\hfill \\ \hfill =\frac{\beta V_0}{2}& \frac{(\beta V_0)^3}{16}+\frac{(\beta V_0)^5}{96}+O((\beta V_0)^7),\hfill \end{array}$$
(24)
yielding for the approximate energy per site

$$\begin{array}{cc}\hfill E=& \frac{1}{L^2}H=\frac{1}{2\beta }+\hfill \\ & +V_0\left(1\frac{\beta V_0}{2}+\frac{(\beta V_0)^3}{16}+O((\beta V_0)^5)\right).\hfill \end{array}$$
(25)
For the roughness, we have to compute

$$\begin{array}{cc}\hfill (h(𝒓)\overline{𝒉})^\mathit{2}_𝑯& =\frac{Z_0}{Z}\times \hfill \\ \hfill \times \underset{n=0}{\overset{\mathrm{}}{}}\frac{(\beta V_0)^n}{n!}& h^2(𝒓)\underset{𝒓^{}}{}\left[\mathrm{cos}𝒉(𝒓^{})\right]^𝒏_{H_0},\hfill \end{array}$$
(26)
assuming that at equilibrium $`\overline{h}=0`$. Neglecting again subdominant terms, we find

$$h^2(𝒓)_𝑯=\frac{𝒁_\mathit{0}}{𝒁}𝒉^\mathit{2}(𝒓)_{𝑯_\mathit{0}},$$
(27)
and hence

$$\begin{array}{cc}\hfill w_H^2=w_{H_0}^2(1& \frac{\beta V_0}{4}+\frac{3(\beta V_0)^2}{64}\hfill \\ \hfill & \frac{19(\beta V_0)^3}{2304}+O((\beta V_0)^4)).\hfill \end{array}$$
(28)
Finally, at low temperatures the height exhibits only small deviations from $`h=0`$, and therefore we can approximate the hamiltonian as

$$\begin{array}{cc}\hfill H& =\underset{i}{}\frac{1}{2}(h(𝒓))^\mathit{2}+𝑽_\mathit{0}[\mathit{1}\mathrm{cos}(𝒉(𝒓))]\hfill \\ & \underset{i}{}\frac{1}{2}(h(𝒓))^\mathit{2}+𝑽_\mathit{0}\frac{𝒉(𝒓)^\mathit{2}}{\mathit{2}},\hfill \end{array}$$
(29)
i.e., $`2N`$ quadratic terms, each one of which, according to the equipartition theorem Huang , contributes with $`T/2`$ to the energy value. Taking into account the global factor $`1/2`$, we conclude that at low temperatures the energy of the OsGM is approximately $`E=T/2`$.

## III Langevin dynamics results

### III.1 Numerical simulation details

We have integrated the Langevin equation (5) corresponding to the hamiltonian (1), on $`L\times L`$ square lattices with periodic boundary conditions, using a stochastic second order Runge-Kutta method simu ; in some cases, we have repeated the simulations with a Heun method maxi , with excellent agreement between both procedures. We are therefore sure that our results are not an artifact of our numerical method, a conclusion further reinforced by the agreement with the theoretical expectations of the previous section as we discuss below. The simulations reported in this paper were carried out with a time step $`\mathrm{\Delta }t=0.01`$ on lattices of sizes $`L=64,128`$ and 256. It is important to stress that we did not perform averages over the quenched noise in the RsGM; however, we checked that the outcome of our simulations did not depend strongly on the realization of the quenched noise or the initial conditions (flat, $`h(𝒓)=\mathit{0}`$; as the substrate, $`h(𝒓)=𝒉^{(\mathit{0})}(𝒓)`$, or random) by repeating several times a number of our simulations. In all cases, simulations consisted of an equilibration time and a measuring period. Eq. (15) predicts that the time needed for equilibration is $`t_\times 500`$ for $`L=128`$ and $`t_\times 2000`$ for $`L=256`$, and hence we used equilibration times of $`5000`$ and $`10000`$ units, respectively; afterwards, we let the system evolve for an equal period, over which we performed thermal averages. Equilibration was ensured in all cases by verifying that the fluctuations of the energy were gaussian and by checking the equality of the specific heat computed from those fluctuations and from derivatives of the mean energy Plischke , as well as by monitoring the evolution in time of the quantities of interest toward a stationary state. As an additional test, we compared the imposed simulation temperature, arising from the noise term, to that measured during the evolution according to the equipartition theorem Huang ; both quantities were always found to agree within a 0.1%. Finally, we did a few very longer runs, whose outcome agreed with that of the shorter runs.

### III.2 Standard RsGM: $`V_0=1`$ and $`ϵ=2\pi `$

We begin by discussing our results for the “canonical” version of both the OsGM and the RsGM, i.e., hamiltonian (1) with $`V_0=1`$, as studied (for the ordered case) in falo ; yo1 . In those works, the roughening temperature was determined by a direct comparison to RG predictions, looking for the temperature at which the height difference correlation function reached a universal (in the RG framework) value, with the result that for the OsGM $`T_R25`$ in our dimensionless units. Remarkably, this is the RG value for $`T_R`$, which makes very tempting to claim that this method indeed yields $`T_R`$ correctly. Monte Carlo simulations of the discrete gaussian model (hamiltonian (1) with $`V_0=0`$ and $`h(𝒓)`$ restricted to integer values) by Shugard et al. Shugard  with the same criterion for locating the transition yielded similar results. However, as RG calculations are perturbative in $`V_0`$ and carried out on the continuum Langevin equation Weeks ; Beijeren ; Barabasi ; Villain2 ; Chui ; Knops , it is not obvious that they apply to a discrete model with $`V_0=1`$, i.e., of the same order as the kinetic term. In view of this, we decided to include the OsGM in this study, both to analyze in detail whether the comparison to the universal RG prediction for the factor is a good tool to find $`T_R`$ and to compare its high temperature phase with that of the RsGM, which should also be of EW type.

The first quantity we discuss, shown in Fig. 1, is the mean energy of both models. As we see, the results are largely independent of the system size, and hence it is unlikely that they are affected by finite size effects. The plot shows that the mean energy of the OsGM reaches the high temperature approximation at $`T_0=16\pm 1`$; on the other hand, the mean energy for the RsGM is never too far from it, although for temperatures lower than $`T_1=4\pm 1`$ the numerical values lie slightly below the high temperature result. At temperatures higher than $`T_0`$ the energies of both models coincide within the accuracy of our simulations. These results suggest that $`T_0`$ could be the roughening temperature, $`T_R`$, for the OsGM and $`T_1`$ the superroughening temperature, $`T_{SR}`$, for the RsGM, because the EW behavior of the mean energy of both models indicates the effective suppresion of the sine term by temperature. The inset in Fig. 1 presents the specific heat, $`C_v`$, of both models, exhibiting a well defined peak in $`C_v`$ for the OsGM with its maximum at temperature $`T=9`$, much lower than $`T_0`$. In contrast, we do not observe any peak for the RsGM; there might be a peak at $`T=3`$, but the evidence is not conclusive. Concerning the peak for the OsGM, we stress again the absence of any finite size effect, consistent with a KT type transition. Furthermore, we believe that this peak is a (Schottky) anomaly Golden  similar to that observed in 2D XY and easy-plane Heisenberg spin models Kawabata  above the KT transition. Recall that when mapping the OsGM to the XY model the temperature of the former maps to the inverse temperature of the latter Chui , hence the observation of the anomaly below the possible transition temperature $`T_0`$. This reinforces our interpretation of $`T_0`$ as the roughening temperature of the OsGM.

The total roughness of both models, shown in Fig. 2, behaves similarly for temperatures higher than $`T_0`$, depending linearly on temperature. In both cases, we see that the slope of the roughness depends on the system size, as predicted by Eq. (16). Whereas the approximation in Eq. (17) expression yields slopes 0.66 and 0.77 for $`L=64`$ and 128, respectively, if we numerically compute the exact result, Eq. (11), the slopes turn out to be 0.71 y 0.82, in excellent agreement with the results of our simulations, 0.71 and 0.83. Below $`T_0`$, the roughness for the OsGM is independent of the size and depends nonlinearly on temperature, whereas above $`T_0`$ we find linear dependence on temperature and clear finite size effects. For the RsGM, the linear behavior extends all the way down to $`T_1`$, and below $`T_1`$ the behavior becomes nonlinear. The slope of the linear region is approximately the same in the first part, from $`T_1`$ to $`T_0`$, and in the second, above $`T_0`$, i.e., the whole linear region is well described by the EW model. This means that above $`T_0`$ the linear model describes accurately the behavior of the OsGM, and hence from now on we identify $`T_0`$ with the roughening temperature $`T_R`$, whereas for the RsGM, the same is true of $`T_1`$ and $`T_{SR}`$. Figure 2 presents also results for the roughness susceptibility, $`\chi _w`$, defined as $`\chi _w^2=[(w^2)^2w^2^2]/T^2`$. For the OsGM, $`\chi _w`$ exhibits a very clear peak at $`T_R`$, and above $`T_R`$ it is the same as for the RsGM; however, this magnitude is very noisy and these results must be taken with caution. In fact, one could identify a peak for the RsGM at $`T_{SR}`$, but different realizations lead to different results, in contrast with the peak for the OsGM, which is the same for all realizations. Below $`T_R`$, the values of $`\chi _w`$ for the OsGM are independent of the system size, whereas above $`T_R`$ they increase with size without any definite scaling.

Figure 3 depicts the height-difference correlation function for the two studied models, and shows that above $`T_R=T_0`$ and $`T_{SR}=T_1`$ they behave as predicted by the linear theory: The slope of the numerical height-difference correlation function is $`0.32`$, indistinguishable from the predicted $`0.318`$ by Eq. (19). In addition, the correlation functions for the OsGM and the RsGM coincide, as shown in the plot for the RsGM. We see that the behavior of the correlation functions is in full agreement with our claims regarding $`T_R`$ and $`T_{SR}`$, and this is further confirmed by the plot of the slope correlation function in Fig. 4. It is important to note that, below $`T_{SR}`$, the behavior of the height difference correlation function for the RsGM is approximately a squared logarithm, as predicted by RG calculations. We postpone discussion of this point to the next section.

Finally we studied another magnitude, namely

$$m_1=\mathrm{cos}\left[h(𝒓)𝒉^{(\mathit{0})}(𝒓)\right].$$
(30)
For the OsGM this is the average computed in the preceding section, whereas for the RsGM it is the average of the cosine of the height referred to the substrate. Figure 5 shows our results: The high temperature approximation, Eq. (24), agrees very accurately with the simulations for temperatures above $`T_R`$ (OsGM) and $`T_{SR}`$ (RsGM). The results for both models are again indistinguishable for temperatures above $`T_R`$. Interestingly, $`m_1`$ is largely independent not only of the system size, but also on the realization of the quenched disorder for the RsGM.

### III.3 Other potential strengths

We now turn to the question of the influence of the potential strength on the RsGM behavior. We have considered two representative values, $`V_0=0.2`$ and $`V_0=5`$, i.e., five times smaller and larger, respectively, than the “canonical” value $`V_0=1`$. The smallest value is close to that considered in Batrouni , $`V_0=0.15`$, and we expect that our results will be comparable to theirs. As before, we begin by discussing the total roughness and the specific heat (see Fig. 6). First of all, for all values of $`V_0`$ there is a temperature above which the roughness value is independent of $`V_0`$ and of the presence or absence of disorder. This means that our identification of this regime with the effective suppresion of any potential effect is indeed correct: Different $`V_0`$ leads only to different transition temperatures. Thus, for the OsGM we find $`T_R^{V_0=0.2}=13\pm 1`$ and $`T_R^{V_0=5}=19\pm 1`$, in agreement with the intuitive expectation that larger potentials need higher temperatures to be suppressed. Aside from this, the general shape of the roughness curve is basically the same for the three values of $`V_0`$. The situation (for the OsGM) is the same as far as the specific heat is concerned: Larger (smaller) $`V_0`$ leads to larger (smaller) anomalies, which are displaced to higher (lower) $`T`$ following the corresponding $`T_R`$. Therefore, we conclude that changing $`V_0`$ does not introduce anything qualitatively new in the OsGM.

The picture for the RsGM is substantially different: Modifying $`V_0`$ does give rise to qualitatively new phenomena. Let us first look at the small $`V_0`$ case. Figure 6 shows that the roughness follows a straight line all the way down to $`T=0`$ (although we cannot exclude that there are nonlinear effects for $`T1`$ with our present resolution). This would suggest that there is no transition in this case, very much like the results of Batrouni and Hwa Batrouni . The upper panel of Fig. 7, where the height difference correlation function is depicted, confirms this interpretation, showing no dependence on temperature in the analyzed range; Fig. 8, for $`m_1`$, agrees with this as well, in so far as the dependence of $`m_1`$ on temperature is well described by the high temperature expansion. In view of this, we can conclude that if there is a transition, it occurs at a temperature smaller than $`T1`$.

Finally, let us consider the large $`V_0`$ case, $`V_0=5`$. The plot of the roughness in Fig. 6 exhibits a striking peak for $`T=5`$, after which the roughness decreases until reaching the high temperature regime (marked by the corresponding OsGM with the same value for $`V_0`$) at $`T=8\pm 1`$. To assess the relevance of this peak in the roughness, we repeated our simulations for $`L=64`$ and performed additional ones for $`L=128`$. The results, collected in Fig. 9, show that the peak is a realization dependent feature. However, in this plot we also see that for temperatures below $`T=5`$ the roughness is roughly independent of the system size, something which we did not observe when $`V_0=1`$ (the lines for $`V_0=1`$ are included in Fig. 9 again for comparison). Hence, even if the peak at $`T=5`$ does not actually exist, that temperature does seem to separate two different regions. In addition, for $`V_0=5`$ the specific heat has a (more smeared) maximum at about the same temperature as that of the roughness, although our data are much noisier and we cannot establish clearly the maximum temperature; the dependence, however, is manifestly non-monotonic. Figure 8 supports our conclusion that $`T_{SR}^{V_0=5}=8\pm 1`$, whereas nothing special is seen as $`m_1`$ goes through $`T=5`$, the roughness maximum. The most intriguing result is the one in Fig. 7 for the height difference correlation function: For $`T5`$, the scaled correlations decrease with temperature but, simultaneously, the correlation length increases, up to $`T=5`$, when it increases beyond the system size. Above that temperature, it follows the same evolution as the $`V_0=1`$ case, finally reaching $`T_{SR}^5=8\pm 1`$. Whereas in this intermediate temperature regime the height difference correlation functions are well described by squared logarithms, Fig. 7 immediately shows that the lowest temperature correlations can by no means be considered squared logarithms. This suggest the presence of a new phase transition at $`T^{}=5\pm 1`$. We will discuss the possible nature of the low temperature phase and the existence of this transition in the next section.

### III.4 Other disorder strengths

We now generalize the RsGM and let $`h^{(0)}(𝒓)`$ to be randomly chosen from a uniform distribution in $`[0,ϵ]`$, with $`0<ϵ<2\pi `$; $`ϵ=2\pi `$ is the case studied in the preceding Section. As representative examples we have considered $`ϵ=0.8`$ and $`ϵ=0.2`$, closer to the RsGM and the OsGM, respectively. From the roughness dependence on temperature, shown in Fig. 10, we see that for both values of $`ϵ`$ the dependence of the roughness on temperature is qualitatively similar to the OsGM, the results for the lower $`ϵ`$ value being practically the same as for the $`ϵ=0`$ case. However, the case $`ϵ=0.8`$ is somewhat different: The low temperature region appears to consist of two straight lines, changing slope at a temperature around $`T^{}=5`$, rather than a nonlinear dependence. By reasoning as above, we identify $`T_{SR}^{ϵ=0.2}T_R=16\pm 1`$ and $`T_{SR}^{ϵ=0.8}=12\pm 1`$, values which are confirmed by the energy behavior (not shown), the height difference (Fig. 11) and slope (not shown) correlation functions, and by the dependence of $`m_1`$ on temperature (not shown).

An interesting question arises from Fig. 11: There is no evidence about the squared logarithmic behavior found for the RsGM and, furthermore, the plots exhibit a finite correlation length below $`T_{SR}`$ for both values of $`ϵ`$. Another intriguing fact is the nonmonotonic dependence of the correlation function on temperature for $`ϵ=0.8`$: From the curve for $`T=1`$, the scaled correlation function decreases up to $`T=5`$; upon further increasing the temperature, the evolution of the curves is very similar to that of the OsGM. This might be connected with the change in slope in the roughness curve mentioned in the preceding paragraph (see Fig. 10), but we have not been able to draw a clearer connection. All this is clear evidence that the behavior of the RsGM is significantly dependent on the disorder strength.

## IV Discussion and conclusions

Let us begin the discussion of the above results by analyzing our findings about the OsGM. Our simulations strongly support that $`T_R=16\pm 1`$ for the OsGM on a lattice with $`V_0=1`$. This is in contrast to the claims in falo ; yo1  that $`T_R=25`$ where a different way of definining the transition, which assumes the validity of the RG approach, was used (see Sec. III B and falo ; yo1 ; Shugard ). Further, the result is also in contrast with the RG prediction itself Weeks ; Beijeren ; Barabasi ; Villain2 ; Chui , which in our units is $`T_R=8\pi `$. However, we believe that the comparison with the linear theory for the EW high temperature phase has a much more physical character while keeping the basic RG ideas, and establishes beyond doubt that for the studied lattices the roughening temperature is $`T_R=16\pm 1`$ for $`V_0=1`$. Another hint in favor of our claim is the finding of the (Schottky) anomaly in the specific heat, which should appear below the transition temperature in view of what occurs for the XY and related models Kawabata . Finally, the fact that we obtain the same results for both the OsGM and the RsGM above $`T_R`$ is clear evidence that the potential is irrelevant (in the RG sense) in that regime and that we have indeed located the transition. Clearly, we cannot exclude the possibility that working on even larger lattices we would find the transition where the RG predicts it, but the absence of any finite size effects even for $`L=256`$ makes this possibility quite unlikely. Another possible reason for the discrepancy is the fact that our simulations are intrinsically discrete in space while RG theories for the OsGM are always applied to the continuum equation; again, very much larger lattices would remove this objection and clarify the effects of discreteness. Aside from that, we have also found that increasing (decreasing) $`V_0`$ increases (decreases) the roughening temperature: In the cases we studied, we found $`T_R^{V_0=0.2}=13\pm 1`$ and $`T_R^{V_0=5}=19\pm 1`$, which is intuitively reasonable as larger potential barriers require larger temperatures for the surface to overcome them. On the other hand, RG calculations are perturbative in $`V_0`$, so one would expect better agreement to the RG prediction for $`V_0=0.2`$, but in fact the agreement is worse in that case.

Let us now turn to the RsGM. In the “canonical” case, $`V_0=1`$, we found a superroughening transition at $`T_{SR}=4\pm 1=T_R/4`$, to be compared to RG predictions that it should occur at $`T_R/2`$. Below $`T_{SR}`$, we have obtained a $`\mathrm{ln}^2`$ dependence of the height difference correlation function, in agreement with RG results. However, we have clearly shown that the superroughening transition temperature depends on $`V_0`$, confirming the earlier report by Batrouni and Hwa Batrouni  on the absence of the transition in Langevin dynamics simulations for $`V_0=0.15`$, and later reports by Rieger Riegercom  and Ruiz-Lorenzo jjfises  who also observed this dependence in Monte Carlo simulations. In the opposite case, $`V_0=5`$, we find that $`T_{SR}^{V_0=5}=8\pm 1`$, considerably higher than the $`V_0=1`$ temperature. This disagrees with the RG predictions of a universal $`T_{SR}`$ independent of $`V_0`$. We believe that the agreement between our results and the previous ones Batrouni ; Riegercom ; jjfises  indeed supports a dependence of $`T_{SR}`$ on $`V_0`$, whose explanation remains an open question as far as RG is concerned.

A second, novel finding arises when considering our numerical results for $`V_0=5`$, which strongly suggest the possibility of two different low temperature phases. In our comments in the preceding paragraph, we took $`T_{SR}^{V_0=5}=8\pm 1`$ interpreting that the superroughening transition implies a change from a $`\mathrm{ln}^2`$ behavior of the height difference correlation function to a $`\mathrm{ln}`$ form (and the rest of EW features). However, the lack of size dependence of the roughness and the specific heat on temperature below $`T^{}=5\pm 1`$, with peaks absent for smaller values of $`V_0`$, raises the possible existence of another phase transition. If one looks at the correlation functions in Fig. 7, it turns out that for temperatures below $`T^{}=4`$ the correlation length is finite, in agreement with the roughness independence on the system size. Whereas the range of correlations above $`T^{}`$, which we believe is infinite, could be a subject of debate as we only have studied sizes up to $`L=256`$, our claim of finite correlation lengths below $`T^{}`$ is difficult to dispute. Further evidence in this regard is shown in Fig. 12, where curves for $`V_0=5`$ at $`T=1`$ are compared for two different system sizes. At this point, it is interesting to recall that in a previous paper us1  it was found that at $`T=0`$, the RsGM with $`V_0=1`$ exhibits a finite correlation length of about 20 lattice units (the reader may find Figs. 2 and 3 of us1  illustrative). Having this in mind, it is not unreasonable to conjecture that there is a $`T^{}`$ for the $`V_0=1`$ case, which could be below $`T=1`$ or close to 1 (the upper curve in Fig. 3 might already show a finite correlation length). We stress that this phase has not been previously reported in works at $`T=0`$: Thus, Zeng et al. Zeng  studied a discrete model but, being different from the Gaussian, their low temperature results can not be compared to ours, and the results of Rieger et al. Rieger1 ; Rieger2  do not allow to conclude anything in this respect. Intuitively, one can expect a finite correlation length phase at low temperatures and large $`V_0`$; in the limit $`V_0\mathrm{}`$, the surface follows the disorder (i.e., $`h(𝒓)=𝒉(𝒓)^{(\mathit{0})}+\mathit{2}𝒏(𝒓)\pi `$ everywhere), but the gradient term smoothens out the lack of correlations of $`h^{(0)}(𝒓)`$, the competition of these two effects yielding a finite correlation length. In a loose sense, this could be interpreted as Anderson localization taking over the coupling between neighboring sites with increasing $`V_0`$. This picture is confirmed by simulations for $`V_0=25`$ (Fig. 12): For such a large value of $`V_0`$ the correlation length is only 1 lattice unit, i.e., correlations reach only nearest neighbors. Clearly, the data presented here is not conclusive, but the conjecture that there are two transitions whose critical temperatures depend on $`V_0`$ is not unreasonable and deserves further consideration.

We have also found that the transition temperature and the correlation functions depend on the disorder strength. This is not unexpected, in so far as the change in the disorder distribution interacts with the periodicity of the sine potential, and therefore it is clear that when $`ϵ=2\pi `$, i.e., in the standard RsGM, we are in a special case. In this respect, our results appear to indicate that the RsGM (with $`ϵ=2\pi `$) is a very specific model, and that its behavior at low temperatures might not be representative of what one would find in an actual experiment, where the disorder can not be so precisely controlled. Another conclusion we may draw from our work is that there might be two classes of behavior at low temperatures for small and large $`ϵ`$: Small $`ϵ`$ models would behave very similarly to the OsGM, whereas large values of $`ϵ`$ would give rise to a more complex phenomenology with, e.g., nonmonotonic behavior of the correlation functions.

As a final conclusion, we remark that the most relevant result of the present work is the determination of the transition temperature from the high temperature phase to the low temperature phase (or phases) of both the OsGM and the RsGM. This poses a number of questions to be addressed either with greater numerical capabilities or with new analytical tools. We believe that the complex phase diagram of the RsGM is being partially unveiled, as our research supports previous findings such as the non universality of the transition temperature. Significantly, once we know where to look for the low temperature phase of the RsGM, we can investigate the nature of that phase (or phases). It is often stated that the RsGM is “glassy” in this regime; however, this assertion has never been really proven nor fully detailed and, furthermore, if the RG picture is qualitatively correct, it has to be recalled that it is a replica-symmetric theory, which implies that the superrough phase would not be glassy in the replica sense jjfises . We have obtained preliminary evidence that there are long-lived metastable states in the low temperature phase of the RsGM unpub , but in view of our present results and their non-universality we will examine this question more carefully in future work. Investigation of the dynamics of the RsGM will also be important; we recall that Batrouni and Hwa Batrouni  did find evidence for a superroughening transition in the dynamics of the model, and hence it would be worth checking their results for larger values of $`V_0`$. We hope that those analyses, along with measurements of nonlinear susceptibilities and of relaxation dynamics, will shed light on this difficult problem. Work along these lines is in progress.

###### Acknowledgements.

We thank Rodolfo Cuerno, Juan Jesús Ruiz Lorenzo, and Raúl Toral for discussions. Work at Leganés was supported by project PB96-0119. Work at Los Alamos was done under the auspices of the U.S. Department of Energy. A.S. and E.M. acknowledge the support by NATO CRG 971090, by DGESIC (Spain) through a “Acción Integrada Hispano-Británica,” and by Los Alamos National Laboratory for a stay at Los Alamos, where part of this work was carried out.


text:
# Magnetic Double Structure for 𝑆={1,1/2} Mixed-Spin Systems

## I Introduction

The Haldane gap system has been one of the most fascinating subjects in condensed matter physics, for which extensive experimental and theoretical investigations have been providing a variety of new interesting phenomena. One of the hot topics is the magnetic double structure in the quasi-one-dimensional (1D) Haldane systems observed for the rare-earth compounds $`R_2\mathrm{BaNiO}_5`$.  In the case with $`R=\mathrm{Y}^{3+}`$, the system has a disordered ground state with the Haldane gap.  If $`R`$ is substituted by other magnetic ions such as $`\mathrm{Nd}^{3+}`$ and $`\mathrm{Pr}^{3+}`$, the $`s=1/2`$ spins on the $`R^{3+}`$ sublattice order magnetically at low temperatures, thus giving rise to the magnetic double structure composed of the Haldane gap excitation and the gapless excitation induced by the long-range order. This system has been theoretically treated as the Haldane chains with a static staggered field which are induced by the ordered $`s=1/2`$ spins on the $`R^{3+}`$ sublattice so far. In a recent paper, however, it has been pointed out that the dynamics of the $`s=1/2`$ spins may be also important, stimulating us to treat the system by a mixed-spin model. Concerning the magnetic double structure, there is another interesting spin system of current interest, i.e. the quantum ferrimagnetic chain composed of two kinds of mixed spins. Several compounds have been already found, which indeed realize the mixed-spin system. This has been stimulating further intensive theoretical studies on this subject. In particular, it has been pointed out that such a ferrimagnetic chain has the double structure for the spin excitation spectrum, which controls the characteristic properties in the ferrimagnetic chain.

The above two subjects, which have been studied experimentally in different contexts, should possess the common interesting physics, because both systems are characterized by the mixture of two kinds of distinct spins. Motivated by these hot topics, we here study the properties of the antiferromagnetic mixed-spin system in detail, for which the $`s=1`$ and $`s=1/2`$ spin chains are stacked alternately. We exploit the two-dimensional (2D) system as the simplest model which possesses the magnetic long-range order at zero temperature. By computing the dispersion relation and the staggered magnetization by means of the Schwinger-boson mean-field theory, we clarify how our system generates the magnetic double structure, which naturally interpolates the above two interesting spin systems. We also show that the results obtained are consistent with the experimental findings for quasi-1D Haldane compounds $`R_2\mathrm{BaNiO}_5`$.

This paper is organized as follows. We introduce the model and then briefly summarize the Schwinger-boson techniques in Sec. II. We discuss the zero-temperature properties in Sec. III, and finally move to the thermodynamic properties in Sec. IV. A brief summary is given in Sec. V.

## II Schwinger-boson Mean-Field Theory

Let us consider a 2D mixed-spin model on the square lattice, which is described by the following Hamiltonian,

$`H`$ $`=`$ $`J_1{\displaystyle \underset{i,j,\eta }{}}\left[𝐒_{2i,2j}^{A_1}𝐒_{2i+\eta ,2j}^{B_1}+𝐒_{2i+\eta ,2j+1}^{B_2}𝐒_{2i,2j+1}^{A_2}\right]`$ (1)
$`+`$ $`{\displaystyle \underset{i,j,\eta }{}}\left[J_2𝐒_{2i,2j}^{A_1}𝐒_{2i,2j+\eta }^{A_2}+J_3𝐒_{2i+1,2j}^{B_1}𝐒_{2i+1,2j+\eta }^{B_2}\right],`$ (2)
where $`𝐒_{i,j}^A`$ ($`𝐒_{i,j}^B`$) is the spin operator at the $`(i,j)`$-th site in the $`(x,y)`$ plane, and $`\eta `$ implies the summation to be taken over nearest-neighbor sites. All the exchange couplings $`J_1,J_2`$ and $`J_3`$ are assumed to be antiferromagnetic. We here focus on the mixed-spin system composed of $`S^A=1/2`$ and $`S^B=1`$ spins, since it is straightforward to generalize the results to the arbitrary-spin case.

In Fig. 1, we have drawn the mixed-spin model schematically. The indices $`A_m`$ and $`B_m`$ label the sublattice in the unit cell. Note that this mixed-spin model is constructed in two ways depending on how we stack the independent spin chains by introducing three kinds of the coupling constants. By taking $`J_3J_1,J_2`$, we can study the characteristic properties of the $`s=1`$ Haldane spin chains coupled with the $`s=1/2`$ gapless spin chains, which may have the relevance to the Haldane gap system in a staggered field observed in the compounds $`R_2\mathrm{BaNiO}_5`$.  On the other hand, by setting the coupling constants $`J_1J_2,J_3`$, we can investigate how the independent ferrimagnetic chains $`(J_2=J_3=0)`$, which have been studied extensively in recent years, are combined to form the 2D system. The advantage of our mixed-spin approach is that we can systematically describe these two interesting systems by continuously varying the parameters $`J_1,J_2`$ and $`J_3`$.

We employ the Schwinger-boson mean-field theory (SBMFT) to study the above mixed-spin model. This method was applied successfully to the 2D Heisenberg model with uniform spins, and then to the bilayer system, the double-exchange system, the ferrimagnetic chain, etc. It is known that the SBMFT can describe the magnetically ordered phase, which is characterized by the condensation of the Schwinger-bosons. In the Schwinger-boson representation, the spin operators are expressed in terms of the boson creation and annihilation operators $`\gamma _\alpha ^{}`$, $`\gamma _\alpha `$, with the Pauli matrix $`𝝈`$ as $`𝐒=\frac{1}{2}\gamma _\alpha ^{}`$$`𝝈`$$`{}_{\alpha \beta }{}^{}\gamma _{\beta }^{}`$ $`(\alpha ,\beta =,)`$. Since the unit cell in our model includes four sites, we necessarily introduce eight kinds of Bose operators $`\gamma _\alpha =a_\alpha ^{(1)},a_\alpha ^{(2)},b_\alpha ^{(1)},b_\alpha ^{(2)}`$ which belong to the $`A_1,A_2,B_1,B_2`$ sublattices, respectively. By imposing the constraint $`\gamma _i^{}\gamma _i+\gamma _i^{}\gamma _i=2S^A`$ or $`2S^B`$ on each site, we can correctly map the original spin system to the boson system. Introducing the Lagrange multipliers $`\lambda _{ij}^A`$ and $`\lambda _{ij}^B`$, the Hamiltonian with the constraints is recast to

$`H`$ $`=`$ $`2{\displaystyle \underset{i,j,\eta }{}}[J_1(A_{ij\eta }^{}A_{ij\eta }+B_{ij\eta }^{}B_{ij\eta })`$ (3)
$`+`$ $`J_2C_{ij\eta }^{}C_{ij\eta }+J_3D_{ij\eta }^{}D_{ij\eta }]`$ (4)
$`+`$ $`{\displaystyle \underset{n=1}{\overset{2}{}}}[{\displaystyle \underset{(ij)A_n}{}}\lambda _{ij}^A({\displaystyle \underset{\sigma }{}}a_{ij\sigma }^{(n)}a_{ij\sigma }^{(n)}2S^A)`$ (5)
$`+`$ $`{\displaystyle \underset{(ij)B_n}{}}\lambda _{ij}^B({\displaystyle \underset{\sigma }{}}b_{ij\sigma }^{(n)}b_{ij\sigma }^{(n)}2S^B)]`$ (6)
$`+`$ $`{\displaystyle \underset{i,j,\eta }{}}[J_1(S^{A_1}S^{B_1}+S^{A_2}S^{B_2})`$ (7)
$`+`$ $`J_2S^{A_1}S^{A_2}+J_3S^{B_1}S^{B_2}],`$ (8)
where four bond operators are introduced as

$`A_{ij\eta }`$ $`=`$ $`a_{2i,2j}^{(1)}b_{2i+\eta ,2j}^{(1)}a_{2i,2j}^{(1)}b_{2i+\eta ,2j}^{(1)},`$ (9)
$`B_{ij\eta }`$ $`=`$ $`a_{2i,2j+1}^{(2)}b_{2i+\eta ,2j+1}^{(2)}a_{2i,2j+1}^{(2)}b_{2i+\eta ,2j+1}^{(2)},`$ (10)
$`C_{ij\eta }`$ $`=`$ $`a_{2i,2j}^{(1)}a_{2i,2j+\eta }^{(2)}a_{2i,2j}^{(1)}a_{2i,2j+\eta }^{(2)},`$ (11)
$`D_{ij\eta }`$ $`=`$ $`b_{2i+1,2j}^{(1)}b_{2i+1,2j+\eta }^{(2)}b_{2i+1,2j}^{(1)}b_{2i+1,2j+\eta }^{(2)}.`$ (12)
We perform a Hartree-Fock decomposition of eq. (8) by taking the thermal average $`<A_{ij\eta }>=A`$, $`<\lambda _{ij}^A>=\lambda _A`$, etc., which means that these values are assumed to be uniform and static.

By diagonalizing the mean-field Hamiltonian via the Bogoliubov transformation, we have

$`H_{\mathrm{MF}}`$ $`=`$ $`{\displaystyle \underset{𝐤\sigma }{}}{\displaystyle \underset{n=1}{\overset{2}{}}}\left[E_𝐤^{(1)}\alpha _{𝐤\sigma }^{(n)}\alpha _{𝐤\sigma }^{(n)}+E_𝐤^{(2)}\beta _{𝐤\sigma }^{(n)}\beta _{𝐤\sigma }^{(n)}\right],`$ (13)
where $`\alpha `$ and $`\beta `$ are the Bose operators for normal modes. The corresponding energy spectrums read

$`E_𝐤^{(1)}`$ $`=`$ $`\sqrt{{\displaystyle \frac{E_0\sqrt{E_1}}{2}}},E_𝐤^{(2)}=\sqrt{{\displaystyle \frac{E_0+\sqrt{E_1}}{2}}},`$ (14)
where $`E_0`$ and $`E_1`$ are given by

$`E_0`$ $`=`$ $`\lambda _A^2+\lambda _B^22d_𝐤^2e_𝐤^2f_𝐤^2,`$ (15)
$`E_1`$ $`=`$ $`(\lambda _A^2\lambda _B^2e_𝐤^2+f_𝐤^2)^2`$ (19)
$`4d_𝐤^2\left((\lambda _A\lambda _B)^2(e_𝐤+f_𝐤)^2\right),`$
with

$`d_𝐤`$ $`=`$ $`2AJ_1\mathrm{cos}k_x,e_𝐤=2CJ_2\mathrm{cos}k_y,`$ (20)
$`f_𝐤`$ $`=`$ $`2DJ_3\mathrm{cos}k_y.`$ (21)
By minimizing the free energy thus obtained at finite temperatures, we end up with the self-consistent equations for $`A=B,C,D,\lambda _A,\lambda _B`$,

$`1+2S^A={\displaystyle \underset{n}{}}{\displaystyle \frac{d𝐤}{\pi ^2}\mathrm{coth}\kappa E_𝐤^{(n)}\frac{E_𝐤^{(n)}}{\lambda _A}},`$ (22)
$`1+2S^B={\displaystyle \underset{n}{}}{\displaystyle \frac{d𝐤}{\pi ^2}\mathrm{coth}\kappa E_𝐤^{(n)}\frac{E_𝐤^{(n)}}{\lambda _B}},`$ (23)
$`8J_1A={\displaystyle \underset{n}{}}{\displaystyle \frac{d𝐤}{\pi ^2}\mathrm{coth}\kappa E_𝐤^{(n)}\frac{E_𝐤^{(n)}}{A}},`$ (24)
$`4J_2C={\displaystyle \underset{n}{}}{\displaystyle \frac{d𝐤}{\pi ^2}\mathrm{coth}\kappa E_𝐤^{(n)}\frac{E_𝐤^{(n)}}{C}},`$ (25)
$`4J_3D={\displaystyle \underset{n}{}}{\displaystyle \frac{d𝐤}{\pi ^2}\mathrm{coth}\kappa E_𝐤^{(n)}\frac{E_𝐤^{(n)}}{D}},`$ (26)
with $`\kappa =1/(2k_BT)`$, where we have assumed that the bond operators which link the $`s=1`$ and $`s=1/2`$ spins take the same mean value, $`A=B`$. This completes our formulation based on the SBMFT. In the following sections, we solve these self-consistent equations to estimate the excitation spectrum and the thermodynamic quantities. Since it is not easy to analytically perform the Bogoliubov transformation, we numerically diagonalize the mean-field Hamiltonian to compute the energy dispersions $`E_𝐤^{(1)}`$ and $`E_𝐤^{(2)}`$.

## III Properties at Zero Temperature

In order to treat the ground state properties, it should be taken into account that the Bogoliubov particles of the $`\alpha `$ branch in eq.(13) may condense at absolute zero, because the excitation energy $`E_𝐤^{(1)}`$ has its minimal value $`E_𝐤^{(1)}=0`$ at $`𝐤=\mathrm{𝟎}`$ while the $`\beta `$ branch has a finite gap even at $`T=0`$. Sarker et al. showed that the long-range order is described by the condensation of the Schwinger-bosons for the ferromagnetic and antiferromagnetic Heisenberg models. This is also the case for our 2D mixed-spin model on a square lattice.

Suppose that the bosons condense at the states of $`\alpha _{}^{(1)}|_{𝐤=\mathrm{𝟎}}`$ and $`\alpha _{}^{(2)}|_{𝐤=\mathrm{𝟎}}`$, by fictitiously applying an infinitesimal external staggered field to the $`A`$ and $`B`$ lattices. The self-consistent equations at $`T=0`$ , which include the Bose condensation, now read,

$`\mathrm{\Gamma }`$ $`=`$ $`{\displaystyle \frac{2}{N^2}}\mathrm{coth}\kappa E_{𝐤=\mathrm{𝟎}}^{(1)}|_\kappa \mathrm{}`$ (27)
$`=`$ $`{\displaystyle \frac{\left[1+2S^A_n\frac{d𝐤}{\pi ^2}\frac{E_𝐤^{(n)}}{\lambda _A}\right]}{E_{𝐤=\mathrm{𝟎}}^{(1)}/\lambda _A}},`$ (28)
$`1+2S^B`$ $`=`$ $`\mathrm{\Gamma }{\displaystyle \frac{E_𝐤^{(1)}}{\lambda _B}}|_{𝐤=\mathrm{𝟎}}+{\displaystyle \underset{n}{}}{\displaystyle \frac{d𝐤}{\pi ^2}\frac{E_𝐤^{(n)}}{\lambda _B}},`$ (29)
$`4A`$ $`=`$ $`\mathrm{\Gamma }{\displaystyle \frac{E_𝐤^{(1)}}{d_𝐤}}|_{𝐤=\mathrm{𝟎}}+{\displaystyle \underset{n}{}}{\displaystyle \frac{d𝐤}{\pi ^2}\mathrm{cos}k_x\frac{E_𝐤^{(n)}}{d_𝐤}},`$ (30)
$`2C`$ $`=`$ $`\mathrm{\Gamma }{\displaystyle \frac{E_𝐤^{(1)}}{e_𝐤}}|_{𝐤=\mathrm{𝟎}}+{\displaystyle \underset{n}{}}{\displaystyle \frac{d𝐤}{\pi ^2}\mathrm{cos}k_y\frac{E_𝐤^{(n)}}{e_𝐤}},`$ (32)
$`2D`$ $`=`$ $`\mathrm{\Gamma }{\displaystyle \frac{E_𝐤^{(1)}}{f_𝐤}}|_{k=0}+{\displaystyle \underset{n}{}}{\displaystyle \frac{d𝐤}{\pi ^2}\mathrm{cos}k_y\frac{E_𝐤^{(n)}}{f_𝐤}}.`$ (34)
We shall solve these equations numerically for given coupling constants $`J_1,J_2`$ and $`J_3`$.

### A Dispersion relation

We start our discussions with the case that may be regarded as the Haldane gap system in a staggered field. When $`J_1=0`$, the 2D system is completely decoupled into the $`s=1`$ massive Haldane chains and the $`s=1/2`$ massless spin chains. When we introduce the interchain couplings among them, the ground state has the long-range order mainly due to the $`s=1/2`$ spins. The important point is that even though we have the long-range order, the Haldane-type gapful excitation still exists. Therefore, as far as the case with small $`J_1`$ and $`J_2`$ is concerned, the system is regarded as the one often called the Haldane gap system in a staggered field. In Fig. 2, we show the dispersion relation calculated for $`J_1=J_2=1/2`$ and $`J_3=1`$.

As seen in this figure, the lower branch is gapless with linear dispersion relation, reflecting the antiferromagnetic long-range order, whereas the upper optical mode is mainly composed of the Haldane-type excitation. It is observed that the dispersion in the $`k_x`$-direction is indeed weak for the optical mode, since the interchain coupling is much smaller than the energy gap for the optical mode, which confirms that the system in Fig. 2 may be regarded as the Haldane gap system in a staggered field. Although we have also performed the calculation for the cases with smaller $`J_1`$ and $`J_2`$, the obtained dispersion of the optical branch becomes almost flat in the $`k_x`$-direction, so that we have not shown them here. Concerning this limit of small $`J_1`$, we here make a brief comment on the validity of the SBMFT. When $`J_1`$ takes too small value, for which the system is almost decoupled into independent $`s=1`$ and $`s=1/2`$ chains, the SBMFT may lead to a pathological result: although the Haldane gap for the $`s=1`$ chain is well described by the SBMFT, it is not the case for the gapless $`s=1/2`$ chains, for which we are left with a gapful phase. Actually, if the value of $`J_1`$ ($`=J_2`$) becomes smaller than $`J_1/J_310^2`$, we encounter a problem that the present self-consistent calculation does not converge, implying that our assumption for the antiferromagnetic ordered state does not hold anymore. Nevertheless, we find that the correct behavior with the antiferromagnetic ground state is still obtained except for this small parameter region.

By increasing the coupling parameters $`J_2,J_3`$ continuously, we naturally enter in the 2D ordered mixed-spin system. Note that although during this process the magnetic double structure is kept unchanged in its typical feature, the nature of the optical mode is gradually changed from the Haldane-gap excitation: i.e. the gapful excitation may be equally contributed from both the $`s=1`$ and $`s=1/2`$ spin sectors. As a reference, we show the dispersion relation for the mixed-spin model with the isotropic bonds $`J_1=J_2=J_3=1`$ in Fig. 3.

If we further decrease the couplings $`J_2`$ and $`J_3`$, the system gradually approaches the quasi-1D ferrimagnetic chains with the periodic arrangement of spins $`1/211/21`$. The dispersion relation obtained in the corresponding parameter region is shown in Fig. 4.

Now the optical mode with a weak dispersion in the $`k_y`$-direction is essentially the same as that found for the ferrimagnetic chain.

In order to clearly observe how the gapless dispersion changes its character in the ferrimagnetic-chain limit, we have shown the low-energy dispersion relation at $`k_y`$=0 in Fig. 5. It is seen that with the decrease of the coupling constants $`J_2`$ and $`J_3`$, the $`k_x`$-linear dependence is gradually changed to the $`k_x^2`$ dependence characteristic of the ferrimagnetic chain except for the small $`k_x`$ region where the 2D antiferromagnetic order still gives rise to the $`k_x`$-linear dependence. At $`J_2=J_3=0`$, the system is reduced to the isolated ferrimagnetic chains, for which the 2D character completely disappears, and thus the gapless dispersion simply follows the $`k_x^2`$ dependence. We note that the ferrimagnetic chains have been experimentally realized in the compounds such as $`\mathrm{NiCu}(\mathrm{C}_2\mathrm{O}_4)_24\mathrm{H}_2\mathrm{O}`$ and $`\mathrm{MnCu}(\mathrm{pba})(\mathrm{H}_2\mathrm{O})_32\mathrm{H}_2\mathrm{O}`$, and have been studied theoretically by many groups. Among others, it has been reported  that this system exhibits the dual properties consistent with our results: the physical quantities are controlled by the $`s=1`$ antiferromagnetic spin chain at high temperatures and by the effective $`s=1/2`$ ferromagnetic spin chain at low temperatures.

The above analysis of the excitation spectrum implies that the quasi-1D Haldane gap system in a staggered field and the quasi-1D weakly coupled ferrimagnetic chains, which have been studied in different contexts experimentally, share common interesting physics inherent in the mixed-spin systems. In particular, it is highly desirable to experimentally observe the magnetic double structure in the excitation spectrum for the ferrimagnetic-chain compounds.

### B Haldane gap in a staggered field

Let us discuss the case of the Haldane gap system in a staggered field in more detail. We here observe how the effective staggered field induced by the $`s=1/2`$ spin chains affects the properties of the Haldane gap, by changing the coupling constants $`J_1`$ and $`J_2`$. For small values of $`J_1`$ and $`J_2`$, we define the effective staggered field on the $`s=1`$ Haldane chain by $`H_{\mathrm{ST}}^{\mathrm{eff}}=J_1<S_z^A>`$, where $`<S_z^A>`$ is the spontaneous staggered magnetization of the $`s=1/2`$ spin chain. We numerically estimate the effective staggered field, and show the obtained results in Fig. 6.

It is seen that the effective staggered field increases monotonously with the increase of $`J_1`$ and $`J_2`$, as should be expected. In a similar way, we also compute the staggered magnetization $`<S_z^B>`$ of the $`s=1`$ spin chains, and as well as the Haldane gap which is defined as the minimum of the excitation energy in the optical branch $`\mathrm{\Delta }=E_{𝐤=(\frac{\pi }{2},0)}^{(2)}`$.

These quantities are shown in Figs. 7 and 8. Following the way used for the experimental analysis,  we have plotted the staggered magnetization (the square of the Haldane gap) as a function of the effective staggered field (the square of the staggered magnetization). It is seen that as the staggered field is increased, the staggered magnetization as well as the Haldane gap $`\mathrm{\Delta }`$ are increased, being consistent with the result pointed out for the Haldane chain system with a static magnetic field.

Since we are dealing with the 2D model at zero temperature, our results may not be directly applied to the experiments. Nevertheless, we can confirm whether the present results are consistent with the experiments for the rare-earth compounds $`R_2\mathrm{BaNiO}_5`$  with $`R=\mathrm{Nd}^{3+}`$ or $`\mathrm{Pr}^{3+}`$. Experimentally, as the temperature is decreased, the system shows the phase transition to the magnetically ordered phase. When the temperature is further decreased, the staggered moment for the $`s=1/2`$ sector develops, and thereby gives rise to the increase of the staggered field, which indeed enhances the magnitude of the Haldane gap. Also, the staggered moment on the Haldane chains increases, as should be expected. These characteristic features are consistent with our results, and in particular, the qualitative behaviors for the staggered magnetization and the Haldane gap shown in Figs. 7 and 8 agree fairly well with experimental findings.

### C Coupled ferrimagnetic chains

Before closing this section, we briefly discuss the case close to the ferrimagnetic chain, which is realized by taking the limit of $`J_1J_2,J_3`$. Even in this one-dimensional limit, the system still exhibits the antiferromagnetic order as far as $`J_2`$ and $`J_3`$ take finite values. As mentioned above, at $`J_2=J_3=0`$, the system shows the ferrimagnetic order, and hence the dispersion relation for the acoustic mode is changed to the quadratic one. In this way, the crossover-like behavior is seen in the gapless mode, whereas much simpler behavior is observed in the gapful mode. In Fig. 9, we display the spin gap $`\mathrm{\Delta }`$ for the optical branch in the ferrimagnetic-chain limit. With decreasing $`J_2`$ and $`J_3`$, the spin gap monotonically decreases, and reaches the value of $`\mathrm{\Delta }`$=1.778 at $`J_2=J_3=0`$, which is very close to that of the quantum Monte Carlo method, $`\mathrm{\Delta }`$=1.767, as already demonstrated by Wu et al. In this way, in the ferrimagnetic-chain limit, the SBMFT may provide the reliable estimates of the physical quantities even at quantitative level. As discussed by Yamamoto and Fukui, this optical mode is mainly composed of excitations in the effective $`s=1`$ antiferromagnetic spin chain, whereas the acoustic mode is given by excitations in the effective $`s=1/2`$ ferromagnetic chain.

## IV Thermodynamic Properties at Finite Temperatures

We now move to the thermodynamic properties at finite temperatures. The self-consistent equations (22)-(26) are solved numerically at finite temperatures to obtain the thermodynamic quantities. We here briefly summarize the results obtained. In the previous section, it has been shown that even for the set of the parameters $`J_1=J_2=1/2,J_3=1`$ ($`J_1=1,J_2=J_3=1/2`$), the system may be approximately described by the Haldane gap system in a staggered field (coupled ferrimagnetic chains), so that we shall show the results for these parameters below.

We start with the effective spin gaps calculated at finite temperatures, which are shown as a function of the temperature in Fig. 10. It is seen that both of the two cases exhibit similar temperature dependence in the spin gap. For each case, there are two distinct spin gaps, corresponding respectively to the optical mode and the acoustic mode, reflecting the double structure for the excitation spectrum. Note that both of two modes should be massive at finite temperatures. It may be physically sensible to regard these spin gaps as the inverse of the correlation lengths. As should be expected, the spin gap for the optical mode increases up to the zero-temperature value with the decrease of the temperature. On the other hand, the spin gap for the acoustic mode decreases, leading to the divergent correlation length which characterizes the magnetically ordered state at $`T=0`$.

The uniform and staggered spin susceptibilities, $`\chi _{\mathrm{uni}}`$ and $`\chi _{\mathrm{stag}}`$, are calculated by using the standard linear-response formulae,

$`\chi _{\mathrm{uni}}={\displaystyle \underset{i,j=1,2}{}}[S_z^{A_i}S_z^{A_j}_{𝐪,\omega }+S_z^{B_i}S_z^{B_j}_{𝐪,\omega }`$ (36)
$`+2S_z^{A_i}S_z^{B_j}_{𝐪,\omega }]|_{𝐪,\omega 0},`$ (37)
(38)
$`\chi _{\mathrm{stag}}=`$ (39)
$`{\displaystyle \underset{i,j=1,2}{}}[(1)^{i+j}(S_z^{A_i}S_z^{A_j}_{𝐪,\omega }+S_z^{B_i}S_z^{B_j}_{𝐪,\omega })`$ (40)
$`+2(1)^{i+j+1}S_z^{A_i}S_z^{B_j}_{𝐪,\omega }]|_{𝐪,\omega 0},`$ (41)
where $`S_z^{A_1}S_z^{A_1}_{𝐪,\omega }`$, etc., are the retarded spin correlation functions.

In Fig. 11, we plot the effective Curie constants $`T\chi _{\mathrm{uni}}`$ and $`T\chi _{\mathrm{stag}}`$ as a function of $`T`$. It is seen that the effective Curie constant for the uniform sector is gradually decreased as the temperature is decreased, while that for the staggered sector diverges, implying that the antiferromagnetic correlation is enhanced at low temperatures.

As seen from the above results, the thermodynamic properties at finite temperatures show quite similar behavior both for the Haldane chain in a staggered field and for the quasi-1D ferrimagnetic chain, reflecting the magnetic double structure inherent in these mixed-spin systems.

## V summary

We have studied the 2D mixed-spin model for which the $`s=1/2`$ and $`s=1`$ spin chains are stacked alternately. This mixed-spin model includes two interesting spin systems addressed recently: the Haldane gap system in a staggered field as well as the ferrimagnetic chain. By calculating the dispersion relation and the thermodynamic quantities by means of the Schwinger-boson mean-field theory, we have discussed the magnetic double structure inherent in our mixed-spin systems. In particular, we have treated systematically the quasi-1D Haldane gap system in a staggered magnetic field and also the mixed-spin chain with the ferrimagnetic ground state. This implies that these two spin systems, which have been studied in different contexts experimentally, should possess interesting physics common to the mixed spin systems. We have also found that the results obtained for the staggered-field effect on the Haldane gap system are qualitatively consistent with the experimental findings in the rare-earth compounds $`R_2\mathrm{BaNiO}_5`$. It remains an interesting problem to evaluate dynamical quantities related to the neutron scattering, etc. Also, it is important to study how the string-order parameter behaves, when the system changes from the Haldane system to the ferrimagnetic chain. These problems are now under consideration.

## acknowledgments

The work is partly supported by a Grant-in-Aid from the Ministry of Education, Science, Sports, and Culture. A. K. is supported by the Japan Society for the Promotion of Science. N. K. wishes to thank K. Ueda for the warm hospitality at ISSP, University of Tokyo.