format: text

text:
# Stochastic unraveling of time-local quantum master equations beyond the Lindblad class

## Abstract

A new method for stochastic unraveling of general time-local quantum master equations (QME) which involve the reduced density operator at time $`t`$ only is proposed. The present kind of jump algorithm enables a numerically efficient treatment of QMEs which are not of Lindblad form. So it opens new large fields of application for stochastic methods. The unraveling can be achieved by allowing for trajectories with negative weight. We present results for the quantum Brownian motion and the Redfield QMEs as test examples. The algorithm can also unravel non-Markovian QMEs when they are in a time-local form like in the time-convolutionless formalism.

Quantum master equations (QMEs) are frequently used to describe time-independent as well as time-dependent phenomena in chemical physics, quantum optics, solid state physics, biological physics, etc. (see Ref. Weiss, 1999 for a number of typical examples). These QMEs describe the time evolution of density matrices which are used in order to represent the mixed nature of the states. Stochastic unraveling is an efficient numerical tool for solving such equations. This method allows one to simulate much larger and more complex systems with many degrees of freedom. It can, for example, be used to accurately describe femtochemical experiments in the liquid phase whose description has been limited until now to models with one or two effective interaction coordinates. In the unraveling scheme one considers an ensemble of stochastic Schrödinger equations (SSEs) which in the limit of a large ensemble resembles the respective QME. The numerical effort scales much more favorably with the size of the basis since one is now dealing with wave functions and not density matrices anymore (for a comparison of direct integrators, see Ref. Kondov et al., 2001). Another aspect of the stochastic methods is the possible physical interpretation of experiments detecting macroscopic fluctuations (e.g. photon counting) in various quantum systems Plenio and Knight (1998). Most of the unraveling schemes Dalibard et al. (1992); Gardiner et al. (1992); Gisin and Percival (1992); Garraway and Knight (1994); Wolfseder and Domcke (1995); Plenio and Knight (1998) have been restricted to QMEs of Lindblad form Lindblad (1975) which ensures that the reduced density matrix (RDM) stays positive semi-definite for all times and all parameters. Nevertheless there are many physical meaningful QMEs which result in positive- or almost positive-definite RDMs although they are not of Lindblad form. The increasing interest in descriptions beyond the Lindblad class such as the quantum Brownian motion Strunz et al. (1999); Yu et al. (1999), the Redfield formalism May and Kühn (2000), non-Markovian schemes Čápek (1994); Meier and Tannor (1999); Breuer et al. (1999), etc. resulted in various efforts to develop new stochastic wave function algorithms.

Strunz et al. Strunz et al. (1999); Yu et al. (1999) developed the non-Markovian Quantum Diffusion Model. In general, this method can also be applied to QMEs in Markov approximation even though they might not preserve positivity (see also Ref. Budini (2000)). A similar approach was also proposed by Gaspard et al. Gaspard and Nagaoka (1999a). Very recently Stockburger and GrabertStockburger and Grabert  developed a method how to exactly represent the RDM of a system coupled to a linear heat bath in terms of SSEs. The numerical properties of this approach need to be explored. Breuer et al. Breuer et al. (1999) extended a scheme which they had used to calculate multi-time correlation functions Breuer et al. (1997) to the unraveling of QMEs. Their technique is based on doubling the Hilbert space. Instead of a single stochastic wave function one has a pair of them Breuer et al. (1999). This scheme conserves Hermiticity of the RDM only on average and not for every single realization. Thus, the deviation from Hermiticity is a quantity with statistical error and one has to perform a huge number of realizations in order to achieve good convergence. Since stability and efficiency are crucial issues for unraveling algorithms we propose in this article an alternative approach which fulfills these criteria.

The aim is to represent, in terms of quantum trajectories, the solution $`\rho (t)`$ of a generalized time-local Hermiticity-conserving QME

$`{\displaystyle \frac{d\rho (t)}{dt}}`$ $`=`$ $`A(t)\rho (t)+\rho (t)A^{}(t)`$ (1)
$`+{\displaystyle \underset{k=1}{\overset{M}{}}}\left\{C_k(t)\rho (t)E_k^{}(t)+E_k(t)\rho (t)C_k^{}(t)\right\}`$
with the total number $`M`$ of dissipative channels and arbitrary operators $`A(t)`$, $`C_k(t)`$, and $`E_k(t)`$. Examples for these operatores are given below. Here we restrict the operators in such a way that the norm of the solution stays conserved. For readability we shall omit the time arguments in the following.

In order to approach the problem let us define a state vector $`(|\psi ,|\varphi )^T`$ spanning a doubled Hilbert space as proposed in Ref. Breuer et al., 1999. Unlike Ref. Breuer et al., 1999 the RDM shall be reproduced by an ensemble average (denoted by overbars) of outer products of the vectors $`|\psi `$ and $`|\varphi `$:

$$\rho =\overline{|\psi \varphi |}+\overline{|\varphi \psi |}.$$
(2)
A particular realization of the stochastic process will be denoted by the pair $`(|\psi ,|\varphi )`$. The averaging is performed over all trajectories possibly including a weighted sum over pure initial states. A vantage of this averaging is the conservation of Hermiticity for every single trajectory in contrast to Ref. Breuer et al. (1999). We note that this small modification improves the numerical efficiency significantly.

For the SSEs let us consider $`2M`$ independent possibly complex noise variables $`\xi _k^i(t)`$. The superscripts denote which of the two terms from the Hermitian pair in the sum in Eq. (1) is taken and subscripts denote the various dissipative channels. All stochastic differentials $`d\xi _k^i(t)`$ are assumed to have zero mean, to be normalized and uncorrelated Gardiner (1985):

$`\overline{d\xi _k^i}=0,\overline{d\xi _k^id\xi _l^j}=\delta _{ij}\delta _{kl}dt.`$ (3)
Next, as an ansatz we construct a SSE which propagates the pair $`(|\psi ,|\varphi )`$

$`d|\psi `$ $`=`$ $`D_1|\psi dt+{\displaystyle \underset{k=1}{\overset{M}{}}}{\displaystyle \underset{i=1}{\overset{2}{}}}S_{1k}^i|\psi d\xi _k^i,`$ (4a)
$`d|\varphi `$ $`=`$ $`D_2|\varphi dt+{\displaystyle \underset{k=1}{\overset{M}{}}}{\displaystyle \underset{i=1}{\overset{2}{}}}S_{2k}^i|\varphi d\xi _k^i.`$ (4b)
The operators $`D_1`$ and $`D_2`$ govern the deterministic and the operators $`S_{jk}^i`$ the stochastic part of the evolution. In general, they may depend on the state vector and explicitly on time. After differentiating Eq. (2), neglecting all terms higher than first order in $`dt`$, and assuming that ensemble averages always factorize van Dorsselaer and Nienhuis (2000) one obtains

$`d\rho `$ $`=`$ $`\left[D_1\overline{|\psi \varphi |}+D_2\overline{|\varphi \psi |}\right]dt`$
$`+{\displaystyle \underset{k=1}{\overset{M}{}}}\left[S_{1k}^1|\overline{\psi \varphi |}S_{2k}^1+S_{2k}^2|\overline{\varphi \psi |}S_{1k}^2\right]dt+h.c.`$
Comparing with Eq. (1) one notes that $`S_{1k}^1`$ has to equal $`S_{2k}^2`$ and $`S_{2k}^1`$ has to equal $`S_{1k}^2`$. Moreover, one can see that $`S_{2k}^1=C_k+\alpha _k^1`$ and $`S_{2k}^2=E_k+\alpha _k^2`$ with $`\alpha _k^1`$ and $`\alpha _k^2`$ being arbitrary scalar functions of $`(|\psi ,|\varphi )^T`$ and possibly of time. Making the latter substitutions in Eq. (Stochastic unraveling of time-local quantum master equations  beyond the Lindblad class) yields the constraint

$`D_1=D_2=A{\displaystyle \underset{k=1}{\overset{M}{}}}\left(\alpha _k^2C_k+\alpha _k^1E_k+\alpha _k^1\alpha _k^2\right).`$ (6)
Any quantum jump method is specified by jump rates $`p_k^i`$ which have to be real scalar functions of $`(|\psi ,|\varphi )`$. If $`n_k^i(t)`$ is the number of jumps in channel $`k`$ and due to term $`i`$ up to time $`t`$, the probability for $`n_k^i(t)`$ to increase by one, i.e. the expectation value of both $`dn_k^i`$ and $`(dn_k^i)^2`$, is equal to $`p_k^idt`$ during the infinitesimal time interval $`dt`$. Thus, the noise variables $`\xi _k^i`$ obeying Eq. (3) are related to $`dn_k^i(t)`$ as van Dorsselaer and Nienhuis (2000)

$`d\xi _k^i={\displaystyle \frac{dn_k^ip_k^idt}{\sqrt{p_k^i}}}e^{i\phi }.`$ (7)
The phase factor $`e^{i\phi }`$ does not change the RDM expressions within each realization and can be set to one. Substituting Eq. (7) into Eq. (4) one finds that $`\alpha _k^i=\sqrt{p_k^i}`$. So the SSEs for our quantum jump method read

$`d|\psi `$ $`=`$ $`\left(A+{\displaystyle \underset{k=1}{\overset{M}{}}}{\displaystyle \frac{p_k^1+p_k^2}{2}}\right)|\psi dt`$
$`+{\displaystyle \underset{k=1}{\overset{M}{}}}\left[\left({\displaystyle \frac{E_k}{\sqrt{p_k^1}}}1\right)dn_k^1+\left({\displaystyle \frac{C_k}{\sqrt{p_k^2}}}1\right)dn_k^2\right]|\psi ,`$
$`d|\varphi `$ $`=`$ $`\left(A+{\displaystyle \underset{k=1}{\overset{M}{}}}{\displaystyle \frac{p_k^1+p_k^2}{2}}\right)|\varphi dt`$
$`+{\displaystyle \underset{k=1}{\overset{M}{}}}\left[\left({\displaystyle \frac{C_k}{\sqrt{p_k^1}}}1\right)dn_k^1+\left({\displaystyle \frac{E_k}{\sqrt{p_k^2}}}1\right)dn_k^2\right]|\varphi .`$
The jump rates $`p_k^1`$ and $`p_k^2`$ still remain free parameters. In the statistical limit their values have no influence on any averaged physical quantity. Nevertheless, it turns out that they can strongly influence the convergence behavior of the jump algorithm, i.e. they determine the statistical error of the observables calculated. A detailed discussion of this influence and utilization of such free parameters can be found in Ref. Felbinger and Wilkens, 1999.

To ensure an efficient scheme with fast convergence one has to require that the norm of every single trajectory is constant in time. Asking for $`\varphi |\varphi `$, $`\psi |\psi `$, etc. being constant in time does not create a stable scheme but the condition of norm preservation of $`|\psi \varphi |+|\varphi \psi |`$

$`\mathrm{Tr}\left\{{\displaystyle \frac{d}{dt}}\left[|\psi \varphi |+|\varphi \psi |\right]\right\}=0`$ (9)
does. Unfortunately, applying this condition does not lead to positive values of the jump rates $`p_k^i`$ for all trajectories at all times. However, since the $`p_k^i`$ are arbitrary real functions, they can be replaced by their absolute values. The price to pay is that we have to introduce an additional weight factor for the trajectories which jumps between one and minus one. In addition, there is a small deviation of the norm from unity because in the regions where the $`p_k^i`$ are replaced by their absolute values norm conservation is no longer guaranteed. But in all our tests this deviation was far below one percent and neither effected numerical stability nor efficiency. The negative weights are actually needed to reconstruct RDMs which are, in general, not positive-semi-definite. If the RDM stays positive-semi-definite during its entire time-evolution the negative weights of some trajectories are not needed, i.e. all trajectories can be normalized to unity and represent physically pure states of the open quantum system. In the examples below the RDM can exhibit negative populations. This unphysical situation could probably be cured by applying an initial slippage to the initial state Gaspard and Nagaoka (1999b); Yu et al. (2000). We note that this physically unreasonable RDMs occur because of unphysical initial states or because the QME is not physically correct or is applied in a parameter region where it is is not valid. Nevertheless an unraveling scheme has to be able to mimic also this unphysical behavior of the QME because in the ensemble average both should fully coincide.

The condition (9) applied to the QME (1) results in the additional constraint

$`A+A^{}={\displaystyle \underset{k=1}{\overset{M}{}}}\left(E_k^{}C_k^{}+C_k^{}E_k^{}\right)`$ (10)
and if applied to the deterministic part of the corresponding SSE (8) yields the total jump rate

$`p={\displaystyle \frac{\varphi |A+A^{}|\psi +\psi |A+A^{}|\varphi }{\varphi |\psi +\psi |\varphi }}.`$ (11)
All partial jump rates can be found subsequently making use of Eqs. (10) and (11):

$`p_k^1`$ $`=`$ $`{\displaystyle \frac{\varphi |C_k^{}E_k^{}|\psi +\psi |E_k^{}C_k^{}|\varphi }{\varphi |\psi +\psi |\varphi }},`$ (12a)
$`p_k^2`$ $`=`$ $`{\displaystyle \frac{\varphi |E_k^{}C_k^{}|\psi +\psi |C_k^{}E_k^{}|\varphi }{\varphi |\psi +\psi |\varphi }}.`$ (12b)
In the rest of this article let us briefly show how the proposed method can be applied to two typical physical problems: the quantum Brownian motion and dissipative electron transfer within Redfield theory. In both cases the systems are described by Markovian QMEs which do not have Lindblad structure. The model of Brownian motion Weiss (1999) describes a particle with mass $`m`$, coordinate $`q`$, momentum $`p`$ and Hamiltonian $`H_\mathrm{S}`$ interacting with a thermal bath. In the high temperature limit of a bath of harmonic oscillators the relevant QME has the form

$`{\displaystyle \frac{d\rho }{dt}}={\displaystyle \frac{i}{\mathrm{}}}[H_\mathrm{S},\rho ]{\displaystyle \frac{i\gamma }{2\mathrm{}}}[q,\{p,\rho \}]{\displaystyle \frac{m\gamma kT}{\mathrm{}^2}}[q,[q,\rho ]]`$ (13)
where $`\gamma `$ is the damping rate. Comparing with Eq. (1) one finds the operators of the jump algorithm ($`M=2`$)

$`E_1`$ $`=`$ $`\sqrt{{\displaystyle \frac{\gamma }{2\mathrm{}}}}q,C_1=i\sqrt{{\displaystyle \frac{\gamma }{2\mathrm{}}}}p,`$ (14a)
$`E_2`$ $`=`$ $`\sqrt{{\displaystyle \frac{m\gamma kT}{\mathrm{}^2}}}q,C_2=E_2,`$ (14b)
$`A`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}H_\mathrm{S}+{\displaystyle \frac{i\gamma }{2\mathrm{}}}qp{\displaystyle \frac{m\gamma kT}{\mathrm{}^2}}qq.`$ (14c)
Modeling the particle as a harmonic oscillator with eigenfrequency $`\omega `$ one can compute the population dynamics depicted in Fig. 1. The initial state of the oscillator is the pure state $`\rho _{33}=1`$. As can be seen, the agreement of the results using our new stochastic method and using a direct integration of the QME is already quite good for one thousand samples.

As a next test for the present quantum jump method we shall demonstrate the stochastic unraveling of the Redfield QME May and Kühn (2000); Kleinekathöfer et al. (2001)

$`\dot{\rho }={\displaystyle \frac{i}{\mathrm{}}}[H_\mathrm{S},\rho ]+{\displaystyle \frac{1}{\mathrm{}^2}}\left\{[\mathrm{\Lambda }\rho ,K]+[K,\rho \mathrm{\Lambda }^{}]\right\}`$ (15)
in which $`\mathrm{\Lambda }`$ is the relaxation operator and $`K`$ the system part of the system-bath interaction May and Kühn (2000). Let us consider a model for electron transfer in which the system includes a single reaction coordinate with the Hamiltonian May and Kühn (2000); Kleinekathöfer et al. (2001)

$`H_\mathrm{S}=H_1|11|+H_2|22|+v_{12}(|12|+|21|)`$ (16)
where $`H_1`$ and $`H_2`$ are the Hamiltonians of two coupled harmonic oscillators with frequency $`\omega `$. We choose a potential configuration in the normal region with no barrier between the two harmonic potentials (change of free energy $`\mathrm{\Delta }E=2\omega `$, reorganization energy $`\lambda =3\omega `$) with inter-center coupling $`v_{12}=\omega `$. The bath is described by a cut-off frequency $`\omega _c=\omega `$ and temperature $`kT=\omega /4`$. The system-bath interaction is characterized by the damping rate $`\mathrm{\Gamma }=\pi \eta /(\mathrm{exp}(1))=\omega /10`$ (see Ref. Kleinekathöfer et al. (2001) for details). After rearrangement of Eq. (15) one can easily identify the operators involved in Eq. (1) ($`M=1`$): $`C_1=K`$, $`E_1=\mathrm{\Lambda }`$, $`A=iH_\mathrm{S}K\mathrm{\Lambda }`$. A Gaussian wave packet located at the donor state $`|1`$ and having energy slightly above the crossing of the harmonic potentials was chosen as initial state. The numerical simulation for about 1000 trajectories provides sufficiently converged and accurate results. Fig. 2 shows the relaxation of the ensemble averaged donor population $`\overline{P_1}=\overline{\psi |11|\varphi }+\overline{\varphi |11|\psi }`$. A widely discussed property of the Redfield equation is that it does not conserve positivity May and Kühn (2000). Although $`\overline{P_1}`$ is always positive the tiny negative fraction in Fig. 3 is evidence for the existence of single realizations with negative $`P_1`$. The simulation of the same system within the so-called diabatic-damping approximation Egorova et al. (2001); Kleinekathöfer et al. (2001) with a Lindblad QME by means of the standard quantum jump method Dalibard et al. (1992); Gardiner et al. (1992); Gisin and Percival (1992); Garraway and Knight (1994); Plenio and Knight (1998) keeps all values of $`P_1`$ well confined between $`0`$ and $`1`$.

To summarize, a new method of stochastic unraveling of QMEs beyond the Lindblad form is proposed and thus new large fields of application for stochastic methods are opened. This progress became possible with the use of the wave–function pair in the doubled Hilbert space and the derivation of stable, almost normalized SSEs. The efficiency is determined by the behavior of the norm of every single trajectory. In this sense the jump rates were used as parameters to influence the efficiency. Negative values for the weight of single trajectories allow for the reconstruction of non positive-semidefinite RDMs if required. The method was successfully tested for a simple electron transfer model and for Brownian motion and should allow for better quantum dynamical simulation of large systems. It can also unravel non-Markovian QMEs when they are in a time-local form like in the time-convolutionless formalism Čápek (1994) or in methods using auxiliary density matrices to include the memory effects Meier and Tannor (1999) as well as post-Markov master equations Yu et al. (2000).


text:
# In-plane dipole coupling anisotropy of a square ferromagnetic Heisenberg monolayer

## I Introduction

The investigation of the magnetic properties of ferromagnetic ultrathin films is a field of intense current interest HeB94 . In this contribution we determine the magnetic in-plane anisotropy of a square two-dimensional (2D) Heisenberg ferromagnet in presence of the magnetic dipole interaction. In 2D magnets the action of magnetic anisotropies is twofold. First, they induce a long range magnetization Kit51 ; Mal76 ; LHR99  with a Curie temperature of the order of the exchange coupling Yab91 , whereas a strictly isotropic 2D Heisenberg magnet does not exhibit any long range magnetic order at finite temperatures (Mermin-Wagner-theorem) Mer66 . Moreover, the long range character of the dipole interaction itself is sufficient to stabilize the magnetization in a 2D magnet. Secondly, the anisotropies determine the easy and hard axes of the magnetization with respect to the lattice frame HeB94 . For a ferromagnetic thin film the dipole coupling prefers an in-plane direction of the magnetization (shape anisotropy, demagnetizing field). Except for very small thicknesses thin films are usually in-plane magnetized. The corresponding in-plane easy axis depends on the lattice symmetry of the film face. For a rectangular (110) monolayer the dipole interaction induces a uniaxial in-plane anisotropy, the easy axis being the in-plane axial direction with the smaller lattice constant, as can easily be shown by calculating the corresponding lattice sums.

On the other hand, for a square monolayer the dipole coupling, although itself not rotational invariant, exhibits a continuous energy degeneracy for classical spins at $`T=0`$. Thus, due to the Mermin-Wagner theorem a long range magnetic order is not expected at finite temperatures. However, as has been shown by Monte Carlo calculations and interacting spin wave theory, a magnetic ordering and a critical temperature exists for dipole coupled spins on a square lattice DMB97 , since the magnetic excitations are not continuously degenerate. In this case a quartic in-plane anisotropy is present, the corresponding easy axes being the edges of the square lattice. In other words, the density of states and thus the entropy depends on the magnetic direction within the lattice. This phenomenon is an example of the ’order-by-disorder’ effect in frustrated magnets Vil79 . In these systems thermal fluctuations or lattice disturbances partly remove frustrations, and a collective magnetic ordering emerges at finite temperatures.

Whereas the existence of a finite magnetization and a quartic anisotropy at finite temperatures for dipole coupled spins on a square lattice has been proven, the strength and the temperature behavior of this anisotropy has not been determined yet. Of particular interest is whether the dipole-coupling-induced quartic anisotropy can be measured, and how this anisotropy compares with the single-ion quartic lattice anisotropy due to the spin-orbit coupling. In this study we will adress these questions with the help of a mean field approximation, taking into account a square spin cluster with different numbers of spins (Oguchi-theory) Ogu55 ; Sma66 , as outlined in Sec.II. Furthermore, we apply a many-body Green’s function theory within the Tyablikov-decoupling (RPA) by considering collective magnetic excitations (spin waves) and interactions between them Tya59 ; Tya67 , which are in particular important for 2D magnetic systems. Our results are summarized and discussed in Sec.III.

## II Theory

The free energy $`F(T,\varphi )`$ is determined within a Heisenberg-type Hamiltonian as a function of the temperature $`T`$ and the in-plane angle $`\varphi `$. A square (001) ferromagnetic monolayer is considered, which is spanned by the $`xz`$\- plane, the $`z`$\- direction refers to $`\varphi =0`$. The isotropic nearest neighbor Heisenberg exchange coupling $`J`$, the magnetic dipole-dipole interaction, and an external magnetic field B between localized Heisenberg spins $`𝐒_i`$ on lattice sites $`i`$ are taken into account:

$``$ $`=`$ $`{\displaystyle \frac{J}{2}}{\displaystyle \underset{\mathrm{}i,j\mathrm{}}{}}𝐒_i𝐒_jg\mu _B{\displaystyle \underset{i}{}}𝐁𝐒_i`$ (1)
$`+{\displaystyle \frac{(g\mu _B)^2}{2}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{i,j}{ij}}{}}{\displaystyle \frac{1}{r_{ij}^5}}\left[r_{ij}^2𝐒_i𝐒_j3\left(𝐫_{ij}𝐒_i\right)\left(𝐫_{ij}𝐒_j\right)\right],`$
with $`g`$ the Landé factor and $`\mu _B`$ the Bohr magneton. The spin quantum number is assumed to be $`S=1/2`$. The lattice vector between sites $`i`$ and $`j`$ is given by $`𝐫_{ij}`$, with $`r_{ij}=|𝐫_{ij}|`$ the distance. All spins are assumed to be aligned parallely (ferromagnetic order), the magnetization $`𝐦_i(T)=𝐦(T)=(0,0,m(T))`$ is directed along the magnetic field $`𝐁=(0,0,B)`$.

To account for a varying angle of the magnetization with respect to the lattice frame, the lattice is rotated by the in-plane angle $`\varphi `$ with respect to the $`z`$\- axis, yielding the rotated lattice vectors:

$$𝐫_{ij}^{}=\left(\begin{array}{c}r_{ij}^z\mathrm{sin}\varphi +r_{ij}^x\mathrm{cos}\varphi \\ 0\\ r_{ij}^z\mathrm{cos}\varphi r_{ij}^x\mathrm{sin}\varphi \end{array}\right).$$
(2)
Clearly, the only source in Eq.(1) for a possible in-plane anisotropy is the double scalar product of the dipole interaction. Within the single-spin (Bragg-Williams-) mean field approximation Sma66  the dipole coupling yields an additional contribution to the molecular field. The resulting single-spin Hamiltonian reads

$`^{(1)}`$ $`=`$ $`\left(qJm(T)+wS(0,0)m(T)+g\mu _BB\right)S_1^z`$ (3)
$`+{\displaystyle \frac{1}{2}}m^2(T)\left(qJ+wS(0,0)\right),`$
with $`q=4`$ the number of nearest neighbors (coordination number), $`w=(g\mu _B)^2/a_0^3`$ the strength of the dipole interaction, $`a_o`$ the lattice constant, $`S(0,0)\mathrm{\hspace{0.17em}4.51681}`$ the corresponding lattice sum BeM69 , and $`m(T)=\mathrm{}S_1^z\mathrm{}`$ the magnetization. Thus, within this approximation the dipole interaction is isotropic for a square monolayer and does not depend on the in-plane angle $`\varphi `$.

Our approach will now be improved by applying the so-called Oguchi method Ogu55 ; Sma66 . A number $`N>1`$ of neighboring spins is considered, the interactions between the $`N`$ spins in this cluster are treated exactly. The remaining lattice is coupled to the spin cluster by a molecular field. To conserve the symmetry of the square lattice, we take into account only square-shaped spin clusters, the three smallest possible clusters are characterized by $`N=4,\mathrm{\hspace{0.33em}9}`$, and 16. Note that the resulting cluster Hamiltonians $`^{(N)}`$ are usually non-diagonal, i.e. contributions proportional to the $`S^x`$\- and $`S^y`$\- spin components may be present. Therefore, as a further approximation we consider only the diagonal elements. This approach guarantees that the magnetization is directed along the external magnetic field. As an example, the Hamiltonian for the four-spin cluster ($`N=4`$) reads:

$`^{(4)}`$ $`=`$ $`(Jw)\left(S_1^zS_2^z+S_1^zS_3^z+S_2^zS_4^z+S_3^zS_4^z\right){\displaystyle \frac{w}{2^{5/2}}}(S_1^zS_4^z+S_2^zS_3^z)`$ (4)
$`3w\left[\mathrm{sin}^2\varphi (S_1^zS_2^z+S_3^zS_4^z)+\mathrm{cos}^2\varphi (S_1^zS_3^z+S_2^zS_4^z)+{\displaystyle \frac{2\mathrm{cos}\varphi \mathrm{sin}\varphi }{2^{5/2}}}(S_1^zS_4^zS_2^zS_3^z)\right]`$
$`\left(2Jm(T)+g\mu _BB\right)\left(S_1^z+S_2^z+S_3^z+S_4^z\right)wm(T)(S_1^z+S_4^z)\left(S(0,0)1{\displaystyle \frac{1}{2^{5/2}}}(1+6\mathrm{cos}\varphi \mathrm{sin}\varphi )\right)`$
$`wm(T)(S_2^z+S_3^z)\left(S(0,0)1{\displaystyle \frac{1}{2^{5/2}}}(16\mathrm{cos}\varphi \mathrm{sin}\varphi )\right)+4Jm^2(T)+2wm^2(T)\left(S(0,0)1{\displaystyle \frac{1}{2^{5/2}}}\right).`$
The partition function for this case is given by:

$$Z^{(4)}(T,B,\varphi )=\underset{S_1^z,S_2^z,S_3^z,S_4^z=S}{\overset{+S}{}}\mathrm{exp}(\beta ^{(4)}),$$
(5)
and the (average) magnetization by:

$`m(T,B,\varphi )`$ $`=`$ $`{\displaystyle \frac{1}{4Z^{(4)}(T,B,\varphi )}}\times `$ (6)
$`{\displaystyle \underset{S_1^z,S_2^z,S_3^z,S_4^z=S}{\overset{+S}{}}}(S_1^z+S_2^z+S_3^z+S_4^z)\mathrm{exp}(\beta ^{(4)}),`$
with $`\beta =1/k_BT`$ and $`k_B`$ the Boltzmann constant. The corresponding expressions for the other investigated sizes are rather lengthy and not shown here. The difference of the free energy $`F(T,B,\varphi )=k_BT\mathrm{ln}Z(T,B,\varphi )`$ between the diagonal ($`\varphi =\pi /4`$) and the axial directions ($`\varphi =0`$) yields the in-plane anisotropy $`𝒦_{4,dip}(T,B)`$, which will be calculated as a function of the temperature $`T`$ and the external magnetic field $`B`$.

It is important to mention that the mean field theory as described in the preceeding subsection does not fulfill the Mermin-Wagner-theorem for the isotropic 2D Heisenberg magnet. In this case the long range magnetic order becomes unstable against collective magnetic excitations with long wavelengths Kit51 ; Mal76 ; Yab91 . Thus, the consideration of these spin waves is very important for the magnetic behavior of a ferromagnetic monolayer. For the determination of the dipole-coupling-induced in-plane anisotropy we apply in addition a many-body Green’s function theory to the Hamiltonian Eq.(1), after the lattice is rotated by the in-plane angle $`\varphi `$, Eq.(2). We consider the Green’s functions $`G_{ij}^\alpha =\mathrm{}\mathrm{}S_i^\alpha ;S_j^{}\mathrm{}\mathrm{}`$, $`\alpha =+,,z`$, which are solved within the Tyablikov-decoupling Tya59 . Since interactions between spin waves are partly taken into account, the magnetic properties can be determined up to the Curie temperature by this method. The respective formalism is described in detail in LHR99 , thus we outline merely necessary extensions in the Appendix. The magnetization $`m(T,B,\varphi )`$, the internal energy $`E(T,B,\varphi )`$, and the free energy $`F(T,B,\varphi )`$ are calculated. For comparison, we determine the corresponding quantities also by considering the Holstein-Primakoff approximation HoP40 , which is valid at low temperatures.

## III Results and Discussion

By application of the described methods we calculate the effective dipole-coupling-induced in-plane magnetic anisotropy $`𝒦_{4,dip}(T,B)`$ for a square ferromagnetic monolayer as a function of the temperature $`T`$ and the applied magnetic field $`B`$. A finite value indicates that this anisotropy is caused by the magnetic dipole coupling, since other sources of magnetic anisotropies are not present in the Hamiltonian Eq.(1). For the strength of the dipole interaction we choose $`(g\mu _B)^2/a_o^3=w=0.01J`$, which is appropriate for $`3d`$\- transition metal ferromagnets. $`𝒦_{4,dip}(T,B)`$ is given in units of the energy difference between the perpendicular and the in-plane magnetization (demagnetizing energy) $`E_{demag}(0)=(3/2)wS(0,0)m^2(0)2\pi wS^2`$ for a ferromagnetic monolayer at $`T=0`$.

The free energy $`F(T,\varphi )`$ exhibits a four-fold symmetry as a function of the in-plane angle $`\varphi `$. The easy magnetic axes are directed along the edges ($`\varphi =0,\pi /2`$, etc.), and the hard axes along the diagonals ($`\varphi =\pi /4,3\pi /4`$, etc.) of the square lattice as obtained in DMB97 . The quartic anisotropy $`𝒦_{4,dip}(T,B)`$ is depicted in Fig.1 as a function of the temperature, calculated by the Oguchi approach. Three different sizes of the spin cluster are considered ($`N=4,\mathrm{\hspace{0.33em}9},\mathrm{\hspace{0.33em}16}`$), as well as different magnetic field strengths $`B`$. The Curie temperature $`T_C(N)`$ decreases with an increasing $`N`$, since additional spin correlations are considered. As can be seen, $`𝒦_{4,dip}(T,B)`$ increases with an increasing temperature, reaching a maximum at $`T_C`$. This temperature behavior is in striking contrast to the corresponding behavior of other anisotropic contributions. Usually the effective anisotropies decrease with an increasing temperature, and vanish at $`T_C`$ for $`B=0`$. $`𝒦_{4,dip}(T,B)`$ exhibits a cusp at $`T_C`$ and decreases for temperatures $`T>T_C`$. By increasing the size $`N`$ of the spin cluster, $`𝒦_{4,dip}(T,B)`$ becomes larger. An application of a magnetic field reduces $`𝒦_{4,dip}(T,B)`$, and the cusp near $`T_C`$ changes to a maximum. For comparison, in Fig.1 we show $`𝒦_{4,dip}(T,B)`$ for $`B=0`$ and $`N=4`$ by considering also non-diagonal elements in the cluster Hamiltonian $`^{(4)}`$, cf. Sec.II. A decrease of $`T_C`$ and an increase of $`𝒦_{4,dip}(T,B)`$ by $`10\%`$ is obtained with respect to the calculations which consider diagonal elements only.

In addition, we have calculated $`𝒦_{4,dip}(T,B)`$ by application of the many-body Green’s function theory as outlined in LHR99  and in the Appendix. For $`w=0.01J`$ the Curie temperature is calculated to be $`k_BT_C^{RPA}=0.373J`$. In Fig.2(a,b) we present the resulting magnetization $`m(T,B)`$ and anisotropy $`𝒦_{4,dip}(T,B)`$ as functions of the temperature and the magnetic field. In accordance with the results obtained from the Oguchi approach, $`𝒦_{4,dip}(T,B)`$ increases with an increasing temperature. Due to the consideration of quantum mechanical spins a finite value of $`𝒦_{4,dip}(T,B)`$ is already present at $`T=0`$, resulting from the zero-point spin motion (quantum fluctuations). Note that the presence of the dipole coupling causes a non-saturated magnetization $`m(T=0)<S=1/2`$ Mal76 ; HoP40 . In Fig.2 also the results as obtained from the Holstein-Primakoff approximation are shown, which neglects interactions between spin waves. The corresponding magnetization $`m(T,B)`$ and anisotropy $`𝒦_{4,dip}(T,B)`$ at low temperatures are close to the results as obtained from the RPA. At elevated temperatures $`TT_C/3`$ the Holstein-Primakoff approximation is no longer valid.

As has already been obtained from Monte Carlo calculations and interacting spin wave theory DMB97 , the dipole-coupling-induced in-plane anisotropy is caused by magnetic fluctuations. This order-by-disorder effect is thus controlled by the magnetic entropy of the system Vil79 . An energetic influence should become apparent at $`T=0`$ already within a mean field approximation. Besides the effects resulting from the zero-point motion, the terms of the dipole interaction dependent on the in-plane angle $`\varphi `$ cancel for a strictly square magnetic lattice, i.e. for a fully saturated ferromagnetic state. At finite temperatures thermal fluctuations disturb this square periodicity, causing an effective quartic in-plane magnetic anisotropy, which reflects the underlying spatial lattice symmetry. The single-spin mean field theory conserves the periodicity of the square magnetic lattice, since all other spin operators are replaced by uniform expectation values. Thus, within this approximation one obtains $`𝒦_{4,dip}(T,B)=0`$ for all temperatures. On the other hand, within the Oguchi approach the interactions between a number of spins are treated exactly, considering a few collective magnetic excitations with short wavelengths ranging over several lattice constants, and resulting in a finite value for $`𝒦_{4,dip}(T,B)`$ for $`T>0`$. The many-body Green’s function theory takes into account spin waves with all possible wavelengths. Since $`𝒦_{4,dip}(T,B)`$ is maximal at elevated temperatures, $`𝒦_{4,dip}(T,B)`$ is caused by magnetic fluctuations in particular with short wavelengths, which become excited in this temperature range.

Within the Oguchi approach we calculate a maximum value of $`𝒦_{4,dip}(T,B)`$ of the order of 0.1% of the demagnetizing energy $`E_{demag}(0)`$. On the other hand, $`𝒦_{4,dip}(T,B)`$ as calculated from the Green’s function method is more than a magnitude larger, namely $`𝒦_{4,dip}(T,B)\mathrm{\hspace{0.17em}0.01}E_{demag}(0)`$ near $`T_C`$. Similarly, the action of a magnetic field $`B`$ on the anisotropy $`𝒦_{4,dip}(T,B)`$ is calculated to be much more pronounced within the Green’s function theory than within the Oguchi approach, cf. Figs.1 and 2. The reason for the strong differences between the results obtained from the two methods is that the former one considers collective magnetic excitations with long wavelengths. These are known to have strong effects on the magnetic properties of ultrathin films Mal76 ; Yab91 . In this case the magnetic field acts merely on spin blocks, i.e. correlated regions of neighboring spins characterized by the short range order parameter Pol75 , resulting in a much larger magnetic response. For instance, the induced magnetization in a ferromagnetic trilayer has been determined to be an order of magnitude larger by consideration of spin waves JBP99 . The collective magnetic excitations are most pronounced for a single magnetic layer Mal76 , the zero-point motion is strongest for the spin number $`S=1/2`$. Note that the theoretical methods applied in the present study are much less demanding than Monte Carlo simulations DMB97 .

The obtained dipole-coupling-induced in-plane anisotropy is small as compared to other (effective) anisotropies, since we have considered interaction strengths appropriate for $`3d`$\- transition metal ferromagnets. A small value for $`𝒦_{4,dip}(T,B)`$ has been conjectured in DMB97 . If a ten times larger ratio $`w/J`$ between the dipole coupling strength and the exchange interaction is assumed, the ratio $`𝒦_{4,dip}(T,B)/E_{demag}(0)`$ increases roughly by the same factor. Since usually the effective anisotropies decrease with an increasing temperature, $`𝒦_{4,dip}(T,B)`$ might become observable at elevated temperatures. For comparison, we show in Figs.1 and 2 also the effective single-ion quartic in-plane anisotropy $`𝒦_{4,so}(T)`$ resulting from the spin-orbit interaction. This quantity is calculated with the help of a thermodynamic perturbation theory Cal66 . By assuming its strength to be $`𝒦_{4,so}(0)=0.01E_{demag}(0)`$, and adapting corresponding Curie temperatures, $`𝒦_{4,so}(T)`$ becomes comparable to $`𝒦_{4,dip}(T,B)`$ at $`T/T_C0.80.9`$. In this temperature range the total quartic in-plane anisotropy should increase again, or exhibit an in-plane magnetic reorientation. As mentioned, this dipole-coupling-induced anisotropy should be more apparent for ferromagnetic (001) thin films with a large $`w/J`$\- ratio.

We like to comment on the fact that a finite anisotropy is obtained for $`T>T_C`$. Usually the effective anisotropy as observed for a collectively ordered ferromagnetic state vanishes above $`T_C`$. However, a vanishing effective anisotropy for $`T>T_C`$ does not indicate that the anisotropy as present in the Heisenberg Hamiltonian disappears. A single (paramagnetic) spin is still subject to e.g. the single-ion uniaxial anisotropy $`K_2`$ also if the net magnetization is zero. The resulting free energy difference between the easy and the hard magnetic directions (’paramagnetic anisotropy’) behaves as $`K_2^2/k_BT`$ for $`K_2k_BT`$, and is rather small if $`K_2`$ is small as compared to the exchange interaction $`Jk_BT_C`$. With regard to the present study, a finite value of $`𝒦_{4,dip}(T,B=0)`$ for $`T>T_C`$ reflects the dipole-coupling-induced anisotropy of a spin block Pol75 .

Note that the free energy $`F(T,B)`$ as obtained from the Green’s function method yields unphysical results for temperatures near and above $`T_C`$. As can be shown Jenup , for large temperatures the free energy $`F(T,B)`$ as calculated by this method does not approach the value $`k_BT\mathrm{ln}2`$, which is the free energy of a single, non-interacting $`S=1/2`$\- spin. This unphysical behavior results in a still increasing $`𝒦_{4,dip}(T,B)`$ for $`T>T_C`$, cf. Fig.2(b). Nevertheless, we expect that the Green’s function approach reflects the correct behavior of $`𝒦_{4,dip}(T,B)`$ for $`T<T_C`$, since its behavior is corroborated by the Oguchi method and the Holstein-Primakoff approach.

We have considered a square monolayer only. A similar effect is expected also for thicker films with a square (001) face, or for a hexagonal (111) thin film. In addition, a corresponding effect may emerge for a three-dimensional cubic lattice, for which the dipole interaction cancels exactly. At elevated temperatures, however, caused by thermal fluctuations the frustration due to this periodicity will be lifted. The dipole coupling then induces a cubic magnetic anisotropy with easy axes directed along e.g. the edges of the cubic lattice. Also is this case an increasing anisotropy with an increasing temperature is expected. To obtain an anisotropic contribution caused by the magnetic dipole interaction via the order-by-disorder effect, we emphasize that the lattice symmetry must be strictly cubic or square. A lattice distortion e.g. due to magnetostrictive effects will probably blur this effect.

P.J.J. gratefully acknowledges the invitation and the hospitality of the I. Institute for Theoretical Physics of the Hamburg University. Fruitful discussions with Prof. K.D. Schotte, Free University of Berlin, are acknowledged.

## Appendix: Many-body Green’s function theory

The following commutator Green’s functions in the frequency space are taken into account: $`G_{ij}^\alpha (\omega )=\mathrm{}\mathrm{}S_i^\alpha ;S_j^{}\mathrm{}\mathrm{}_\omega `$, $`\alpha =+,,z`$. By considering the Hamiltonian Eq.(1), these Green’s functions are solved in the usual way by the equation of motion LHR99 ; Tya59 ; Tya67 . The higher order Green’s functions are approximated by the Tyablikov-decoupling (RPA) Tya59 , e.g. $`\mathrm{}\mathrm{}S_k^zS_i^+;S_j^{}\mathrm{}\mathrm{}\mathrm{}S^z\mathrm{}\mathrm{}\mathrm{}S_i^+;S_j^{}\mathrm{}\mathrm{}`$, resulting in a system of linear equations for the $`G_{ij}^\alpha (\omega )`$. By application of a Fourier transformation into the 2D wave vector space $`𝐤_{}𝐤`$, one obtains $`G^z(𝐤,\omega )=0`$, if $`\mathrm{}S^+\mathrm{}=\mathrm{}S^{}\mathrm{}=0`$. The remaining two Green’s functions $`G^\pm (𝐤,\omega )`$ are easily evaluated LHR99 :

$`G^+(𝐤,\omega )`$ $`=`$ $`2m(T){\displaystyle \frac{\omega +a(𝐤)}{\omega ^2\epsilon ^2(𝐤)}},`$ (7)
$`G^{}(𝐤,\omega )`$ $`=`$ $`2m(T){\displaystyle \frac{b(𝐤)}{\omega ^2\epsilon ^2(𝐤)}},`$ (8)
with $`m(T)=\mathrm{}S^z\mathrm{}`$. The magnon dispersion relation $`\epsilon (𝐤)`$ is given by

$`\epsilon ^2(𝐤)`$ $`=`$ $`a^2(𝐤)b^2(𝐤),`$ (9)
$`a(𝐤)`$ $`=`$ $`g\mu _BB+m(T)\left(J(0)J(𝐤)\right)`$
$`+wm(T)[S(0,0)+S(k_z,k_x)(1{\displaystyle \frac{3}{2}}\mathrm{cos}^2\varphi )`$
$`+S(k_x,k_z)(1{\displaystyle \frac{3}{2}}\mathrm{sin}^2\varphi )+3T(k_z,k_x)\mathrm{cos}\varphi \mathrm{sin}\varphi ],`$
$`b(𝐤)`$ $`=`$ $`{\displaystyle \frac{3}{2}}wm(T)[S(k_z,k_x)\mathrm{cos}^2\varphi +S(k_x,k_z)\mathrm{sin}^2\varphi `$ (11)
$``$ $`2T(k_z,k_x)\mathrm{cos}\varphi \mathrm{sin}\varphi ].`$
$`J(𝐤)=2J(\mathrm{cos}k_xa_0+\mathrm{cos}k_za_0)`$ is the Fourier transform of the exchange interaction, with $`a_o`$ the lattice constant. The oscillating lattice sums BeM69  are defined by

$`S(k_x,k_z)`$ $`=`$ $`{\displaystyle \underset{l,n=\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{l^2}{(l^2+n^2)^{5/2}}}\mathrm{exp}(ik_xa_0lik_za_0n),`$
$`T(k_z,k_x)`$ $`=`$ $`{\displaystyle \underset{l,n=\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{ln}{(l^2+n^2)^{5/2}}}\mathrm{exp}(ik_za_0lik_xa_0n),`$
where the terms with $`l=n=0`$ have to be omitted. The following expressions are valid for the spin quantum number $`S=1/2`$. The magnetization $`m(T)`$ is obtained from the spectral theorem Tya67 :

$`m(T,\varphi )`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}\mathrm{}S^{}S^+\mathrm{}(𝐤)`$ (13)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{m(T)}{N}}{\displaystyle \underset{𝐤}{}}\left({\displaystyle \frac{a(𝐤)}{\epsilon (𝐤)}}\mathrm{coth}x1\right),`$
with $`x=\epsilon (𝐤)/2k_BT`$. The summation extends over the first Brillouin zone, $`N`$ denotes the number of $`𝐤`$-points. By calculating the expectation value $`\mathrm{}S_i^zS_j^z\mathrm{}`$ with the help of $`G^+(𝐤,\omega )`$ Tya67 , we obtain for the internal energy $`E(T,\varphi )`$ per spin:

$`E(T,\varphi )`$ $`=`$ $`E_0+{\displaystyle \frac{m(T,\varphi )}{2N}}{\displaystyle \underset{𝐤}{}}(e_1(𝐤)({\displaystyle \frac{a(𝐤)}{\epsilon (𝐤)}}\mathrm{coth}x1)`$ (14)
$`+(e_2(𝐤){\displaystyle \frac{b(𝐤)}{\epsilon (𝐤)}}+\epsilon (𝐤))\mathrm{coth}xa(𝐤)+b(𝐤)),`$
with the denotations

$`E_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}g\mu _BB{\displaystyle \frac{1}{8}}\left(J(0)+wS(0,0)\right),`$ (15)
$`e_1(𝐤)`$ $`=`$ $`g\mu _BB+{\displaystyle \frac{1}{2}}\left(J(0)J(𝐤)\right)`$ (16)
$`+{\displaystyle \frac{w}{2}}\left(S(0,0)+S(k_z,k_x)+S(k_x,k_z)\right),`$
$`e_2(𝐤)`$ $`=`$ $`{\displaystyle \frac{1}{2}}J(𝐤){\displaystyle \frac{w}{2}}\left(S(k_z,k_x)+S(k_x,k_z)\right).`$ (17)
Finally, the free energy $`F(T,\varphi )`$ is calculated from a temperature integral over $`E(T,\varphi )`$ Tya67 . The Holstein-Primakoff approximation HoP40  is obtained by replacing the magnetization $`m(T)`$ on the right side of Eq.(13) by its saturation value $`S=1/2`$.


text:
# Spectrum of screening masses in the (3+1)D SU(2) pure gauge theory near the critical temperature Talk presented by A. Papa.

## 1 INTRODUCTION

It is well known that both the 3D Ising model

$$S=\beta \underset{<n,m>}{}s_ns_m$$

and the 3D lattice regularized $`\varphi ^4`$ theory

$$S=\beta \underset{<n,m>}{}\varphi _n\varphi _m+\underset{n}{}\varphi _n^2+\lambda \underset{n}{}(\varphi _n^21)^2$$

undergo a second order phase transition at $`\beta =\beta _c`$ from a phase where the Z(2) symmetry is broken ($`\beta >\beta _c`$) to a phase when the Z(2) symmetry is restored ($`\beta <\beta _c`$). The order parameter of this transition is the expectation value of the fundamental field, i.e. the magnetization.

The spectrum of massive excitations in the broken symmetry phase of 3D Ising model and $`\varphi ^4`$ theory can be determined from the analysis of wall-wall correlations of spin operators. By performing Monte Carlo simulations near $`\beta _c`$, it has been possible to detect three $`0^+`$ states . Results for their mass ratios agree in the two models: $`m_2/m_1=1.83(3)`$ and $`m_3/m_1=2.45(10)`$ ($`m_1`$ is the fundamental mass). Since the 3D Ising model and the 3D (lattice) $`\varphi ^4`$ theory belong to the same universality class, this result shows that universality is valid far beyond critical indices.

The first excitation ($`m_2`$) is below the two particle production threshold, it represents therefore a non-perturbative state, identified as a bound state of the fundamental excitation . The second mass state ($`m_3`$) could be a fictitious one which mimics the two-particle cut.

According to duality, the broken phase of the 3D Ising model can be put in correspondence with the confined phase of the 3D Z(2) gauge theory. This implies that mass ratios in the glueball spectrum in the confined phase of 3D Z(2) gauge theory should agree with mass ratios in the broken phase of 3D (spin) Ising model. This has been verified in Ref.  using Monte Carlo results  for the $`0^+`$ glueball masses in the 3D Z(2) gauge theory: $`m_2/m_1=1.88(2)`$ and $`m_3/m_1=2.59(4)`$.

According to universality, the broken phase of 3D Ising model corresponds to the deconfined phase of the (3+1)D SU(2) pure gauge theory at finite temperature , with the magnetization of the spin model mapped into the Polyakov loop, $`P(x,y,z)\text{Tr}_{n_4=1}^{N_t}U_4(x,y,z,an_4)`$. This implies that the same mass ratios observed in the broken phase of 3D Ising should appear in the spectrum of the slopes of the exponential decays of wall-wall correlations of Polyakov loop operators, i.e. in the spectrum of screening masses, in the deconfined phase of (3+1)D SU(2) at finite temperature near criticality. In this work we present a numerical test of this prediction.

## 2 NUMERICAL RESULTS

A naive way to determine $`0^+`$ screening masses is to study the connected wall-wall correlation function

$$G(|z_1z_2|)\overline{P}(z_1)\overline{P}(z_2)\overline{P}(z_1)\overline{P}(z_2)$$

$$=A(e^{m|z_1z_2|}+e^{m(N_z|z_1z_2|)})$$

$$+B(e^{m^{}|z_1z_2|}+e^{m^{}(N_z|z_1z_2|)})+\mathrm{}.$$

Here $`\overline{P}(z)1/(N_xN_y)_{x,y}P(x,y,z)`$ is the wall-average of Polyakov loops, $`m`$ is the lowest mass in the $`0^+`$ channel, $`m^{}`$ the first excited level in the $`0^+`$ channel and the dots represent the contribution from higher mass excitations and from multi-particle cuts. Since $`\overline{P}(z)0`$ in the high temperature phase, in principle also a $`z`$-independent term in $`G(|z_1z_2|)`$ of the form $`C_0\mathrm{exp}[mN_z]`$ should be included (see ). However, such term is evidently sub-dominant for $`m^{}<2m`$, as we expect to be the case here.

We performed Monte Carlo simulations using the overheat-bath updating algorithm  with Kennedy-Pendleton improvement . We studied the naive wall-wall correlation at $`\beta =2.33`$ on a $`18^2\times 36\times 4`$ lattice (statistics 1M) and at $`\beta =2.36`$ on a $`18^3\times 4`$ (statistics 1.5M). These values of $`\beta `$ are slightly above the critical value $`\beta _c=2.29895(10)`$ determined in Ref. . We determined the effective mass defined as

$$G(z)=C(e^{m_{\text{eff}}z}+e^{m_{\text{eff}}(N_zz)})$$

and found in both cases that data for $`m_{\text{eff}}`$ as a function of $`z`$ reach a plateau value corresponding to the fundamental mass: $`m=0.3660(20)`$ at $`\beta =2.33`$ and $`m=0.5356(26)`$ at $`\beta =2.36`$. The deviations from the plateau at the smallest values of $`z`$ can be attributed to lattice artifacts and to the possible effect of other physical states. In order to single out the latter effect, we rescaled the values of $`m_{\text{eff}}`$ and $`z`$ by $`m1/\xi `$ and $`\xi `$, respectively, and put together on the same plot data from $`\beta =2.33`$ and $`\beta =2.36`$. As shown in Fig. 1, data from the two different $`\beta `$’s fall on the same curve before reaching the plateau, indicating the physical origin of these effects.

In order to detect this excited state, we adopted the so-called “variational” method . The method includes the following steps: define a basis of (wall-averaged) operators $`\{O_\alpha \}`$, compute the connected cross-correlation matrix among them

$$C_{\alpha \beta }(z)=O_\alpha (z)O_\beta (0)O_\alpha (z)O_\beta (0),$$

diagonalize $`C_{\alpha \beta }(z)`$ to obtain the eigenvalues

$$\lambda _i(z)e^{m_iz}+e^{m_i(N_zz)}$$

and finally extract the masses $`m_i`$. In practical simulations not more than two (sometimes three) leading eigenvalues allow to extract a signal. The effectiveness of the method relies on the choice of a “good” set of operators; it is convenient to define operators living on different length scales (for instance, by use of recursive definitions). Moreover, by defining operators with non-trivial transformation under spatial rotation, it is possible to look for states with non-zero angular momentum.

As a first set of operators we considered the one which can be built adapting to the present case the recursive procedure used in Ref. :

$$P^{(0)}(x,y,z)=P(x,y,z),$$

$$P^{(n+1)}(x,y,z)=\text{sign}(u)\left[(1w)|u|+wy\right],$$

$`u`$ $`=`$ $`{\displaystyle \frac{1}{4}}(P^{(n)}(xa,y,z)+P^{(n)}(x,ya,z)`$
$`+`$ $`P^{(n)}(x+a,y,z)+P^{(n)}(x,y+a,z)),`$
with $`w=0.1`$ and $`y=P(x,y,z)`$. We considered a set of six such (wall-averaged) operators corresponding to the smoothing steps 0, 3, 6, 9, 12 and 15. The effective masses extracted by the variational method applied to these set of operators are shown in Fig. 2 which shows a clear plateau at a value $`m(0^+)=0.3629(50)`$, in agreement with the naive determination of the fundamental mass. Data coming from the next-to-leading eigenvalue and corresponding to the first excited level are quite noisy. Nevertheless, if we take as “plateau” value the first value of the effective mass which agrees at 1$`\sigma `$ level with the preceding and the following value, we get $`m^{}(0^+)=0.684(48)`$ and a mass ratio $`m^{}(0^+)/m(0^+)=1.89(16)`$ in good agreement with the prediction from universality. Further analysis has shown that this excited state is coupled with the operators with larger scale in the set we considered. Indeed, taking into account only two operators (smoothing level 0 and 15) the same result is obtained. Instead, performing the analysis taking only the first four operators (smoothing levels 0, 3, 6, 9) gives $`m(0^+)=0.3615(40)`$ and $`m^{}(0^+)=0.848(28)`$, with a mass ratio $`m^{}(0^+)/m(0^+)=2.35(10)`$ which is compatible with the value of $`m_3/m_1`$ in the 3D Ising model and could be the analog of the two-particle cut in the same model.

Finally, in order to look for states with non-zero momentum, we considered the following set of $`2^+`$ operators, inspired by Ref. ,

$`\overline{P}^{(n)}(z)`$ $`=`$ $`{\displaystyle \frac{1}{N_xN_y}}{\displaystyle \underset{x,y}{}}[P(x,y,z)P(x+na,y,z)`$
$``$ $`P(x,y,z)P(x,y+na,z)],`$
with $`n=1,\mathrm{},5`$. The variational method gave for the fundamental $`2^+`$ mass the value $`m(2^+)=1.23(14)`$. The ratio of this value to the fundamental mass in the $`0^+`$ channel (from the naive determination) $`m(0^+)=0.3660(20)`$ gives $`m(2^+)/m(0^+)=3.36(40)`$. This result can be compared with the corresponding mass ratio in the glueball spectrum of the 3D gauge Ising model which, by duality, is connected to the spectrum of massive excitations in the broken phase of the 3D (spin) Ising model. In the 3D gauge model we have $`m(2^+)/m(0^+)=2.59(4)`$ and $`m^{}(2^+)/m(0^+)=3.23(7)`$, so our result for the lowest screening mass in the $`2^+`$ channel seems to be in agreement with the first excited mass $`m^{}(2^+)`$, more than with the fundamental $`m(2^+)`$.