text:
# Grain boundary motion in layered phases

## I Introduction

This paper addresses the motion of a grain boundary separating two regions of parallel convective rolls, as can be observed in a Rayleigh-Bénard convection cell of large aspect ratio. Each semi-infinite region (or grain) is comprised of roll patterns of same wavenumber $`q_0`$, but the corresponding wavevectors are mutually perpendicular. Our main focus in this paper is the relationship between local roll curvature and grain boundary motion. Although the model equations used in our analysis are appropriate for Rayleigh-Bénard convection near onset , we expect that the qualitative features of our findings also hold in others systems that exhibit layered phases, such as, for example, lamellar phases in weakly segregated block copolymers .

Consider a disordered system brought into a layered phase, e.g. by a temperature quench in the case of a diblock copolymer, or by a change in Rayleigh number in a Rayleigh-Bénard convection cell. Domains comprised of rolls (in Rayleigh-Bénard convection) or lamellae (in block copolymers near a symmetric mixture) quickly form that have a well defined characteristic wavenumber $`q_0`$ (in the case of Rayleigh-Bénard convection near threshold, for example, $`q_0`$ lies on the marginal stability boundary against a zig-zag instability ). Due to translational and rotational invariance, the spontaneous relaxation from the initial disordered state leads in practice to a large number of such domains, and a sufficiently large system remains isotropic macroscopically. Such a configuration also contains a large density of defects, such as grain boundaries, disclination and dislocations.

Defect motion in two dimensional layered phases has been studied extensively, especially in connection with the evolution of convective rolls in Rayleigh-Bénard cells . A primary question is how defect velocities are related to features of the background surrounding them, such as roll periodicity. It is well known, for example, that dislocations climb , or that grain boundary motion between two domains with straight rolls provides a wavelength selection mechanism . On the other hand, less attention has been paid to the relationship between defect motion and roll curvature . We study in this paper a particular situation that appears to be prevalent (even if idealized) during the formation and evolution of layered structures: curvature induced motion of a grain boundary separating two semi-infinite ordered domains. We first show that if the rolls of one of the domains are periodically modulated along their transverse direction, the perturbation decays in time with a rate proportional to $`q^4`$, where $`q`$ is the wavenumber of the modulation. In addition, a novel feature is associated with this relaxation: the average position of the boundary does not remain stationary, but undergoes a net displacement such that the size of the region with straight unperturbed rolls increases. Analytic calculations based on the amplitude equation formalism show that this motion is a nonlinear effect, that it occurs in a time scale much shorter than the linear relaxation of the curved rolls. The total distance traveled by the boundary scales as $`ϵ^{1/2}`$, where $`ϵ`$ is the reduced Rayleigh number. Consequently, the grain boundary can travel large distances, even for very small initial perturbations.

The analysis of grain boundary motion and relaxation, while interesting in its own right, is also expected to contribute to our understanding of the coarsening of an ensemble of grains. Linear analyses of boundary or defect motion have already been used in the past to predict coarsening exponents . The method reproduces known exponents in the case of models $`A`$ and $`B`$ in the lexicon of Hohenberg and Halperin , and has also been used to predict coarsening exponents for O(N) vector models  (see also ). Coarsening laws for layered phases, on the other hand, remain largely unexplained .

In Section II, we describe the configuration of the grain boundary studied, as well as the model equations. Section III presents a linear stability analysis of long wavelength transverse modulations near the boundary. The numerical method used to validate our solutions is outlined in Section IV, as well as the motivation for the nonlinear analysis given in Section V. Finally, in Section VI, we present conclusions, and briefly discuss possible implications of our results on coarsening rates.

## II Grain boundary configuration and governing equations

The base or reference state is a stationary and planar grain boundary separating two semi-infinite domains $`A`$ and $`B`$ of straight convective rolls. Both sets of rolls have the same wavelength $`\lambda _0=2\pi /q_0`$. It is in general expected from previous studies that such grain boundaries are stable against low wavenumber perturbations . We wish to study here the decay of a perturbation of wavenumber $`qq_0`$ applied in the direction transverse to the rolls of one of the domains, referred to as to domain $`A`$ (see Fig. 4a). In this article, we consider for all our analytic calculations the particular case in which the two sets of rolls are oriented at right angles relative to each other in the way depicted on Fig. 4b. In this case, the rolls of domain $`A`$ are parallel to the grain boundary itself. We expect that a $`90^{}`$ grain boundary is the boundary of lowest energy, and hence that it is the most common in an extended system that evolves spontaneously. However, our analysis can be generalized to cases such as shown on Fig. 4a, in which $`A`$ and $`B`$ rolls have arbitrary angles with respect to the grain boundary, and we do not anticipate qualitative changes in our conclusions.

Our focus is on transverse modulations since this is the simplest kind of perturbation that induces roll curvature without any average roll compression (compression is associated with longitudinal perturbations along the roll normal). Our study differs from that of ref.  in which a grain boundary was perturbed with a modulation of wave vector having components both transverse and longitudinal to the rolls. As was shown in ref., the longitudinal component of the perturbation dominated the relaxation. We analyze here the relaxation of a pure transverse (or curvature) mode, and obtain results that are qualitatively different from those of ref. .

Let an order parameter $`\psi `$ represent, for example, the vertical velocity in the mid-plane of a convecting fluid layer. Just above onset, when the reduced Rayleigh number $`ϵ=(RR_c)/R_c1`$ ($`R_c`$ is the critical Rayleigh number for instability), $`\psi `$ can be expanded as the superposition of two waves with slowly varying complex amplitudes $`A`$ and $`B`$ ,

$$\psi (x,y,t)=\frac{1}{2}[A(x,y,t)e^{iq_0x}+B(x,y,t)e^{iq_0y}+c.c];A|_{x=\mathrm{}}=0,B|_{x=+\mathrm{}}=0.$$
(1)
The amplitudes satisfy a set of two coupled Ginzburg-Landau equations that can be derived from the equation for $`\psi `$ (for instance the Swift-Hohenberg model ) by multi scale analysis ,

$$\frac{A}{t}=\frac{\delta F}{\delta \overline{A}},\frac{B}{t}=\frac{\delta F}{\delta \overline{B}},$$
(2)
where $`\delta F/\delta \overline{A}`$ is the functional derivative with respect to the complex conjugate of $`A`$, and $`F`$ is a Lyapunov functional,

$`F`$ $`=`$ $`{\displaystyle }d\stackrel{}{r}={\displaystyle }d\stackrel{}{r}\{ϵ(|A|^2+|B|^2)+{\displaystyle \frac{g}{2}}(|A|^4+|B|^4)+g_{}|A|^2|B|^2+`$ (4)
$`\xi _0^2|(_x{\displaystyle \frac{i}{2q_0}}_y^2)A|^2+\xi _0^2|(_y{\displaystyle \frac{i}{2q_0}}_x^2)B|^2\}.`$
The coherence length $`\xi _0`$ is of the order of $`q_0^1`$, $`g`$ and $`g_{}`$ are interaction coefficients, and the time scale factor has been set to unity. These three parameters depend on the particular model equation considered for the original field $`\psi `$.

The nontrivial stationary solutions $`\{A_0,B_0\}`$ of the governing set of equations (2)-(4) describe a planar boundary, and depend only on the coordinate $`x`$ normal to the boundary. They are defined by,

$`0`$ $`=`$ $`ϵA_0+\xi _0^2_x^2A_0gA_0^3g_{}B_0^2A_0`$ (5)
$`0`$ $`=`$ $`ϵB_0{\displaystyle \frac{\xi _0^2}{4q_0^2}}_x^4B_0gB_0^3g_{}A_0^2B_0.`$ (6)
This system of equations was extensively studied in refs. . It is a planar grain boundary of thickness proportional to $`\xi _0/\sqrt{ϵ}`$. Note that the system of equations is not invariant under permutation of $`A_0`$ and $`B_0`$. The amplitude $`A_0`$ of the rolls parallel to the interface vanishes as $`\mathrm{exp}(x\sqrt{ϵ}/\xi _0)`$ when $`x\mathrm{}`$, and saturates to $`(ϵ/g)^{1/2}\mathrm{tanh}(x\sqrt{ϵ}/\xi _0)`$ when $`x+\mathrm{}`$. The behavior of the amplitude of the rolls perpendicular to the interface is slightly different: $`B_0(x)(ϵ/g)^{1/2}\mathrm{exp}(x\sqrt{ϵ}/\xi _0)`$ when $`x\mathrm{}`$, and there is a point $`x^{}`$ such that $`B_0(x^{})=0`$. It is customary to define the location of the grain boundary by the point $`x^{}`$. To a good approximation, $`B_0(x)0`$ for $`xx^{}`$.

We emphasize that a stationary solution only exists when the wavenumber of the solution equals the marginal wavenumber $`q_0`$, and that at first order in $`ϵ`$, the solutions $`A_0`$ and $`B_0`$ only depend on the slow spatial variable $`ϵ^{1/2}x`$, but not on the fast scale $`x`$. Hence, the location of the boundary and the phase of the rolls are independent .

## III Linear stability analysis of a 90 grain boundary

We linearize Eqs. (2) around $`\{A_0,B_0\}`$, and assume the following perturbed solutions (note that the form of the perturbation explicitly assumes that the boundary does not undergo a net average displacement),

$`A(x,y,t)`$ $`=`$ $`(A_0(x)+\stackrel{~}{a})e^{i\stackrel{~}{\varphi }},`$ (7)
$`B(x,y,t)`$ $`=`$ $`B_0(x)+\stackrel{~}{b}.`$ (8)
The phase $`\stackrel{~}{\varphi }`$ represents a transverse distortion, as shown in Fig. 4b, and is of the general form

$$\stackrel{~}{\varphi }=\varphi (x)\mathrm{cos}(qy)e^{\sigma t},\stackrel{~}{\varphi }1.$$
(9)
The real fields $`\stackrel{~}{a}`$ and $`\stackrel{~}{b}`$ are amplitude corrections and similarly read,

$`\stackrel{~}{a}=a(x)\mathrm{cos}(qy)e^{\sigma t}`$ (10)
$`\stackrel{~}{b}=b(x)\mathrm{cos}(qy)e^{\sigma t}.`$ (11)
In the above expansion, we have neglected the imaginary part of $`B`$. $`m(B)`$ can written as $`d(x)\mathrm{sin}(qy)e^{\sigma t}`$; yet, because of the fourth order $`x`$-derivative in the equation satisfied by $`B`$ (the second relation in Eq. (2)), $`d(x)`$ is of order $`ϵ`$ lower than the real part function $`b(x)`$, and hence negligible. Inserting Eqs. (7)-(8) into Eqs. (2), we find,

$`\sigma a`$ $`=`$ $`ϵa+\xi _0^2[_x^2a(q^4/4q_0^2)a(q^2/q_0)_x(\varphi A_0)](3gA_0^2+g_{}B_0^2)a2g_{}A_0B_0b`$ (12)
$`\sigma A_0\varphi `$ $`=`$ $`\xi _0^2[A_0^1_x(A_0^2_x\varphi )(q^4/4q_0^2)A_0\varphi +(q^2/q_0)_xa]`$ (13)
$`\sigma b`$ $`=`$ $`ϵb\xi _0^2q^2b(3gB_0^2+g_{}A_0^2)b2g_{}B_0A_0a.`$ (14)
We further expand the solutions in power series of the small parameter $`q^2`$, and define, $`a=a_0+q^2a_2+q^4a_4+..`$; $`\varphi =\varphi _0+q^2\varphi _2+q^4\varphi _4+..`$; $`b=b_0+q^2b_2+q^4b_4+..`$. Since the limit $`q0`$ corresponds to a uniform translation of the interface, $`\sigma =q^2\sigma _2+q^4\sigma _4+..`$. We now analyze the resulting equations order by order in $`q^2`$.

At order $`q^0`$, Eq. (13) reduces to $`A_0^2_x\varphi _0=\mathrm{cst}`$. Since $`A_0(x)^2`$ diverges exponentially at $`\mathrm{}`$ and $`_x\varphi `$ must remain finite on the whole interval $`[\mathrm{},\mathrm{}]`$, the only admissible solution is $`_x\varphi _0=0`$. Hence $`\varphi _0=`$cst is a free parameter that represents the magnitude of the initial phase modulation. With the notation of Fig. 4b, $`\varphi _0`$ is simply related to the magnitude of the roll deformation $`\delta x_0`$ through $`\varphi _0=q_0\delta x_0`$ ($`\delta x_0\lambda _0`$). At order $`q^0`$, Eqs.(12) and (14) can be written as

$$H_0\left(\begin{array}{c}a_0\hfill \\ b_0\hfill \end{array}\right)=0.$$

If we take the $`x`$-derivative of Eqs. (5) and (6) we recover the above equation, with $`a_0`$ and $`b_0`$ replaced by $`_xA_0`$ and $`_xB_0`$. Hence the solutions at this order are:

$$a_0(x)=\alpha _0_xA_0,b_0(x)=\alpha _0_xB_0,$$
(15)
where $`\alpha _0`$ is a constant to be determined from the solvability condition at the next order.

At order $`q^2`$, the solutions $`a_2`$ and $`b_2`$ satisfy,

$$H_0\left(\begin{array}{c}a_2\hfill \\ b_2\hfill \end{array}\right)=\left(\begin{array}{c}(\alpha _0\sigma _2+\varphi _0\xi _0^2/q_0)_xA_0\hfill \\ \alpha _0(\sigma _2+\xi _0^2)_xB_0\hfill \end{array}\right).$$
(16)
The solvability condition at order $`q^2`$ requires that the right hand side vector be orthogonal to the kernel of the adjoint of $`H_0`$. Since $`H_0`$ is hermitian, and $`H_0(_xA_0,_xB_0)=0`$, we find,

$$(\alpha _0\sigma _2+\varphi _0\xi _0^2/q_0)_{\mathrm{}}^{\mathrm{}}(_xA_0)^2𝑑x+\alpha _0(\sigma _2+\xi _0^2)_{\mathrm{}}^{\mathrm{}}(_xB_0)^2𝑑x=0.$$
(17)
On the other hand, at order $`q^2`$ one can replace $`\sigma \varphi `$ by $`\sigma _2\varphi _0`$, and $`a`$ by $`\alpha _0_xA_0`$ in Eq. (13). After multiplying Eq. (13) by $`A_0`$, and integrating over $`x`$, one finds

$$A_0^2_x\varphi _2=\frac{\sigma _2\varphi _0}{\xi _0^2}_{\mathrm{}}^xA_0^2(u)𝑑u\frac{\alpha _0}{q_0}\left[A_0(x)_xA_0(x)_{\mathrm{}}^x(_uA_0)^2𝑑u\right]+C,$$
(18)
where $`C`$ is a constant of integration. However, the condition that $`|_x\varphi _2|<\mathrm{}`$ when $`x\mathrm{}`$ requires $`C=0`$. Furthermore, in the limit $`x+\mathrm{}`$, $`A_0(ϵ/g)^{1/2}`$, hence $`^xA_0^2(u)𝑑ux`$. The gradient of $`\varphi _2`$ remains finite only if $`\sigma _2=0.`$ Therefore Eq. (17) with $`\sigma _2=0`$ now yields the value of the solvability constant $`\alpha _0`$,

$$\alpha _0=\delta x_0\left(_{\mathrm{}}^{\mathrm{}}(_xA_0)^2𝑑x\right)\left(_{\mathrm{}}^{\mathrm{}}(_xB_0)^2𝑑x\right)^1,$$
(19)
where the definition $`\varphi _0=q_0\delta x_0`$ has also been used.

To summarize our results up to this point, the perturbed amplitudes at order $`q^0`$ are given by,

$`A`$ $`=`$ $`\left[A_0(x)+\alpha _0\mathrm{cos}(qy)e^{\sigma t}_xA_0(x)\right]e^{i\stackrel{~}{\varphi }}A_0\left(x+\alpha _0\mathrm{cos}(qy)e^{\sigma t}\right)e^{iq_0\delta x_0\mathrm{cos}(qy)e^{\sigma t}},`$ (20)
$`B`$ $``$ $`B_0\left(x+\alpha _0\mathrm{cos}(qy)e^{\sigma t}\right),`$ (21)
where we have used Eqs. (7)-(8) and (15). Equations (20) and (21) show that the amplitude moduli $`|A|`$ and $`|B|`$ of the weakly modulated rolls at this order equal the unmodulated profiles with a simple change in the coordinate origin. By analogy with the planar case , we define the location of the modulated boundary by the set of points $`\{x_g,y\}`$ such that $`B(x_g,y)=0`$. Eq. (21) indicates that the location of the boundary is given by the curve

$$x_g(y)=x^{}\alpha _0\mathrm{cos}(qy)e^{\sigma t},$$

with $`B_0(x^{})=0`$. It is a local translation relative to the planar boundary by $`\alpha _0\mathrm{cos}(qy)e^{\sigma t}`$ along the $`x`$-axis. Note, on the other hand, that the lines of constant phase of $`A`$ are given by $`x=cst\delta x_0\mathrm{cos}(qy)e^{\sigma t}`$ instead, so that $`\delta x_0\mathrm{cos}(qy)e^{\sigma t}`$ represents the local deformation of the straight rolls. As shown by Eq. (19), boundary and roll deformations are not independent: the lengths $`\alpha _0`$ and $`\delta x_0`$ are related to each other through the non-uniform profiles $`A_0`$ and $`B_0`$. A quite non-intuitive result is the opposite directions of the related translations, evidenced by the minus sign in Eq. (19). Consider for instance a point where the lines of constant phase of $`A`$ are displaced towards the region $`B`$; the actual position of the boundary, however, is displaced towards region $`A`$. This effect can be seen more easily on Figure 4 (see also Section IV).

At order $`q^4`$ integration of Eq. (13) leads to,

$$A_0^2_x\varphi _4=_{\mathrm{}}^x\left(\frac{\sigma _4}{\xi _0^2}+\frac{1}{4q_0^2}\right)\varphi _0A_0^2(u)𝑑u\frac{1}{q_0}_{\mathrm{}}^xA_0(u)_ua_2du,$$

where the second order result $`\sigma _2=0`$ has been used. The second integral in the right hand side is finite, but the first one diverges as $`x`$ when $`x+\mathrm{}`$, except if $`\sigma _4/\xi _0^2+1/4q_0^2=0`$. We therefore conclude that, at leading order in $`q`$

$$\sigma =\frac{\xi _0^2}{4q_0^2}q^4.$$
(22)
This is one of the central results of this section: the modulated boundary is expected to relax exponentially with a rate proportional to $`q^4`$.

The linear analysis shows that the interface is always stable with respect to long wavelength perturbations. This is in agreement with the case studied in ref. , although the authors concluded that the decay rate $`\sigma q^2`$ instead. Such a behavior was found because the wave vector of the modulation considered in their work had a nonzero component in the longitudinal direction. We note that the decay rate of Eq. (22) is in fact identical to that of a single wave ($`A`$-rolls only in Figure 4) . According to Eq. (22), the relaxation of the boundary is completely determined by the relaxation of the $`A`$ rolls far from it.

Unfortunately, in spite of the fact that the main predictions of the perturbation analysis agree well with a numerical solution of the Swift-Hohenberg equation, the linear analysis just discussed is not uniformly valid (we will further elaborate on this point in Section V). Consider Eq. (18) for $`\varphi _2`$ in the limit $`x+\mathrm{}`$. The leading behavior is given by

$$\varphi _2\frac{\alpha _0}{q_0A_0^2(\mathrm{})}\left(_{\mathrm{}}^{\mathrm{}}(_uA_0)^2𝑑u\right)x.$$
(23)
Therefore the phase perturbation $`\varphi \varphi _0+q^2\varphi _2`$, which is assumed to be small in the derivation of Eq. (13), diverges as $`x`$ for large $`x`$ through the second order correction $`\varphi _2`$ at any finite time $`t`$. As a consequence, from Eq. (9), the $`y`$-component of the wavevector of the roll pattern is,

$$_y\stackrel{~}{\varphi }=q\mathrm{sin}(qy)e^{\sigma t}(\varphi _0+q^2\varphi _2),$$

and also diverges with $`x`$. This implies that, far enough from the boundary, the $`y`$-component of the wave vector becomes larger than the $`x`$ component, given by $`q_0+_x\stackrel{~}{\varphi }q_0`$. This solution represents an increasingly zig-zagging wave, which is of course not physical. As it is apparent from Eq. (23), the singularity arises because of the presence of a non-uniform amplitude of the pattern near the boundary $`_xA_00`$. In fact, we think it is possible that the linearized system with non-constant coefficients (12)-(14) does not have bounded eigenfunctions.

It appears that the failure of the linear analysis is related to an additional important feature of boundary motion, not taken into account in the previous calculation: a net translation of the average position of the boundary towards the $`A`$-region, as shown by the numerical calculations presented in the following Section. After we present the numerical evidence, we argue in Section V that this feature cannot be obtained at the linear level in the amplitude of the perturbation, and that it is in fact a singular contribution in the limit $`ϵ0`$.

## IV Numerical solution of the Swift-Hohenberg equation

In order to test the predictions of the previous Section, and to obtain further insight into grain boundary motion, we have numerically solved the Swift-Hohenberg model of Rayleigh-Bénard convection,

$$\frac{\psi }{t}=ϵ\psi \frac{\xi _0^2}{4q_0^2}(q_0^2+^2)^2\psi \psi ^3.$$
(24)
We consider the evolution from an initial condition that corresponds to the geometry shown in Fig. 4b. We note that solving Eq. (24) in the case of curved rolls is much simpler than solving the two dimensional Ginzburg-Landau equations (2).

Equation (24) is solved numerically with a pseudo-spectral method. Further details on the algorithm and the time integration scheme can be found in ref. . The stability of the algorithm allows relatively large values of the time step, which is fixed to 0.4 in the dimensionless time units of Eq. (24). Equation (24) is then discretized on a square grid of size $`512\times 512`$, and occasionally $`256\times 256`$. Spatial discretization is such that there are $`8`$ grid nodes per wavelength $`\lambda _0`$. In all the following numerical examples, we have chosen $`\xi _0=2q_0^1`$, which is close to the value that corresponds to the stress free boundary conditions in the original system of fluid mechanical equations, from which the Swift-Hohenberg equation is derived as an approximation. We use periodic boundary conditions in both directions, and hence the initial condition comprises two symmetric, well separated grain boundaries located at $`x=1/4`$ and $`x=3/4`$, in units of the system size. We have verified that the numerical solutions of the Swift-Hohenberg equation for the stationary straight grain boundary coincides with those obtained by directly solving the amplitude equations (5)-(6), with the appropriate parameters that follow from the multi-scale analysis of the Swift-Hohenberg equation ($`g=3/4`$ and $`g_{}=3/2`$). The modulated initial condition is implemented with the help of the one dimensional solution, as

$$\psi (x,y,t=0)=A_0(x)\mathrm{cos}[q_0x+q_0\delta x_0\mathrm{cos}(qy)]+B_0(x)\mathrm{cos}[q_0y].$$
(25)
A typical configuration is shown in Figure 4, taken at the intermediate time $`t=1600`$, with $`ϵ=0.04`$, $`q=q_0/16`$ and a system size of $`256\times 256`$. Note that this solution is consistent with our earlier assumption of neglecting $`m(B)`$ in the perturbation equations (Eq. (8)), since the $`B`$ rolls remain straight near the grain boundary despite of the transverse modulation of the $`A`$ rolls. The opposition in phase of the roll profile and of the boundary location, given in Eq. (19), is also well reproduced numerically Yet, this effect is easy to observe only for relatively large values of $`\delta x_0`$, as shown on Fig. 4.

Figure 4 shows our numerical results for the inverse decay rate $`\sigma ^1`$ as a function of $`q`$. These results are obtained for a system of size $`512\times 512`$, so that the number of rolls that separate the two grain boundaries is twice larger than in Figure 4. We checked that the planar grain boundaries did not appreciably interact, and therefore remained stationary. The value of $`\sigma `$ was determined from an exponential fit to the decay of the phase of a given roll, for a small initial $`\delta x_0`$. We first compute numerically $`\sigma `$ in the absence of any grain boundary, i.e. with the single wave $`A`$. We find very good agreement with Eq. (22), given that there are no adjustable parameters in the theoretical curve. When grain boundaries are present, the numerical results (square symbols) are somewhat higher, but still compare well with Eq. (22). They remain closer to the law $`\sigma q^4`$, rather than to, say, $`\sigma q^5`$, and certainly than to $`\sigma q^2`$. The small discrepancy observed is probably due to the finite size of the system: although the two grain boundaries are separated by roughly $`30`$ rolls, they slightly interact during relaxation. As a check, we have computed $`\sigma (q=3q_0/32)`$ again in a system of size $`1024\times 1024`$, with same value of $`\lambda _0`$. The difference between that numerical result and the theoretical curve is reduced by a factor of $`2.5`$ compared with a system size of $`512\times 512`$.

The numerical solution reveals an additional feature which was not included in our linear analysis: while the boundary modulation is relaxing, there is net motion of the grain boundary towards the $`A`$ roll region. This motion is a pure amplitude mode, as the lines of constant roll phase remain stationary on average. Therefore, after a long time, the modulated rolls have completely relaxed but the size of the region of $`A`$ rolls has decreased, and the grain boundary is located at a distance $`d_{\mathrm{}}`$ away from its initial average location. The distance $`d_{\mathrm{}}`$ exceeds the initial deformation $`\delta x_0`$, and it can be even significantly larger than the roll wavelength $`\lambda _0`$. In the configuration shown in Fig. 4 (system of size $`256\times 256`$), the magnitude of the perturbation was deliberately chosen relatively large ($`\delta x_0=\lambda _0`$), so that, at late times, the modulated region completely disappears by annihilation of the two grain boundaries.

We now report our results for $`d_{\mathrm{}}`$ in a bigger system of size $`512\times 512`$, and with $`q=3q_0/32`$. Figure 4 shows our results for the total traveled distance $`d_{\mathrm{}}`$ as a function of $`\delta x_0`$. The complete relaxation of the perturbation requires times around $`10^5`$, and the location of the (nearly) flat boundary is again defined to be the point $`x^{}`$ at which $`B(x^{})=0`$ . This point if obtained from the field $`\psi `$ through the relation $`B(x)=[\psi (x,y=\lambda _0)\psi (x,y=3\lambda _0/2)]/2`$ for a flat interface (see Eq. (1)). We observe in the figure that the total traveled distance increases nonlinearly with $`\delta x_0`$. The following Section gives an interpretation of this feature, and further discussion of the results.

## V Defect motion through energy relaxation

We study in this section a more general form for the perturbed amplitudes than those given in Eqs. (20)-(21). We assume a perturbed solution of the form,

$`A`$ $`=`$ $`A_0\left(x\kappa \delta x(t)\mathrm{cos}(qy)d(t)\right)e^{iq_0\delta x(t)\mathrm{cos}(qy)}`$ (26)
$`B`$ $`=`$ $`B_0\left(x\kappa \delta x(t)\mathrm{cos}(qy)d(t)\right)`$ (27)
where

$$\delta x(t)=\delta x_0e^{\sigma t},\sigma =\frac{\xi _0^2}{4q_0^2}q^4,\mathrm{and}\kappa =(_xA_0)^2𝑑x/(_xB_0)^2𝑑x.$$
(28)
The quantity $`d(t)`$ represents the net distance traveled by the grain boundary at time $`t`$. A linear stability analysis based on (26)-(27) along the lines of Section III leads to $`d(t)=0`$ for an interface initially located at the origin. As we show below, the net displacement depends nonlinearly on the initial modulation, and hence cannot be obtained from a linearized set of equations. We use here a different method that has been widely used to study of defect motion in potential systems . It is based on the identity,

$$\frac{dF}{dt}=2𝑑\stackrel{}{r}\left(\left|\frac{A}{t}\right|^2+\left|\frac{B}{t}\right|^2\right),$$
(29)
directly derived from Eqs. (2)-(4).

The left hand side of Eq. (29) can be calculated approximately by substituting Eqs. (26) and (27) into Eq. (4), and taking the time derivative. After some algebra, we find,

$$\frac{dF}{dt}=𝑑\stackrel{}{r}\left\{\dot{d}_x_0+\frac{\xi _0^2q^4}{4}\mathrm{cos}^2(qy)\left[2\delta x\delta \dot{x}A_0^2(u)\delta x^2\dot{d}_xA_0^2(u)\right]\right\},$$
(30)
up to order $`\delta x^2`$. $`_0`$ is the free energy density of a planar boundary, and $`u=x\kappa \delta x(t)\mathrm{cos}(qy)d(t)`$. The term involving $`_0`$ vanishes since $`_0(x=+\mathrm{})=_0(x=\mathrm{})`$. In Eq. (30), we have neglected the contributions from the terms $`_yA_0(u)`$ and $`_yB_0(u)`$; they are of order $`ϵ^{1/2}`$ smaller than $`A_0(u)_ye^{i\stackrel{~}{\varphi }}`$ because of the slow variations of the amplitudes compared with the roll periodicity. In addition, we have used the approximation $`_tA_0(u)\dot{d}_xA_0(u)`$ by neglecting $`\kappa \delta \dot{x}(t)\mathrm{cos}(qy)_xA_0(u)`$: the contributions proportional to this last term vanish with $`\mathrm{cos}(qy)𝑑y`$ at leading order in $`ϵ`$ when spatial integration over $`y`$ is performed in Eq. (30). For the same reason, all terms proportional to $`\delta x`$ in Eq. (30) do not contribute.

By using similar considerations, the right hand side of Eq. (29) at leading order in $`ϵ`$ and $`\delta x(t)`$ reads,

$$2𝑑\stackrel{}{r}\left\{\dot{d}^2[(_xA_0(u))^2+(_xB_0(u))^2]+(\delta \dot{x})^2q_0^2\mathrm{cos}^2(qy)A_0^2(u)\right\}.$$
(31)
By combining Eqs. (30) and (31), and by using Eq. (28) for $`\delta x(t)`$, we find that all terms that do not involve $`\dot{d}`$ cancel, and obtain

$$\dot{d}(t)=\left(\frac{\xi _0^2}{4q_0^2}q^4\right)\frac{ϵ}{4g}\frac{[q_0\delta x(t)]^2}{_{\mathrm{}}^{\mathrm{}}𝑑x\left[(_xA_0)^2+(_xB_0)^2\right]},$$
(32)
where we have used the identities $`A_0^2(+\mathrm{})=ϵ/g`$ and $`_0^L𝑑y\mathrm{cos}^2(qy)/L=1/2`$, where $`L`$ is the length of the grain boundary. Equation (32) shows that $`\dot{d}>0`$; i.e., in the configuration shown in Fig. 4 the motion of the boundary is such that the modulated rolls are progressively replaced by the advancing perpendicular straight rolls. Furthermore, and contrary to the usual motion of kinks, the velocity is not constant but decreases as $`\delta x(t)^2`$. Since at long times $`\delta x(t)0`$, the boundary eventually stops. In Appendix A, we show that Eq. (32) generalizes the equation describing the uniform motion of a grain boundary through a time independent modulated background. The method presented in Appendix A also provides a simpler alternative way to extend the result (32) to the case in which the transversally modulated rolls are not parallel to the grain boundary.

Equations (28) and (32) can now be used to calculate the total distance traveled by the grain boundary $`d_{\mathrm{}}`$, allowing further comparison with the numerical results of Section IV. We find,

$$d_{\mathrm{}}=\frac{ϵ}{8g}\frac{(q_0\delta x_0)^2}{_{\mathrm{}}^{\mathrm{}}𝑑x\left[(_xA_0)^2+(_xB_0)^2\right]}.$$
(33)
Note that, in contrast to $`\dot{d}`$, $`d_{\mathrm{}}`$ does not depend on the wavenumber of the modulation $`q`$, for small $`q`$. Since $`A_0(x)=\sqrt{ϵ/g}f(\sqrt{ϵ}x/\xi _0)`$ at leading order in $`ϵ`$, and if a similar form is assumed for $`B_0`$, it is easy to see from Eqs. (32) and (33) that $`\dot{d}`$ and $`d_{\mathrm{}}`$ are proportional to $`ϵ^{1/2}`$ . Hence the displacement of the boundary diverges close to onset, regardless of the initial amplitude of the modulation. Eqs. (32) and (33) also show that the magnitude of this displacement depends quadratically on the perturbation $`\delta x_0`$, explaining why it can not be derived from the linear analysis presented in Section III. Finally, we note that the velocity $`\dot{d}(t)`$ is proportional to $`ϵ^{1/2}\delta x(t)^2q^4`$, at leading order in $`ϵ`$, $`\delta x(t)`$ and $`q`$.

The mechanism for net grain boundary motion can be understood as the follows: Relaxation requires energy decrease while straight rolls have the same energy regardless of their orientation (both base states A and B, are in principle degenerate). However, when the boundary is modulated in the way depicted in Fig. 4, only the rolls parallel to the boundary are distorted at lowest order in $`ϵ`$, whereas the perpendicular rolls remain straight. Therefore, the energy of the modulated configuration can decrease by both linear relaxation of the initially curved rolls and by nonlinear net displacement of the grain boundary, the effect of which is to replace curved rolls by straight ones. Close to onset, relaxation through the nonlinear mode is dominant, and even a small perturbation $`\delta x_0`$ can induce a large boundary displacement $`d_{\mathrm{}}\delta x_0`$. Of course, for fixed $`ϵ`$, $`d_{\mathrm{}}`$ goes to zero with $`\delta x_0`$. It is precisely the observation that the amplitude of this nonlinear mode diverges at $`ϵ0`$ that suggests that the linearized problem discussed in Section III may not have any bounded eigenfunctions.

At fixed $`ϵ`$, it is instructive to obtain the magnitude of $`\delta x_0`$ above which the velocity of the boundary is much larger than the phase velocity of the rolls; i.e.,

$$\dot{d}(t=0)|\delta \dot{x}(t=0)|.$$
(34)
In this limit, boundary motion occurs at approximately constant roll phase. Given that $`\delta \dot{x}(0)=\sigma \delta x_0`$, and by using Eq. (32) for $`\dot{d}(0)`$, Eq. (34) leads to

$$\delta x_0\left(\frac{\lambda _0}{2\pi }\right)^2\frac{4g}{ϵ}_{\mathrm{}}^{\mathrm{}}𝑑x\left[(_xA_0)^2+(_xB_0)^2\right].$$

With the numerical solutions obtained for $`A_0(x)`$ and $`B_0(x)`$ from the coupled Ginzburg-Landau equations, we find that $`\delta x_0\mathrm{5.0\; 10}^2\lambda _0`$ for $`ϵ=0.04`$, and $`\delta x_0\mathrm{7.4\; 10}^2\lambda _0`$ for $`ϵ=0.1`$. Therefore we conclude that nonlinear motion dominates linear relaxation even for relatively small initial perturbations. We recall that our perturbative calculation is valid if $`\delta x_0<\lambda _0`$.

We have verified that $`d_{\mathrm{}}`$ does not appreciably depend on the wavelength of the modulation, as shown in Eq. (33). Figure 4 shows our analytical and numerical results of the total distance traveled by the grain boundary as a function of $`\delta x_0`$. The theoretical curve has no adjustable parameters, and is based on the numerical resolution of Eqs. (5)-(6) with the values $`g=3/4`$ and $`g_{}=3/2`$ that correspond to the Swift-Hohenberg equation. We performed two series of runs, for $`ϵ=0.04`$ and $`ϵ=0.1`$. At $`ϵ=0.04`$, Eq. (33) is in quantitative agreement with the numerical data for small modulations. The discrepancy observed at large modulation indicates the breakdown of the linear relaxation of $`\delta x(t)`$: perturbations start decaying at a faster rate than (22), so that Eq. (33) is over-estimating the traveled distance at short times. At $`ϵ=0.1`$, we expect a smaller displacement $`d_{\mathrm{}}`$, for fixed $`\delta x_0`$, than in the case $`ϵ=0.04`$. This is what we observe. However, the numerical results show that $`d_{\mathrm{}}`$ does not vary continuously with $`\delta x_0`$, but rather it increases discontinuously. As it is apparent in Figure 4 (crosses), $`d_{\mathrm{}}`$ takes now only multiple values of half the wavelength of the pattern. A lower threshold is needed to displace the grain boundary, which we estimate to be $`\delta x_c\lambda _0/7`$. These observations cannot be explained by the present analysis, and reveal non-adiabatic effects that go beyond the amplitude equation formalism. Such effects have been previously observed in other studies of fronts propagation . The underlying cause is that with increasing $`ϵ`$ the variations of $`A`$ in the slow variable $`ϵ^{1/2}x`$ do not decouple with the fast variable describing the underlying periodic structure. In this case, the solution for the amplitudes depends on the position of the rolls.

## VI Conclusions

We have found an additional mechanism for grain boundary motion in layered phases which is driven by roll curvature in the vicinity of the boundary. We have analytically studied the motion of an isolated grain boundary within the amplitude equations formalism, and numerically by direct solution of the Swift-Hohenberg equation. Transverse perturbations of the rolls near the grain boundary decay exponentially in time with a rate proportional to $`q^4`$, where $`q`$ is the wavenumber of the perturbation. In addition, we have found a nonlinear mode involving net translation of the average location of the boundary. On the basis of functional minimization, this mode is preferred to simple linear relaxation near onset. The net velocity of the grain boundary has been computed, and is given in Eq. (32). This result can be interpreted as giving the velocity as the ratio between a time dependent external force imposed by curved rolls, and a drag term or mobility that depends on the local amplitude profile of the grain boundary. Close to onset, the drag term goes to zero as $`ϵ^{3/2}`$, and makes the interface velocity (and the total traveled distance) diverge as $`ϵ^{1/2}`$. We have argued that this motion cannot be explained by linear analysis.

The precise relationship between these results and the related problem of coarsening of layered phases requires further study. However, we can offer a few comments based on the results presented in this paper. If the temporal evolution of a disordered configuration were controlled by the relaxation of an ensemble of grain boundaries moving through a background of curved rolls, then Eq. (32) could be used to infer the coarsening rate. If $`v`$ is a characteristic grain boundary speed, and $`\kappa `$ is the characteristic curvature of the rolls ahead of the boundary, then Eq. (32) leads to $`v\kappa ^2`$. If self-similarity holds during the coarsening of the structure, this last relationship would imply a coarsening law $`l(t)t^{1/3}`$, where $`l(t)`$ is any measure of the linear scale of the structure. This result for the coarsening law disagrees with previous literature on the subject .

Although the configuration studied in this paper is idealized, is is conceivable that it describes one among possibly many competing mechanisms during coarsening. A numerical solution of the two dimensional Swift-Hohenberg equation close to onset shows that curved rolls are essentially immobile due to topological constraints (mostly disclinations), whereas grain boundaries move over large distances. It is likely that the motion of the latter is in part due to a background of curved rolls with a characteristic curvature which is set by the spatial distribution of disclinations. Disclinations, in turn, can be eliminated by grain boundary motion. If the characteristic length scales associated with both defects (grain boundary perimeter and disclination separation) are proportional to each other, as required by self-similarity, then $`t^{1/3}`$ would be a possible contribution to the coarsening law. Whether this mechanism would dominate the asymptotic behavior as $`l\mathrm{}`$ is being currently investigated.

###### Acknowledgements.

This research has been supported by the U.S. Department of Energy, contract No. DE-FG05-95ER14566, and also in part by the Supercomputer Computations Research Institute, which is partially funded by the U.S. Department of Energy, contract No. DE-FC05-85ER25000.

## A Grain boundary motion without phase relaxation

Equation (32) can be derived in an different way, which is similar to the calculation of the velocity of a climbing dislocation . Suppose that a set of rolls is modulated along the transverse direction, with wave vector $`\stackrel{}{q}`$, and is rigidly held with a time independent modulation amplitude $`\delta x_f`$. Hence, the system is allowed to evolve by rigid translation alone, and no roll (or boundary) relaxation is allowed. We then introduce the perturbations (26)-(27) with $`\delta x(t)=\delta x_f`$. Equation (29) can be simply recast as

$$\dot{d}𝑑y\{(x=\mathrm{},y)(x=+\mathrm{},y)\}=2\dot{d}^2𝑑\stackrel{}{r}\left[|_xA|^2+|_xB|^2\right],$$
(A1)
where $`\dot{d}`$ is the constant velocity of the grain boundary. The left-hand side of Eq. (A1) represents the free energy lost per unit time by the system as the grain boundary moves. The straight rolls, of lower free energy, replace curved ones. The difference $`dF`$ is computed that corresponds to increasing the area of the domain occupied by $`B`$ straight rolls by an amount $`\dot{d}dt`$, and by decreasing by the same amount the area of the domain occupied by $`A`$ curved rolls. The free energy gain is easily computed by using Eq. (4) with $`B(\mathrm{})=(ϵ/g)^{1/2}`$, $`A(\mathrm{})=0`$, and $`B(+\mathrm{})=0`$, $`A(+\mathrm{})=(ϵ/g)^{1/2}e^{iq_0\delta x_f\mathrm{cos}(qy)}`$. In the right-hand side of Eq. (A1), $`|_xA(B)|`$ is replaced by $`|_xA_0(B_0)|`$ at leading order. Eq. (A1) yields Eq. (32) straightforwardly, where $`\delta x(t)`$ must be replaced by $`\delta x_f`$. By using Eq. (A1), the numerator of Eq. (32) can now be interpreted as a (time dependent) external force acting on the line defect, the denominator being the drag term or friction coefficient. It is interesting to note that the more rigorous analysis of the combined effects of roll relaxation and front advance that leads to Eq. (32) gives the same expression as when roll relaxation is omitted. Beside being simpler, an additional advantage of the energy argument presented here is that it can be applied to a configuration where the set of rolls $`A`$ and $`B`$ are not parallel nor perpendicular to the grain boundary. In this case, the external force $`𝑑y\{(x=\mathrm{},y)(x=+\mathrm{},y)\}`$ remains unchanged, by rotational invariance. Only the drag term would probably differ, because of the different grain boundary profile. However, we expect that the scaling behavior $`\dot{d}\delta x^2q^4`$ still holds.

Fig. 1a

Fig. 1b

Fig. 2


format: text

text:
# Convergent variational calculation of positronium-hydrogen-atom scattering lengths

## Abstract

We present a convergent variational basis-set calculational scheme for elastic scattering of positronium atom by hydrogen atom in S wave. Highly correlated trial functions with appropriate symmetry are needed for achieving convergence. We report convergent results for scattering lengths in atomic units for both singlet ($`=3.49\pm 0.20`$) and triplet ($`=2.46\pm 0.10`$) states.

PACS Number(s): 34.90.+q, 36.10.Dr

Lately, there has been interest in the experimental  and theoretical  studies of ortho positronium (Ps) atom scattering by different neutral atomic and molecular targets. The Ps-H system is theoretically the most simple and fundamental and a complete understanding of this system is necessary before a venture to more complex targets . There have been R-matrix , close-coupling (CC)  and model-potential  calculations for Ps-H scattering. Here we present a convergent variational basis-set calculational scheme for low-energy Ps-H scattering in S wave below the lowest Ps-excitation threshold at 5.1 eV. Using this method we report numerical results for scattering length of electronic singlet and triplet states.

A recent study based on a regularized nonlocal electron-exchange model potential  yielded low-energy (total) cross sections in agreement with experiment for Ps scattering by He , Ne , Ar  and H2 . For the Ps-H system the model-potential results for S-wave singlet binding and resonance energies are in agreement with accurate variational estimates . It would be interesting to see if the model-potential result for the Ps-H singlet scattering length agrees with the present convergent calculation.

Because of the existence of two identical fermions (electrons) in the Ps-H system, one needs to antisymmetrize the full wave function before attempting a solution of the scattering problem. The position vectors of the electrons $``$ $`𝐫_\mathrm{𝟏}`$ (Ps) and $`𝐫_\mathrm{𝟐}`$ (H) $``$ and positron ($`𝐱`$) with respect to (w.r.t.) the massive proton at the origin are as shown in figure 1. We also use the position vectors $`𝐬_j=(𝐱+𝐫_j)/2`$, $`\rho _j=𝐱𝐫_j`$, $`j=1,2`$, $`𝐫_{12}=𝐫_1𝐫_2`$. The fully antisymmetric state $`\psi _𝐤^𝒜`$ of Ps-H scattering is given by $`|\psi _𝐤^𝒜=𝒜_1|\psi _𝐤^1=(1\pm P_{12})|\psi _𝐤^1=|\psi _𝐤^1\pm |\psi _𝐤^2`$ where $`𝐤`$ is the incident momentum, the antisymmetrizer $`𝒜_1`$ is $`(1+P_{12})`$ for the singlet state and $`(1P_{12})`$ for the triplet state with $`P_{12}`$ the permutation operator of electrons 1 and 2. The function $`\psi _𝐤^1`$ refers to the Ps-H wave function with electron 1 forming the Ps as in figure 1 and $`\psi _𝐤^2`$ refers to the same with the two electrons interchanged.

The full Ps-H Hamiltonian $`H`$ can be broken in two convenient forms as follows $`H=H_1+V_1=H_2+V_2`$ where $`H_1`$ includes the full kinetic energy and intracluster interaction of H and Ps for the arrangement shown in figure 1 and $`V_1`$ is the sum of the intercluster interaction between H and Ps in the same configuration, $`H_2`$ and $`V_2`$ refer to the same quantities with the two electrons interchanged:

$$V_1=[\frac{1}{x}\frac{1}{r_1}+\frac{1}{r_{12}}\frac{1}{\rho _2}],V_2=P_{12}V_1=[\frac{1}{x}\frac{1}{r_2}+\frac{1}{r_{12}}\frac{1}{\rho _1}].$$
(1)
The fully antisymmetric state satisfies the Lippmann-Schwinger equation

$$|\psi _𝐤^𝒜=|\varphi _𝐤^1+G_1V_1|\psi _𝐤^𝒜.$$
(2)
where the channel Green’s function $`G_1(E+i0H_1)^1`$ and the incident wave $`|\varphi _𝐤^1`$ satisfies $`(EH_1)|\varphi _𝐤^1=0.`$ The incident Ps energy $`E=6.8k^2`$ eV. We are using atomic units (au) in which $`a_0=e=m=\mathrm{}=1`$, where $`e`$ ($`m`$) is the electronic charge (mass) and $`a_0`$ the Bohr radius. Using the definition of the antisymmetrized state we rewrite (2) as

$$|\psi _𝐤^1=|\varphi _𝐤^1+G_1M_1|\psi _𝐤^1$$
(3)
$$M_1=V_1𝒜_1+(EH_1)(1𝒜_1)𝒜_1V_1+(1𝒜_1)(EH_1).$$
(4)
The properly symmetrized transition matrix for elastic scattering is defined by $`\varphi _𝐤^1|T^𝒜|\varphi _𝐤^1=\varphi _𝐤^1|V_1|\psi _𝐤^𝒜=\varphi _𝐤^1|V_1𝒜_1|\psi _𝐤^1=\psi _𝐤^1|𝒜_1V_1|\varphi _𝐤^1`$ . A basis-set calculational scheme for the transition matrix can be obtained from the following expression

$`\varphi _𝐤^1|T^𝒜|\varphi _𝐤^1=\psi _𝐤^1|𝒜_1V_1|\varphi _𝐤^1+\varphi _𝐤^1|𝒜_1V_1|\psi _𝐤^1\psi _𝐤^1|𝒜_1V_1M_1G_1𝒜_1V_1|\psi _𝐤^1.`$ (5)
Using (3), it can be verified that (5) is an identity if exact scattering wave fumctions $`\psi _𝐤^1`$ are used. If approximate wave functions are used, (5) is stationary w.r.t. small variations of $`|\psi _𝐤^1`$ but not with $`\psi _𝐤^1|`$. This one-sided variational property emerges because of the lack of symmetry of the formulation in the presence of explicit antisymmetrization operator $`𝒜_1`$. However, this variational property can be used to formulate a basis-set calculational scheme with the following trial functions

$$|\psi _𝐤^1=\underset{n=1}{\overset{N}{}}a_n|f_n,\psi _𝐤^1|=\underset{m=1}{\overset{N}{}}b_mf_m|.$$
(6)
Substituting (6) into (5) and using this variational property w.r.t. $`|\psi _𝐤^1`$ we obtain

$$\psi _𝐤^1|=\underset{m=1}{\overset{N}{}}\varphi _𝐤^1|𝒜_1V_1|f_nD_{nm}f_m|$$
(7)
$$(D^1)_{mn}=f_m|𝒜_1V_1[𝒜_1V_1+(1𝒜_1)(EH_1)]G_1𝒜_1V_1|f_n.$$
(8)
Using the variational form (7) and definition $`\varphi _𝐤^1|T^𝒜|\varphi _𝐤^1=\psi _𝐤^1|𝒜_1V_1|\varphi _𝐤^1`$ we obtain the following basis-set calculational scheme for the transition matrix

$$\varphi _𝐤^1|T^𝒜|\varphi _𝐤^1=\underset{m,n=1}{\overset{N}{}}\varphi _𝐤^1|𝒜_1V_1|f_nD_{nm}f_m|𝒜_1V_1|\varphi _𝐤^1.$$
(9)
(8) and (9) are also valid for the K matrix and in partial-wave form where the momentum-space integration over the Green’s function $`G_1`$ should be performed with the principal-value prescription.

In the calculation, the basis functions are taken in the following form

$$f_m(𝐫_2,\rho _1,𝐬_1)=\phi (r_2)\eta (\rho _1)e^{\delta _mr_2\alpha _m\rho _1\beta _ms_1\gamma _m(\rho _2+r_{12})\mu _m(x+r_1)}\frac{\mathrm{sin}(ks_1)}{ks_1}$$
(10)
where $`\phi (r)=\mathrm{exp}(r)/\sqrt{\pi }`$ and $`\eta (\rho )=\mathrm{exp}(0.5\rho )/\sqrt{8\pi }`$ represent the H(1s) and Ps(1s) wave functions, respectively. For elastic scattering the direct Born amplitude is zero and the exchange correlation dominates scattering. To be consistent with this, the direct terms in the form factors $`f_m|𝒜_1V_1|\varphi _𝐤^1`$ and $`\varphi _𝐤^1|𝒜_1V_1|f_n`$ are zero with the above choice of correlations in the basis functions via $`\gamma _m`$ and $`\mu _m`$. This property follows as the above function is invariant w.r.t. the interchange of $`𝐱`$ and $`𝐫_1`$ whereas the remaining part of the integrand in the direct terms changes sign under this transformation. A proper choice of the correlation parameters $`\gamma _m`$ and $`\mu _m`$ is crucial for obtaining good convergence.

In the following we specialize to the K-matrix formulation in S wave at zero energy, when $`\mathrm{sin}(ks_1)/(ks_1)=1`$ in (10). The useful matrix elements of the present approach are explicitly written as

$$\varphi _p^1|𝒜_1V_1|f_n=\pm \frac{1}{2\pi }\phi (𝐫_1)\eta (\rho _2)\frac{\mathrm{sin}ps_2}{ps_2}[V_1]f_n(𝐫_2,\rho _1,𝐬_1)𝑑𝐫_2𝑑\rho _1𝑑𝐬_1$$
(11)
$$f_m|𝒜_1V_1|\varphi _p^1=\pm \frac{1}{2\pi }f_m(𝐫_1,\rho _2,𝐬_2)[V_1]\phi (𝐫_2)\eta (\rho _1)\frac{\mathrm{sin}ps_1}{ps_1}𝑑𝐫_2𝑑\rho _1𝑑𝐬_1$$
(12)
$$f_m|𝒜_1V_1|f_n=\pm \frac{1}{4\pi }f_m(𝐫_1,\rho _2,𝐬_2)[V_1]f_n(𝐫_2,\rho _1,𝐬_1)𝑑𝐫_2𝑑\rho _1𝑑𝐬_1$$
(13)
$`f_m|M_1G_1𝒜_1V_1|f_n{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑pf_m|𝒜_1V_1|\varphi _p^1\varphi _p^1|𝒜_1V_1|f_n`$ (14)
where the so called off-shell term $`(1𝒜_1)(EH_1)`$ has been neglected for numerical simplification in this calculation. This term is expected to contribute to refinement over the present calculation. In this convention the on-shell K-matrix element at zero energy is the scattering length: $`a=\varphi _0^1|K^𝒜|\varphi _0^1.`$

All the matrix elements above can be evaluated by a method presented in . We describe it in the following for $`\varphi _p^1|𝒜_1V_1|f_n`$ of (11). By a transformation of variables from $`(𝐫_2,\rho _1,𝐬_1)`$ to $`(𝐬_1,𝐬_2,𝐱)`$ with Jacobian $`2^6`$ and separating the radial and angular integrations, the form factor (11) is given by

$`\varphi _p^1|𝒜_1V_1|f_n`$ $`=`$ $`\pm {\displaystyle \frac{2^6}{16\pi ^3}}{\displaystyle _0^{\mathrm{}}}s_2^2𝑑s_2{\displaystyle \frac{\mathrm{sin}(ps_2)}{ps_2}}{\displaystyle _0^{\mathrm{}}}s_1^2𝑑s_1e^{\beta _ns_1}{\displaystyle _0^{\mathrm{}}}x^2𝑑xe^{\mu _nx}`$ (15)
$`\times `$ $`{\displaystyle e^{(ar_1+b\rho _1/2)}e^{(cr_2+d\rho _2/2)}e^{\gamma _nr_{12}}[V_1]𝑑\widehat{s}_1𝑑\widehat{s}_2𝑑\widehat{x}}`$ (16)
where $`a=1+\mu _n`$, $`b=2\alpha _n+1`$, $`c=1+\delta _n`$ and $`d=2\gamma _n+1`$. Recalling that $`𝐫_j=2𝐬_j𝐱`$, $`𝐫_{12}=2(𝐬_1𝐬_2)`$, $`\rho _j=2(𝐱𝐬_j),j=1,2`$, we employ the following expansions of the exponentials in (16)

$`e^{a|2𝐬𝐱|b|𝐱𝐬|}={\displaystyle \frac{4\pi }{sx}}{\displaystyle \underset{lm}{}}G_l^{(a,b)}(s,x)Y_{lm}^{}(\widehat{s})Y_{lm}(\widehat{x})`$ (17)
$`{\displaystyle \frac{e^{a|2𝐬𝐱|b|𝐱𝐬|}}{|2𝐬𝐱|}}={\displaystyle \frac{4\pi }{sx}}{\displaystyle \underset{lm}{}}J_l^{(a,b)}(s,x)Y_{lm}^{}(\widehat{s})Y_{lm}(\widehat{x})`$ (18)
$`{\displaystyle \frac{e^{a|2𝐬𝐱|b|𝐱𝐬|}}{|𝐬𝐱|}}={\displaystyle \frac{4\pi }{sx}}{\displaystyle \underset{lm}{}}K_l^{(a,b)}(s,x)Y_{lm}^{}(\widehat{s})Y_{lm}(\widehat{x})`$ (19)
$`{\displaystyle \frac{e^{a|𝐬_\mathrm{𝟏}𝐬_\mathrm{𝟐}|}}{|𝐬_\mathrm{𝟏}𝐬_\mathrm{𝟐}|}}={\displaystyle \frac{4\pi }{s_1s_2}}{\displaystyle \underset{lm}{}}A_l^{(a)}(s_1,s_2)Y_{lm}^{}(\widehat{s}_1)Y_{lm}(\widehat{s}_2)`$ (20)
$`e^{a|𝐬_\mathrm{𝟏}𝐬_\mathrm{𝟐}|}={\displaystyle \frac{4\pi }{s_1s_2}}{\displaystyle \underset{lm}{}}B_l^{(a)}(s_1,s_2)Y_{lm}^{}(\widehat{s}_1)Y_{lm}(\widehat{s}_2)`$ (21)
where the $`Y_{lm}`$’s are the usual spherical harmonics. Using (17) $``$ (21) in (16) we get

$`\varphi _p^1|𝒜_1V_1|f_n`$ $`=`$ $`\pm 2^8{\displaystyle _0^{\mathrm{}}}e^{\beta _ns_1}𝑑s_1{\displaystyle _0^{\mathrm{}}}𝑑s_2{\displaystyle \frac{\mathrm{sin}(ps_2)}{ps_2}}{\displaystyle _0^{\mathrm{}}}𝑑xe^{\mu _nx}{\displaystyle \underset{l=0}{\overset{L}{}}}(2l+1)`$ (22)
$`\times `$ $`[{\displaystyle \frac{1}{x}}G_l^{(a,b)}(s_1,x)G_l^{(c,d)}(s_2,x)B_l^{(2\gamma _n)}(s_1,s_2)J_l^{(a,b)}(s_1,x)`$ (23)
$`\times `$ $`G_l^{(c,d)}(s_2,x)B_l^{(2\gamma _n)}(s_1,s_2)+{\displaystyle \frac{1}{2}}G_l^{(a,b)}(s_1,x)G_l^{(c,d)}(s_2,x)`$ (24)
$`\times `$ $`A_l^{(2\gamma _n)}(s_1,s_2){\displaystyle \frac{1}{2}}G_l^{(a,b)}(s_1,x)K_l^{(c,d)}(s_2,x)B_l^{(2\gamma _n)}(s_1,s_2)].`$ (25)
where the $`l`$-sum is truncated to $`l=L`$. This evaluation avoids complicated angular integrations involving $`𝐬_1`$, $`𝐬_2`$ and $`𝐱`$. These integrals take a simple form requiring straightforward numerical computation of certain radial integrals only, which must, however, be carried out carefully. The functions $`G_l`$, $`J_l`$, $`K_l`$ etc. are easily calculated using (17) $``$ (21):

$$G_l^{(a,b)}(s,x)=\frac{sx}{2}_1^1𝑑tP_l(t)e^{a|2𝐬𝐱|b|𝐱𝐬|}$$
(26)
where $`P_l(t)`$ is the usual Legendre polynomial and $`t`$ is the cosine of the angle between $`𝐬`$ and $`𝐱`$. The integrals (12) and (13) can be evaluated similarly. For example

$`f_m|𝒜_1V_1|\varphi _p^1`$ $`=`$ $`\pm 2^8{\displaystyle _0^{\mathrm{}}}e^{\beta _ms_2}𝑑s_2{\displaystyle _0^{\mathrm{}}}𝑑s_1{\displaystyle \frac{\mathrm{sin}(ps_1)}{ps_1}}{\displaystyle _0^{\mathrm{}}}𝑑xe^{\mu _mx}{\displaystyle \underset{l=0}{\overset{L}{}}}(2l+1)`$ (27)
$`\times `$ $`[{\displaystyle \frac{1}{x}}G_l^{(c,d)}(s_1,x)G_l^{(a,b)}(s_2,x)B_l^{(2\gamma _m)}(s_1,s_2)J_l^{(c,d)}(s_1,x)`$ (28)
$`\times `$ $`G_l^{(a,b)}(s_2,x)B_l^{(2\gamma _m)}(s_1,s_2)+{\displaystyle \frac{1}{2}}G_l^{(c,d)}(s_1,x)G_l^{(a,b)}(s_2,x)`$ (29)
$`\times `$ $`A_l^{(2\gamma _m)}(s_1,s_2){\displaystyle \frac{1}{2}}G_l^{(c,d)}(s_1,x)K_l^{(a,b)}(s_2,x)B_l^{(2\gamma _m)}(s_1,s_2)]`$ (30)
$`f_m|𝒜_1V_1|f_n`$ $`=`$ $`\pm 2^7{\displaystyle _0^{\mathrm{}}}e^{\beta _ns_1}𝑑s_1{\displaystyle _0^{\mathrm{}}}𝑑s_2e^{\beta _ms_2}{\displaystyle _0^{\mathrm{}}}𝑑xe^{(\mu _n+\mu _m)x}{\displaystyle \underset{l=0}{\overset{L}{}}}(2l+1)`$ (31)
$`\times `$ $`[{\displaystyle \frac{1}{x}}G_l^{(e,f)}(s_1,x)G_l^{(g,h)}(s_2,x)B_l^{(2\gamma _{mn})}(s_1,s_2)J_l^{(e,f)}(s_1,x)`$ (32)
$`\times `$ $`G_l^{(g,h)}(s_2,x)B_l^{(2\gamma _{mn})}(s_1,s_2)+{\displaystyle \frac{1}{2}}G_l^{(e,f)}(s_1,x)G_l^{(g,h)}(s_2,x)`$ (33)
$`\times `$ $`A_l^{(2\gamma _{mn})}(s_1,s_2){\displaystyle \frac{1}{2}}G_l^{(e,f)}(s_1,x)K_l^{(g,h)}(s_2,x)B_l^{(2\gamma _{mn})}(s_1,s_2)]`$ (34)
where $`e=1+\delta _m+\mu _n`$, $`f=2\alpha _n+2\gamma _m+1`$, $`g=1+\delta _n+\mu _m`$, $`h=2\alpha _m+2\gamma _n+1`$ and $`\gamma _{mn}=\gamma _m+\gamma _n`$.

We tested the convergence of the integrals by varying the number of integration points in the $`x`$, $`s_1`$ and $`s_2`$ integrals in (22), (27) and (31) and the $`t`$ integral in (26). The evaluation of (22), (27) and (31) essentially involves four-dimensional integration, which is performed with caution. The $`x`$ integration was relatively easy and 20 Gauss-Legendre quadrature points appropriately distributed between 0 and 16 were enough for convergence. In the evaluation of integrals of type (26) 40 Gauss-Legendre quadrature points were sufficient for adequate convergence. The convergence in the numerical integration over $`s_1`$ and $`s_2`$ was achieved with 300 Gauss-Legendre quadrature points between 0 and 12. The maximum value of $`l`$ in the sum in (22), (27) and (31), $`L`$, is taken to be 6 which is sufficient for obtaining the convergence with the partial-wave expansions (17) $``$ (21).

Table 1: Singlet ($`a_s`$) and triplet ($`a_t`$) Ps-H scattering lengths for different $`L`$ and $`N`$.

|  | $`L=`$ | 0 | $`L=`$ | 2 | $`L=`$ | 4 | $`L=`$ | 6 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`N`$ | $`a_t`$ | $`a_s`$ | $`a_t`$ | $`a_s`$ | $`a_t`$ | $`a_s`$ | $`a_t`$ | $`a_s`$ |
| 6 | 4.00 | 3.61 | 3.66 | 4.00 | 3.07 | 3.75 | 2.92 | 3.74 |
| 7 | 3.98 | 4.11 | 3.25 | 3.88 | 2.99 | 4.06 | 2.85 | 4.00 |
| 8 | 3.93 | 4.12 | 3.35 | 3.96 | 2.66 | 3.72 | 2.54 | 3.73 |
| 9 | 3.97 | 4.15 | 3.42 | 3.83 | 2.65 | 3.92 | 2.55 | 4.06 |
| 10 | 3.97 | 4.22 | 3.42 | 3.92 | 2.57 | 4.42 | 2.50 | 3.45 |
| 11 | 4.01 | 3.82 | 3.44 | 3.91 | 3.55 | 3.75 | 2.67 | 3.80 |
| 12 |  |  |  |  | 2.54 | 3.72 | 2.46 | 3.73 |
| 13 |  |  |  |  | 2.48 | 3.47 | 2.46 | 3.49 |

In a numerical calculation a judicial choice of the parameters in (10) is needed for rapid convergence. As the variational method does not provide a bound on the result, the method could converge to a wrong scattering length if an inappropriate (incomplete) basis set is choosen. After some experimentation we find that to obtain proper convergence the parameters $`\delta _n`$ and $`\alpha _n`$ should be taken to have both positive and negative values, $`\gamma _n`$ and $`\mu _n`$ should include values close to unity and $`\beta _n`$ should have progressively increasing values till about 1.5. The results reported in this work are obtained with the following parameters for the functions $`f_n,n=1,\mathrm{},13`$: $`\{\delta _n,\alpha _n,\beta _n,\gamma _n,\mu _n\}\{0.5,0.25,0.3,0.01,0.02\},`$ $`\{0.5,0.25,0.5,0.04,0.02\},`$ $`\{0.5,0.25,0.7,0.03,0.06\},`$ $`\{0.2,0.1,0.6,0.2,0.2\},`$ $`\{0.1,0.1,0.8,0.25,0.25\},`$ $`\{0.2,0.2,0.6,0.35,0.35\},\{0.1,0.1,0.7,0.4,0.4\},\{0.15,0.2,0.8,0.5,0.5\},`$ $`\{0.12,0.12,1,0.7,0.7\},`$ $`\{0.2,0.2,1.2,0.9,0.9\},`$ $`\{0.1,0.2,1.3,1,0.7\},`$ $`\{0.2,0.1,1.4,0.7,1\},`$ $`\{0.3,0.15,1.5,1,1\}`$

In table 1 we show the convergence pattern of the present calculation w.r.t. the number of partial waves $`L`$ and basis functions $`N`$ used in the calculation. The convergence is smooth with increasing $`L`$. However, as the present calculation does not produce a bound on the result, the convergence is not monotonic with increasing $`N`$. The lack of an upper bound in this calculation is clearly revealed in table 1 where the results do not decrease monotonically as the number of terms in the trial wave function is increased. The unbounded nature of the results is consistent with the noted oscillation of the scattering lengths as $`N`$ increases. The oscillation is larger in the singlet state where it is more difficult to obtain convergence. Although, the final result for the largest $`N`$ and $`L`$ is supposed to be the most accurate, it is not known whether this result is larger or smaller than the exact one. Also, an estimate of error of this result is not known. It is difficult to provide a quantitative measure of convergence. However, from the noted fluctuation of the results for large $`N`$ for $`L=6`$ we believe the error in the triplet scattering length to be less than 0.10 and in the singlet scattering length to be less than 0.20. The final results of the present calculation are those for $`N=13`$ and $`L=6`$ with the above estimate of error: $`a_s=3.49\pm 0.20`$ au and $`a_t=2.46\pm 0.10`$ au.

The maximum number of functions ($`N=13`$) used in this calculation is also pretty small, compared to those used in different Kohn-type variational calculations for electron-hydrogen ($`N=56`$)  and positron-hydrogen ($`N286`$) scattering. Because of the explicit appearance of the Green’s function, the present basis-set approach is similar to the Schwinger variational method. Using the Schwinger method, convergent results for electron-hydrogen  and positron-hydrogen  scattering have been obtained with a relatively small basis set ($`N10`$). These suggest a more rapid convergence in these problems with a Schwinger-type method.

Now we compare the present results with those of other calculations. While a static-exchange calculation by Hara and Fraser  yielded $`a_s=7.28`$ au and $`a_t=2.48`$ au a model potential calculation by Drachman and Houston  produced $`a_s=5.33`$ au and $`a_t=2.36`$ au. A 22-state R-matrix calculation by Campbell et al.  and a six-state CC calculation by Sinha et al.  yielded $`a_s=5.20`$ au, $`a_t=2.45`$ au and $`a_s=5.90`$ au, $`a_t=2.32`$ au, respectively. These triplet scattering lengths are in qualitative agreement with the present result: $`a_t=2.46\pm 0.10`$. However, the singlet scattering lengths of these calculations have not yet converged. The disagreement in the singlet channel shows that it is more difficult to get the converged result in the attractive singlet channel than in the repulsive triplet channel. This is consistent with the common wisdom that a scattering model at low energies is more sensitive to the detail of an effective attractive interaction than to that of a repulsive interaction. The nonconvergence of these results for the singlet scattering length was conjectured before . Using a correlation between the S-wave scattering length and binding energy for the Ps-H system, the value $`a_s=3.5`$ au was predicted in that study  in excellent agreement with the present result: $`a_s=3.49\pm 0.20`$.

To summarize, we have formulated a convergent basis-set calculational scheme for S-wave Ps-H elastic scattering below the lowest inelastic threshold using a variational expression for the transition matrix. We illustrate the method numerically by calculating the singlet and triplet scattering lengths: $`a_s=3.49\pm 0.20`$ au and $`a_t=2.46\pm 0.10`$ au.

The work is supported in part by the Conselho Nacional de Desenvolvimento - Científico e Tecnológico, Fundação de Amparo à Pesquisa do Estado de São Paulo, and Financiadora de Estudos e Projetos of Brazil.

Figure Caption:

1. Different position vectors for the Ps-H system w.r.t the proton (p) at the origin in arrangement 1 with electron 1 forming the Ps.


text:
# Modification of the standard model for the lanthanides

## Abstract

We show that incorporation of strong electron correlations into the Kohn-Sham scheme of band structure calculations leads to a modification of the standard model of the lanthanides and that this procedure removes the existing discrepancy between theory and experiment concerning the ground state properties. Within the picture suggested, part of the upper Hubbard $`f`$-band is occupied due to conduction band-$`f`$-mixing interaction (that is renormalized due to correlations) and this contributes to the cohesive energy of the crystal. The lower Hubbard band has zero width and describes fermionic excitations in the shell of localized $`f`$-s. Fully self-consistent calculations (with respect to both charge density and many-electron population numbers of the $`f`$-shell) of the equilibrium volume $`V_0`$ and the bulk modulus of selected lanthanides have been performed and a good agreement is obtained.

Published as: Sol. State Commun. 115, 7-12 (2000).
For strongly correlated systems, state-of-the-art band structure methods fail to give an accurate description. The reason for the failure is well-known: the strong intra-shell Hubbard repulsion, $`U`$, is underestimated in the calculations. If the $`f`$-$`f`$ interaction is sufficiently strong, the $`f`$-shell forms a localized multiplet (Russel-Saunders coupling) and can accurately be treated as atomic like. This is the basic idea behind the so-called standard model of the lanthanides and many of the physical properties seem consistent with it, for example the Curie-Weiss behavior of the susceptibility , the gross features of the equilibrium volumes , and the structural properties . When using ab initio calculations employing the local density approximation (LDA) and treating the $`f`$-electrons as core states, a too large equilibrium volume, $`V_0`$, is obtained (for all other elements, LDA underestimates the equilibrium volume). On the other hand, when the 4$`f`$’s are treated as itinerant electrons, $`V_0`$ is much too low compared to experiment and the observed localized moment (Curie-Weiss law) is absent. Calculations by Delin et al.  show that treating the 4$`f`$-electrons as core states is a very good approximation for the late lanthanides, but for the lighter lanthanides the disagreement between theory and experiment is gradually increasing when one goes from heavy to the lighter lanthanides. Thus, even though the standard model is essentially correct, it has to be modified slightly to better describe the cohesive properties of the lanthanides. Such a modification requires a method combining an ab initio band structure calculations (AIBSC) and a many-body approach. A variety of such methods are developed. The ideas of orbital polarization, LDA+U method, Hartree-Fock type of approximations, method of phase shift and more advanced methods, like three-body Faddeev equations, have in some cases been successfully applied to systems with strong electron correlations (SEC). New approaches: on the basis of the Gutzwiller wave function, dynamical mean field and the LDA++-approach have also been used for the description of different SEC systems. In many cases the AIBSC-based density of electron states is used as input to a model calculation using a non-self-consistent fixed lattice parameter and Hubbard $`U`$. However, since a large value of Hubbard repulsion restricts the available electron phase space, wave functions and, correspondingly, the matrix elements of mixing and hopping, changes. Recently an attempt to describe SEC by means of slave-boson method (i.e. Hubbard $`U=\mathrm{}`$) build into the scheme of linear muffin-tin orbital method (LMTO) was made. The author found that the results are unsatisfactory and had to use a scheme different from the standard slave-boson self-consistency schemes. However, in all these attempts the influence of SEC on the cohesive properties has not been questioned yet. Our approach to improve on the standard model is based on the following physical picture: the Hubbard repulsion is sufficiently large to form large energy gaps between the $`(n1)`$-, $`n`$\- and $`(n+1)`$-electron configurations not only in the heavy but also in the light lanthanides . Near the atomic limit all $`f`$-electrons are in the strong coupling regime and hopping and mixing interactions lead to the formation of a lower Hubbard sub-bands (LHB, $`|n1|n`$ transitions) and a upper sub-band (UHB, $`|n|n+1`$ transitions). In our case, however, the lower transitions are substantially below the bottom of the conduction bands (Ce is an exception) and, hence, they do not experience any mixing interaction. Hubbard $`U`$ increases with atomic number, due to the more localized nature of the 4$`f`$’s, thus the LHB goes down in energy while the UHB goes up and therefore it experiences smaller mixing and is populated with fewer states. These many-electron local excitations correspond to localized core electrons within the one-electron picture. The upper transitions are in the energy region where the delocalized electrons do exist. Therefore, the overlap and mixing interaction may delocalize these transitions, forming an UHB, filled only by a small amount, $`\eta `$, of electrons. This picture corresponds to the multi-orbital Hubbard-Anderson model. The total occupation number can be written as $`n_f=n+\eta `$, where $`n`$ is an integer, which is determined by the nuclear charge valence, and $`0<\eta <1`$. Thus, most of the $`f`$-occupation is in the LHB. This is consistent with the observed Curie-Weiss behavior of the high-temperature (above Tc) susceptibility. In the lanthanides the width of the appropriate $`f`$-sub band is much smaller than the Hubbard band gap. This allows one to formulate an effective Schrödinger (or Dirac) equation (SE) for two different energy regions. As is known from Hubbard’s models, the single-electron spectral weight at $`U=0`$ is split by the Hubbard repulsion $`U`$ to a set of weights, corresponding to an electron transition from one ion state to another. These weights renormalize the matrix elements of the hopping and mixing interaction, which in turn, cause a narrowing of the bands due to correlations. The theory is developed for a non-orthogonal basis set, as for instance used in the linear muffin-tin orbital method in the atomic sphere approximation (LMTO-ASA). Here we apply the formalism only in the lowest approximation which requires less complicated changes for methods used in self consistent ab initio calculation. Compared to conventional considerations of the periodical Anderson model (PAM) which do not include structure-dependent vertex corrections we solve the Schrödinger equation for the electrons in fields generated by the nuclei and electrons in every iteration, with the charge density obtained from a system of equations for Greens functions (GF) and spectral weights, thus renormalizing self-consistently the matrix elements of the PAM. On the other hand compared to the normal LDA approach we solve self-consistently additional system of equations for many-electron spectral weights and obtain a bandwidth which is reduced by correlations.

Technically our method, briefly described, consists of the following steps:
1. The total many-electron Hamiltonian is written in terms of an LMTO basis set.
2. The $`f`$-electron operators are transformed into a Hubbard representation.
3. Analytical calculations are performed for the electron Greens functions for the corresponding PAM in the LMTO representation in the Hubbard-I approximation (strong-coupling regime).
4. Renormalizing factors to the $`f`$-$`(spd)`$-blocks of the LMTO Hamiltonian and overlap matrices are introduced which are derived from a comparison of the frequency and overlap matrixes. They arise due to strong correlations in the system for Greens functions.
5. Self-consistent, ab initio calculations of the system with SEC are performed.
In practice, the main effect of these steps consists of the separation of the $`f`$-electron system into two subsystems and the correlation driven narrowing of the $`f`$-bands in the conventional LMTO Hamiltonian. This, as will be demonstrated below, leads to an improved description of the bonding of the lanthanide metals. In order to describe our theory in more detail, the secondary quantized full many-electron Hamiltonian, $``$, is written in terms of an LMTO-basis set. An unperturbed Hamiltonian and a perturbation is formulated as $`=(T+V_{ne}+V_{LDA})+(H_{ee}V_{LDA})=_0+^{}.`$ The main effect of the $`f`$-part single site correction from Coulomb interaction,

$$_U^{}=\frac{1}{2}U_{m_1m_2m_3m_4}f_{m_1\sigma }^{}f_{m_2\sigma ^{}}^{}f_{m_3\sigma ^{}}f_{m_4\sigma }v_{l=3,\sigma }^{(LDA)}f_{m_l\sigma }^{}f_{m_l\sigma },$$
(1)
is to localize part of the $`f`$-spectra, similar to a so called LDA+U solution. Let us denote all non-$`f`$-operators in the LMTO representation as $`c_{𝐤L}`$. Then the LMTO part of the Hamiltonian, $`H_{LL^{}}(𝐤)a_{𝐤L}^{}a_{𝐤L^{}}`$, can be written as a sum of $`s,p,d,`$ and $`f`$-electrons, and their mixing interaction,

$`_0={\displaystyle \underset{𝐧m_l\sigma }{}}ϵ_f^0f_{𝐧m_l\sigma }^{}f_{𝐧m_l\sigma }+{\displaystyle \underset{𝐤L,L^{}}{}}H_{L,L^{}}(𝐤)c_{𝐤L}^{}c_{𝐤L^{}}+`$ (2)
$`{\displaystyle \underset{\mathrm{𝐧𝐦},m_l\sigma m_l^{}}{}}t_{\mathrm{𝐧𝐦}}^{m_lm_l^{}}f_{𝐧m_l\sigma }^{}f_{𝐦m_l^{}\sigma }+{\displaystyle \underset{𝐤,𝐧,L,m_l\sigma }{}}[H_{L,m_l\sigma }(𝐤)e^{i\mathrm{𝐤𝐑}_𝐧}c_{𝐤L\sigma }^{}f_{𝐧m_l\sigma }+h.c.],`$ (3)
where $`t_{\mathrm{𝐧𝐦}}^{m_lm_l^{}}=_𝐤e^{i𝐤(𝐑_𝐧𝐑_𝐦)}[H_{m_l\sigma ,m_l^{}\sigma }(𝐤)\delta _{m_l,m_l^{}}ϵ_f^0]`$ is the hopping matrix element between $`f`$-orbitals and $`ϵ_f^0=_𝐤H_{m_l\sigma ,m_l\sigma }(𝐤)`$. Since the LMTO-ASA basis is non-orthogonal with an overlap matrix $`\chi _{𝐧L}|\chi _{𝐦L^{}}=𝒪_{𝐧L,𝐦L^{}}`$, the anticommutator is $`\{a_{𝐧L},a_{𝐦L^{}}^{}\}=(𝒪^1)_{𝐧L,𝐦L^{}}`$. Therefore, the equation for the bare fermion Greens function, $`G_{𝐧L,𝐦L^{}}(\omega )T[a_{𝐧L}(t)a_{𝐦L^{}}^{}(0)]_\omega ,`$ for this Hamiltonian is

$$[\omega 𝒪_{𝐧L,𝐦_\mathrm{𝟏}L_1}H_{𝐧L,𝐦_\mathrm{𝟏}L_1}]G_{𝐦_\mathrm{𝟏},L_1,𝐦L^{}}(\omega )=\delta _{𝐧,𝐦}\delta _{L,L^{}}.$$
(4)
Thus, we find (by construction) in the square brackets the conventional secular problem for the LMTO-ASA method (or any LDA based method). Introducing a full set of states in an orbital representation $`|\mathrm{\Gamma }_0)=|0,|\mathrm{\Gamma }_1)=|\gamma ,|\mathrm{\Gamma }_2)=|\gamma ,\gamma ^{},\mathrm{}`$ ($`\gamma `$ denoting a single particle $`f`$-state, $`m_l\sigma `$) allows to express any single-site operator $`\widehat{𝒜}`$ in terms of $`X`$-operators, $`\widehat{𝒜}=_{\mathrm{\Gamma },\mathrm{\Gamma }^{}}\mathrm{\Gamma }|\widehat{𝒜}|\mathrm{\Gamma }^{}X^{\mathrm{\Gamma },\mathrm{\Gamma }^{}}`$, $`X^{pq}=|pq|`$ and where $`X^{\mathrm{\Gamma },\mathrm{\Gamma }}N_\mathrm{\Gamma }`$ is the occupation number for the state $`\mathrm{\Gamma }`$. Now we can express the zero $`f`$-Greens functions, for the LMTO problem, in terms of $`X`$-Greens functions as $`G_\gamma ^{(0)}(\omega )=_a|(f_\gamma )_a|^2P_aD_a(\omega ),`$ where $`D_a(\omega )=\frac{1}{\omega \mathrm{\Delta }_a},`$ and an $`X`$-operator GF $`𝒢^{a\overline{a}}=P^aD_a`$. Here $`a=(\mathrm{\Gamma }_n,\mathrm{\Gamma }_{n+1})`$ ($`\overline{a}=(\mathrm{\Gamma }_{n+1},\mathrm{\Gamma }_n)`$) are possible Hubbard transitions, $`P_{\mathrm{\Gamma }_n,\mathrm{\Gamma }_{n+1}}=N_{\mathrm{\Gamma }_n}+N_{\mathrm{\Gamma }_{n+1}}`$, $`\mathrm{\Delta }_{\mathrm{\Gamma }_{n+1},\mathrm{\Gamma }_n}=E_{\mathrm{\Gamma }_{n+1}}E_{\mathrm{\Gamma }_n}\mu `$, $`E_{\mathrm{\Gamma }_n}nϵ_f^0+Un(n1)/2`$ and $`\mu `$ is the chemical potential.

The $`f`$-electron operators can be expanded in the Hamiltonian $``$ into $`X`$-operators for the set of many-electron states described above. Then we arrive at an Anderson like model

$`={\displaystyle \underset{𝐧\mathrm{\Gamma }}{}}E_\mathrm{\Gamma }X_𝐧^{\mathrm{\Gamma }\mathrm{\Gamma }}+{\displaystyle \underset{𝐤,L,L^{}}{}}H_{L,L^{}}(𝐤)c_{𝐤L}^{}c_{𝐤L^{}}+`$ (5)
$`{\displaystyle \underset{𝐤,𝐧,\mu ,a}{}}[H_{L,\mu }(𝐤)(f_\mu )_ae^{i\mathrm{𝐤𝐑}_𝐧}c_{𝐤L}^{}X_𝐧^a+h.c.].+{\displaystyle \underset{𝐧,𝐦,a,b}{}}t_{\mathrm{𝐧𝐦}}^{ab}X_𝐧^aX_𝐦^b.`$ (6)
Here $`t_{\mathrm{𝐧𝐦}}^{ab}=_{\mu ,\nu }(f_{𝐧\mu }^{})^a(f_{𝐦\nu }^{})^bt_{\mathrm{𝐧𝐦}}^{\mu \nu }`$. For the derivation of the equations for the Greens functions, $`𝒢^{(cX)}T[c_{𝐧L}(t)X_𝐧^{M\mathrm{\Gamma }}(0)]`$ and $`𝒢^{(XX)}T[X_𝐧^{\mathrm{\Gamma }M}(t)X_𝐧^{M\mathrm{\Gamma }}(0)]`$, (which will be compared to the equation for the regular Kohn-Sham (LMTO) problem) we need the anticommutators $`\{c,X\}`$, we find $`\{c_{𝐧L},X_𝐧^{}^{\mathrm{\Gamma }^{}\mathrm{\Gamma }}\}=_{𝐧_\mathrm{𝟏}\mu \mathrm{\Gamma }_1\mathrm{\Gamma }_2}𝒪_{𝐧L,𝐧_\mathrm{𝟏}\mu }^1(f_\mu )_{\mathrm{\Gamma }_1,\mathrm{\Gamma }_2}\{X_{𝐧_\mathrm{𝟏}}^{\mathrm{\Gamma }_1,\mathrm{\Gamma }_2},X_𝐧^{}^{\mathrm{\Gamma }^{},\mathrm{\Gamma }}\}.`$ Physically, this implies that the $`f`$-part of the conduction-band wavefunction, coming from other sites, experiences strong correlations on the reference ion. The Hubbard-I as well as mean field approximation correspond to $`\{X_𝐧^{\mathrm{\Gamma }_1,\mathrm{\Gamma }_2},X_𝐧^{\mathrm{\Gamma }^{},\mathrm{\Gamma }}\}P_{\mathrm{\Gamma }^{},\mathrm{\Gamma }}\delta _{\mathrm{\Gamma }^{},\mathrm{\Gamma }_2}\delta _{\mathrm{\Gamma },\mathrm{\Gamma }_1}`$. Introducing renormalized ”Fermi” operators $`\stackrel{~}{f}^{\mathrm{\Gamma }^{},\mathrm{\Gamma }}=X^{\mathrm{\Gamma }^{},\mathrm{\Gamma }}/\sqrt{P_{\mathrm{\Gamma }^{},\mathrm{\Gamma }}}.`$ We see that the matrix form of the GF, $`𝒢^{(ij)}`$, with $`i,j=c,X`$, coincides with the GF of the Kohn-Sham equation, $`G`$, with the only difference being $`A_{L,\mu }\stackrel{~}{A}_{L,\mu }=_{\mathrm{\Gamma }^{},\mathrm{\Gamma }}A_{L,\mu }(f_\mu )_{\mathrm{\Gamma }^{},\mathrm{\Gamma }}\sqrt{P_{\mathrm{\Gamma }^{},\mathrm{\Gamma }}}`$, where $`A=H`$ or $`𝒪`$. Thus the most essential, technical part of our theory leads to an effective secular equation which, when diagonalized gives eigenvalues with many-body corrections.

As an example let us consider praseodymium metal, which has two localized $`f`$-electrons. Therefore, $`|ϵ_f^0|`$ has a value which gives an energy minimum for $`n=2`$. For briefness we will use the following notation for the $`f`$-orbitals: $`\gamma =(m_l=3,)1,(m_l=2,)2`$ and $`\gamma =\nu `$ for $`\gamma 1,2`$. $`\mathrm{\Gamma }_{(12)}(12)`$, $`\mathrm{\Gamma }_{(12\nu )}(12\nu )`$. Let us now consider a simple polarized solution, where the $`f`$-orbitals with $`\gamma =1,2`$ are fully occupied, while the rest of the $`f`$-electrons occupy the $`f^3`$-states, $`|\mathrm{\Gamma }_3=|1,2,\nu `$. When mixing is absent, $`n_f=2`$ and $`N_{(12)}^0=1`$ whereas all other population numbers, $`N_\mathrm{\Gamma }^0`$, are zero. Then $`𝒢_{\gamma =1}^0=D_{2(12)},𝒢_{\gamma =2}^0=D_{1(12)}`$ and $`𝒢_{\gamma =\nu }^0=D_{(12)(12\nu )}.`$ We denote $`\mathrm{\Delta }_{(12)2}=\mathrm{\Delta }_{(12)1}\mathrm{\Delta }_1,`$ and $`\mathrm{\Delta }_{(12\nu )(12)}\mathrm{\Delta }_2`$. In the limiting case, $`\mathrm{\Delta }_2,\mathrm{\Delta }_3\mathrm{}`$, no $`f`$’s are present in the conduction band, and therefore this corresponds exactly to the limit of the standard model. The $`f_\nu `$-bands are slightly above the Fermi energy, and mixing interaction transfers $`f`$-character into the $`s,p,d`$-electron states. In a general case it is, of course, impossible to write equations for $`N_\mathrm{\Gamma }`$ in terms of only orbital GF’s, $`𝒢_\gamma `$, because equations for higher correlation functions are needed. However, if for simplicity we assume that these 12 $`\nu `$-bands are occupied symmetrically, i.e., $`N_{(12\nu )}=(1/12)_\nu ^{}N_{(12\nu ^{})}`$, we immediately find

$$\{\begin{array}{c}n_f=2N_{(12)}+3_\nu N_{(12\nu )}=2+\eta ,\hfill \\ N_{(12)}+_\nu N_{(12\nu )}=1,\hfill \end{array}$$
(7)
i.e., $`N_{(12\nu )}=\eta /12`$ and $`N_{(12)}=1\eta `$ (for the other rare-earth elements a slight modification of these expressions are needed). Thus, renormalization of the population numbers leads to the form of $`f`$-locators (single-site, $`𝐤`$-independent part of the GF)

$$𝒢_1=𝒢_2=\frac{1\eta }{\omega \mathrm{\Delta }_1}+\frac{\eta }{\omega \mathrm{\Delta }_2}\text{ and }𝒢_\nu =\frac{111\eta /12}{\omega \mathrm{\Delta }_2}+\frac{11\eta /12}{\omega \mathrm{\Delta }_3},$$
(8)
where $`\mathrm{\Delta }_3\mathrm{\Delta }_{(12\nu \nu ^{})(12\nu ^{})}`$. The center of the non-renormalized LDA-$`f`$-band is in $`\omega =\mathrm{\Delta }_2`$ and $`𝒢_1^{LDA}=𝒢_2^{LDA}=𝒢_\nu ^{LDA}=1/(\omega \mathrm{\Delta }_2).`$ Let us ignore for the moment the non-orthogonality. For a diagonal effective hopping $`t_\gamma (𝐤,\omega )`$, of any origin, we have $`G_\gamma =[(G_{0,\gamma }^{(at)})^1t_\gamma ]^1`$. Therefore, the expressions for the local, $`f`$-part, of the self-energies are $`\mathrm{\Sigma }_\gamma =(G_\gamma ^{LDA})^1U_\gamma ^{}(1\delta _{\gamma ,\gamma ^{}})n_\gamma ^{}(G_{0,\gamma }^{(at)})^1`$. Since all poles in $`𝒢_\gamma `$ are well separated, one can formulate two effective Schrödinger equations near $`\omega \mathrm{\Delta }_1`$ and $`\omega \mathrm{\Delta }_2`$. Renormalization factors $`(1(\frac{\mathrm{\Sigma }}{\omega })_{\omega =\mathrm{\Delta }_i})`$ are, of course, given by spectral weights (numerators) of the GF’s in these poles. Since $`\eta `$ is small the bandwidths of the upper Hubbard sub bands for $`\gamma =1,2`$ are also small and, being above the Fermi level, they are empty. Within the scenario considered the lower sub bands have weight slightly less than one, $`(1\eta )`$. For the bands $`\nu `$, $`\mathrm{\Delta }_3`$ is far above the Fermi energy. Thus the upper sub bands are empty and the spectral weight in the pole $`\omega =\mathrm{\Delta }_2`$ is $`(1\frac{11\eta }{12})`$. Therefore, the effective Schrödinger equation for bands $`1`$ and $`2`$ at $`\omega \mathrm{\Delta }_1=\mathrm{\Delta }_2U`$ gives $`\omega (\mathrm{\Delta }_2U)(1\eta )t_1(𝐤,\omega )=0`$, whereas for the bands $`\nu `$ at $`\omega \mathrm{\Delta }_2`$ we find $`\omega \mathrm{\Delta }_2(111\eta /12)t_\nu (𝐤,\omega )=0`$.
Since in the case of LDA calculations with delocalized $`f`$-electrons the center of the $`f`$-band is $`\mathrm{\Delta }_2`$ we conclude that the potential for the 1st and 2nd orbitals is shifted down by $`U`$ (similar to the LDA+U - result, and self interaction correction (SIC)-theory  but without factors $`n_\mu `$). The next conclusion is that the described factors should be to introduced into the LMTO overlap and Hamiltonian matrices. Thus, we calculate the charge density and a self consistent potential, corresponding to the given $`\eta `$, where $`\eta `$ comes from the equation of self-consistency

$$\eta =\frac{1}{\pi }\underset{\omega ,\nu ,𝐤}{}f(\omega )Im𝒢_\nu ^{(12,12\nu ),(12\nu ,12)}(\omega +i\delta ,𝐤),$$
(9)
since the parameter $`P_{\mathrm{\Gamma }^{},\mathrm{\Gamma }}=1\frac{11}{12}\eta `$ renormalizes $``$ and $`𝒪`$. Here $`𝒢_\nu ^{(12,12\nu ),(12\nu ,12)}`$ is obtained from the transformed secular matrix for the normal LMTO-GF’s, $`(\stackrel{~}{𝒪}E\stackrel{~}{})^1`$, via $`\stackrel{~}{f}`$. For Pr we find the following picture. Two of the f-orbitals are located at an energy -$`U`$ lower than the remaining $`f`$-states and below the bottom of the conduction bands. Therefore, we treat those $`f`$-states as core states removing the somewhat arbitrary choice of a value for the Hubbard $`U`$. Compared to LDA+U and SIC methods, in our approach the spd-f hybridization is reduced by the renormalizing factor $`\sqrt{P}`$. The theory is generalized for the other lanthanides straightforwardly.

Before giving numerical examples of cohesive properties (e.g. equilibrium volume) obtained from our method we describe how the total energy was calculated. As is well-known, the equation for the exchange-correlation potential $`v_{xc}`$ derived by Sham connects it with the exact self-energy of the electron system. Making use of a strong-coupling (SC) perturbation theory (PT) developed by us we have performed the analysis of the contributions to the self-energy and found that a there exist a one-to-one correspondence between the graphs for the self-energy in the standard weak-coupling (WC) PT (in e.g. the random phase approximation) and the sequence of graphs in SCPT. Using the facts that: a) Sham’s equation should produce a $`v_{xc}`$ corresponding to the particular SE choice; b) established correspondence between the WC and SC PT graphs; c) the statement of Kotani, that the static random phase approximation reproduces well the DFT-LDA calculation with the standard choice of $`v_{xc}`$, we come to the conclusion that, within the approximation chosen, the original LDA functional form for $`v_{xc}`$ can be used. However, the charge density is calculated from renormalized fermions which differ from the standard fermions by the renormalizing factors $`\sqrt{P}`$. Hence the charge density,$`\rho _P(r)`$, also depends on $`\sqrt{P}`$. The mixing and hopping in the effective Kohn-Sham equation should be renormalized by the same $`P`$ as described above and a part of the renormalized $`f`$-electrons should be described as localized (lower transitions ) and the rest (coming from the upper transitions ) as delocalized. The separation is caused by self-interaction correction automatically generated in PT. Finally, the equation of self-consistency for parameter $`P`$ (Eqn 9) should be added.

We have hence calculated the total energy using the regular LDA total energy expression, using the electron density $`\rho _P(r)`$, i.e. $`E_{LDA}(\rho _P(r)`$, where the parameter $`P`$ which comes from electron correlations, is found self-consistently from additional equation and regulates distribution of spectral weight between lower and upper parts of energy.

In Fig. 1A we show the equilibrium volumes of some selected lanthanides, from experiment, uncorrected Kohn-Sham calculations (LMTO-ASA) and from our corrected theory. Note that the many-body corrections are presented together with shifted data due to an incomplete basis set used in the LMTO-ASA calculations. This shift was calculated by comparing (LDA) volumes obtained from a full-potential-(FP)-LMTO, multi-basis set calculation with (LDA) volumes from LMTO-ASA calculations, $`\mathrm{\Delta }R_{sws}=R_0^{FP}R_0^{LMTOASA}`$. An overall good agreement is obtained, and for the heavier lanthanides the renormalization correction becomes vanishingly small (as seen for Tm). The only parameter of our model, the position of the UHB with respect to the Fermi level, was taken from bremsstrahlung isochromat spectra (BIS) . In principle this value can be found from a super cell calculation changing the $`f`$-occupancy at one site. For cerium, $`\mathrm{\Delta }_{n,n+1}`$ is $``$3.5 eV above the Fermi level, and shifting this to 3.5$`\pm `$ 1 eV, only affected the equilibrium volume with $`\pm `$0.5%. For samarium the position of the UHB is $``$0.6 eV, and changes in the position with $`\pm `$0.1 eV, changed the volume with $`\pm `$0.8%. For the other elements changes were of the same order. Therefore, we conclude that the results are rather insensitive with respect to the position of the UHB. Fig. 1B shows the bulk modulus. Note that this value is quite sensitive to the fitting procedure (pressure-volume fit) and to the structure, thus, it should be considered only as an approximate value.

The physical picture obtained is consistent with previous theoretical work on Fermi surfaces  where a better agreement between theory and experiment was obtained when including, ad-hoc, the 4f states as valence electrons even though part of these 4f states were treated as occupied core states. It is also consistent with results obtained from self interaction corrected density functionals on Pr metal. Unlike the method of orbital polarization the present theory is a true many-body correction, and thus presents a different physical approach. The decisive parameter in our case (unlike to the LDA+U method) is not the absolute value of the Hubbard $`U`$ but the position of the transition $`E(\mathrm{\Gamma }_{n+1})E(\mathrm{\Gamma }_n)`$ with respect to Fermi energy. We have taken this parameter from experimental data  and this led us to an improvement in the calculated cohesive properties. First loop corrections, affecting the position of the levels, are small in our case . Fluctuations coming from charge excitations are expected to alter the situation at non-zero temperature, however, this remains an open question.

We have shown that many-body corrections to the standard model improves the description of the cohesive properties of the lanthanides, and explains the observed discrepancy for the equilibrium volumes for the light lanthanides. The only input to the present theory is the position of the upper Hubbard transition, $`\mathrm{\Delta }_{\mathrm{\Gamma }+1,\mathrm{\Gamma }}`$, which was taken from bremsstrahlung isochromat spectra (BIS) data. We found that a moderate change in energy of this transition gives only a minor change in the equilibrium volume. The calculated volumes are in good agreement with the experimental ones, which means that the excited states contributes to the ground state properties. The presented theory is designed to deal with strongly correlated electron systems and is integrated in an electronic structure method, and made fully self-consistent. It describes an additional contribution to the chemical bonding from a fraction $`f`$-electrons hybridized with the conduction electrons. Compared to for instance the LDA+U, orbital polarization or SIC methods the present theory involves a decrease of mixing due to strong correlations, causing the $`f`$-bandwidth to be reduced.

We are thankful for financial support from the Swedish Natural Science Research Council. The critical reading of our manuscript by Prof. Helmut Eschrig is greatly appreciated.