text:
# Area-preserving dynamics of a long slender finger by curvature: a test case for the globally conserved phase ordering

## I INTRODUCTION

This work is motivated by the recent developments in phase-ordering theory. Phase ordering is the growth of order from disorder via domain growth and coarsening. As phase ordering systems are strongly nonlinear and highly disordered, their theoretical description remains challenging . To get insight, it is often useful to consider the coarsening dynamics of simple test objects. One such object is a spherical droplet of the “minority” phase. It serves as the building block of simplified phase-ordering theories for locally conserved (bulk-diffusion-controlled)  and globally conserved (interface-controlled)  systems.

Being too a simple object, a spherical droplet may not give sufficient insight. For example, it’s dynamics (shape-preserving shrinking or expansion) look the same in both locally, and globally conserved systems. On the other hand, there are important differences in the phase-ordering behavior of locally and globally conserved systems. In both cases, an ensemble of droplets exhibits Ostwald ripening, and the corresponding mean-field theories, due to Lifshitz and Slyozov  and Wagner , respectively, predict dynamic scaling behavior of the droplet size distribution (with different dynamic exponents). However, for finite volume fractions, the difference in the type of conservation law leads to a different role of correlations. Correlations between neighboring droplets are much more important in the locally-conserved systems , than in the globally-conserved ones . Even bigger differences have been found in the phase ordering dynamics of locally and globally conserved systems with long-range correlations. Coarsening of fractal clusters shows dynamic scale invariance and “normal” scaling in the case of global conservation , and breakdown of scale invariance and anomalous scaling in locally conserved systems .

It is natural to attribute these differences to a basic difference in the character of transport: in globally conserved systems the transport is uninhibited by Laplacian screening effects, typical for locally conserved diffusion-controlled systems. To study this basic difference in a simple setting, the dynamics of a single long and slender finger is often considered. It has been found recently  that, in the locally-conserved, diffusion-controlled system, the finger acquires a dumbbell shape and shows non-trivial dynamic scalings for the finger length and the “ball” size, while the initial finger width remains almost constant until a late stage of the dynamics. Furthermore, a finger-shaped domain (a “needle”) served as a test object in still another conserved coarsening system: the one controlled by edge diffusion . Finger-like objects appear naturally in two-dimensional simulations of dewetting , etc. In this work we investigate the globally conserved dynamics of a long slender finger. An additional motivation for studying the finger dynamics is a recent observation that, at a late stage of coarsening of fractal clusters, the cluster morphology shows long branches, or fingers, in both locally , and globally  conserved systems.

Here is an outline of the rest of the paper. In Section 2 we briefly review the Ginzburg-Landau equation with a global conservation law (a phase field model for globally-conserved phase ordering) and its general sharp interface formulation in two dimensions. Under a certain condition (that will be elucidated) this formulation is reduced to the area-preserving dynamics by curvature. In Section 3 we formulate a moving boundary problem for the finger dynamics and develop a perturbation theory which enables us to obtain an intermediate asymptotic solution for the finger shape. In Section 4 we return to the globally-conserved Ginzburg-Landau equation and investigate the finger dynamics numerically. Section 5 includes a brief discussion.

## II globally conserved phase ordering: a phase field model and its sharp-interface limit

Globally-conserved phase ordering can be described by a simple phase field model . This model represents the dynamics in terms of a simple gradient descent

$$\frac{u}{t}=\frac{\delta F}{\delta u}$$
(1)
with the Ginzburg-Landau free energy functional

$$F[u]=\left[(1/2)(u)^2+V(u)+Hu\right]𝑑𝐫$$
(2)
and a double-well potential $`V(u)=(1/4)(1u^2)^2`$. The effective “magnetic field” $`H=H(t)`$ varies in time so as to impose global conservation law

$$<u(𝐫,t)>=\frac{1}{L_xL_y}u(𝐫,t)𝑑𝐫=\text{const}.$$
(3)
Here $`u(𝐫,t)`$ is the order parameter field, $`L_x1`$ and $`L_y1`$ are the linear dimensions of the system (a two-dimensional rectangular box), and the integration is carried out over the entire box. Eqs. (1)-(3) yield the nonlocal Ginzburg-Landau equation

$$\frac{u}{t}=^2u+uu^3<uu^3>.$$
(4)
Either no-flux, or periodic boundary conditions are assumed.

For a wide class of initial conditions this coarsening system segregates, at late times, into large domains of ”phase 1” and ”phase 2” separated by thin domain walls (which width is of order unity) . Correspondingly, a sharp interface theory can be developed for these late times . At these times the magnetic field $`H(t)`$ is already small, $`H(t)1`$, and slowly varies with time. The phase field in the phases $`1`$ and $`2`$ is almost uniform and rapidly adjusts to the current value of $`H(t)`$, so that $`u=1H(t)/2`$ and $`1H(t)/2`$, respectively. The normal velocity of the interface is given by

$$v_n(s,t)=\kappa (s,t)\frac{3}{\sqrt{2}}H(t),$$
(5)
where $`\kappa `$ is the local curvature and $`s`$ is the coordinate along the interface. The positive sign of $`v_n`$ corresponds to the interface moving towards phase 1, while $`\kappa `$ is positive when the interface is convex toward phase 2.

The dynamics of $`H(t)`$ is determined by

$$\dot{H}(t)=\frac{4\mathrm{\Lambda }(t)}{L_xL_y}\left(\overline{\kappa (s,t)}\frac{3}{\sqrt{2}}H(t)\right),$$
(6)
where

$$\overline{\kappa (s,t)}=\frac{1}{\mathrm{\Lambda }(t)}\kappa (s,t)𝑑s,$$
(7)
and $`\mathrm{\Lambda }(t)`$ is the total interface length. Equations (5) and (6) make a complete set and provide a general sharp-interface formulation to this problem. In some cases this formulation can be further simplified . Let us compute the rate of change of the domain area, $`\dot{A}`$:

$$\dot{A}(t)=v_n(s,t)𝑑s=\mathrm{\Lambda }(t)\left[\overline{\kappa (s,t)}\frac{3}{\sqrt{2}}H(t)\right].$$
(8)
If the two terms in the right hand side of Eq. (8) balance each other,

$$H(t)\frac{\sqrt{2}}{3}\overline{\kappa (s,t)},$$
(9)
the domain area remains approximately constant. Then, using Eqs. (5) and (9), we obtain

$$v_n(s,t)=\kappa (s,t)\overline{\kappa (s,t)}.$$
(10)
Model (10) is known as area-preserving flow by curvature . Dynamics (10) reduce the interface length of the system , leaving the areas of each of the phases constant. The only stable steady state of Eq. (10) (not attached to the boundaries of the system) is a single perfectly circular domain . Under what condition is model (10) a good approximation to the more general sharp-interface model (5) and (6)? Consider Eq. (6) which includes the same combination $`\overline{\kappa }(3/\sqrt{2})H(t)`$. Let us treat the term in the right hand side of Eq. (6) perturbatively, putting $`H(t)=H^{(0)}(t)+h(t)`$, where $`H^{(0)}(t)=(\sqrt{2}/3)\overline{\kappa }`$ is the leading term, and $`h(t)`$ is a subleading term. Keeping terms up to order $`h(t)`$ we obtain:

$$h(t)=\frac{\sqrt{2}L_xL_y}{12\mathrm{\Lambda }}\dot{H}^{(0)}(t)=\frac{L_xL_y}{18\mathrm{\Lambda }}\dot{\overline{\kappa }}$$
(11)
If there is only one (single-connected) domain, $`\overline{\kappa }=2\pi /\mathrm{\Lambda }`$. Therefore,

$$h(t)=\frac{\pi L_xL_y}{9}\frac{\dot{\mathrm{\Lambda }}}{\mathrm{\Lambda }^3}.$$
(12)
Requiring $`|h(t)||H^{(0)}(t)|`$, we obtain the “area-preservation” criterion

$$\frac{L_xL_y|\dot{\mathrm{\Lambda }}|}{\mathrm{\Lambda }^2}1.$$
(13)
In Section 3 we will investigate the area-preserving dynamics of a long slender finger, and verify criterion (13) a posteriori.

## III Area-preserving finger dynamics: a perturbation theory

Let the initial condition for dynamics (10) represent a long and narrow rectangular bar with length $`2a_0`$ and width $`2\mathrm{\Delta }_0`$. We will see that, with time, the bar evolves into a finger-shaped object. We place the origin of a Cartesian coordinate system in the finger’s center and denote by $`x`$ the coordinate along the finger and by $`y=y(x,t)`$ the instantaneous location of the finger boundary. Because of the symmetry with respect to the $`x`$-axis, we will be interested only in the upper boundary dynamics: $`y(x,t)0`$. In terms of $`y(x,t)`$ the local curvature is given by :

$$\kappa (s,t)=\frac{y_{xx}}{(1+y_x^2)^{3/2}},$$

where $`y_xy(x,t)/x`$. The average curvature $`\overline{\kappa (s,t)}`$ is

$$\overline{\kappa (s,t)}=\frac{2\pi }{\mathrm{\Lambda }}=\frac{\pi }{2_0^{a(t)}(1+y_x^2)^{1/2}𝑑x},$$

where $`a(t)`$ is the time-dependent position of the finger tip. Since

$$\frac{y}{t}=(1+y_x^2)^{1/2}v_n,$$

dynamics (10) imply the following evolution equation for $`y(x,t)`$:

$$\frac{y}{t}=\frac{y_{xx}}{1+y_x^2}+\frac{\pi (1+y_x^2)^{1/2}}{2_0^{a(t)}(1+y_x^2)^{1/2}𝑑x}.$$
(14)
$`a(t)`$ is defined by the boundary conditions $`y(x=a,t)=0`$ and $`y_x(x=a,t)=\mathrm{}`$. An additional boundary condition follows from the symmetry with respect to the $`y`$-axis: $`y_x(x=0,t)=0`$. Equation (14), a nonlinear integro-differential equation in partial derivatives, describes a free-boundary problem. It is easy to check that the dynamics (14) preserve the finger area $`A`$.

It can be proved  that any convex domain evolving by Eq. (10) \[and therefore by Eq. (14)\] finally becomes a circle with the same area $`A`$. We are interested in the finger dynamics at times much smaller than those needed for approaching this equilibrium (but still large enough so that details of the initial shape are forgotten). The corresponding strong double inequality will be presented at the end of this Section.

As long as the ratio of the finger width to its length is small, we can use it as the small parameter of our theory. Indeed, near the finger tip the first (local) term in the right hand side of Eq. (14) is of the order of the inverse finger width, while the second (nonlocal) term is of the order of the inverse finger length. Therefore, the nonlocal term can be treated perturbatively. The unperturbed equation,

$$\frac{y}{t}=\frac{y_{xx}}{1+y_x^2},$$
(15)
was obtained more than 40 years ago by Mullins in the context of motion of grain boundaries . It corresponds to a non-conserved motion by curvature: $`v_n(s,t)=\kappa `$. We notice immediately that the perturbed problem, Eq. (14), possesses an integral of motion (the finger area) which is absent in the corresponding unperturbed problem. This situation is unusual: more often than not perturbations destroy integrals of motion pertinent to the unperturbed system. This unusual property will be exploited in the following, as we will require area conservation in each order of the perturbation expansion.

Mullins  obtained a one-parameter family of traveling-wave solutions of Eq. (15) describing a half-infinite constant-width finger retreating by its tip’s curvature:

$$y=Y(x,t)=\frac{1}{c}\mathrm{arccos}\left[\mathrm{exp}\left(c\xi \right)\right],$$
(16)
where $`\xi =a(t)x`$, $`a(t)=a_0ct`$ is the time-dependent position of the finger tip, and the constant speed of retreat $`c`$ is the parameter of this family of solutions. The finger’s half-width $`\mathrm{\Delta }=\pi /(2c)`$ is constant. We will show that this solution (with two important modifications) can serve as the leading (or zero) order approximation of our perturbation theory for the area-preserving finger dynamics. First, we introduce a (slow) time-dependence in the parameter $`c`$. Second, we “stick together” two identical (very long but finite) half-fingers. The resulting ansatz is the following:

$$Y(x,t)=\frac{1}{c(t)}\mathrm{arccos}\left[\mathrm{exp}\left(c(t)\xi \right)\right],$$
(17)
where $`\xi =a(t)|x|`$, and $`a(t)>0`$ is yet unknown. The finger half-width $`\mathrm{\Delta }(t)=\pi /2c(t)`$ is now time-dependent; its growth with time makes up, in the area conservation, for the finger shortening. As a function of $`x`$, ansatz (17) is continuous everywhere. Its $`x`$-derivative is not continuous at $`x=0`$. However, as far as $`ϵ(t)=\mathrm{\Delta }(t)/a(t)`$ is very small: $`ϵ(t)1`$, the $`x`$-derivatives at $`x=0`$ from left and right are of order of $`\mathrm{exp}(1/ϵ)`$, that is exponentially small. This and other exponentially small effects will be neglected throughout the paper. Because of the symmetry with respect to the $`y`$-axis, we will consider only the region $`0xa(t)`$. The perturbation expansion for $`y(x,t)`$ is

$$y(x,t)=Y(x,t)+\delta y(x,t)+\mathrm{},$$
(18)
where $`\delta y(x,t)`$ is the subleading (or first order) term. Inserting Eq. (17) into Eq. (14) and keeping only the zero order terms, we obtain

$$\dot{a}=c=\pi /(2\mathrm{\Delta }).$$
(19)
The evolution equation for $`\mathrm{\Delta }(t)`$ follows from the requirement of area conservation in the zero order. Calculating the area under the $`Y(x,t)`$ profile, we obtain:

$$_0^{a(t)}Y(x,t)𝑑x=a\mathrm{\Delta }\frac{2\mathrm{\Delta }^2\mathrm{ln}2}{\pi }+O(\mathrm{exp}(1/ϵ)).$$
(20)
One can see that the finger area in the zero order is simply $`A=4a\mathrm{\Delta }`$. Demanding $`\dot{A}=0`$ and using Eq. (19), we obtain:

$$\dot{\mathrm{\Delta }}=\pi /(2a).$$
(21)
We will see later that the same equation follows from the analysis of the first order correction. The solution of the zero order equations (19) and (21) is

$$a^{(0)}(t)=a_0\mathrm{exp}(\pi t/2a_0\mathrm{\Delta }_0),$$
(22)
and

$$\mathrm{\Delta }^{(0)}(t)=\mathrm{\Delta }_0\mathrm{exp}(\pi t/2a_0\mathrm{\Delta }_0).$$
(23)
Therefore, in the zero order, the finger shortens (and its width grows) exponentially with time. The characteristic growth time is of order of the finger area. The zero order solution yields a criterion for the validity of the perturbation theory. The aspect ratio of the evolving finger should be very small which leads to

$$ta_0\mathrm{\Delta }_0\mathrm{ln}\frac{a_0}{\mathrm{\Delta }_0},$$
(24)
This condition sets an upper limit for the times for which our perturbation expansion is valid.

We now turn to calculating the small correction $`\delta y(x,t)`$. Since the nonlocal term of Eq. (14) is already of order $`ϵ`$, we evaluate it on the zero-order solution $`Y(x,t)`$ which yields $`\pi (2a)^1(1+Y_x^2)^{1/2}`$. Introduce a new variable $`u`$:

$$u(x,t)=\mathrm{exp}(c(t)\xi ).$$
(25)
Substituting Eq. (18) into Eq. (14) and linearizing, we get a linear partial differential equation for $`\delta y(x,t)`$:

$`\delta y_t`$ $`=`$ $`(1u^2)\delta y_{xx}2cu^2\delta y_x`$ (27)
$`{\displaystyle \frac{(c+\dot{a})u}{(1u^2)^{1/2}}}+{\displaystyle \frac{\dot{c}}{c^2}}q(u)+{\displaystyle \frac{\pi }{2a(1u^2)^{1/2}}},`$
where

$$q(u)=\frac{u\mathrm{ln}u}{(1u^2)^{1/2}}+\mathrm{arccos}u.$$
(28)
It is convenient to go over to new independent variables $`\xi `$ and $`t`$. We define $`\delta y(x,t)=f(\xi ,t)`$ and obtain $`\delta y_t=f_t+\dot{a}f_\xi `$. As the $`t`$-dependence in the new variables is slow, the term $`f_t`$ can be neglected in this order of the perturbation scheme. Therefore, we are left with an ordinary differential equation (the slow time enters as a parameter):

$`(1u^2)f_{\xi \xi }+c(2u^2+1)f_\xi `$ $`=`$ $`{\displaystyle \frac{(c+\dot{a})u}{(1u^2)^{1/2}}}{\displaystyle \frac{\dot{c}}{c^2}}q(u)`$ (30)
$`{\displaystyle \frac{\pi }{2a(1u^2)^{1/2}}}.`$
Going over from $`\xi `$ to $`u`$ and defining $`f(\xi )=g(u)`$ we obtain

$`(1u^2)g_{uu}3ug_u`$ $`=`$ $`{\displaystyle \frac{(c+\dot{a})}{c^2}}{\displaystyle \frac{1}{u(1u^2)^{1/2}}}{\displaystyle \frac{\dot{c}}{c^4}}{\displaystyle \frac{q(u)}{u^2}}`$ (32)
$`{\displaystyle \frac{\pi }{2ac^2u^2(1u^2)^{1/2}}}.`$
Define $`\mathrm{\Phi }(u)=g_u(u)`$. The general solution of the homogeneous equation for $`\mathrm{\Phi }(u)`$ is $`C\mathrm{\Phi }_0(u)`$, where

$$\mathrm{\Phi }_0(u)=(1u^2)^{3/2},$$
(33)
and $`C`$ is an arbitrary constant. Therefore, we look for the general solution of Eq. (32) in the form $`\mathrm{\Phi }(u)=C(u)\mathrm{\Phi }_0(u)`$, where $`C(u)`$ is an unknown function. Substituting this into Eq. (32) and integrating, we obtain

$$C(u)=\frac{(c+\dot{a})}{c^2}\mathrm{ln}u\frac{\dot{c}}{c^4}s(u)+\frac{\pi }{2ac^2u}+\alpha (t),$$
(34)
where

$`s(u)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}^2u+{\displaystyle \frac{1}{2}}\mathrm{arccos}^2u\mathrm{ln}u`$ (36)
$`{\displaystyle \frac{(1u^2)^{1/2}\mathrm{arccos}u}{u}},`$
and $`\alpha (t)`$ is an arbitrary function of time. Going back to $`g(u)`$, we can write the first order correction to the finger shape as

$$g(u)=_0^uC(u)\mathrm{\Phi }_0(u)𝑑u+\beta (t),$$
(37)
where $`\beta (t)`$ is another arbitrary function of time. The boundary conditions for $`g(u)`$ are $`g(0)=0`$ and $`g(1)=0`$. The first boundary condition gives $`\beta (t)=0`$, while the second boundary condition should be used to find an equation for $`\dot{a}`$ in the first order. Before doing this, an analysis of the integrals in Eq. (37) should be made. To prevent divergence of $`g(u)`$ at $`u=1`$ we must require $`\alpha (t)=\pi /(2ca)`$. Similarly, to prevent divergence of $`g(u)`$ at $`u=0`$ we must require $`\dot{c}/c^2=1/a`$. The latter equation is equivalent to Eq. (21) obtained in the leading order from the requirement of area conservation.

The equation for $`\dot{a}`$ in the first order of the perturbation theory is obtained from the boundary condition $`g(1)=0`$ which requires calculating the integral in Eq. (37) from $`u=0`$ to $`u=1`$. This calculation yields

$$\dot{a}=c+\frac{3\mathrm{ln}2}{a}=\frac{\pi }{2\mathrm{\Delta }}+\frac{3\mathrm{ln}2}{a}.$$
(38)
It is now possible to obtain an explicit expression for $`g(u)`$. An integration from $`0`$ to $`u`$ in Eq. (37) gives:

$`g(u)`$ $`=`$ $`{\displaystyle \frac{1}{c^2a}}[w_1(u)\mathrm{arcsin}u\mathrm{ln}\left({\displaystyle \frac{8u}{e}}\right)+`$ (40)
$`{\displaystyle \frac{\pi }{2}}\mathrm{ln}\left({\displaystyle \frac{2}{e\left[1+(1u^2)^{1/2}\right]}}\right)+w_2(u)],`$
where

$$w_1(u)=\frac{2u\mathrm{ln}u\mathrm{ln}(8u^{1/2}/e)+u\mathrm{arccos}^2u+\pi (1u)}{2(1u^2)^{1/2}},$$
(41)
and

$$w_2(u)=2_0^u\frac{\mathrm{arcsin}u}{u}𝑑u,$$
(42)
In the absence of a conservation law, we would have to go to the second order of the perturbation scheme to obtain the equation for $`\dot{c}`$ in the first order. In our area-preserving perturbation scheme this equation can be obtained already in the first order. Notice that the first-order correction obeys the following conservation law:

$$\frac{d}{dt}_0^{a(t)}\delta y(x,t)𝑑x=\dot{a}[f(\xi =a,t)f(\xi =0,t)]=0,$$
(43)
where the last equality results from the boundary conditions. This relation combined with the requirement that the finger area is preserved leads to

$$\frac{d}{dt}_0^{a(t)}Y(x,t)𝑑x=0.$$
(44)
It should be stressed that functions $`a(t)`$ and $`\mathrm{\Delta }(t)`$ entering Eq. (44) now include first order corrections. Now, using Eqs. (20), (38) and (44) and keeping terms only up to the first order in $`ϵ`$, we obtain

$$\dot{\mathrm{\Delta }}=\frac{\pi }{2a}\frac{\mathrm{\Delta }\mathrm{ln}2}{a^2}.$$
(45)
Equations (38) and (45) are the equations for $`\dot{a}`$ and for $`\dot{\mathrm{\Delta }}`$ in the first order of the perturbation scheme. This set of equations is autonomous and therefore integrable. In the framework of the perturbation theory it is more consistent to solve the equations perturbatively by linearization. To do this, we exploit the conservation law (44) which, together with Eq. (20), allows us to express $`\mathrm{\Delta }(t)`$ as a function of $`a(t)`$ or vice versa. The expression obtained is then linearized with respect to $`\mathrm{\Delta }/a`$ and inserted into Eq. (38) (or (45)). This yields the following equations:

$$\dot{a}=\frac{\pi a}{2a_0\mathrm{\Delta }_0}+\frac{4\mathrm{ln}2}{a}\frac{a\mathrm{ln}2}{a_0^2},$$
(46)
$$\dot{\mathrm{\Delta }}=\frac{\pi \mathrm{\Delta }}{2a_0\mathrm{\Delta }_0}+\frac{\mathrm{\Delta }\mathrm{ln}2}{a_0^2}\frac{\mathrm{\Delta }\mathrm{ln}2}{a^2}\frac{\mathrm{\Delta }^3\mathrm{ln}2}{a_0^2\mathrm{\Delta }_0^2}.$$
(47)
We are looking for solutions in the form of $`a(t)=a^{(0)}(t)+a^{(1)}(t)`$ and $`\mathrm{\Delta }(t)=\mathrm{\Delta }^{(0)}(t)+\mathrm{\Delta }^{(1)}(t)`$, where $`a^{(1)}(t)`$ and $`\mathrm{\Delta }^{(1)}(t)`$ are first order corrections. Solving the resulting linear equations, we obtain:

$`a^{(1)}(t)`$ $`=`$ $`{\displaystyle \frac{8\mathrm{\Delta }_0\mathrm{ln}2}{\pi }}\mathrm{sinh}\left({\displaystyle \frac{\pi t}{2a_0\mathrm{\Delta }_0}}\right)`$ (49)
$`{\displaystyle \frac{t\mathrm{ln}2}{a_0}}\mathrm{exp}\left({\displaystyle \frac{\pi t}{2a_0\mathrm{\Delta }_0}}\right),`$
$`\mathrm{\Delta }^{(1)}(t)`$ $`=`$ $`{\displaystyle \frac{4\mathrm{\Delta }_0^2\mathrm{ln}2}{\pi a_0}}\mathrm{exp}\left({\displaystyle \frac{\pi t}{a_0\mathrm{\Delta }_0}}\right)\mathrm{sinh}\left({\displaystyle \frac{\pi t}{2a_0\mathrm{\Delta }_0}}\right)+`$ (51)
$`{\displaystyle \frac{\mathrm{\Delta }_0t\mathrm{ln}2}{a_0^2}}\mathrm{exp}\left({\displaystyle \frac{\pi t}{2a_0\mathrm{\Delta }_0}}\right).`$
We checked that when solving Eqs.(38) and (45) numerically, the results differ from the approximate analytical results (49) and (51) by less than $`1`$%. Finally, requiring $`a^{(0)}(t)a^{(1)}(t)`$ and $`\mathrm{\Delta }^{(0)}(t)\mathrm{\Delta }^{(1)}(t)`$ and using (22), (23), (49), and (51), one can see that the linearization procedure is valid for times obeying the same strong inequality (24). Therefore, the linearization procedure is consistent with the perturbation theory.

The validity of our perturbative solution is limited by not too long times, see inequality (24). On the other hand, one should wait until irrelevant details of the initial bar shape are forgotten, and the solution approaches $`Y(x,t)`$. Typically, it takes time needed for the bar tip to pass the distance of order of the characteristic tip size $`1/c\mathrm{\Delta }`$. The tip speed is of order $`1/\mathrm{\Delta }`$, so the “waiting time” of the theory can be estimated as $`\mathrm{\Delta }_0^2`$. Our perturbative solution represents, therefore, an intermediate asymptotics, valid for

$$\mathrm{\Delta }_0^2ta_0\mathrm{\Delta }_0\mathrm{ln}\frac{a_0}{\mathrm{\Delta }_0}.$$
(52)
Finally, we should check the area-preservation criterion (13). In the leading order we have: $`\mathrm{\Lambda }(t)a^{(0)}(t)`$ and $`\dot{\mathrm{\Lambda }}(t)\dot{a}^{(0)}(t)1/\mathrm{\Delta }^{(0)}(t)`$. Using (22) and (23) in (13) we obtain:

$$ta_0\mathrm{\Delta }_0\mathrm{ln}(\stackrel{~}{f}a_0),$$
(53)
where $`\stackrel{~}{f}=(a_0\mathrm{\Delta }_0)/(L_xL_y)`$ is the area fraction occupied by the finger.

## IV numerical solution

In order to test the predictions of the perturbation theory we performed numerical simulations with the nonlocal Ginzburg-Landau equation (4) (with no-flux boundary conditions), taking a long and narrow rectangular bar as the initial condition. We used an explicit Euler integration scheme to advance the solution in time, and second order central differences to discretize the Laplace operator. A mesh size of $`\mathrm{\Delta }x=\mathrm{\Delta }y=0.5`$ was found sufficient for an accurate resolution of the interface. A time step of $`\mathrm{\Delta }t=0.05`$ was required for numerical stability. The bar width $`2\mathrm{\Delta }_0=50`$ was chosen to guarantee area conservation with a good accuracy \[see Eq. (53)\]. The simulation was carried out up to a time $`t_f=910^4`$ which was long enough to distinguish between an exponential and linear dynamic behavior for $`a(t)`$ and $`\mathrm{\Delta }(t)`$. The bar length $`2a_0=810^3`$ was chosen, so that the parameter $`ϵ(t)=\mathrm{\Delta }(t)/a(t)`$ was sufficiently small even at the end of simulation (when it reached $`0.1`$). The bar was placed in the center of a rectangular box $`8020\times 220`$. Because of symmetry with respect to the $`x`$ and $`y`$ axes only the quadrant $`x>0,y>0`$ was actually simulated. We checked that the finger area was conserved with an accuracy better than $`1`$%.

The time dependence of the finger half-length $`a(t)`$ obtained from the numerical simulations is shown in Fig. 1. The same graph also shows the theoretical values of $`a(t)`$ in the zeroth order \[$`a^{(0)}(t)`$, Eq. (22)\], and first order \[$`a^{(0)}(t)+a^{(1)}(t)`$, Eq. (49)\] of the perturbation theory. A good agreement between the numerical and analytical results is obtained already in the zeroth order. This agreement is improved further by the first order correction. For example, at the final time of the simulation $`t_f=910^4`$, there is a $`9`$% difference between the numerical value of $`a(t)`$ and $`a^{(0)}(t)`$, and only a $`1.4`$% difference between the numerical value of $`a(t)`$ and $`a^{(0)}(t)+a^{(1)}(t)`$.

The numerically found half-width of the finger $`\mathrm{\Delta }(t)`$ is shown in Fig. 2. Its comparison with the zero-order values $`\mathrm{\Delta }^{(0)}(t)`$ \[Eq. (23)\] and with the first-order values $`\mathrm{\Delta }^{(0)}(t)+\mathrm{\Delta }^{(1)}(t)`$ \[Eq. (51)\] is also shown. Again, there is a good agreement already in the zeroth order, and this agreement is further improved by the first order correction. At the final time of the simulation, the numerical value of $`\mathrm{\Delta }(t)`$ differs by $`6`$% from $`\mathrm{\Delta }^{(0)}(t)`$, and by less than $`2`$% from $`\mathrm{\Delta }^{(0)}(t)+\mathrm{\Delta }^{(1)}(t)`$.

Fig. 3 shows the finger shape $`y(x,t)`$ found numerically, at the final time of the simulation $`t_f=910^4`$ (notice the different scales in the $`x`$ and $`y`$ axes). Also shown are the following zero- and first order theoretical results evaluated at $`t=t_f`$:

$$Y^{(0)}(x,t)=\frac{1}{c^{(0)}(t)}\mathrm{arccos}[\mathrm{exp}(c^{(0)}(t)\xi ^{(0)})],$$
(54)
where $`\xi ^{(0)}=a^{(0)}(t)x`$, and

$$y_s(x,t)=\frac{1}{c_s(t)}\mathrm{arccos}[\mathrm{exp}(c_s(t)\xi _s)]+\delta y(x,t),$$
(55)
where $`\xi _s=a_s(t)x`$, $`a_s=a^{(0)}+a^{(1)}`$, and $`c_s=c^{(0)}+c^{(1)}`$. In addition, shown is the function $`Y_{num}(x,t)`$ defined by:

$$Y_{num}(x,t)=\frac{1}{c_{num}(t)}\mathrm{arccos}[\mathrm{exp}(c_{num}(t)\xi _{num})],$$
(56)
where $`\xi _{num}=a_{num}(t)x`$, and $`a_{num}(t)`$ and $`c_{num}(t)`$ are the numerical values of $`a(t)`$ and $`c(t)`$, respectively. The function $`Y_{num}(x,t)`$ is also evaluated at $`t=t_f`$. The comparison with the zero-order prediction $`Y^{(0)}(x,t_f)`$ shows a good agreement even at this relatively late time, when $`\mathrm{\Delta }/a0.1`$. The agreement is further improved by the first-order theory. The comparison with $`Y_{num}(x,t_f)`$ shows an excellent agreement, implying that the subleading term $`\delta y(x,t)`$ is very small compared to the leading term $`Y(x,t)`$. Fig. 4 which shows the theoretical result for $`\delta y(x,t)`$ at $`t=t_f`$, calculated by using Eqs. (40)-(42), confirms the later conclusion. We should add that a very good agreement was found between the theoretical $`\delta y(x,t)`$ and the numerical $`\delta y(x,t)`$, which is obtained by subtracting $`Y_{num}(x,t)`$ from the numerical result for $`y(x,t)`$.

## V CONCLUSIONS

We investigated the dynamics of a long slender finger-shaped domain undergoing globally-conserved (interface-controlled) coarsening. We worked in the parameter region where these dynamics can be reduced to area-preserving motion by curvature. The nonlinear moving boundary problem for the finger shape was solved by using an asymptotic perturbation theory which employed, as the zero-order solution, a straightforward modification of the Mullins solution for a finger retreating by the curvature of its tip. Both the leading, and the subleading terms of the solution were calculated analytically and verified by a numerical solution of the corresponding phase-field equation: the Ginzburg-Landau equation with a global conservation law.

The coarsening dynamics of the finger look quite differently from those observed in the case of a locally conserved (diffusion-controlled) system . The finger preserves its simple shape, while its characteristic length and width are changing on the same time scale. Therefore, the global character of transport (uninhibited by Laplacian screening effects characteristic of the locally conserved system) manifests itself already in the simple setting of finger dynamics.

## ACKNOWLEDGEMENTS

This work was supported in part by the Israel Science Foundation, administered by the Israel Academy of Sciences and Humanities.


format: text

text:
# Vertical Density Matrix Algorithm: A Higher-Dimensional Numerical Renormalization Scheme based on the Tensor Product State Ansatz

## I Introduction

Since the density matrix renormalization group (DMRG) method was invented by S.R.White,  the method has been applied to various problems in one-dimensional (1D) quantum systems and two-dimensional (2D) classical systems.  Such a great success of the DMRG has been stimulating us to extend the algorithm to the one which can handle higher-dimensional systems, principally 2D quantum systems and 3D classical systems.

We should recall that, in the DMRG, the matrix-product structure of the wavefunction of the target states (usually the ground state or the maximum-eigenvalue eigenstate) is essential.  From this point of view, the tensor product state (TPS) which is a natural higher-dimensional generalization of the matrix product state, should play a key role in the “higher-dimensional DMRG”. A simple but non-trivial example of the TPS is the ground state of the 2D valence-bond-solid (VBS) type quantum spin systems, where the wavefunction is expressed as a product of local finite-dimensional tensors, with all the tensor indices being contracted.  As for 3D classical statistical systems, the maximum-eigenvalue eigenstate of the layer-to-layer transfer matrix of the 3D classical system can exactly be represented as the TPS, if we allow the tensor dimension to be infinite. We should note that we can reduce the calculation of the expectation value of the TPS to a statistical average in a lower-dimensional classical system; $`(D+1)`$-dimensional classical (or $`D`$-dimensional quantum) problem reduces to a $`D`$-dimensional classical statistical problem. In fact, the properties of the 2D VBS model have been studied in terms of a 2D vertex model associated with the TPS wavefunction.

When developing the TPS formulation, the most important step is the determination (optimization) of the local tensor in the TPS. In Refs. , direct variational formulations are employed for the optimization of the local tensors. As compared with these “direct” method, our novel algorithm given in the present paper for 3D classical system is more like the original DMRG. The local tensor is updated by the “block-spin basis transformation” along the vertical direction. Since this transformation is constructed in terms of the density matrix made along the “vertical direction” (transfer direction of the transfer matrix), we call this algorithm the vertical density matrix algorithm (VDMA). We apply the VDMA to the 3D Ising model and discuss its efficiency. We also report the application of the VDMA to the 2D transverse field Ising model, with help of the Suzuki-Trotter transformation .

This paper is organized as follows: In Sec.II, we briefly explain the VDMA for 3D classical spin systems, taking the 3D Ising model as an example. In Sec.III we show the numerical result for the 3D Ising model and the 2D transverse field Ising model. The last section is devoted to the conclusion.

## II Method

Let us consider the 3D Ising model on the simple cubic lattice of the size $`N\times N\times 2L`$ in X, Y, and Z directions. Suppose that $`L`$ and $`N`$ are sufficiently large, and the neighboring Ising spins $`\sigma `$ and $`\sigma ^{}`$ have ferromagnetic interaction $`J\sigma \sigma ^{}`$. Then the Boltzmann weight for the unit cube is written as

$`W\left(\begin{array}{cccc}\overline{\sigma }_{ij}& \overline{\sigma }_{i^{}j}& \overline{\sigma }_{i^{}j^{}}& \overline{\sigma }_{ij^{}}\\ \sigma _{ij}& \sigma _{i^{}j}& \sigma _{i^{}j^{}}& \sigma _{ij^{}}\end{array}\right)\mathrm{exp}[{\displaystyle \frac{K}{2}}(\sigma _{ij}\sigma _{i^{}j}+\sigma _{i^{}j}\sigma _{i^{}j^{}}+\sigma _{i^{}j^{}}\sigma _{ij^{}}+\sigma _{ij^{}}\sigma _{ij}+\overline{\sigma }_{ij}\overline{\sigma }_{i^{}j}+\overline{\sigma }_{i^{}j}\overline{\sigma }_{i^{}j^{}}`$ (3)
$`+\overline{\sigma }_{i^{}j^{}}\overline{\sigma }_{ij^{}}+\overline{\sigma }_{ij^{}}\overline{\sigma }_{ij}+\sigma _{ij}\overline{\sigma }_{ij}+\sigma _{i^{}j}\overline{\sigma }_{i^{}j}+\sigma _{i^{}j^{}}\overline{\sigma }_{i^{}j^{}}+\sigma _{ij^{}}\overline{\sigma }_{ij^{}})],`$ (4)
where $`i^{}=i+1`$, $`j^{}=j+1`$, and $`K=J/T`$. The locations of the spin variables are shown in Fig. 1.

For the book keeping, let us here introduce some notations for spin variables. Write the configuration of the four spins surrounding a plaquette in the $`XY`$ plane as

$$\{\sigma _{ij}\}=(\sigma _{ij}\sigma _{i^{}j}\sigma _{i^{}j^{}}\sigma _{ij^{}}),$$
(5)
where the position of the plaquette can be labeled by the index $`ij`$. Then the Boltzmann weight is simply written as $`W\left\{\begin{array}{c}\overline{\sigma }_{ij}\\ \sigma _{ij}\end{array}\right\}`$. Also for a spin layer in the $`XY`$ plane, we denote the configuration of the $`N\times N`$ spins as,

$$[\sigma ]\left[\begin{array}{ccc}\sigma _{11}& \mathrm{}& \sigma _{1N}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \sigma _{N1}& \mathrm{}& \sigma _{NN}\end{array}\right].$$
(6)
Using these notations, the transfer matrix $`T`$ from a layer $`[\sigma ]`$ to the next layer $`[\overline{\sigma }]`$ is written as:

$$T[\overline{\sigma }|\sigma ]=\underset{i+j=even}{}W\left\{\begin{array}{c}\overline{\sigma }_{ij}\\ \sigma _{ij}\end{array}\right\}.$$
(7)
In the product of Eq. (7), the spin variables are shared by the adjacent cubes in the diagonal direction and thus the Boltzmann weights make the checkerboard pattern in the $`XY`$-plane.(see Fig.1-(b))

Our goal is to evaluate the maximum-eigenvalue $`\lambda _{\mathrm{max}}`$ of the transfer matrix $`T`$ and the corresponding eigenvector $`|\psi _{\mathrm{max}}`$, using the TPS representation of the eigenvector. In order to analyze the TPS structure of $`|\psi _{\mathrm{max}}`$, let us consider the power method briefly, which is the simplest but powerful technique to calculate $`\lambda _{\mathrm{max}}`$ and $`|\psi _{\mathrm{max}}`$. Define the vector $`|\psi _L`$ with

$$|\psi _L=T^L|\psi _0,$$
(8)
where $`|\psi _0`$ is an ‘initial’ vector that is not orthogonal to $`|\psi _{\mathrm{max}}`$. Then the maximum-eigenvalue eigenvector $`|\psi _{\mathrm{max}}`$ is obtained as

$$|\psi _{\mathrm{max}}=\mathrm{const}\times \underset{L\mathrm{}}{lim}|\psi _L,$$
(9)
where the details of $`|\psi _0`$ are not important.

In Fig. 2, we show the graphical representation of $`|\psi _L`$, where we can see the structure of $`|\psi _L`$ more clearly. As is seen in this figure, the $`L`$ unit cubes are piling vertically up to the surface, and then the product of the vertically-arranged Boltzmann weights can be regarded as a local tensor; We define the local tensor $`A_L`$ at the $`ij`$-plaquette as,

$$A_L\left\{\begin{array}{c}\sigma _{ij}^L\\ \xi _{ij}^L\end{array}\right\}=W\left\{\begin{array}{c}\sigma _{ij}^L\\ \sigma _{ij}^{L1}\end{array}\right\}W\left\{\begin{array}{c}\sigma _{ij}^{L1}\\ \sigma _{ij}^{L2}\end{array}\right\}\mathrm{}W\left\{\begin{array}{c}\sigma _{ij}^1\\ \sigma _{ij}^0\end{array}\right\},$$
(10)
where $`\{\sigma _{ij}^L\}`$ is the spin configuration at the $`ij`$-plaquette on the “surface layer”, and the auxiliary variable $`\xi _{ij}^L(\sigma _{ij}^0\sigma _{ij}^1\sigma _{ij}^2\mathrm{}\sigma _{ij}^{L1})`$ denotes the configuration for the spins under the surface. Using the local tensors defined above, $`|\psi _L`$ is represented as a TPS:

$$|\psi _L=\underset{[\xi _L]}{}\underset{i+j=even}{}A_L\left\{\begin{array}{c}\sigma _{ij}^L\\ \xi _{ij}^L\end{array}\right\}.$$
(11)
Taking the limit $`L\mathrm{}`$, we have the maximum-eigenvalue eigenvector $`|\psi _{\mathrm{max}}`$, which is now represented as the product of the local tensor $`A_{\mathrm{}}`$.

From the practical view point of the usual power method, the eigenvector is improved with the relation $`|\psi _{L+1}=T|\psi _L`$ iteratively. In the framework of the TPS, this power method is reformulated as the recursion relation for the local tensor:

$$A_{L+1}\left\{\begin{array}{c}\sigma _{ij}^{L+1}\\ \xi _{ij}^{L+1}\end{array}\right\}=W\left\{\begin{array}{c}\sigma _{ij}^{L+1}\\ \sigma _{ij}^L\end{array}\right\}A_L\left\{\begin{array}{c}\sigma _{ij}^L\\ \xi _{ij}^L\end{array}\right\},$$
(12)
with which we can carry out the iterative calculation, until $`A_L`$ gives a good approximation of $`A_{\mathrm{}}`$. However it is generally difficult to store the tensor $`A_L`$ in the computer memory for a sufficiently large $`L`$, because if $`\xi _{ij}^L`$ has $`M`$ states, then the extended auxiliary variable $`\xi ^{L+1}`$ has $`2M`$ states.

In order to restrict the number of the auxiliary variable, we now import the idea of the DMRG into the TPS. The essence of the DMRG is that the increased number of states for $`\xi _{ij}^{L+1}`$ can be reduced, by using the “projection operator” generated from the “density matrix”. In the present case, the appropriate density matrix should be constructed for the spin variables $`(\sigma _{ij}^{L+1},\xi _{ij}^{L+1})`$ in the vertical direction (we thus call this density matrix as “the vertical density matrix”). Introducing the “transposed” local tensor:

$$\overline{A}\left\{\begin{array}{c}\eta _{ij}\\ \sigma _{ij}\end{array}\right\}=W\left\{\begin{array}{c}\overline{\sigma }_{ij}^0\\ \overline{\sigma }_{ij}^1\end{array}\right\}W\left\{\begin{array}{c}\overline{\sigma }_{ij}^1\\ \overline{\sigma }_{ij}^2\end{array}\right\}\mathrm{}W\left\{\begin{array}{c}\overline{\sigma }_{ij}^{L1}\\ \sigma _{ij}\end{array}\right\},$$
(13)
the explicit form of the vertical density matrix is defined as (See Fig 3.),

$`\rho _{kl}(\sigma _{kl}^{}\xi _{kl}^{}|\sigma _{kl}\xi _{kl})`$ $`=`$ $`\underset{[\sigma ],[\xi ],[\eta ]}{{\displaystyle }^{^{}}}\left[\underset{i+j=even}{{\displaystyle }^{^{}}}\overline{A}\left\{\begin{array}{c}\eta _{ij}\\ \sigma _{ij}\end{array}\right\}A\left\{\begin{array}{c}\sigma _{ij}\\ \xi _{ij}\end{array}\right\}\right]\overline{A}\left\{\begin{array}{c}\eta _{k^{\prime \prime }l^{\prime \prime }}\\ \stackrel{ˇ}{\sigma }_{k^{\prime \prime }l^{\prime \prime }}\end{array}\right\}A\left\{\begin{array}{c}\stackrel{ˇ}{\sigma }_{k^{\prime \prime }l^{\prime \prime }}\\ \stackrel{ˇ}{\xi }_{k^{\prime \prime }l^{\prime \prime }}\end{array}\right\}`$ (14)
$`\times `$ $`\overline{A}\left\{\begin{array}{c}\eta _{kl}\\ \sigma _{kl}\end{array}\right\}A\left\{\begin{array}{c}\sigma _{kl}\\ \xi _{kl}\end{array}\right\},`$ (15)
where $`^{^{}}`$ denotes the configuration sum for the all spin variables except $`\sigma _{kl}^{(^{})},\xi _{kl}^{(^{})}`$, and $`^{^{}}`$ means the product for the site index $`(i,j)`$ except at $`(i,j)=(k,l),(k^{\prime \prime },l^{\prime \prime })`$ with $`k^{\prime \prime }=k1`$ and $`l^{\prime \prime }=l1`$. In Eq. (15), we have also used the notation for the “checked spins”:

$`\left\{\begin{array}{c}\stackrel{ˇ}{\sigma }_{k^{\prime \prime }l^{\prime \prime }}\\ \stackrel{ˇ}{\xi }_{k^{\prime \prime }l^{\prime \prime }}\end{array}\right\}\left(\begin{array}{cccc}\sigma _{k^{\prime \prime }l^{\prime \prime }}& \sigma _{kl^{\prime \prime }}& \sigma _{kl}^{}& \sigma _{k^{\prime \prime }l}\\ \xi _{k^{\prime \prime }l^{\prime \prime }}& \xi _{kl^{\prime \prime }}& \xi _{kl}^{}& \xi _{k^{\prime \prime }l}\end{array}\right).`$

Further, we have omitted the sub(super)script $`L`$ if it is apparent.

For the case of the isotropic 3D Ising model, the vertical density matrix $`\rho _{kl}`$ becomes independent of the site index $`kl`$ in the thermodynamic limit. Thus we write $`\rho _{kl}`$ simply as $`\rho `$ hereafter. Moreover it should be noted that, for the isotropic case, the local tensors $`A`$ and $`\overline{A}`$ satisfies the relation:

$$\overline{A}\left(\begin{array}{cccc}\eta _{ij}& \eta _{i^{}j}& \eta _{i^{}j^{}}& \eta _{ij^{}}\\ \sigma _{ij}& \sigma _{i^{}j}& \sigma _{i^{}j^{}}& \sigma _{ij^{}}\end{array}\right)=A\left(\begin{array}{cccc}\sigma _{ij}& \sigma _{ij^{}}& \sigma _{i^{}j^{}}& \sigma _{i^{}j}\\ \eta _{ij}& \eta _{ij^{}}& \eta _{i^{}j^{}}& \eta _{i^{}j}\end{array}\right).$$
(16)
Following the spirit of the DMRG, we diagonalize $`\rho `$ to have the eigenvalues $`w_{\stackrel{~}{\xi }}`$ in the decreasing order $`w_1w_2\mathrm{}(0)`$:

$$\underset{\sigma \xi }{}\rho (\sigma ^{}\xi ^{}|\sigma \xi )U(\sigma \xi |\stackrel{~}{\xi })=U(\sigma ^{}\xi ^{}|\stackrel{~}{\xi })w_{\stackrel{~}{\xi }},$$
(17)
where $`U(\sigma \xi |\stackrel{~}{\xi })`$ is the eigenvector for $`w_\xi `$. By taking $`U(\sigma \xi |\stackrel{~}{\xi })`$ with $`\stackrel{~}{\xi }1,\mathrm{},M`$, we construct the projection operator $`U`$, which is a $`2M\times M`$ rectangular matrix. Operating $`U`$ to the spin variables on each “edge” of $`A`$, we then make the renormalized local tensor $`\stackrel{~}{A}`$:

$`\stackrel{~}{A}\left\{\begin{array}{c}\tau _{ij}\\ \stackrel{~}{\xi }_{ij}\end{array}\right\}`$ $``$ $`{\displaystyle \underset{\left\{\sigma _{ij}\right\}\left\{\xi _{ij}\right\}}{}}A\left(\begin{array}{cccc}\tau _{ij}& \tau _{i^{}j}& \tau _{i^{}j^{}}& \tau _{ij^{}}\\ \sigma _{ij}& \sigma _{i^{}j}& \sigma _{i^{}j^{}}& \sigma _{ij^{}}\\ \xi _{ij}& \xi _{i^{}j}& \xi _{i^{}j^{}}& \xi _{ij^{}}\end{array}\right)U(\sigma _{ij}\xi _{ij}|\stackrel{~}{\xi }_{ij})U(\sigma _{i^{}j}\xi _{i^{}j}|\stackrel{~}{\xi }_{i^{}j})U(\sigma _{i^{}j^{}}\xi _{i^{}j^{}}|\stackrel{~}{\xi }_{i^{}j^{}})`$ (21)
$`\times `$ $`U(\sigma _{ij^{}}\xi _{ij^{}}|\stackrel{~}{\xi }_{ij^{}}).`$ (22)
Using Eqs. (12) and (22) recursively, we can now calculate the effective local tensor $`\stackrel{~}{A}_{\mathrm{}}`$. However we encounter another problem in this process; it is also a numerically heavy problem to compute the vertical density matrix with Eq. (15).

Let us next explain how to overcome the difficulty in calculating the vertical density matrix. The key idea is that we can consider Eq. (15) as a kind of 2D classical spin system with a point defect. To see it, we here define a new tensor

$`G\left\{\begin{array}{c}\eta _{ij}\\ \sigma _{ij}\\ \xi _{ij}\end{array}\right\}\overline{A}\left\{\begin{array}{c}\eta _{ij}\\ \sigma _{ij}\end{array}\right\}A\left\{\begin{array}{c}\sigma _{ij}\\ \xi _{ij}\end{array}\right\},`$ (23)
which is graphically represented in Fig. 4.

Regarding the spin variable along the $`z`$ axis as a $`2M^2`$ state single spin:

$$𝝈_{ij}=\left(\begin{array}{c}\eta _{ij}\\ \sigma _{ij}\\ \xi _{ij}\end{array}\right),$$
(24)
we can see that $`G\left\{𝝈_{ij}\right\}`$ becomes the Boltzmann weight for the 2D effective classical model. In fact, the vertical density matrix $`\rho `$ can be expressed as the density matrix for the point defect in the 2D classical model:

$`\rho _{kl}(\sigma _{kl}^{}\xi _{kl}^{}|\sigma _{kl}\xi _{kl})=\underset{[𝝈]}{{\displaystyle }^{^{}}}\left[\underset{i+j=even}{{\displaystyle }^{^{}}}G\left\{𝝈_{ij}\right\}\right]G\left\{\stackrel{ˇ}{𝝈}_{k^{\prime \prime }l^{\prime \prime }}\right\}G\left\{𝝈_{kl}\right\}`$ (25)
where the meaning of the prime at the summation and the product is the same as Eq. (15), and the “checked spin” is given by

$`\left\{\stackrel{ˇ}{𝝈}_{k^{\prime \prime }l^{\prime \prime }}\right\}=\left\{\begin{array}{c}\eta _{k^{\prime \prime }l^{\prime \prime }}\\ \stackrel{ˇ}{\sigma }_{k^{\prime \prime }l^{\prime \prime }}\\ \stackrel{ˇ}{\xi }_{k^{\prime \prime }l^{\prime \prime }}\end{array}\right\}.`$

We calculate each component of the vertical density matrix (25) as a statistical average of a special local observable in 2D. To this end we apply the corner transfer matrix renormalization group (CTMRG)  which has been known to be quite efficient for 2D classical statistical systems.

Thus we have obtained a closed algorithm to calculate the local tensor $`A`$ with Eqs. (12) and (22), assisted by the CTMRG for the 2D effective classical model. We summarize the numerical procedure as follows:

* For the free-boundary condition in the $`Z`$ direction, define the initial local tensor $`A_1`$ as

$$A_1\left\{\begin{array}{c}\sigma _{ij}\\ \sigma _{ij}^{}\end{array}\right\}=W\left\{\begin{array}{c}\sigma _{ij}\\ \sigma _{ij}^{}\end{array}\right\}.$$
(26)
For the Ferro boundary condition, define $`A_1`$ as

$$A_1\left\{\begin{array}{c}\sigma _{ij}\\ \sigma _{ij}^{}\end{array}\right\}=W\left\{\begin{array}{c}\sigma _{ij}\\ \sigma _{ij}^{}\end{array}\right\}\times W\left\{\begin{array}{c}\sigma _{ij}^{}\\ +1\end{array}\right\},$$
(27)
where $`\{+1\}`$ means that the four Ising spins are aligned upward.
* Define the effective Boltzmann weight of the 2D classical spin system with Eq. (23).
* Perform the CTMRG calculation for the 2D effective classical spin system with the Boltzmann weight $`G`$ and obtain the vertical density matrix $`\rho `$.
* Diagonalize $`\rho `$ and construct the projection operator $`U`$.
* Renormalize the local tensor $`A`$ with Eq. (22).
* Return (b), until the local tensor $`A`$ is converged.

In this VDMA calculation, the accuracy is determined by the number of retained basis $`M`$ for the auxiliary variables $`\xi `$ and $`\eta `$, and $`m`$ for the CTMRG calculation in 2D classical system. We can check the convergence of the computed quantities with respect to $`M`$ and $`m`$.

## III Results

### A The 3D Ising model

Fig. 5 shows the spontaneous magnetization $`\sigma `$ calculated by using the VDMA. For comparison we also show the results of the 3D-version of the Kramers-Wannier (KW) approximation  and the Talapov and Blöte’s Monte Carlo results . For each $`M`$, we have carried out the VDMA calculations for $`m=4`$, $`8`$, $`12`$, and $`16`$, where we can find the good convergence with respect to $`m`$. We observe that the convergence with respect to $`M`$ is also sufficient in the off-critical region. Near the critical point, however, the magnetization becomes smaller as $`M`$ is increased, implying that large $`M`$ is needed for calculation in the critical region.

### B The 2D transverse field Ising model

Let us next consider the VDMA for the 2D transverse field Ising (TFI) model on the square-lattice, which is one of the fundamental quantum spin models in 2D. The Hamiltonian of the 2D TFI model is given by

$$H=J\underset{i,j}{}\sigma _i^z\sigma _j^z\mathrm{\Gamma }\underset{i}{}\sigma _i^x,$$
(28)
where $`J(>0)`$ is the ferromagnetic coupling constant, $`\sigma _i^\alpha `$ ($`\alpha =x,y,`$or $`z`$) are the Pauli matrices, and $`i,j`$ denotes the nearest-neighbor pairs on the square lattice. The transverse field $`\mathrm{\Gamma }`$ induces the quantum fluctuation into the system. At zero temperature, the TFI model exhibits the quantum phase transition at $`\mathrm{\Gamma }=\mathrm{\Gamma }_c`$, below which the spontaneous magnetization $`\sigma `$ behaves as $`\sigma (\mathrm{\Gamma }_c\mathrm{\Gamma })^\beta `$. The critical field $`\mathrm{\Gamma }_c`$ and the critical exponent $`\beta `$ have been estimated as $`\mathrm{\Gamma }_c=3.06`$ and $`\beta =0.31`$ by the quantum Monte Carlo (QMC) simulation.  In the following we consider the VDMA for the 2D TFI model at zero temperature, and show the results of $`\sigma `$.

As was described in the previous section, the VDMA is formulated for the 3D classical systems. In order to apply the VDMA to the 2D TFI model, we map the model to the 3D anisotropic Ising model by using the Suzuki-Trotter transformation . The partition function $`Z`$ of the 2D TFI model is obtained as the limit of the 3D anisotropic Ising model:

$`Z=\underset{L\mathrm{}}{lim}\mathrm{Tr}\mathrm{exp}\left[K_h{\displaystyle \underset{\tau =1}{\overset{L}{}}}{\displaystyle \underset{i,j}{}}\sigma _{i,\tau }\sigma _{j,\tau }+K_v{\displaystyle \underset{\tau =1}{\overset{L}{}}}{\displaystyle \underset{i}{}}\sigma _{i,\tau }\sigma _{i,\tau +1}\right],`$ (29)
where $`\sigma _{i,\tau }`$ is the Ising variable at the position $`i`$ and imaginary time $`\tau `$. The effective couplings $`K_h`$ and $`K_v`$ in Eq. (29) are given by

$`K_h=ϵJ,`$ (30)
$`K_v={\displaystyle \frac{1}{2}}\mathrm{log}[\mathrm{tanh}(ϵ\mathrm{\Gamma })],`$ (31)
$`ϵ=1/(TL),`$ (32)
where the subscripts $`h`$ and $`v`$ denote the horizontal ($`X`$-$`Y`$) direction and the vertical (Trotter) one, respectively. We can perform the VDMA calculation for this anisotropic 3D Ising model.

We should make a comment about the boundary condition, before proceeding to details. As was seen in the previous section, the open boundary condition is assumed in the VDMA. However, the periodic boundary condition is imposed along the Trotter direction in Eq. (29). As far as the zero-temperature properties are concerned, the boundary condition is inessential due to the double limit $`T0`$ and $`L0`$, allowing us to apply the VDMA to Eq. (29).

For a fixed value of $`ϵ`$, we calculate the magnetization $`\sigma (ϵ)`$ with the VDMA for the infinite volume. After obtaining $`\sigma (ϵ)`$ for various $`ϵ`$, we take the $`ϵ0`$ limit by extrapolation. In the actual calculation, we have observed the following $`ϵ`$-dependence

$$\sigma (ϵ)=\sigma (0)+\mathrm{const}\times ϵ^2,$$
(33)
which we adopt for the extrapolation formula.

In Fig. 6, we show the $`\sigma `$ .vs. $`ϵ^2`$ plots at $`\mathrm{\Gamma }=3.0`$ as an example, which are obtained with the VDMA of the numbers of retained bases $`(M,m)=(2,8)`$ and $`(3,18)`$, where the convergence with respect to $`m`$ is rapid. In the region $`\mathrm{\Gamma }>3.2`$, $`\sigma `$ converges to 0 smoothly.

In Fig. 7, we show the $`\sigma `$-$`\mathrm{\Gamma }`$ curve and the results of the series expansion for comparison. In $`\mathrm{\Gamma }<2.6`$, the VDMA results are sufficiently reliable, where the good convergence about $`m`$ and $`M`$ can be seen. We can further see the good agreement with the series expansion in the small field region ($`\mathrm{\Gamma }<2.0`$). In the vicinity of the critical point, however, the calculated magnetization exhibits the $`M`$ dependence. The roughly estimated critical field from the VDMA calculation is about 3.2, which is $`4\%`$ larger than the QMC one.

## IV Conclusion

In this paper we have constructed a novel higher-dimensional numerical renormalization algorithm which utilizes the natural tensor-product form of maximum-eigenvalue eigenstate $`|\psi _{\mathrm{max}}`$ of the transfer matrix. In our algorithm, called VDMA (vertical density matrix algorithm), the local tensor forming $`|\psi _{\mathrm{max}}`$ is iteratively updated using the vertical density matrix. We have successfully applied the VDMA to 3D Ising model. The VDMA can also be applied to 2D quantum systems using the Suzuki-Trotter transformation to 3D classical statistical systems, which we have demonstrated for the 2D transverse field Ising model. Application of the VDMA to other 2D quantum systems such as the Heisenberg model is an important subject of study, which we are now undertaking.

## Acknowledgments

Y. H. and K. O. was supported by the Japan Society for the Promotion of Science (JSPS). This work was partially supported by the “Research for the Future” program from the JSPS (JSPS-RFTF97P00201) and by the Grant-in-Aid for Scientific Research from Ministry of Education, Science, Sports and Culture (No. 12640393 and No. 11640376).


text:
# Abstract

## Abstract

Below a (comparatively large) measure of non-Hermiticity $`Z=Z_0^{(crit)}>0`$ of a $`𝒫𝒯`$ symmetrically complexified square well, bound states are constructed non-numerically. All their energies prove real and continuous in the (Hermitian) limit $`Z0`$. Beyond the threshold $`Z_0^{(crit)}`$ (and, in general, beyond $`Z_m^{(crit)}`$ at $`m=0,1,\mathrm{}`$) the lowest two real energies (i.e., $`E_{2m}`$ and $`E_{2m+1}`$) are shown to merge and disappear.

PACS 03.65.Ge, 03.65.Fd

## 1 Introduction

An interest of physicists in complex potentials $`V(x)=\mathrm{Re}V(x)+i\mathrm{Im}V(x)`$ with the generalized symmetry property

$$\mathrm{Re}V(x)=\mathrm{Re}V(x),\mathrm{Im}V(x)=\mathrm{Im}V(x)$$
(1)
dates back to the perturbative study of imaginary cubic anharmonic oscillators

$$V(x)=\omega x^2+i\lambda x^3$$
(2)
by Caliceti et al . In the early nineties, an increase of this interest  was inspired by the role of the imaginary cubic anharmonicity in field theory where eq. (1) mimics the fundamental parity times time-reversal (= $`𝒫𝒯`$) symmetry of phenomenological Lagrangians . Under certain circumstances, the $`𝒫𝒯`$ symmetric and Hermitian models can even prove mathematically equivalent .

In the late nineties, Bender and Boettcher  analysed the one-parametric family of the power-law $`\omega 0`$ models

$$V(x)=x^2(ix)^\delta $$
(3)
by the quasi-classical and purely numerical means. On this basis they conjectured that all the $`\delta >0`$ bound-state energies $`E_n=E_n(\delta )`$ form a real and discrete, smooth continuation of the well known harmonic-oscillator $`\delta =0`$ spectrum. An extension of this study inspired them later to apply the conjecture (formulated, originally, by Bessis at $`\delta =1`$ ) to a still broader class of interactions. Within the resulting “generalized” (so called $`𝒫𝒯`$ symmetric) quantum mechanics , there appears a growing number of interesting studies, promoting the ideas of supersymmetry , exact semiclassical techniques , functional analysis  and perturbation theory . All of them provide a consistent picture of a theory with certain (not quite well understood) limitations. Even in the above “guiding” example (3) the $`𝒫𝒯`$ symmetry breaks down spontaneously at $`\delta <0`$ .

Phenomenological appeal of the $`𝒫𝒯`$ symmetric modifications need not even stop before the “sacred” quantum electrodynamics . Still, the formalism abounds in open questions and its mathematical foundations are mostly conjectures. This is the reason why the attention quickly spreads from the phenomenologically oriented eq. (3) to its exactly solvable alternatives. In this direction, the partially solvable extension of the $`\delta =2`$ quartic case  and the discovery of the exact solvability of the certain version of the $`\delta \mathrm{}`$ limiting case  were the first and encouraging successes. They were followed by the $`𝒫𝒯`$ symmetric regularization of the more-dimensional harmonic oscillator  and by the formulation and solution of many further shape-invariant models -.

Once we move beyond the domain of analytic potentials, numerical studies provide significantly less encouraging results . This is the reason why the “most elementary” square well

$$\begin{array}{c}\mathrm{Re}V(x)=0,x(1,1)\\ \mathrm{Im}V(x)=Z,x(1,0)\\ \mathrm{Im}V(x)=Z,x(0,1)\end{array}$$
(4)
has always been $`𝒫𝒯`$ symmetrized just in the (various, non-equivalent) $`\delta \mathrm{}`$ limits of eq. (3.

## 2 Solution

Presumably, the reason for absence of the “forgotten” $`𝒫𝒯`$ symmetric square well (4) in the literature lies in an ambiguity of its continuation beyond the discontinuities at $`x=\pm 1`$. We are going to treat this problem simply by imposing the Dirichlet boundary conditions on our (complex) wave functions,

$$\psi (\pm 1)=0.$$
(5)
Having made the latter decision the explicit computations are really elementary. Their essence lies in the easy construction of the general right and left solutions

$$\begin{array}{c}\psi _+=c_+e^{\kappa x}+d_+e^{\kappa x}\\ \psi _{}=c_{}e^{\kappa ^{}x}+d_{}e^{\kappa ^{}x}\end{array}$$
(6)
with the complex $`\kappa `$ and its conjugate $`\kappa ^{}`$. We re-parametrize the (by assumption, real) energies $`E=t^2s^2`$ and the measure of non-Hermiticity $`Z=2st`$. This gives the complex exponents $`\kappa =sit`$ since $`\kappa ^2=V(x)E`$ is also complex (and constant for all $`x(0,1)`$). Our solutions (6) are made explicit, for every real $`E`$ and $`Z`$, when we use the inverse formulae

$$t=\frac{1}{\sqrt{2}}\left(E+\sqrt{E^2+Z^2}\right)^{1/2},$$

$$s=\frac{Z}{\sqrt{2}}\left(E+\sqrt{E^2+Z^2}\right)^{1/2}.$$

In the light of the $`𝒫𝒯`$ symmetry of our Hamiltonian $`H=𝒫𝒯H𝒫𝒯`$, the product $`𝒫𝒯\psi (x)\psi ^{}(x)`$ will satisfy the same Schrödinger equation as $`\psi (x)`$. Hence, in the origin, we are permitted to normalize our bound states in a $`𝒫𝒯`$ symmetric way,

$$\psi _+(0)=\psi _{}(0)=1,_x\psi _+(0)=_x\psi _{}(0)=iA.$$

These conditions contain a free real parameter $`A`$ and are equivalent to the matching of wave functions,

$$\begin{array}{c}c_\pm =1d_\pm ,\\ \kappa (12d_+)=\kappa ^{}(12d_{})=iA.\end{array}$$

This implies that we know all the coefficients in eq. (6),

$$d_+=\frac{1}{2}i\frac{A}{2\kappa },d_{}=\frac{1}{2}i\frac{A}{2\kappa ^{}}.$$

It is easy to satisfy the external boundary conditions (5) and reduce them to the elementary prescription

$$iA=\kappa \mathrm{coth}\kappa .$$
(7)
In terms of the real parameters $`s`$ and $`t`$ this represents a system of two algebraic equations,

$$A=\frac{s}{\mathrm{tan}t}+t\mathrm{tanh}s=s\mathrm{tan}t+\frac{t}{\mathrm{tanh}s}.$$
(8)
Its first part defines the (necessarily, real) value of $`A=A(s,t)`$. A re-arrangement of the second algebraic equation is elementary and gives the rule

$$2t\mathrm{sin}2t=2s\mathrm{sinh}2s.$$
(9)
As long as both its sides depend on the mere absolute values of the respective variables $`t`$ and $`s`$, we may pick up $`t0`$, insert the definition of $`s=Z/2t`$ and solve this equation numerically.

## 3 Discussion

By construction, the spectrum of the real energies $`E=E_n`$ is defined in terms of the roots $`t_n`$ of eq. (9), $`E_n=t_n^2s_n^2`$, $`s_n=Z/2t_n`$, $`n=0,1,\mathrm{}`$. After we re-scale $`tT=2t/\pi `$, we immediately see that in the limit $`Z0`$ this spectrum degenerates to the known Hermitian one, with $`T_n=n+1`$ etc. From our equation (9) it is obvious that at the very small $`Z`$ (and, hence, $`s0`$), the change of the above roots $`T_n`$ (and, of course, energies) remains very small as well.

At the larger (and, say, positive) $`Z`$, the analysis of our quantization condition (9) becomes significantly simplified by the re-scalings $`sS_0=2s\mathrm{sinh}2s`$ and $`S_0S`$ such that, say, $`S_0=4\mathrm{sinh}^2S`$. Both these steps represent a one-to-one mapping exhibiting the strict leading-order quasi-linearity $`sS`$ achieved at both ends of our half-axis, i.e., for $`s0`$ as well as for $`s\mathrm{}`$. Moreover, this introduces just a minimal deformation of the scale of the coordinate $`s`$ (practically invisible on a picture) and replaces equation (9) by a new one,

$$\pi T\mathrm{sin}\pi T=4\mathrm{sinh}^2S.$$
(10)
This new relation is explicitly solvable. The resulting analytic and $`Z`$independent formula $`S=X(T)`$ defines a curve in our new $`ST`$ plane. Its $`T<14`$ part is displayed in Figure 1.

In order to determine the separate roots $`T_n`$ (and, hence, the spectrum of energies), it remains for us to recollect the $`Z`$dependent, hyperbolic constraint $`s=Z/2t`$. After we translate its form in our new variables, $`S=Y(Z,T)`$, we discover that it has a monotonous, hyperbolic shape which depends on $`Z`$. Figure 1 offers a few samples. We may conclude that each point of the intersection of our two curves $`X(T)`$ and $`Y(Z,T)`$ determines a $`Z`$dependent root $`T=T(Z)`$ and, hence, a real energy.

The $`Z`$ dependence of the roots $`T(Z)`$ is, in general, smooth. A non-perturbative effect is only encountered at certain critical values $`Z=Z^{(crit)}`$. In the vicinity of these points, the two leftmost roots $`T_0(Z)<T_1(Z)`$, $`Z<Z^{(crit)}`$ merge in a single, doubly degenerate real root $`T_0(Z^{(crit)})=T_1(Z^{(crit)})`$. Figure 1 gives two examples and shows that even the smallest one of these critical values is already quite large ($`Z_0^{(crit)}4.48`$). Below this bound we may summarize that

* our complexified square well generates the infinite set of real energies;
* the roots $`T_n`$ which define these energies are almost equidistant, especially at the higher $`n`$;
* as expected, the standard square-well-type behaviour of the spectrum is reproduced at the small $`Z`$ and for the highly excited states.

In the strongly non-Hermitian domain, the lowest doublets of states subsequently disappear. Presumably, their energies dissolve in conjugate pairs in complex plane. The $`𝒫𝒯`$ symmetry of their wave functions breaks down. The related solutions become “unphysical” and have to be omitted in a way paralleling the similar “disappearance of states” at $`\delta <0`$ in the model (3) of ref. . The resulting sudden upward jump of the ground-state energy definitely enters the list of paradoxes, emerging in the other exactly solvable models. Their present list already involves the unavoided level crossings in the harmonic and Coulomb oscillators , the anomalously large excitations in the $`𝒫𝒯`$ symmetrized but bounded Rosen-Morse field , some unexpected manifestations of the strong singularities in the Pösch-Teller and Eckart models , and a spontaneous re-ordering of levels in the Morse asymmetric well .

In this context, the role of the complexified square well is as exceptional as in the Hermitian limit. We may expect that its generalizations with more points of discontinuity will remain tractable analytically. This could be of a significant help, say, within perturbation theory . In the more pragmatic numerical setting, one should recollect that the so called Prüfer transformation  (i.e., nothing but a square-well-inspired use of the locally exponential solutions) found a firm place in the standard computer software . Last but not least, one has to keep in mind that in the Hermitian quantum mechanics the use of the locally constant forces could also clarify the various manifestations of the pertaining Sturm Liouville theory . An appropriate $`𝒫𝒯`$ symmetrization of this theory is expected to be quite a difficult task . In such a direction, also the knowledge of our present solutions could mediate a further progress, hopefully, in the near future.

## Acknowledgement

Partially supported by the GA AS grant Nr. A 104 8004.