text:
# 1 Introduction

## 1 Introduction

The concept of universality is one of the fundamental principles in the theory of critical phenomena. According to the universality hypothesis, the long-range behavior of a statistical system near a phase transition is governed by fluctuations of order parameters and is essentially determined by symmetries, no matter how complicated the underlying microscopic dynamics is. For systems with $`_2`$ symmetry universality leads to the description in terms of the one-component $`\varphi ^4`$ field theory (see e.g. Refs. ).

It is known that an order parameter in three-dimensional systems does not fluctuate strongly and the mean-field approximation (Landau theory) works reasonably well. The systematic methods that start from the mean-field approximation and perturbatively take into account fluctuative corrections give very accurate predictions for universal quantities such as critical exponents. However, a detailed numerical study of correlation functions in the scaling regime of various systems with $`_2`$ symmetry revealed deviations from the mean-field approximation which are not too large numerically, but at the same time have no qualitative explanation within perturbation theory. These deviations arise in the hierarchy of length scales in the phase with broken $`_2`$ symmetry.

The only macroscopic scale in the problem, of course, is the correlation length $`\xi `$. A generic correlator decreases as $`e^{r/\xi }`$ at infinity, but there are some correlation functions that decrease more rapidly. Composite operators in the Ising model whose correlation functions decay faster than $`e^{r/\xi }`$ were constructed in  and include, as typical examples, operators of non-zero angular momentum. What are the length scales associated with these operators? The answer is simple in the Gaussian approximation, when correlators factorize: a correlation function of composite operators that factorizes on $`n`$ irreducible correlators behaves as $`(e^{r/\xi })^n`$. The characteristic decay rate of any correlation function is therefore an integer multiple of the inverse correlation length.

The fact that fluctuations are not exactly Gaussian should not change this conclusion. The ratio of the correlation length to any other length scale is an integer to any order in perturbation theory! This can be easily understood from the field theory perspective. The behavior of a correlation function at infinity is determined by the singularities of its Fourier transform in the complex momentum plane. These singularities are associated with particles described by the quantum field. The only singularities that arise in perturbative $`\varphi ^4`$ theory are the single-particle pole associated with the elementary quantum and the multiparticle cuts. Correlation functions of operators that couple to the single-particle state decay as $`e^{\mu r}`$, where $`\mu `$ is the inverse correlation length, or the particle mass. Correlation functions of operators that have zero overlap with the one-particle state (because of the angular momentum conservation, for instance) do not have poles in the complex momentum plane. Their only singularities are multiparticle branch cuts that start at $`p^2=n^2\mu ^2`$. Consequently, such correlators decay as $`e^{n\mu r}`$ with some integer $`n`$.

The numerical simulations clearly show that the above picture does not hold in the phase with broken $`_2`$ symmetry. The firmest evidence comes from the existence of a spin zero state with mass just below the two-particle threshold: $`M1.8\mu `$, strongly supported by numerical results . This and other non-perturbative states give visible contributions to the quantities which are sensitive to distances of order of the correlation length. For instance, the spin-spin correlator in the Ising model near the critical point cannot be fitted by a simple pole-plus-cut structure suggested by perturbation theory and definitely contains an extra pole below the cut .

The simulations show that dimensionless ratios like $`M/\mu `$ become constant when approaching the critical temperature. Moreover they have the same value within statistical errors in two different realizations of the universality class: the Ising model itself and the lattice regularized $`\varphi ^4`$ theory. These two facts allow us to be confident about the universal character of these effects. Specifically, we expect the same spectrum to appear in all realizations of the universality class, as long as one moves away from the critical point by perturbing with a $`_2`$-symmetric operator towards the region where the symmetry is spontaneously broken.

To resolve the aforementioned contradiction, we proposed to interpret the non-perturbative states in the broken symmetry phase of $`\varphi ^4`$ theory as bound states of elementary quanta of the scalar field . We shall develop this idea further in this paper.

Let us briefly recollect some facts about bound states. Two-particle bound states show up as poles in the Fourier transform of correlation functions below the branch cut associated with the two-particle threshold. In three dimensions, the threshold singularity is logarithmic, and, though the two-particle cut comes from one-loop Feynman diagrams and thus is accompanied by a coupling constant, the smallness of the coupling is compensated by the large logarithm. In the vicinity of the threshold, perturbation theory breaks down. Fortunately, the large logarithms can be resummed. This resummation substantially changes an analytic structure of correlation functions near the threshold and produces a pole that is invisible in perturbation theory. This pole is associated with a two-particle bound state. Its position is determined by the Bethe-Salpeter (BS) equation, which can be solved order by order in perturbation theory.

In this work we analyze the BS equation for $`3D`$ $`\varphi ^4`$ theory in the broken symmetry phase. The analysis shows that a bound state of the elementary quanta indeed exists. We compute its binding energy at the next-to-leading order. While the value we obtain at the leading order is in very good agreement with the Monte Carlo data, the series appear to be badly behaved, so that a reliable numerical estimate of the subleading corrections to the binding energy would require the knowledge of a large number of terms in the series.

The $`3D`$ $`_2`$ gauge theory is related to the Ising model by an exact duality transformation. Since the duality transformation maps the confining phase of the gauge theory into the broken symmetry phase of the spin model, all of our considerations apply to the glueball spectrum of the gauge theory as well. Indeed, we shall prove that duality implies the exact equality between the spectra of the spin model in the broken symmetry phase and the gauge theory in the confining phase for all values of the coupling.

This correspondence suggests a novel interpretation of the glueballs of the gauge theory as bound states of the dual spin model, which in turn allows one to understand the qualitative features of the glueball spectrum, and in particular its peculiar angular momentum dependence. Bound states with different values of the angular momentum are easily seen to appear in the low-temperature expansion of the $`3D`$ Ising model. Simple group-theoretic arguments determine how many elementary quanta are necessary to construct a bound state of any given angular momentum and parity. Assuming that the binding energy is always much smaller than the mass of the elementary components, and that this description of the spectrum in terms of bound states is conserved when one moves from the deep low-temperature region to the scaling region, one ends up with a prediction for the qualitative features of the glueball spectrum in $`_2`$ gauge theory, and in particular for the angular momentum dependence of the various masses. This prediction is in very good agreement with the numerical data for the glueball spectrum, thus lending support to the interpretation of glueballs as bound states of the dual spin model.

The rest of the paper is organized as follows: in Sec. 2 we will present the detailed computations for the bound state energy in the framework of the BS equation for $`3D`$ $`\varphi ^4`$ field theory; in Sec. 3 we study the dependence of the binding energy on the magnetic field. Sec. 4 contains a proof of the exact correspondence between the spectrum of the $`3D`$ Ising model in the broken symmetry phase and the $`_2`$ gauge theory in the confining phase. In Sec. 5 we describe in detail the argument, already presented in , that allows us to explain the angular momentum dependence of the glueball spectrum of the $`_2`$ gauge theory using the dual bound state picture. In Sec. 6 we discuss obtained results and outline possible future developments.

## 2 Bound states in $`3D`$ $`\varphi ^4`$: BS equation

We consider $`3D`$ $`\varphi ^4`$ theory in the broken symmetry phase. The Euclidean action is

$$S=d^3x\left[\frac{1}{2}_\mu \varphi _\mu \varphi +\frac{g}{4!}(\varphi ^2v^2)^2\right].$$
(2.1)
The field $`\varphi `$ acquires an expectation value and the perturbative expansion is performed around the stable vacuum $`\varphi =v`$. Extracting the fluctuating part,

$$\varphi =v+\phi ,$$
(2.2)
and setting

$$m^2=\frac{gv^2}{3}$$
(2.3)
we write the action as

$$S=d^3x\left[\frac{1}{2}_\mu \phi _\mu \phi +\frac{1}{2}m^2\phi ^2+\frac{1}{3!}\sqrt{3gm^2}\phi ^3+\frac{1}{4!}g\phi ^4\right].$$
(2.4)
The quartic coupling $`g`$ has a dimension of mass, so the dimensionless expansion parameter is the ratio $`g/m`$.

If the coupling is weak, the bound state lies very close to the threshold. We shall see that the binding energy is exponentially small in $`g/m`$. While the exponent can be easily found , computing the pre-exponential factor requires sufficiently complicated next-to-leading order calculations. The calculations are to certain extent similar to those of Ref.  where the BS equation was solved to the first few orders in perturbation theory for (1+1)-dimensional scalar field theory with the cosine interaction.

### 2.1 Leading order solution

The two-particle bound state manifests itself as a pole in the four-point function. The pole arises after summation of the ladder diagrams (Fig. 1a). Each ladder diagram develops a kinematic singularity when the momentum flowing through the diagram is near the two-particle threshold: $`p^24m^2`$, since the momentum integration in each loop then contains a region where two of the intermediate states simultaneously go on-shell and the propagators in the loop integral simultaneously have poles. This singularity produces large logarithms that compensate for powers of the coupling, so near the threshold all ladder diagrams are of the same order. The vertices in Fig. 1a correspond to the sum of the two-particle irreducible (2PI) diagrams, that is, the diagrams that cannot be split into disconnected pieces by cutting two internal lines whose total momentum is $`p`$. 2PI diagrams are analytic at the two-particle threshold in the $`s`$-channel.

The position of the pole is determined by the BS equation, graphically represented in Fig. 1b. The BS equation is the homogeneous part of the Dyson equation for ladder diagrams :

$$\mathrm{\Gamma }(p_1,p_2)=\frac{d^3q}{(2\pi )^3}\mathrm{\Gamma }(q,pq)G(q)K(q,pq;p_1,p_2)G(pq),p=p_1+p_2.$$
(2.5)
Here, $`G`$ is the propagator and $`K`$ is the BS kernel, the aforementioned sum of 2PI diagrams. The BS equation is a linear eigenvalue problem for the wave function $`\mathrm{\Gamma }`$ and the bound state mass:

$$M^2p^2.$$
(2.6)
We take:

$$M=(2\mathrm{\Delta })m,$$
(2.7)
where $`\mathrm{\Delta }`$ is small.

To the leading order in $`g/m`$, the BS kernel is given by the sum of three diagrams in Fig. 2. The diagram $`a`$ contributes $`g/2`$ to the kernel:

$$K_a=\frac{g}{2}.$$

The diagram $`b`$ gives

$$K_b=\frac{3gm^2}{2}\frac{1}{m^2M^2}=\frac{g}{2}+O(\mathrm{\Delta }).$$

The contribution of the diagrams $`c`$ and $`d`$ depends on the momentum transfer in the $`t`$ and $`u`$ channels. However, the momentum transfer is very small when external momenta and both of the propagators in (2.5) are on shell: simple kinematic arguments show that the typical momentum in the $`t`$ or $`u`$ channel propagators is of order $`\mathrm{\Delta }m`$. To the first approximation, we can just set it to zero:

$$K_c=\frac{3gm^2}{2}\frac{1}{(p_1q)^2+m^2}\frac{3g}{2},$$

$$K_d=\frac{3gm^2}{2}\frac{1}{(p_2q)^2+m^2}\frac{3g}{2}.$$

Collecting all three terms together, we get:

$$K=2g.$$
(2.8)
The kernel corresponds to a local four-point interaction in this approximation. The effective quartic coupling appears to be negative. This happens because the exchange diagram $`c`$ dominates over the pure four-point vertex $`a`$ and the annihilation graph $`b`$, so the interaction is attractive. No matter how weak the interaction is, it will bind elementary quanta into bound states.

Since the binding energy is small, the non-relativistic approximation should be valid at small $`g/m`$. The non-relativistic limit of the BS equation with the positive local kernel is the two-body Schrödinger equation with attractive $`\delta `$-function potential. This quantum-mechanical problem has a bound state at any value of the coupling with exponentially small binding energy . An approach based on the non-relativistic approximation was elaborated in . This approach is simple and physically transparent, but it only allows us to find the binding energy up to a numerical factor. Here we follow a more systematic route based on the full BS equation.

The BS equation with the local kernel Eq. (2.8) corresponds to summing an infinite series of bubble diagrams, and can be written as

$$1=2g\frac{d^3q}{(2\pi )^3}\frac{1}{(q^2+m^2)[(pq)^2+m^2]}.$$
(2.9)
Evaluating the integral,

$`{\displaystyle \frac{d^3q}{(2\pi )^3}\frac{1}{(q^2+m^2)[(pq)^2+m^2]}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{arcsin}\left(\sqrt{\frac{p^2}{4m^2+p^2}}\right)}{4\pi \sqrt{p^2}}}={\displaystyle \frac{\mathrm{ln}\left(\frac{2m+M}{2mM}\right)}{8\pi M}}`$ (2.10)
$`=`$ $`{\displaystyle \frac{\mathrm{ln}\left(\frac{4}{\mathrm{\Delta }}\right)}{16\pi m}}+O(\mathrm{\Delta }\mathrm{ln}\mathrm{\Delta }),`$
we obtain:

$$1=2g\frac{\mathrm{ln}\left(\frac{4}{\mathrm{\Delta }}\right)}{16\pi m}.$$
(2.11)
This equation predicts an exponentially small binding energy:

$$\mathrm{\Delta }=\mathrm{const}\mathrm{exp}\left(\frac{8\pi m}{g}\right),$$
(2.12)
in accord with the result obtained in the non-relativistic approximation ***The coupling $`\lambda `$ used in  differs from $`g`$ by a factor of 24..

The normalization of the binding energy cannot be deduced from Eq. (2.11). We cannot keep the constant term in (2.10), because it has the same magnitude as logarithmic terms in the next order of perturbation theory, which are proportional to $`g^2/m^2\mathrm{ln}\mathrm{\Delta }g/m`$. There are other sources of corrections of order $`g/m`$, as well. We need to compute all of them to fix the normalization of the binding energy.

### 2.2 Corrections to propagator

One-loop corrections to the propagator are given by three diagrams in Fig. 3. They shift the pole of the propagator to

$$\mu ^2=m^2\frac{3gm}{16\pi }\mathrm{ln}3g\mathrm{\Lambda },$$
(2.13)
where by $`\mathrm{\Lambda }`$ we denote the linearly divergent integral

$$\mathrm{\Lambda }\frac{d^3q}{(2\pi )^3}\frac{1}{q^2+m^2}.$$
(2.14)
Expanding the third diagram in Fig. 3 near the mass shell, we get:

$$G^1(k)=(p^2+\mu ^2)\left[1+\left(3\frac{9}{4}\mathrm{ln}3\right)\frac{g}{24\pi m}\right]+O\left((p^2+\mu ^2)^2\right)$$
(2.15)
and

$$G(k)=\frac{1}{p^2+\mu ^2}\left[1\left(3\frac{9}{4}\mathrm{ln}3\right)\frac{g}{24\pi m}\right]+O(1).$$
(2.16)
In evaluating the integral (2.10), the corrections to the propagator will result in the replacement of the bare mass $`m`$ by the renormalized (physical) mass $`\mu `$ and in the multiplication of the whole integrand by the wave function renormalization factor.

It is natural to express everything in terms of the physical mass gap (2.13). In particular, we set:

$$M=(2\mathrm{\Delta })\mu .$$
(2.17)
This definition affects the BS kernel, since the contribution of the diagram $`b`$ in Fig. 2 depends on how exactly we define $`M`$:

$`K_b`$ $`=`$ $`{\displaystyle \frac{3gm^2}{2}}{\displaystyle \frac{1}{m^2M^2}}={\displaystyle \frac{3g}{2}}{\displaystyle \frac{m^2}{\mu ^2M^2}}\left[1\left({\displaystyle \frac{3gm}{16\pi }}\mathrm{ln}3+g\mathrm{\Lambda }\right){\displaystyle \frac{1}{\mu ^2M^2}}\right]`$ (2.18)
$`=`$ $`{\displaystyle \frac{g}{2}}{\displaystyle \frac{g^2}{8\pi m}}\mathrm{ln}3{\displaystyle \frac{2g^2\mathrm{\Lambda }}{3m^2}}+O(\mathrm{\Delta })+O(g^3).`$
This differs from the result based on our previous definition by an amount

$$\delta K_b=\frac{g^2\mathrm{ln}3}{8\pi m}\frac{2g^2\mathrm{\Lambda }}{3m^2}.$$
(2.19)
### 2.3 One-loop corrections to the kernel

The one-loop corrections to the kernel of the BS equations are given by thirteen diagrams in Fig. 4. Since we are interested only in the logarithmically enhanced terms in the BS equation, all these diagrams can be evaluated on shell and at zero momentum transfer in the $`t`$ and $`u`$ channels.

The linearly divergent diagrams 1-3 sum up to zero:

$$\delta K_{13}=0.$$
(2.20)
The diagrams 4-7 cancel the divergence in the diagram $`b`$, the second term in Eq. (2.19):

$$\delta K_{47}=\frac{2g^2\mathrm{\Lambda }}{3m^2}.$$
(2.21)
All linear divergences eventually have canceled. It is satisfying to see that the net one-loop contribution to the kernel is finite, as it should be.

The remaining diagrams give:

$$\delta K_{813}=(96+1+124+2)\frac{g^2}{16\pi m}=7\frac{g^2}{8\pi m}.$$
(2.22)
Collecting together the contributions (2.19)–(2.22) we get the one-loop correction to the BS kernel:

$$\delta K=(7\mathrm{ln}3)\frac{g^2}{8\pi m}.$$
(2.23)
### 2.4 Momentum dependence

Finally, we consider the most subtle corrections that come from the momentum dependence of the kernel and of the wave function. Those were neglected in the leading order calculation because the logarithmic enhancement of the bubble diagram near threshold comes from the region of integration in (2.5) where all momenta are almost on shell. This allowed us to forget about the momentum dependence of the wave function and to set it to one. In fact, we can always require that the wave function is equal to the one on shell. This can be regarded as a normalization condition. But to compute non-logarithmic terms we need to know the complete off-shell wave function.

Once the on-shell condition on the external momenta $`p_1,p_2`$ is relaxed, we can no longer neglect the momentum transfer in the diagrams $`c`$ and $`d`$ in Fig. 2. The BS equation then gives:

$`\mathrm{\Gamma }(p_1,p_2)`$ $`=`$ $`{\displaystyle \frac{d^3q}{(2\pi )^3}\frac{1}{(q^2+m^2)[(pq)^2+m^2]}\left[\frac{3gm^2}{(qp_1)^2+m^2}g\right]}`$ (2.24)
$`=`$ $`{\displaystyle \frac{g}{8\pi m}}\mathrm{ln}\left({\displaystyle \frac{1}{\mathrm{\Delta }}}\right)\left({\displaystyle \frac{3m^2}{4m^2+p_1^2+p_2^2}}{\displaystyle \frac{1}{2}}\right)`$
$`=`$ $`{\displaystyle \frac{3m^2}{4m^2+p_1^2+p_2^2}}{\displaystyle \frac{1}{2}}+O\left({\displaystyle \frac{g}{m}}\right).`$
### 2.5 Next-to-leading order BS equation

With all corrections taken into account, the BS equation reads:

$`1`$ $`=`$ $`\left[1\left(3{\displaystyle \frac{9}{4}}\mathrm{ln}3\right){\displaystyle \frac{g}{24\pi \mu }}\right]^2`$ (2.25)
$`\times {\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}[{\displaystyle \frac{3\mu ^2}{4\mu ^2+p^2+(pq)^2}}{\displaystyle \frac{1}{2}}]{\displaystyle \frac{1}{(q^2+\mu ^2)[(pq)^2+\mu ^2]}}`$
$`\times \left[{\displaystyle \frac{3g\mu ^2}{(qp_1)^2+m^2}}g+(7\mathrm{ln}3){\displaystyle \frac{g^2}{8\pi \mu }}\right],`$
where external momenta $`p_1`$ and $`p_2`$ are on shell: $`p_{1,2}^2=\mu ^2`$. We need to retain terms of order $`g/\mu `$ and $`(g^2/\mu ^2)\mathrm{ln}\mathrm{\Delta }`$ on the right hand side. For this reason, $`O(g/m)`$ corrections to the wave function can be neglected, since they vanish on shell and therefore do not lead to the logarithmic enhancement. The loop integrals encountered in (2.25) are listed in Appendix and in (2.10). Using these results, we get after a little of algebra:

$$\mathrm{log}\mathrm{\Delta }=\frac{8\pi \mu }{g}+2\mathrm{log}2\frac{1}{2}\mathrm{log}3+O\left(\frac{g}{\mu }\right)$$
(2.26)
that isNote that the pre-exponential factor quoted in the published version of Ref.  is not correct.

$$\mathrm{\Delta }=\frac{4}{\sqrt{3}}\mathrm{exp}\left(\frac{8\pi \mu }{g}\right).$$
(2.27)
### 2.6 Comparison with numerical results

To compare Eq. (2.26) with numerical results one has to know a value for the dimensionless coupling $`\widehat{u}g/\mu `$. The safest thing to do is to relate $`\widehat{u}`$ to the renormalized coupling $`u_R`$ defined e.g. as in Ref., whose critical value is precisely known from Monte Carlo simulations

$$u_R=14.3(1)$$
(2.28)
in terms of $`u_R`$ we have, using results from

$$\widehat{u}=u_R\left[1+u_R\left(\frac{15}{128\pi }+\frac{3\mathrm{log}3}{32\pi }\right)\right]$$
(2.29)
so that

$$\mathrm{log}\mathrm{\Delta }=\frac{8\pi }{u_R}+\left(\frac{15}{16}+\frac{\mathrm{log}3}{4}+2\mathrm{log}2\right)+O(u_R)$$
(2.30)
Therefore the leading-order result is

$$\mathrm{log}\mathrm{\Delta }=\frac{8\pi }{u_R}=1.76(2)$$
(2.31)
corresponding to

$$\mathrm{\Delta }=0.172(3).$$
(2.32)
This result is compatible with the assumption of weak binding energy that underlies the Bethe-Salpeter formalism, and is in a very good agreement with the numerical result,

$$\mathrm{\Delta }_{MC}=0.17(3).$$
(2.33)
However, the parameter of expansion in perturbation theory is comparatively large, and the NLO correction is not small: the $`O(1)`$ contribution to $`\mathrm{log}\mathrm{\Delta }`$ is $`2.6`$, that is actually greater than the LO contribution. Therefore a reliable estimate of the subleading corrections to the binding energy seems to require the computation of a large number of terms in perturbative series, so as to make a resummation possible, which is typical for application of $`3d`$ $`\varphi ^4`$ theory to critical phenomena.

## 3 Dependence on the magnetic field

It is possible to study the dependence of the bound state masses on the magnetic field. Because the magnetic field couples to the order parameter, we know how to introduce it in the effective description in terms of $`\varphi ^4`$ theory. We can take the magnetic field into account by adding a linear term to the action:

$$S=d^3x\left[\frac{1}{2}_\mu \varphi _\mu \varphi +\frac{g}{4!}(\varphi ^2v^2)^2h\varphi \right].$$
(3.34)
The magnetic field shifts the minimum of the potential, which now is determined by the equation

$$\frac{g}{6}\varphi _0(\varphi _0^2v^2)=h.$$
(3.35)
Expanding the order parameter around $`\varphi _0`$,

$$\varphi =\varphi _0+\phi ,$$
(3.36)
we get for the action:

$$S=d^3x\left[\frac{1}{2}_\mu \phi _\mu \phi +\frac{1}{2}m_h^2\phi ^2+\frac{1}{3!}g_3\phi ^3+\frac{1}{4!}g\phi ^4\right]+S_0,$$
(3.37)
where the mass and the trilinear coupling are

$$m_h^2=\frac{g}{6}(3\varphi _0^2v^2),$$
(3.38)
$$g_3=g\varphi _0.$$
(3.39)
The sum of the diagrams in Fig. 2 gives the tree-level BS kernel:

$$K_h=\frac{g_3^2}{m_h^2}+\frac{g_3^2}{2}\frac{1}{m_h^24m_h^2}\frac{g}{2}=g\frac{7\varphi _0^2+v^2}{6\varphi _0^22v^2}.$$
(3.40)
Repeating the same steps as we used to calculate the binding energy at zero magnetic field, we get

$$\mathrm{\Delta }=\mathrm{const}\mathrm{exp}\left(\frac{16\pi m_h}{K_h}\right),$$
(3.41)
where $`\mathrm{\Delta }`$ is the dimensionless binding energy: $`M_h=(2\mathrm{\Delta })m_h`$.

It is useful to introduce the dimensionless variable

$$\xi =\frac{\varphi _0}{v},$$
(3.42)
which satisfies the equation

$$\xi ^3\xi =\frac{6}{gv^3}h=\frac{2}{m^3}\sqrt{\frac{g}{3}}h.$$
(3.43)
The binding energy is expressed in terms of $`\xi `$ and $`m`$, the mass at zero magnetic field, as

$$\mathrm{\Delta }=\mathrm{const}\mathrm{exp}\left(\frac{8\pi m}{g}\frac{(6\xi ^22)^{3/2}}{7\xi ^2+1}\right).$$
(3.44)
It is easy to see that $`\mathrm{\Delta }`$ is a decreasing function of $`h`$. Thus the magnetic field loosens the binding of composite states and shifts their masses closer to the threshold.

## 4 Duality and the spectrum of Ising systems

The $`3D`$ Ising model and $`_2`$ gauge theory are related by an exact duality transformation. The broken symmetry phase of the spin model is mapped into the confining phase of the gauge theory. The purpose of this section is to show that, in this phase, the duality relationship implies an exact coincidence of the spectra of the two theories. The proof relies on the existence of a non-zero interface tension and therefore applies exclusively to the broken symmetry phase of the spin model.

There are several books and reviews which discuss duality in spin systems and in particular in Ising-type models (see e.g. Ref. ). However the duality transformation is usually treated in the thermodynamic limit only, while to study the spectrum of the transfer matrix it is necessary to extend the analysis to finite lattices. This is possible as long as one carefully takes into account all possible boundary conditions.

The novel feature of our approach, which greatly simplifies the whole analysis, is the use of the Transfer Matrix (TM) formalism. Let us first recall the definitions of the lattice models we are concerned with, and the notion of duality in the thermodynamic limit.

* The spin Ising model

The Ising model is defined by the action

$$S_{spin}=\beta _{spin}\underset{n,m}{}s_ns_m,$$
(4.45)
where the field variable $`s_n`$ takes the values $`1`$ and $`+1`$; $`n(n_0,n_1,n_2)`$ labels the sites of a simple cubic lattice of size $`L_0`$, $`L_1`$ and $`L_2`$ in the three directions. The notation $`n,m`$ in Eq. (4.45) indicates that the sum is taken over pairs of nearest neighbor sites only.
* The gauge Ising model

The building blocks of the $`_2`$ gauge model are the link variables $`g_{n;\mu }\{1,1\}`$, which play the role of gauge fields. Denoting by $`\mu `$ the direction of the link, the action is

$$S_{gauge}=\beta \underset{n,\mu <\nu }{}g_{n;\mu \nu }$$
(4.46)
where $`g_{n;\mu \nu }`$ are the plaquette variables, defined by

$$g_{n;\mu \nu }=g_{n;\mu }g_{n+\mu ;\nu }g_{n+\nu ;\mu }g_{n;\nu }.$$
(4.47)
This action is invariant under local $`_2`$ gauge transformations defined as follows: one chooses arbitrarily a subset of the sites of the lattice and changes signs of all variables defined on the links which end on these sites (if two neighboring sites belong to the chosen subset the link that joins them is changed twice, that is not changed). It is immediately clear that the plaquette values, like any other product of links along a closed path, are invariant under this gauge transformation. For more details about this model, see e.g. Ref. .
* Duality

There is an exact duality transformation which relates the Ising model and the $`_2`$ gauge model. This transformation is known as Kramers–Wannier duality. It relates the partition functions of the two models evaluated at two different values of the coupling constants:

$`Z_{gauge}(\beta )Z_{spin}(\stackrel{~}{\beta })`$
$`\stackrel{~}{\beta }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{log}\left[\mathrm{tanh}(\beta )\right],`$ (4.48)
where $`\stackrel{~}{\beta }`$ will be denoted as the “dual coupling” in what follows.

It is easy to see that low values of $`\beta `$ are mapped into high values of $`\stackrel{~}{\beta }`$ and vice versa. Thus the confining region of the gauge theory is mapped into the broken symmetry phase of the spin model. In particular the end points of these two phases, the deconfinement transition and the magnetization transition, are mapped into each other.

An important feature of the dual transformation on a lattice of finite size is that it does not conserve the boundary conditions (BC). In the thermodynamic limit this fact becomes irrelevant, but on lattices of finite extent it cannot be neglected. In particular the $`_2`$ gauge model with periodic BC in all directions is mapped by duality into the Ising spin model with fluctuating BC, so that the partition function is given by the sum of the partition functions with all possible choices of periodic ($`p`$) or antiperiodic ($`a`$) BC:

$$Z_{gauge,ppp}(\beta )=\frac{1}{2}c^V(Z_{spin,ppp}(\stackrel{~}{\beta })+Z_{spin,ppa}(\stackrel{~}{\beta })+\mathrm{}+Z_{spin,aaa}(\stackrel{~}{\beta })),$$
(4.49)
where $`VL_0L_1L_2`$ is the volume of the lattice, and $`c`$ is a constant which can be easily evaluated, but is irrelevant for our purposes.

This result is discussed in full generality (for a generic lattice geometry and symmetry group) in . In the particular case of the Ising model on a cubic lattice it can be easily obtained by a direct implementation of the duality transformation.

Let us now consider the duality transformation in the framework of the transfer matrix approach. The $`0`$ direction will be our “time”. If we choose periodic BC in the $`0`$ direction without specifying the BC in the $`1`$ and $`2`$ directions we obtain:

$$Z_{spin,pxy}=\text{tr}(𝐓_{\mathrm{𝐱𝐲}})^{L_0},$$
(4.50)
where $`𝐓_{xy}`$ denotes the transfer matrix of the model in which $`x,y\{a,p\}`$ BC are chosen in the 1 and 2 directions respectively.

The antiperiodic BC in the $`0`$ direction can be obtained by acting with a spin–flip operator $`𝐏`$, which changes the sign of all spins in a given time slice. Thus we may write

$$Z_{spin,axy}=\text{tr}𝐏(𝐓_{xy})^{L_0}.$$
(4.51)
Since the operators $`𝐓_{xy}`$ and $`𝐏`$ commute, they have a common set of eigenfunctions. Let us denote the eigenvalues of $`𝐓_{xy}`$ by $`\lambda _{xy,i}`$ and those of $`𝐏`$ by $`p_i`$. The possible values of $`p_i`$ are $`1`$ and $`1`$. States that are symmetric in the magnetization have $`p_i=1`$ and those that are antisymmetric have $`p_i=1`$.

Thus we can write

$$Z_{spin,pxy}+Z_{spin,axy}=\text{tr}(1+𝐏)(𝐓_{xy})^{L_0}=2\underset{i,sym}{}\lambda _{xy,i}^{L_0},$$
(4.52)
where $`_{i,sym}`$ means that the sum is restricted to the states that are symmetric in the magnetization (See Ref.  for further details on this type of construction).

Using this result we can write the duality relation in terms of the transfer matrix eigenvalues:

$$\underset{j}{}\gamma _j^{L_0}=c^{L_0L_1L_2}\left[\underset{i,sym}{}\lambda _{pp,i}^{L_0}+\underset{i,sym}{}\lambda _{pa,i}^{L_0}+\underset{i,sym}{}\lambda _{ap,i}^{L_0}+\underset{i,sym}{}\lambda _{aa,i}^{L_0}\right],$$
(4.53)
where $`\gamma _j`$ are the eigenvalues of the transfer matrix of the gauge system.

It is instructive to show explicitly that the sums on the two sides of Eq. (4.53) have the same number of terms. On the spin side (right hand side of Eq. (4.53)) we have $`2^{L_1L_21}`$ terms for each sector. In fact, the transfer matrix is a $`2^{L_1L_2}\times 2^{L_1L_2}`$ matrix, but only half of the eigenvalues fulfil the symmetry requirement. Thus taking together all the four sectors we end up with a sum of $`2^{L_1L_2+1}`$ terms. On the gauge side we have $`2L_1L_2`$ link variables, but we may fix the gauge. The maximum number of gauge variables that we may fix, in order to preserve the periodic BC, is exactly $`L_1L_21`$. This is easy to see: we may fix, say, all the links in the $`1`$ direction except those in the last row, where we may fix all the links in the $`2`$ direction, except the last one. Thus we end up with only $`L_1L_2+1`$ degrees of freedom left. Hence the corresponding transfer matrix has $`2^{L_1L_2+1}`$ eigenvalues, that is exactly the same number as for the spin model.

Since for finite values of $`L_1`$ and $`L_2`$ the two sums contain a finite number of terms and since Eq. (4.53) holds for any integer $`L_0`$, each $`\gamma _j`$ has to have an exact (up to an overall factor $`c^{L_1L_2}`$) counterpart in the left hand side of the equation. The overall factor cancels when one takes the ratios of all the eigenvalues to the lowest one to obtain the physical spectrum, and so the latter coincides in the two models.

However our goal is to compare the spectra of the gauge theory and the spin model both with periodic BC in all directions, while Eq. (4.53) involves a sum over many different choices of the spin model BC. Therefore we must find a relationship between the $`\lambda _{pp,i}`$ eigenvalues and those belonging to the other three sectors. This is easily done by noticing that in the low temperature phase of the Ising model with antiperiodic boundary conditions at least one interface has to be created in the system. Therefore at leading order (neglecting the cases in which more than one interface appear):

$$Z_{ppa}\mathrm{exp}(\sigma L_0L_1)Z_{ppp},$$
(4.54)
where $`\sigma `$ is the interface tension.

For the transfer matrix eigenvalues this gives

$$\lambda _{pp,max}\mathrm{exp}(\sigma L_1)\lambda _{pa,max}$$
(4.55)
where, with the notation $`\lambda _{pp,max}`$ we denote the low lying states of the spectrum Notice that this last argument does not hold in the $`2D`$ case, where the anti-periodic eigenvalues are only suppressed by a constant, independent of the system size..

So we can conclude that for sufficiently large $`L_1`$ and $`L_2`$, the low lying spectrum of the $`Z_2`$ gauge theory in the confining phase with periodic boundary conditions coincides with the symmetric sector of the Ising spin model spectrum in the low temperature phase with periodic boundary conditions only. (This last argument was already presented in ).

## 5 The dual bound state picture of glueballs in $`_2`$ gauge theory

In this section we show that the interpretation of the spectrum of three-dimensional models in the universality class of the $`3D`$ Ising model as bound states of the elementary quanta provides an analytical tool to explain the qualitative features of the glueball spectrum of $`_2`$ gauge theory, and especially its angular momentum dependence.

Consider first the spin Ising model in the low-temperature region. Here the existence of bound states can be immediately inferred from the diagrammatics of the low-temperature expansion. This was shown in  for the $`4D`$ case, where however all bound states are expected to disappear in the continuum limit due to triviality.

The mechanism is very simple, and is best explained in the transfer matrix formalism. As discussed in , to the leading order in the low-temperature expansion of the transfer matrix, all time slices are forced to have the same configuration of spins, and the the eigenvectors are given by all the possible spin configurations in a single time slice. The eigenvalues of the Hamiltonian (that is minus the logarithm of the transfer matrix) are given by the energy of each configuration.

The ground state then corresponds to the configuration in which all the spins on any time slice are parallel; the first excited state corresponds to one flipped spin, which requires four bonds to be frustrated. To proceed, one has to flip two spins: if these are chosen in non-neighboring sites, the energy is just twice the one of the first excitation, since eight bonds will be frustrated. However one can flip two neighboring spins at the cost of frustrating six bonds only: the corresponding state is a bound state of the fundamental excitation, with mass just below the two-particle threshold.

One can also construct states of given angular momentum and parity by choosing linear combinations of time slice configurations with nontrivial transformation properties under spatial rotations and parity reflections. This can be done systematically by using standard group theory results, namely the theory of representations of the dihedral group $`D^4`$, which is the relevant group on a square lattice. This analysis is performed in Ref. , to which we refer for details.

The important point is that any given angular momentum requires a certain minimum number $`n_c`$ of spin flips. This means that a bound state of such angular momentum will be composed of at least $`n_c`$ elementary quanta. In Tab. 1 we have reported the values of $`n_c`$ for the various values of the angular momentum. For example, Fig. 5 shows the simplest lattice operator corresponding to angular momentum $`J=2`$ and parity $`P=1`$, which requires $`n_c=3`$ elementary quanta, while Fig. 6 shows how states with quantum numbers $`J=0,2`$ and $`P=+,`$ can be constructed with four elementary quanta. If we assume that the binding energy is always much smaller that the common mass of the elementary constituents of the bound states, we end up with a prediction for the angular momentum dependence of the glueball spectrum of the $`_2`$ gauge theory in the strong coupling regime.

It is not at all guaranteed, of course, that the spectrum will retain these same features when going from the strong coupling regime to the scaling region: however the data from Monte Carlo simulations of the gauge theory show that this is actually the case. The third column of Tab. 1 contains the Monte Carlo results for the glueball masses (from Ref. ), normalized to the mass of the lowest glueball. This shows that our dual bound state picture indeed explains the peculiar qualitative features of the glueball spectrum.

In particular, some characteristic approximate degeneracies in the spectrum find a natural explanation in this picture: for example the states $`0^{+,}`$ and $`2^\pm `$ are nearly degenerate because they are both bound states of $`n_c=3`$ elementary constituents. Bound states of $`n_c=4`$ constituents give rise to the $`0^{}`$ and $`2^{\pm ,}`$ states, again explaining their near degeneracy, and the same applies to $`0^,`$ and $`1^\pm `$ when interpreted as $`n_c=5`$ bound states§§§The $`2^+`$ state could be realized with $`n_c=2`$, but a general theorem (see Ref. ) forces all states with $`J0`$ to be degenerate in parity.. We do not expect these degeneracies to be exact, since there is no reason to expect the binding energies to be exactly the same for bound states of different angular momentum.

Once the minimum number $`n_c`$ of constituents for a given value of the angular momentum has been reached, it is easy to see that states with the same angular momentum can be constructed out of any number $`n>n_c`$ of constituents. This suggests that the approximate degeneracies just discussed appears between pairs of states only because, in general, only two states can be detected with numerical methods in each channel. It is very likely that the degeneracies actually group together larger and larger sets of states as the number of constituents is increased. For example, we might conjecture that the states $`0^,`$ and $`1^\pm `$ ($`n_c=5`$) are nearly degenerate with the (as yet undetected) states $`0^{+,}`$, $`2^{\pm ,}`$. This degeneracy pattern is a prediction of the dual bound state picture of glueballs.

## 6 Summary and further developments

The broken symmetry phase of $`3D`$ models with $`_2`$ symmetry shows a rich spectrum of massive excitations, contrary to naive expectations based on the perturbative field-theoretical treatment. The non-perturbative states in the spectrum can be interpreted as bound composites of the elementary particle excitation. When applied to $`3D`$ $`_2`$ gauge theory, this approach gives a nice explanation of the angular momentum dependence of the glueball spectrum. In particular, an observed rather peculiar degeneracy pattern naturally arises in the bound-state interpretation.

It would be interesting to consider the contribution of the bound states to universal amplitude ratios. The obvious way to do that is to use the ladder approximation for correlation functions, which incorporates all diagrams that blow up near the two-particle threshold. It would be also interesting to understand if there are bound states in the systems with $`N>1`$ order parameters. The physics is completely different in this case because the dominating forces are long-range due to the Goldstone bosons of the broken $`O(N)`$ symmetry.

### Acknowledgements

K.Z. would like to thank P. Simon for discussions. The work of K.Z. was supported by NSERC of Canada, by Pacific Institute for the Mathematical Sciences and in part by RFBR grant 98-01-00327 and RFBR grant 00-15-96557 for the promotion of scientific schools.

## Appendix A

$$\frac{d^3q}{(2\pi )^3}\frac{1}{[4\mu ^2+q^2+(pq)^2](q^2+m^2)[(pq)^2+m^2]}=\frac{1}{32\pi \mu ^3}\mathrm{ln}\left(\frac{4}{9\mathrm{\Delta }}\right).$$
(A.1)
$$\frac{d^3q}{(2\pi )^3}\frac{1}{(q^2+m^2)[(pq)^2+m^2][(qp_1)^2+m^2]}=\frac{1}{16\pi \mu ^3}\mathrm{ln}\left(\frac{4}{9\mathrm{\Delta }}\right).$$
(A.2)
$`{\displaystyle \frac{d^3q}{(2\pi )^3}}`$ $`{\displaystyle \frac{1}{[4\mu ^2+q^2+(pq)^2](q^2+m^2)[(pq)^2+m^2][(qp_1)^2+m^2]}}`$ (A.3)
$`={\displaystyle \frac{1}{32\pi \mu ^5}}\left[\mathrm{ln}\left({\displaystyle \frac{4}{9\mathrm{\Delta }}}\right){\displaystyle \frac{2}{3}}\right].`$
All integrals are calculated on shell: $`p_{1,2}^2=\mu ^2`$, and near the two-particle threshold: $`p^2(p_1+p_2)^2=M^2=(2\mathrm{\Delta })^2\mu ^2`$.


format: text

text:
# Phase-Field Formulation for Quantitative Modeling of Alloy Solidification

\[

## Abstract

A phase-field formulation is introduced to simulate quantitatively microstructural pattern formation in alloys. The thin-interface limit of this formulation yields a much less stringent restriction on the choice of interface thickness than previous formulations and permits to eliminate non-equilibrium effects at the interface. Dendrite growth simulations with vanishing solid diffusivity show that both the interface evolution and the solute profile in the solid are well resolved.

\]

The phase-field approach has emerged as a method of choice to simulate microstructural evolution during solidification . This approach has the well-known advantage that it avoids to explicitly track a sharp boundary by smearing the interface region over some thickness $`W`$. However, it has the disadvantage that it is hard to use quantitatively. This is because it is often computationally too stringent to choose $`W`$ small enough to resolve the desired sharp-interface limit of the phase-field model, even on computers of today. This is especially true for small growth rates where the scale of the microstructure is typically several orders of magnitude larger than the microscopic capillary length $`d_0`$.

Progress has recently been made to overcome this difficulty by using a “thin-interface” analysis of the phase-field model  where $`W`$ is assumed small compared to the scale of the pattern, but not smaller than $`d_0`$. For the standard symmetric model (with equal thermal conductivities in the solid and liquid), Karma and Rappel (KR) have shown that this thin-interface limit yields two essential benefits . Firstly, it maps onto the standard set of sharp-interface equations that one obtains in the classical sharp-interface limit where $`W/d_00`$, but yields a much less stringent restriction on $`W/d_0`$ that renders the computations tractable. Secondly, it makes it possible to eliminate interface kinetic effects by a specific choice of phase-field model parameters.

These two properties combined have made this thin-interface limit ideally suited to model dendritic growth in pure materials quantitatively at low undercooling when used in conjunction with efficient numerical algorithms . However, how to extend these results in a useful way to the more general case where the two phases do not have symmetrical properties has remained an open challenge. Using a distinct thin-interface analysis, Almgren showed that the two-sided model with unequal thermal conductivities maps onto a modified set of sharp-interface equations that is plagued by finite interface thickness effects . These include (i) a temperature jump across the interface, (ii) an interface stretching correction to the heat conservation condition at the interface (Stefan condition) also found previously in , and (iii) a surface diffusion correction to the same condition.

The same effects plague the thin-interface limit of existing phase-field models of alloy solidification  and make them inadequate to model quantitatively the low growth rate regime of experimental relevance. In alloys, (i) translates into a chemical potential jump across the interface associated with the well-known effect of solute trapping , and (ii) and (iii) modify the mass conservation condition at the interface. Here, I present a phase-field formulation of alloy solidification that makes it possible to eliminate *simultaneously* all three effects. Furthermore, I demonstrate that it yields the same computational benefits as the thin-interface limit of the symmetric model  for dendritic growth.

For clarity of exposition, I first discuss the thin-interface limit of a simpler model that describes an idealized binary alloy with parallel liquidus and solidus slopes, and which reduces to the standard Hele-Shaw flow problem  in its Laplace limit. I then consider a realistic dilute alloy model with unequal slopes that reduces to the former model in the limit where the partition coefficient $`k1`$. The equations of the first model are

$`\tau _t\varphi =W^2^2\varphi _\varphi \left[f(\varphi )+\lambda g(\varphi )u\right],`$ (1)
$`_tc+\stackrel{}{}\stackrel{}{j}=0,`$ (2)
where $`c`$ is the concentration defined as the mole fraction of B in a binary alloy of A and B,

$$\stackrel{}{j}=\mathrm{\Delta }c_0Dq(\varphi )\stackrel{}{}ua_tW\mathrm{\Delta }c_0_t\varphi \stackrel{}{}\varphi /|\stackrel{}{}\varphi |,$$
(3)
and

$$uc/\mathrm{\Delta }c_0+h(\varphi )/2(c_{s0}+c_{l0})/(2\mathrm{\Delta }c_0)$$
(4)
is a dimensionless measure of the departure of the chemical potential from its equilibrium value with $`u=0`$ in equilibrium. $`c_{s0}`$ ($`c_{l0}`$) are the equilibrium concentrations in the solid (liquid) at a fixed temperature $`T_0`$, $`\mathrm{\Delta }c_0c_{l0}c_{s0}`$, $`f(\varphi )=\varphi ^2/2+\varphi ^4/4`$ is a double well function with the minima $`\varphi =\pm 1`$ corresponding to the solid ($`+1`$) and liquid ($`1`$), $`g(\varphi )`$ is an odd function of $`\varphi `$ with $`g(\pm 1)=\pm 1`$ and vanishing first and second derivatives $`g^{}(\pm 1)=g^{\prime \prime }(\pm 1)=0`$, and $`h(\varphi )`$ is an odd function of $`\varphi `$ with $`h(\pm 1)=\pm 1`$ that can be chosen independently of $`g(\varphi )`$ for non-variational dynamics .

The above model reduces to the symmetric model for the choice $`q(\varphi )=1`$ and $`a_t=0`$. In this case, the thin-interface limit of this model maps onto the Stefan problem defined by: $`_tu=D^2u`$ in both phases, the Stefan condition $`V=D(_nu|^+_nu|^{})`$, where $`V`$ is the interface velocity and $`_nu|^\pm `$ is the normal gradient of $`u`$ on the solid $`()`$ and liquid side $`(+)`$ of the interface, and the velocity-dependent Gibbs-Thomson condition

$$u=d_0\kappa \beta V,$$
(5)
where

$$d_0=a_1W/\lambda ,\mathrm{and}\beta =a_1\left[\tau /(W\lambda )a_2W/D\right].$$
(6)
The expressions for the coefficients $`a_1`$ and $`a_2`$ are identical to those derived by KR  and interface kinetics can be eliminated ($`\beta =0`$) by choosing $`\lambda =D\tau /(a_2W^2)`$.

Next, for the alloy case, $`q(\varphi )`$ must now be chosen to vary from $`q(1)=1`$ in the liquid to $`q(+1)=D_{\mathrm{solid}}/D`$ in the solid. I consider explicitly the one-sided limit where $`D_{\mathrm{solid}}/D0`$, but the results also extend to the more realistic case where $`D_{\mathrm{solid}}/D1`$. The essential new term that yields the desired thin interface limit is the *anti-trapping* mass current that corresponds to the second term on the right-hand-side of equation (3), and is only non-vanishing in the diffuse interface region. It produces a solute flux from the solid to the liquid along the direction normal to the interface that counter-balances the trapping current associated with the jump of chemical potential across the interface: $`\mathrm{\Delta }uu^+u^{}`$, where $`u^\pm `$ correspond to the liquid ($`+`$) and solid ($``$) sides of the interface. Thus the anti-trapping current makes it possible to eliminate this jump while still leaving enough freedom to choose the other functions in the model to eliminate the corrections to the mass conservation condition. Repeating the analyses of KR  and Almgren , I obtain that $`\mathrm{\Delta }u`$ vanishes if $`F^+=F^{}`$, where

$$F^\pm =_0^\pm \mathrm{}𝑑\eta \left[p(\varphi _0(\eta ))p(1)\right],$$
(7)
and

$$p(\varphi _0)=\left(h(\varphi _0)1+a_t\sqrt{2}(1\varphi _0^2)\right)/q(\varphi _0).$$
(8)
In the above definitions, $`\varphi _0(\eta )=\mathrm{tanh}(\eta /\sqrt{2})`$ is the equilibrium phase-field profile, where $`\eta `$ is a coordinate that runs normal to the interface scaled by $`W`$, and I have used the identity $`_\eta \varphi _0=(1\varphi _0^2)/\sqrt{2}`$. Next, the mass conservation condition has the form

$$V=D_nu|^+c_1W\kappa Vc_2WD_s^2u$$
(9)
where $`c_1H^+H^{}`$, $`c_2Q^+Q^{}`$, and

$`H^\pm `$ $`=`$ $`{\displaystyle _0^\pm \mathrm{}}𝑑\eta \left[h(\varphi _0(\eta ))h(1)\right],`$ (10)
$`Q^\pm `$ $`=`$ $`{\displaystyle _0^\pm \mathrm{}}𝑑\eta \left[q(\varphi _0(\eta ))q(1)\right].`$ (11)
The second and third term on the right-hand-side of equation (9) represent the solute redistribution due to stretching the interface and by diffusion along its arclength $`s`$, respectively. These two terms appear, equivalently, at two successive orders in the thin interface limit considered by KR, or both at second order in the distinct thin-interface limit considered by Almgren .

In summary, one is left with three conditions to satisfy: (i) $`F^+=F^{}`$ (chemical potential jump), (ii) $`H^+=H^{}`$ (stretching), and (iii) $`Q^+=Q^{}`$ (surface diffusion). It is actually possible to satisfy two of these three conditions simultaneously within the standard phase-field formulation without the anti-trapping current. For example, with $`a_t=0`$ in equation (3), the choice $`h(\varphi )=1(1\varphi )^2/2`$ and $`q(\varphi )=(1\varphi )/2`$ satisfy (i) and (iii), but not (ii), and it is also possible to satisfy (i) and (ii) but not (iii). However, test simulations show that with either $`c_1`$ or $`c_2`$ non-vanishing, the morphological stability of the interface that sets the initial scale of the pattern becomes significantly altered for computationally tractable choices of $`W`$ , which is also easy to check analytically by repeating the standard Mullins-Sekerka analysis with the modified boundary condition (9). Similarly, it can be shown that solvability theory predictions for the dendrite tip become modified. Consequently, all three conditions must be satisfied to lift the restriction on $`W`$.

With the anti-trapping current present, one is now free to make the simplest choices $`h(\varphi )=\varphi `$ and $`q(\varphi )=(1\varphi )/2`$ that satisfy $`H^+=H^{}`$ and $`Q^+=Q^{}`$, respectively. By choosing $`a_t=1/(2\sqrt{2})`$, one can then reduce the function $`p(\varphi _0)`$ to the simple form $`p(\varphi _0)=\varphi _01`$ that satisfies $`F^+=F^{}`$ (and a non-vanishing amount of trapping can also be obtained by varying $`a`$). Remarkably, with these choices, the thin interface limit of the one-sided model produces a velocity-dependent Gibbs-Thomson condition that is *identical* to the one of the symmetric model. Consequently the expressions for $`a_1`$ and $`a_2`$ that determine $`d_0`$ and $`\beta `$ are the same as in Ref. : $`a_1=I/J`$ and $`a_2=(K+JF)/(2I)`$ where $`I=_{\mathrm{}}^+\mathrm{}𝑑\eta (_\eta \varphi _0)^2`$, $`J=g(1)g(+1)`$, $`FF^\pm =\sqrt{2}\mathrm{ln}2`$, and

$$K=_{\mathrm{}}^+\mathrm{}𝑑\eta _\eta \varphi _0(\eta )g^{}(\varphi _0(\eta ))_0^\eta 𝑑\xi \varphi _0(\xi ),$$
(12)
where I have defined $`g^{}(\varphi )_\varphi g(\varphi )`$. It also follows that the standard Hele-Shaw flow problem with an infinite viscosity contrast  can be simulated by taking the limit $`_tu0`$ of the present phase-field model.

Consider now the standard one-sided dilute alloy model defined by the set of equations

$`_tc`$ $`=`$ $`D^2c,`$ (13)
$`c_l(1k)V`$ $`=`$ $`D_nc|^+,`$ (14)
$`c_l/c_l^0`$ $`=`$ $`1(1k)d_0\kappa `$ (15)
where $`d_0=\gamma T_M/[L|m|(1k)c_l^0]`$ is the chemical capillary length, $`T_M`$ is the melting temperature, $`L`$ is the latent heat of melting, $`m`$ is the liquidus slope, $`k=c_s/c_l`$ is the partition coefficient where $`c_l`$ ($`c_s`$) is the concentration on the liquid (solid) side of the interface, and solidification is again assumed to take place isothermally.

To construct a thin interface limit that maps onto the free-boundary problem (13)-(15), I follow the same procedure of adding a local anti-trapping current in order to eliminate the jump of chemical potential together with the other terms. The equations are

$$\tau _t\varphi =W^2^2\varphi _\varphi \left[f(\varphi )+\frac{\lambda }{1k}g(\varphi )(e^u1)\right],$$
(16)
with the same continuity relation (2) as before,

$`\stackrel{}{j}`$ $`=`$ $`Dcq(\varphi )\stackrel{}{}ua_tWc_l^0(1k)e^u_t\varphi \stackrel{}{}\varphi /|\stackrel{}{}\varphi |,`$ (17)
$`u`$ $`=`$ $`\mathrm{ln}\left[{\displaystyle \frac{2c/c_l^0}{1+k(1k)h(\varphi )}}\right],`$ (18)
and where $`g(\varphi )`$, $`h(\varphi )`$, and $`q(\varphi )`$ obey the same limits at $`\varphi =\pm 1`$. $`u`$ is again a dimensionless measure of the departure of the chemical potential $`\mu `$ from equilibrium (namely here $`uv_0(\mu \mu _E)/(RT_0)`$ where $`R`$ is the rare gas constant and $`v_0`$ is the molar volume assumed to be constant). The logarithmic dependence of $`u`$ on $`c`$ is related to the entropy of mixing in the free-energy density as in previous models . The main difference here is again the addition of the anti-trapping mass current. Also, the present thin-interface analysis differs from the analysis of Kim et al.  that does not consider interface stretching and surface diffusion, and assumes that $`q(\varphi )`$ is constant in the interface region.

The condition for eliminating surface diffusion becomes now $`Z^+=Z^{}`$, where $`Z^\pm =_0^\pm \mathrm{}𝑑\eta \left[M(\eta )M(\pm \mathrm{})\right]`$ where $`M(\eta )q(\varphi _0(\eta ))c_0(\eta )`$, and

$$c_0(\eta )=c_l^0\left[1+k(1k)h(\varphi _0(\eta ))\right]/2$$
(19)
is the equilibrium concentration profile across a stationary interface. Therefore choosing

$$q(\varphi )=\frac{1\varphi }{1+k(1k)h(\varphi )}$$
(20)
satifies this condition. The choices $`h(\varphi )=\varphi `$ and $`a_t=1/(2\sqrt{2})`$ then make the jump of $`u`$ and the interface stretching term vanish up to finite corrections that turn out to be small in computations and will be discussed elsewhere. Consequently, the conditions for $`c`$ on the two sides of the interface have the desired form, $`c_l/c_l^0=1(1k)d_0\kappa (1k)\beta V`$ and $`c_s=kc_l`$ where the expressions for $`d_0`$, $`\beta `$, $`a_1`$, and $`a_2`$ are again identical to those for the symmetric model quoted earlier here. Therefore, $`\beta `$ can again be made to vanish. The dilute alloy model is easily shown to reduce to the parallel slope model in the limit $`k1`$ by making the change of variable $`U=(e^u1)/(1k)`$. A small solid diffusivity can also be modeled by adding $`(1+\varphi )D_{\mathrm{solid}}/(2D)`$ to the right-hand-side of equation (20), which has a negligible effect on the thin-interface limit in typical alloys where solid diffusion is substitutional ($`D_{\mathrm{solid}}/D10^4`$) .

Finally, the model is straightforward to extend to directional solidification with the standard frozen temperature approximation $`T(z)=T_0+G(zV_pt)`$ that yields (for $`\beta =0`$) the interface condition

$$c_l/c_l^0=1(1k)d_0\kappa (1k)(zV_pt)/l_T,$$
(21)
where $`V_p`$ is the pulling velocity of the sample along the $`z`$-axis, $`l_T=|m|(1k)c_l^0/G`$ is the thermal length , $`G`$ is the temperature gradient, and $`c_l^0=c_{\mathrm{}}/k`$ where $`c_{\mathrm{}}c(z=+\mathrm{})`$. It simply suffices to substitute $`e^u`$ by $`e^u+(1k)(zV_pt)/l_T`$ in equation (16).

The convergence of the model was examined by carrying out two-dimensional simulations of isothermal dendritic growth. These simulations are directly analogous to the ones carried out previously to test the thin interface limit of the symmetric model . Crystalline anisotropy was included by generalizing equation (16) to a standard anisotropic form

$`\tau (\mathrm{\Theta })_t\varphi =_\varphi \left[f(\varphi )+{\displaystyle \frac{\lambda }{1k}}g(\varphi )(e^u1)\right]`$ (22)
$`+\stackrel{}{}(W(\mathrm{\Theta })^2\stackrel{}{}\varphi )_x\left(W(\mathrm{\Theta })W^{}(\mathrm{\Theta })_y\varphi \right)`$ (23)
$`+_y\left(W(\mathrm{\Theta })W^{}(\mathrm{\Theta })_x\varphi \right),`$ (24)
where $`\mathrm{\Theta }\mathrm{arctan}(_y\varphi /_x\varphi )`$ is the angle between the direction normal to the phase-field interface and the $`x`$-axis. As a result, the anisotropic form of $`d_0`$ and $`\beta `$ become $`d_0(\theta )=a_1[W(\theta )+W^{\prime \prime }(\theta )]/\lambda `$ and $`\beta (\theta )=a_1\left[\tau (\theta )/(W(\theta )\lambda )a_2W(\theta )/D\right]`$ in the interface condition where $`a_1=0.8839`$, and $`a_2=0.6267`$ for the common choice $`_\varphi g(\varphi )=(1\varphi ^2)^2`$. Furthermore, I chose the standard form of four-fold anisotropy $`W(\theta )=Wa_s(\theta )`$ with $`a_s(\theta )=1+ϵ_4\mathrm{cos}4\theta `$ and made $`\beta (\theta )`$ vanish by letting $`\tau (\theta )=\tau a_s(\theta )^2`$ and $`\lambda =D\tau /(a_2W^2)`$.

I simulated equations (2) and (22) with $`\stackrel{}{j}`$ and $`u`$ defined by equations (17) and (18), respectively. I compare the results of the present model with $`a_t=1/(2\sqrt{2})`$ and $`q(\varphi )`$ defined by equation (20), to the more standard choice (i.e. similar to previous models) that has no anti-trapping current ($`a_t=0`$) and uses the simplest scaled diffusivity function $`q(\varphi )=(1\varphi )/2`$; $`h(\varphi )=\varphi `$ and $`_\varphi g(\varphi )=(1\varphi ^2)^2`$ in both models. The former model has the desired thin-interface limit with local equilibrium at the interface, whereas the latter has both a chemical potential jump and surface diffusion. I used a simple finite-difference Euler method with $`\mathrm{\Delta }x=0.4`$ and $`\mathrm{\Delta }t=0.008`$, $`W=\tau =1`$, $`ϵ_4=0.02`$, $`k=0.15`$, and the scaled supersaturation $`\mathrm{\Omega }=(c_l^0c_{\mathrm{}})/[c_l^0(1k)]=0.55`$ where $`c_{\mathrm{}}`$ is the initial alloy concentration. In all simulations, the initial condition consisted of a circular seed of radius $`r=22d_0`$, $`u=\mathrm{ln}(1(1k)\mathrm{\Omega })`$, and $`c`$ defined by equation (18) that varies from $`c_{\mathrm{}}`$ in liquid to $`kc_{\mathrm{}}`$ in solid.

The dimensionless dendrite tip velocity $`Vd_0/D`$ is plotted vs the dimensionless time $`tD/d_0^2`$ for the two models and two different ratios of $`d_0/W`$ in Fig. 1. Since the simulation time scales $`(d_0/W)^5`$, the runs with $`d_0/W`$ twice smaller are $`32`$ times shorter. Furthermore, I have plotted in Fig. 2 the scaled concentration $`c_s(x)/c_l^0`$ in the solid vs the scaled position $`x/d_0`$ along the central dendrite axis for the two different models. Plots in Figs. 1-2 must superimpose when results are converged. These plots show, as expected, that the present model is well converged in this range of $`d_0/W`$ that is comparable to the one studied in the symmetric model , whereas the standard model is not. This is especially true for the microsegregation profile that is still far from being converged in the latter model, even for the largest ratio $`d_0/W=0.544`$. In contrast, this profile is already well converged for a twice smaller ratio in the present model; it agrees, self-consistently, within a few percent with the Gibbs-Thomson relation $`c_s(x)/c_l^0=k[1(1k)d_0/\rho ]`$ where $`\rho `$ is the dendrite tip radius in the simulation. It will be shown elsewhere that this dramatic difference of convergence for microsegregation is due to the fact that the amount of solute trapped $`p\mathrm{ln}p`$ for small velocity ($`p=WV/D1`$) in the standard model.

The present results demonstrate that the phase-field method can be successfully extended to model quantitatively microstructural pattern formation in alloys with a realistic solid diffusivity. For this important application, it is potentially more advantageous than the level set method  since it does not require the explicit computation of the interface velocity. These results also revive the hope to extend the phase-field method to model accurately a wide range of other interfacial patterns with a strong asymmetry between phases.

This research is supported by U.S. DOE Grant No. DE-FG02-92ER45471 and benefited from computer time allocation at NERSC and NU-ASCC. I thank Roger Folch and Mathis Plapp for valuable discussions.