text:
# Nonperturbative Renormalization Flow and Essential Scaling for the Kosterlitz-Thouless Transition

## Abstract

The Kosterlitz-Thouless phase transition is described by the nonperturbative renormalization flow of the two dimensional $`\phi ^4`$-model. The observation of essential scaling demonstrates that the flow equation incorporates nonperturbative effects which have previously found an alternative description in terms of vortices. The duality between the linear and nonlinear $`\sigma `$-model gives a unified description of the long distance behaviour for $`O(N)`$-models in arbitrary dimension $`d`$. We compute critical exponents in first order in the derivative expansion.

The Kosterlitz-Thouless phase transition  may describe the critical behaviour of various two dimensional systems. It poses a challenge to our theoretical understanding due to several uncommon features. The low temperature phase exhibits a massless Goldstone boson like excitation despite the fact that the global $`U(1)`$ symmetry is not spontaneously broken by a standard order parameter. In this phase the critical exponents depend on the temperature. In the high temperature phase the approach to the transition is not governed by critical exponents but rather by essential scaling. Modern nonperturbative renormalization group methods have described the characteristic properties of the low temperature phase qualitatively . In this letter we extend this analysis to the essential scaling in the high temperature phase and provide for a quantitatively accurate discussion of the low temperature phase. We present a universal renormalization group description of models with $`O(N)`$ symmetry in arbritrary dimension.

We employ the concept of the effective average action $`\mathrm{\Gamma }_k`$ , which equals the effective action $`\mathrm{\Gamma }`$ apart from the fact that in the former only fluctuations with momenta larger than $`k`$ are included. For large momentum scales $`\mathrm{\Gamma }_k`$ therefore equals some appropriately regularized microscopic action while for $`k=0`$ we have $`\mathrm{\Gamma }_0=\mathrm{\Gamma }`$. Thus $`\mathrm{\Gamma }_k`$ interpolates between microscopic and macroscopic scales. Since $`\mathrm{\Gamma }`$ is the generating functional for the 1PI correlation functions it specifies directly the quantities of interest like the correlation length $`\xi =m_R^1`$. The flow of $`\mathrm{\Gamma }_k`$ is obeys an exact renormalization group equation ,

$$_t\mathrm{\Gamma }_k[\phi ]=\frac{1}{2}\text{Tr}\left((\mathrm{\Gamma }_k^{(2)}[\phi ]+R_k)^1_tR_k\right).$$
(1)
Here $`\mathrm{\Gamma }_k^{(2)}`$ is the second functional derivative, $`_t`$ denotes the logarithmic derivative $`k/k`$ and the trace (in momentum space) reads $`\text{Tr}=_{a=1}^Nd^dq/(2\pi )^d`$. The cutoff $`R_k(q^2)`$ suppresses the low momentum modes. We use a cutoff of the form

$$R_k(q^2)=Z_kq^2r(q^2/k^2)=\frac{Z_kq^2}{\mathrm{exp}(q^2/k^2)1},$$
(2)
where the wave function renormalization $`Z_k`$ will be fixed later. In order to solve equation (1) numerically one has to truncate the most general form of $`\mathrm{\Gamma }_k`$. We introduce dimensionless, renormalized fields, $`\stackrel{~}{\phi }_a=Z_k^{1/2}k^{(2d)/2}\phi _a`$, $`\stackrel{~}{\rho }=(1/2)\stackrel{~}{\phi }_a\stackrel{~}{\phi }_a`$ and parametrize $`\mathrm{\Gamma }_k`$ in first order in a derivative expansion by

$`\mathrm{\Gamma }_k`$ $`=`$ $`{\displaystyle }d^dx\{k^du_k(\stackrel{~}{\rho })+{\displaystyle \frac{1}{2}}k^{d2}z_k(\stackrel{~}{\rho })_\mu \stackrel{~}{\phi }_a_\mu \stackrel{~}{\phi }_a`$ (4)
$`+{\displaystyle \frac{1}{4}}k^{d2}\stackrel{~}{y}_k(\stackrel{~}{\rho })_\mu \stackrel{~}{\rho }_\mu \stackrel{~}{\rho }+O(^4)\}.`$
The flow of $`\mathrm{\Gamma }_k`$ is then given by the flow of the quantities $`u`$, $`z`$ and $`\stackrel{~}{y}`$ which are dimensionless and depend only on the $`O(N)`$-invariant $`\stackrel{~}{\rho }`$ and on $`k`$. We denote by $`\kappa `$ the running minimum of the potential $`u_k(\stackrel{~}{\rho })`$ and fix $`Z_k`$ by requiring $`z_k(\kappa )=1`$.

The partial differential equations

$`_tu`$ $`=`$ $`du+(d2+\eta )\stackrel{~}{\rho }u^{}`$ (6)
$`+2v_d(N1)l_0^d(w,z,\eta )+2v_dl_0^d(\stackrel{~}{w},\stackrel{~}{z},\eta ),`$
$`_tz`$ $`=`$ $`\eta z+\stackrel{~}{\rho }z^{}(d2+\eta )`$ (15)
$`(4v_d/d)\stackrel{~}{\rho }^1\{.m_{2,0}^d(w,z,\eta )`$
$`2m_{1,1}^d(w,\stackrel{~}{w},z,\stackrel{~}{z},\eta )+m_{0,2}^d(\stackrel{~}{w},\stackrel{~}{z},\eta ).\}`$
$`2v_d(\stackrel{~}{z}z)\stackrel{~}{\rho }^1\{.l_1^d(\stackrel{~}{w},\stackrel{~}{z},\eta )`$
$`(2/d)(\stackrel{~}{z}z)l_2^{d+2}(\stackrel{~}{w},\stackrel{~}{z},\eta ).\}`$
$`2v_dz^{}\{.(N1)l_1^d(w,z,\eta )`$
$`(8/d)n_{1,1}^d(w,\stackrel{~}{w},z,\stackrel{~}{z},\eta )`$
$`+\left(5+2z^{\prime \prime }\stackrel{~}{\rho }/z^{}\right)l_1^d(\stackrel{~}{w},\stackrel{~}{z},\eta )`$
$`(4/d)z^{}\stackrel{~}{\rho }l_{1,1}^{d+2}(w,\stackrel{~}{w},z,\stackrel{~}{z},\eta ).\},`$
$`_t\stackrel{~}{z}`$ $`=`$ $`\eta \stackrel{~}{z}+\stackrel{~}{\rho }\stackrel{~}{z}^{}(d2+\eta )`$ (29)
$`2v_d(\stackrel{~}{z}^{}+2\stackrel{~}{\rho }\stackrel{~}{z}^{\prime \prime })l_1^d(\stackrel{~}{w},\stackrel{~}{z},\eta )`$
$`+8v_d\stackrel{~}{\rho }\stackrel{~}{z}^{}(3u^{\prime \prime }+2\stackrel{~}{\rho }u^{\prime \prime \prime })l_2^d(\stackrel{~}{w},\stackrel{~}{z},\eta )`$
$`+4v_d\left(2+1/d\right)\stackrel{~}{\rho }(\stackrel{~}{z}^{})^2l_2^{d+2}(\stackrel{~}{w},\stackrel{~}{z},\eta )`$
$`(8/d)v_d\stackrel{~}{\rho }(3u^{\prime \prime }+2\stackrel{~}{\rho }u^{\prime \prime \prime })^2\stackrel{~}{m}_4^d(\stackrel{~}{w},\stackrel{~}{z},\eta )`$
$`(16/d)v_d\stackrel{~}{\rho }\stackrel{~}{z}^{}(3u^{\prime \prime }+2\stackrel{~}{\rho }u^{\prime \prime \prime })\stackrel{~}{m}_4^{d+2}(\stackrel{~}{w},\stackrel{~}{z},\eta )`$
$`(8/d)v_d\stackrel{~}{\rho }(\stackrel{~}{z}^{})^2\stackrel{~}{m}_4^{d+4}(\stackrel{~}{w},\stackrel{~}{z},\eta )`$
$`+(N1)v_d\{.`$
$`2\left(\stackrel{~}{z}^{}\stackrel{~}{\rho }^1(\stackrel{~}{z}z)\right)l_1^d(w,z,\eta )`$
$`(8/d)\stackrel{~}{\rho }(u^{\prime \prime })^2m_4^d(w,z,\eta )`$
$`(16/d)\stackrel{~}{\rho }u^{\prime \prime }z^{}m_4^{d+2}(w,z,\eta )`$
$`(8/d)\stackrel{~}{\rho }(z^{})^2m_4^{d+4}(w,z,\eta )`$
$`+4(\stackrel{~}{z}z)u^{\prime \prime }l_2^d(w,z,\eta )`$
$`.+4(z^{}(\stackrel{~}{z}z)+(1/d)\stackrel{~}{\rho }(z^{})^2)l_2^{d+2}(w,z,\eta )\}`$
invoke the “threshold functions” $`l_{n_1,n_2}^d`$, $`m_{n_1,n_2}^d`$, $`\stackrel{~}{m}_{n_1,n_2}^d`$ and $`n_{n_1,n_2}^d`$ defined by the integral

$$\frac{1}{2}\underset{0}{\overset{\mathrm{}}{}}𝑑yy^{\frac{d}{2}1}\stackrel{~}{_t}\left\{\frac{X}{(p(y)+w)^{n_1}(\stackrel{~}{p}(y)+\stackrel{~}{w})^{n_2}}\right\},$$
(30)
with $`X=1,y(_yp)^2,y(_y\stackrel{~}{p})^2,y_yp`$ for $`l,m,\stackrel{~}{m},n`$ respectively. We have defined $`u^{}=u/\stackrel{~}{\rho }`$ and we use the shorthands $`v_d^1=2^{d+1}\pi ^{d/2}\mathrm{\Gamma }(d/2),w=u^{},\stackrel{~}{w}=u^{}+2\stackrel{~}{\rho }u^{\prime \prime },\stackrel{~}{z}=z+\stackrel{~}{\rho }\stackrel{~}{y},p(y)=y(z+r(y))`$, and $`\stackrel{~}{p}(y)=y(\stackrel{~}{z}+r(y)).`$ The derivative $`\stackrel{~}{}_t`$ only acts on the $`k`$-dependence of the cutoff $`R_k`$, i.e. $`\stackrel{~}{}_tp(y)=y\left(\eta r(y)+2y_yr(y)\right)`$ and we note that $`\stackrel{~}{}_t_yp=_y\stackrel{~}{}_tp`$. Finally we abbreviate $`l_{n,0}^d=l_n^d`$ etc., where $`l_0^d`$ is defined by the rule $`(p+w)^n\mathrm{log}(p+w)`$. The expression for the anomalous dimension $`\eta =_t\mathrm{ln}Z_k`$ can be obtained from the identity $`_tz_k(\kappa _k)0`$. For $`N=1`$ one has $`\stackrel{~}{y}=0`$, $`z=\stackrel{~}{z}`$.

These equations are valid in arbitrary dimension. In the very simple approximation $`z=\stackrel{~}{z}=1`$, $`u=\frac{1}{2}\lambda (\stackrel{~}{\rho }\kappa )^2`$ these equations already give a correct qualitative picture for $`O(N)`$ symmetric models in arbitrary dimension . We show that the present version leads to quantitatively accurate results. This is done by a numerical solution with initial values specified at a microscopic scale $`k=\mathrm{\Lambda }`$. For $`\kappa 1`$ the evolution is dominated by the $`N1`$ Goldstone modes ($`N>1`$). More precisely, the threshold functions at the minimum vanish rapidly for $`\stackrel{~}{w}=2\kappa u^{\prime \prime }(\kappa )1`$. For $`N>2`$ the coupling of the nonlinear $`\sigma `$\- model for the Goldstone bosons is given by $`\kappa ^1`$.

We concentrate first on $`d=2`$, where the universality of the $`\beta `$-function for the nonlinear coupling implies for $`\kappa 1`$ an asymptotic form ($`N2`$)

$$_t\kappa =\beta _\kappa =\frac{N2}{4\pi }+\frac{N2}{16\pi ^2}\kappa ^1+𝒪(\kappa ^2).$$
(31)
In the linear description one can easily obtain an equation for $`\kappa `$ by using $`_tu_k(\kappa _k)=0`$ together with equation (6). By evaluating the above equations for large $`\kappa `$ it is possible to compare with (31). Previously it was found  that in a much simpler truncation one already obtains the correct lowest order in the above expansion. In order to reproduce the exact two loop result one has, however, to go even beyond the truncation (4).

From equation (31) we expect that $`\kappa `$ will run only marginally at large $`\kappa `$. As a consequence the flow of the action follows a single trajectory for large $`t`$ and can be characterized by a single scale. Notice that the perturbative $`\beta `$-function vanishes for $`N=2`$ since the Goldstone bosons are not interacting in the abelian case. Thus for large $`\kappa `$ one expects a line of fixed points which can be parameterized by $`\kappa `$. This fact plays a major role in the discussion of the Kosterlitz-Thouless transition below. It is responsible for the temperature dependence of the critical exponents. We have evaluated $`\beta _\kappa =\beta _\kappa ^{(1)}(N2)/(4\pi )+\beta _\kappa ^{(2)}\kappa ^1(N2)/(16\pi ^2)+\mathrm{}`$ numerically from the solution of eq. (6) - (29) and extracted the expansion coefficients for large $`\kappa `$ (see table I).

For the nonabelian nonlinear sigma model in $`d=2`$, $`N>2`$ there exists an exact expression  for the ratio of the renormalized mass $`m_R`$ and the scale $`\mathrm{\Lambda }_{\overline{MS}}`$ which characterizes the two loop running coupling in the $`\overline{MS}`$ scheme by dimensional transmutation. The flow equation (1) together with a choice of the cutoff $`R_k`$ and the initial conditions also defines a renormalization scheme. The corresponding parameter $`\mathrm{\Lambda }_{ERGE}`$ specifies the two loop perturbative value of the running coupling $`\kappa ^1`$ similar to $`\mathrm{\Lambda }_{\overline{MS}}`$ in the $`\overline{MS}`$ scheme. The numerical solution of the flow equation permits us to compute $`m_R/\mathrm{\Lambda }_{ERGE}`$. (Two loop accuracy would be needed for a quantitative determination of $`\mathrm{\Lambda }_{ERGE}/\mathrm{\Lambda }_{\overline{MS}}`$.) In table I we compare our results with the exact value of $`m_R/\mathrm{\Lambda }_{\overline{MS}}`$. We also report the ratio $`m_R/k_s`$ with $`k_s`$ defined by $`\kappa (k_s)=0`$.

The abelian case, $`N=2`$, is known to exhibit a special kind of phase transition which is usually described in terms of vortices . The characteristics of this transition are a massive high temperature phase and a low temperature phase with divergent correlation length but zero magnetization. The anomalous dimension $`\eta `$ depends on $`T`$ below $`T_c`$ and is zero above. It takes the exact value $`\eta _{}=0.25`$ at the transition. The most distinguishing feature is essential scaling for the temperature dependence of $`m_R`$ just above $`T_c`$,

$$m_Re^{\frac{b}{(TT_c)^\zeta }},\zeta =\frac{1}{2}.$$
(32)
Our aim is to draw a picture of this transition from the flow equations without ever examinig vortex configurations.

We have already mentioned the existence of a line of fixed points for large values of $`\kappa `$, which is relevant for the low temperature phase. The contribution of a massless (Goldstone) boson in the renormalization group equation ($`w(k)=0`$) is responsible for the finite value of $`\eta `$. This in turn drives the expectation value of the unrenormalized field to zero (even for a nonvanishing renormalized expectation value $`\kappa `$),

$$\rho _0=\kappa /Z_k\kappa (k/\mathrm{\Lambda })^\eta .$$
(33)
We can observe this line of fixed points to a good approximation (cf. fig. 1), although the vanishing of $`\beta _\kappa `$ is not exact (we find $`\beta _\kappa 410^5\kappa ^1+𝒪(\kappa ^2)`$). The line of fixed points ends at a phase transition which corresponds to a microscopic parameter $`\kappa _\mathrm{\Lambda }`$ as an initial value of the flow. In order to verify essential scaling we have to examine the flow for values of $`\kappa _\mathrm{\Lambda }`$ just below that point, $`\kappa _\mathrm{\Lambda }=\kappa _\mathrm{\Lambda }+\delta \kappa _\mathrm{\Lambda }`$, $`\delta \kappa _\mathrm{\Lambda }(TT_c)`$. Then $`\kappa _k`$ crosses zero at the scale $`k_s`$ and we find the mass by continuing the flow in the symmetric regime (minimum at $`\kappa =0`$). In figure 2 we plot $`(\mathrm{ln}(m_R/\mathrm{\Lambda }))^2`$ against $`\kappa _\mathrm{\Lambda }`$ and find excellent agreement with the straight line (32).

How does $`\beta _\kappa `$ have to look like in order to yield essential scaling? Since there is only one independent scale near the transition, one expects $`m_R(T)=C_sk_s(T)`$, where $`k_s`$ denotes the scale at which $`\kappa `$ vanishes, i.e. $`\kappa (k_s,T)=0`$. For $`\kappa `$ close to and below $`\kappa _{}`$ we parameterize $`\beta _\kappa `$ (this approximation is not valid for very small $`\kappa `$)

$$\beta _\kappa =\frac{1}{\nu }(\kappa _{}\kappa )^{\zeta +1}.$$
(34)
For conventional scaling one expects $`\zeta =0`$ and the correlation length exponent is given by $`\nu `$. Integrating equation (34) yields for $`\zeta 0`$, $`\delta \kappa =\kappa \kappa _{}`$:

$$\mathrm{ln}(k/\mathrm{\Lambda })=\frac{\nu }{\zeta }\left(\frac{1}{(\delta \kappa )^\zeta }\frac{1}{(\delta \kappa _\mathrm{\Lambda })^\zeta }\right).$$
(35)
For $`k=k_s`$ the first term $`(\delta \kappa )^\zeta `$ is small and independent of $`T`$ (since $`\delta \kappa (k_s)=\kappa _{}`$) and equation (35) yields the essential scaling relation (32) for $`\zeta =1/2`$. Usually, the microscopic theory is such that one does not start immediately in the vicinity of the critical point and the approximation (34) is not valid for $`k\mathrm{\Lambda }`$. However, if one is near the critical temperature the trajectories will stay close to the critical one, $`\kappa _c(t)`$, with $`\kappa _c(0)=\kappa _\mathrm{\Lambda }`$. This critical trajectory converges rapidly to its asymptotic value $`\kappa _{}`$ and $`\beta _\kappa `$ gets close to eq. (34) at some scale $`\mathrm{\Lambda }^{}<\mathrm{\Lambda }`$. As a result, one may use equation (35) only in its range of validity ($`k<\mathrm{\Lambda }^{}`$) and observe that $`\delta \kappa _\mathrm{\Lambda }^{}`$ is also proportional to $`T_cT`$. The numerical verification of (34) is quite satisfactory: Fitting our data yields $`\kappa _{}=0.248`$, $`\zeta =0.502`$ and $`\nu ^1=2.54`$. The uncertainty for $`\zeta `$ is approximately $`\pm 0.05`$. The numerical values of $`\beta _\kappa `$ and the approximation (34) are shown in figure 3.

One can use the information from figure 2 or 3 in order to determine $`\kappa _{}`$ and therefore $`\eta _{}=\eta (\kappa _{})`$, the anomalous dimension at the transition. We plot $`\eta `$ against $`\kappa `$ in figure 4. One reads off $`\eta _{}=0.287\pm 0.007`$ where the error reflects the two methods used to compute $`\kappa _{}`$ and does not include the truncation error. For $`\kappa _\mathrm{\Lambda }>\kappa _\mathrm{\Lambda }`$ or $`T<T_c`$ the running of $`\kappa `$ essentially stops after a short “initial running” towards the line of fixed points. One can infer from figure 4 the temperature dependence of the critical exponent $`\eta `$ for the low temperature phase. In summary all the relevant characteristic features of the Kosterlitz-Thouless transition are visible within our approach.

We end this letter by reporting the values of the critical exponents obtained in our approximation (6) - (29) for the “standard” second order phase transitions for $`d=2`$, $`N=0,1`$ and $`d=3`$, $`N0`$. In tables II, III they are compared with exact results or “averages” (only for the simplicity of the display!) of results from various other methods .

The agreement is very satisfactory! We also characterize in tables IV, V the scaling solution relevant for the second order transition by quoting $`\kappa _{},\lambda _{}=u_{}^{\prime \prime }(\kappa _{}),u_3=u_{}^{\prime \prime \prime }(\kappa _{})`$ as well as $`z_{}^{}(\kappa _{})`$ and $`\stackrel{~}{z}_{}(\kappa _{})`$. We conclude that the first order in the derivative expansion of the exact flow equation for the effective average action gives a quantitatively accurate picture of all phase transitions of scalar models in the $`O(N)`$ universality class for arbitrary dimension $`2d4`$.


format: text

text:
# Optimizing Schrödinger functionals using Sobolev gradients: Applications to Quantum Mechanics and Nonlinear Optics

## 1 Introduction

The observation of Nature reveals that many unforced continuous systems tend to accommodate into stationary configurations, in which the distributions of mass, charge, velocity, etc, do not change throughout time. In the language of mathematical modeling all configurations are represented by points of a certain space of functions, $`\psi (𝐱)W`$, while the tendency of the system to lie into any of these states is given by a functional, $`E(\psi ):W`$, the energy, whose minima are precisely those stationary “states”. For this reason it is possible to see many physical problems written as variational principles of the type “find $`\psi W`$ such that $`E(\psi ):W`$ achieves a minimum on $`W`$”. In most situations the functional to be minimized has a dependence on $`\psi `$ of the form

$$E(\psi )=f(\psi (𝐱),\psi (𝐱))d^nx.$$
(1)
However, a complete analytical description of the minima of the functional $`E(\psi )`$ is usually not possible. In this paper we will introduce several techniques for performing this study numerically, focusing on the minimization of $`E(\psi )`$ subject to physical constraints.

From a practical point of view, this problem is similar to that of finding the minima of a real function defined over a finite-dimensional space, such as $`^n`$. First, a definition of derivative of the functional, $`E(\psi )`$, must be chosen. If the domain of the functional $`W`$ is equipped with some scalar product we may use the Frechet derivative which is given by a first order expansion of the functional around a function $`\psi `$

$$E(\psi +\delta )=E(\psi )+\delta ,E(\psi )+E(\psi ),\delta +O\left(\delta ^2\right).$$
(2)
The *critical points,* $`\psi _c`$, are defined as the points where the first order variation of the functional vanishes for any perturbation $`\delta `$. That is, the derivative vanishes in a weak sense

$$\delta ,E(\psi _c)=0,\delta .$$
(3)
Just like in the finite-dimensional case it is possible to show that any minimum of the functional must be also a critical point. Thus a common approach is to solve Eq. (3) and verify a posteriori which solutions are actually minima of the functional. If $`W=L^2(^n)`$ this procedure gives us the well known Euler-Lagrange equations of the problem, which is a partial derivatives equation (PDE)

$$\frac{f}{\overline{\psi }}\frac{f}{\overline{\psi }}=0.$$
(4)
However, we have no guarantee to reduce the complexity of the problem, as it is by no means trivial to solve Eq. (4). Also we are likely to obtain more solutions than we actually need, since not only minima, but also maxima and saddle node points will satisfy the Lagrange equations.

To avoid these problems some other methods are used which aim at finding the minima of the functional directly, constructing minimizing sequences, $`\{\psi _i\}`$, whose limit is a minimum of the functional: $`\psi _c=lim_i\mathrm{}\psi _i`$. These methods will be discussed in the following sections.

The outline of this paper is as follows. In §2 we recall the definition of Sobolev gradients as given in . We derive a formal solution to the problem of finding these gradients which is based on the inversion of a positive hermitian operator. In §3 we derive an explicit expression for the Sobolev gradients in the trigonometric Fourier basis and comment on its implementation using Fast Fourier Transforms (FFT).

In §4 and §5 we apply the previous tools to two physical problems. Using descent techniques with Sobolev gradients over Fourier spaces we will find the ground states of a Bose-Einstein condensate in a rotating magnetic trap and the excited states for coupled laser beams propagating through a nonlinear medium. Both physical systems are modeled by nonlinear equations of Schrödinger type and present difficulties when traditional minimization techniques are used. We comment on the great improvements that are achieved using Sobolev gradients. Finally in §6 we summarize our results and offer some conclusions.

## 2 Sobolev gradients

### 2.1 Direct solutions of the variational problem

There are two traditional approaches to the problem of finding the minima of a functional using a discrete basis. The first one expands the unknown solution using a Fourier basis $`\{\varphi _k\}`$, $`\psi =_kc_k\varphi _k,`$ and defines a new functional over the finite–dimensional space

$$E(\{c_k\})E\left(c_k\varphi _k\right).$$
(5)
The functional is then minimized using methods which are well known from the the domain of finite–dimensional problems, e.g. Newton’s method or nonlinear conjugate gradient.

This procedure is quite straightforward and there is a huge amount of literature and tools which can be immediately applied to (5). However, for some types of problems one has to deal with highly nonlinear algebraic equations with many of terms on each equation, something which is computationally too expensive to work with.

The second approach involves what is known as *descent techniques*. The idea is to manipulate the original functional (1) building an analytic equation for a minimizing trajectory in the target space $`W`$. This equation is then discretized and solved on a suitable basis. In the *continuous steepest descent* version, the trajectory $`\psi (t):W`$ is continuous and defined by a PDE which involves the gradient of the functional (2)

$$\frac{\psi }{t}=E.$$
(6)
The *discrete steepest descent* technique is computationally cheaper since instead of requiring an integrator for the PDE it constructs a discrete succession of estimates to the minimum, $`\{\psi _{k+1}=\psi _k+\lambda _kE(\psi _k)\}`$, by locally minimizing $`E(\psi _k+\lambda E)`$ with respect to the real parameter $`\lambda `$.

In this paper we will only deal with descent techniques. The first and most important reason is that we will work with the definition of $`E(\psi _k)`$ trying to improve its convergence. As it was already shown in , this work pays off: a good choice of the gradient improves convergence by several orders of magnitude. The second motivation is that by focusing on the gradient our algorithms will be essentially independent on the descent method, which leaves space for further improvement. For instance one might apply these techniques to a nonlinear conjugate gradient method —which dynamically adjusts the search direction, $`d_kE(\psi _k)`$, with an estimate that takes into account the history of the evolution—. Finally, by working with Eq. (6) or its discrete version we avoid the complex nonlinearities that arise in other methods, and we will have a more ample choice of Fourier basis to work with.

### 2.2 Ordinary gradients

It is customary in the literature to work in spaces which are equipped with an $`L^2`$ scalar product and its corresponding norm

$`\psi ,\varphi _{L^2}`$ $``$ $`{\displaystyle \overline{\psi }\varphi },`$ (7)
$`\psi _{L^2}^2`$ $``$ $`{\displaystyle |\psi |^2}.`$ (8)
If one does so and works with Eq. (2) then the formal definition of the gradient is Lagrange’s one

$$E(\psi )=\frac{E}{\overline{\psi }}\frac{E}{\left(\overline{\psi }\right)}.$$
(9)
We will refer to this definition as the “ordinary” gradient, to distinguish it from the different definitions that we will derive below.

### 2.3 Sobolev gradients

Following the ideas from  we will move our problem to a different space which is the Sobolev space of functions, such that $`\psi `$ and its derivatives, $`\psi `$, have a well defined $`L^2`$-norm:

$$^1\{\psi /\psi ,\psi L^2\}.$$
(10)
This Sobolev space will also be equiped with a scalar product and a norm

$`\psi ,\varphi `$ $``$ $`{\displaystyle \left[\overline{\psi }(𝐱)\varphi (𝐱)+\overline{\psi }(𝐱)\varphi (𝐱)\right]d^nx},`$ (11)
$`\psi ^2`$ $``$ $`{\displaystyle \left[|\psi |^2+|\psi |^2\right]d^nx}.`$ (12)
To obtain a new explicit expression for the gradient of the functional in the Sobolev space we will follow a less rigorous derivation than in . Performing a first order expansion of our functional around a trial state $`\psi `$ we obtain

$$E(\psi +\epsilon \delta )=E(\psi )+\epsilon \overline{(\delta ,\delta )}\left(\begin{array}{c}\frac{E}{\overline{\psi }}\\ \frac{E}{\psi }\end{array}\right)+c.c.+O(\epsilon ^2).$$
(13)
We have to turn this expression into something like Eq. (2). This means that we have to find some $`\varphi `$ such that

$$\left[\overline{\delta }\frac{E}{\overline{\psi }}+\overline{\delta }\frac{E}{(\overline{\psi })}\right]+c.c=\left[\overline{\delta }\varphi +\overline{\delta }\varphi \right]+c.c.$$
(14)
If we integrate by parts and impose that this equality be satisfied for all perturbations, $`\delta `$, the problem has a formal solution which is given by a Lagrange equation

$$\left(1\mathrm{}\right)\varphi =\frac{E}{\overline{\psi }}\frac{E}{(\overline{\psi })}.$$
(15)
In consequence, our formal expression for the Sobolev gradient of $`E(\psi )`$ finally reads

$$_SE\left(1\mathrm{}\right)^1E.$$
(16)
Here $`_SE`$ stands for the Sobolev gradient, $`E`$ is the ordinary one, and $`(1\mathrm{})^1`$ represents the inverse of a linear and strictly positive definite operator.

## 3 Sobolev gradients on Discrete Fourier spaces

In the rest of the paper we will work with functions which are defined over a $`d`$-dimensional rectangular volume $`\mathrm{\Omega }\{x\mathrm{\Pi }_i[a_i,b_i]\}`$ with side lengths given by $`L_i=b_ia_i`$. We customarily define an orthogonal set of basis functions, $`\varphi _n=e^{ik_nx}`$, over $`\mathrm{\Omega }`$ where $`\left\{k_n=2\pi (\frac{n_1}{L_1},\mathrm{},\frac{n_d}{L_d}),n_iZ\right\}`$

It is well-known that it is possible to expand any continuous function $`f(x)`$ with periodic boundary conditions using this basis

$$f(x)=\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}\widehat{f}_n\varphi _n(x)$$
(17)
where

$$\widehat{f}_n=\frac{1}{V}_\mathrm{\Omega }\overline{\varphi }_n(x)f(x).$$
(18)
In the previous formula $`V=\mathrm{\Pi }_iL_i`$ is the volume of $`\mathrm{\Omega }`$ and arises because of the lack of normalization of the basis functions, a common practice which saves some computation time.

To discretized the problem we will work within the set of functions sampled over a set of evenly spaced points from $`\mathrm{\Omega }`$, $`\{x_n=(n_1h_1,\mathrm{},n_dh_d),n_i=0,\mathrm{},N_i1\}`$. Here $`h_i`$ represents the spacing along the $`i`$-th dimension, $`n`$ is a vector of non–negative integers and we denote a sampled function by an index, as in $`f_nf(x_n)`$.

Due to this discretization our previous Fourier basis is now redundant. It can be reduced to a finite subset of functions which may represent any sampled function. These functions are given by $`\left\{k_n=2\pi (\frac{n_1}{L_1},\mathrm{},\frac{n_d}{L_d}),n_i=M_i+1,\mathrm{},M_i\right\},`$ where $`M_i=[N_i/2]`$ is the integer quotient of $`N_i`$ divided by two. In this basis a sampled function is given by an expansion which reads

$$f_m=\underset{n}{}\widehat{f}_n\varphi _n(x_m),$$
(19)
and which makes use of the same coefficients as Eq. (18).

The advantage of the finite Fourier basis over other approaches is that it provides an approximant of any function whose error is of order $`O(L_i/N_i)^p`$ where $`p`$ is the maximum differentiability of the sampled function. Furthermore, there is a numerically efficient method known as Fast Fourier Transform which allows one to compute the Fourier coefficients up from the sampled function, $`f_n\widehat{f}_m`$, and vice versa.

Solving Eq. (16) numerically in a discrete Fourier basis is simple. Let us say that we have computed the ordinary gradient and that its sampled version has some Fourier coefficients

$$E(𝐱_n)=\widehat{e}_m\varphi _m(x_n).$$
(20)
and let us assume that there exists a certain solution to Eq. (16) and that it has another discrete Fourier expansion

$$_sE(𝐱_n)=\widehat{s}_m\varphi _m(x_n).$$
(21)
Then by virtue of Eq. (16)

$$\widehat{s}_m=\frac{\widehat{e}_m}{1+k_m^2}.$$
(22)
That is, in the sampled space the Sobolev gradient represents a preconditioning of the ordinary gradient, such that the most oscillating modes are more attenuated. Furthermore, due to this very simple expression computing the Sobolev preconditioning is computationally cheap and involves only minor changes to existing computer codes based on Fourier transforms.

## 4 Applications to Quantum Mechanics

### 4.1 The problem

In this section we will apply the Sobolev gradients technique to a timely problem from Quantum Physics. The system that we will study is a dilute gas of bosonic atoms which are cooled down to ultralow temperatures at which their dynamics become synchronized. When the temperature is low enough, the gas or “condensate” may be described using a single complex wave function, $`\psi (𝐱,t)`$ which is ruled by the so called Gross-Pitaevskii equation, a type of Nonlinear Schrödinger equation that for the Bose gas in a rotating trap with adimensional units reads

$$i_t\psi (𝐱,t)=\left[\frac{1}{2}\mathrm{}+V(𝐱)+g|\psi (𝐱,t)|^2\mathrm{\Omega }L_z\right]\psi (𝐱,t).$$
(23)
Here $`g^+,\mathrm{\Omega }`$, $`L_z=i\left(x_1_2x_2_1\right)`$ is an hermitian operator whose expected value $`L_z=\overline{\psi }L_z\psi ,`$ represents the *angular momentum* of the condensate along an axis of the trap.

There is conserved quantity associated to equation (23) which is called the *energy functional* of the condensate

$$E(\psi )=\frac{1}{2}\left\{|\psi |^2+\overline{\psi }\left[V(𝐱)+\frac{1}{2}g|\psi |^2\mathrm{\Omega }L_z\right]\psi \right\}d^nx.$$
(24)
Our objective in this part of the work will be to *find the solutions which are the minima of the energy subject to a restriction of the $`L^2`$-norm*

$$|\psi |^2N.$$
(25)
The particular value of $`N`$ is imposed by the experimental conditions. Furthermore, without this restriction the absolute minimum of the energy is always reached at the trivial solution $`\psi =0`$.

In Quantum Mechanics the variational formulation of the problem is traditionally converted into a Lagrange equation (4), which is nothing but the Gross-Pitaveskii equation for so called stationary states. In short the word “stationary” refers to solutions of the type

$$\psi (𝐱,t)=e^{i\mu t}\varphi _\mu (𝐱).$$
(26)
The pair $`\{\mu ,\varphi _\mu (𝐱)\}`$ satisfies a nonlinear eigenvalue problem

$$\mu \varphi _\mu (𝐱)=\left[\frac{1}{2}\mathrm{}+V(𝐱)+g|\varphi _\mu (𝐱)|^2\mathrm{\Omega }L_z\right]\varphi _\mu (𝐱).$$
(27)
Due to the difficulty of solving directly the problem (27) we will try to find a direct solution to the variational problem.

### 4.2 Numerical methods: imaginary time evolution

We can search the minima of (24) using descent techniques modified to account for the restriction on the norm (25). The first way to do this is to use a version of the continuous steepest descent which is known as *imaginary time evolution* \[imaginary-time\]. We can summarize this method with the following set of equations

$`\nu (𝐱,\tau )`$ $`=`$ $`\sqrt{{\displaystyle \frac{N}{\sigma _{L^2}^2}}}\sigma (𝐱,\tau ),`$ (28)
$`{\displaystyle \frac{\sigma }{\tau }}(𝐱,\tau )`$ $`=`$ $`E(\nu ),`$ (29)
$`E(\nu )`$ $`=`$ $`\left[{\displaystyle \frac{1}{2}}\mathrm{}+V(𝐱)+g|\nu |^2\mathrm{\Omega }L_z\right]\nu .`$ (30)
Here we see that $`\nu (𝐱,\tau )`$ evolves continuously maintaining a fixed $`L^2`$–norm $`N`$ and following the direction of decreasing energy given by $`E.`$ Indeed it is easy to show that $`\frac{}{\tau }\left[E(\nu (𝐱,\tau ))\right]0`$. Hence, the limit given by

$$\varphi (𝐱)=\underset{\tau \mathrm{}}{lim}\nu (𝐱,\tau )$$
(31)
is at least a critical point of the energy, if not a minimum111As it is common with these local minimization procedures, it remains the problem that iterations may be trapped in a critical point which is not a minimum. Linear stability analysis may then be applied to check the validity of the solution.

From a practical point of view, the simplest way to find the minimizer using imaginary time evolution is to repeatedly integrate Eq. (29) for a very short time, $`\mathrm{\Delta }t`$, apply Eq. (28) for the newly found $`\sigma (𝐱,t+\mathrm{\Delta }t)`$ and use the new estimate for $`\nu `$ to redefine $`\sigma \nu `$ and repeat the procedure until convergence. This way one avoids the problem that according to Eq.(29) the norm of $`\sigma `$ may grow indefinitely. The same consideration applies to the remaining methods that we will present here.

Although the method from Eq. (29) was derived using an ordinary gradient, nothing prevents us from applying our Sobolev preconditioning and all results should be still valid. If we do so, our new equations are

$`\nu (𝐱,\tau )`$ $`=`$ $`\sqrt{{\displaystyle \frac{N}{\sigma _{L^2}^2}}}\sigma (𝐱,\tau ),`$ (32)
$`{\displaystyle \frac{\sigma }{\tau }}(𝐱,\tau )`$ $`=`$ $`_SE(\nu ),`$ (33)
$`_SE(\nu )`$ $`=`$ $`\left(1\mathrm{}\right)^1\left[{\displaystyle \frac{1}{2}}\mathrm{}+V(𝐱)+g|\nu |^2\mathrm{\Omega }L_z\right]\nu ,`$ (34)
and the critical point is still given by Eq. (31).

### 4.3 Numerical methods: minimization of the free energy

While the imaginary time evolution is easy to understand and to implement, the fact that it performs the descent over a path of functions with a certain norm makes it too restricted and sometimes too slow. A different approach is to define a new functional called the *free energy* with a Lagrange multiplier that takes care of the restriction on the norm. In Quantum Mechanics this free energy is usually defined as

$$F_{QM}(\psi )=E(\psi )\mu N(\psi ),$$
(35)
because it preserves the linear form of the equations. Since our equations are already nonlinear we define a free energy functional more conveniently as

$$F(\psi )=E(\psi )+\frac{1}{2}\left(N(\psi )\lambda \right)^2.$$
(36)
First and most important it is not difficult to show that any absolute or relative minimum of $`F(\psi )`$ is also a minimum of $`E(\psi )`$ subject to Eq. (25) with a nonlinear eigenvalue given by $`\mu =N(\psi )\lambda `$.

Secondly, as we will prove in Appendix A, $`F(\psi )`$ must have at least one finite norm minimum, something which cannot be easily assured for $`F_{QM}`$.

The practical advantage of our new functional consists in that fixing $`\lambda `$ and $`\mathrm{\Omega }`$ we can perform a continuous descent over the whole domain of $`F(\psi )`$ without renormalizing the solution on each iteration —i.e., the search space is larger. The new equation that we must integrate is thus

$$\frac{\nu }{\tau }(𝐱,\tau )=F(\nu )=\left[\frac{1}{2}\mathrm{}+V(𝐱)+g|\nu |^2\mathrm{\Omega }L_z\right]\nu .$$
(37)
Using a Sobolev preconditioning Eq. (37) becomes

$$\frac{\nu }{\tau }(𝐱,\tau )=_SF(\nu )=\left(1\mathrm{}\right)^1\left[\frac{1}{2}\mathrm{}+V(𝐱)+g|\nu |^2\mathrm{\Omega }L_z\right]\nu .$$
(38)
In both cases the actual solution is still given by the limit of Eq. (31).

### 4.4 Numerical results

Up to this point we have shown four different numerical methods, two of them incorporating a Sobolev preconditioning \[Eqs. (33) and (38)\] and two without it \[Eqs. (29) and (37)\]. We have compared the efficiencies of these methods for several test situations. The details of our study are as follows.

All methods have been implemented using the discrete Fourier basis mentioned above. This applies both to the calculation of derivatives and to the application of the Sobolev preconditioning. We restricted our simulations to two-dimensional problems with radially symmetric potential, $`V(𝐱)=\frac{1}{2}|𝐱|^2`$, over a grid of 128$`\times `$128 points. Practical experience with more complex problems shows that our results generalize to higher dimensionality and denser grids.

Due to the requirements of the method, both variants of imaginary time evolution \[Eq. (29) and (33)\] are integrated using a Runge-Kutta-Fehlberg, where the tolerance is adapted as the solution converges to its target. On the other hand, instead of performing a continuous descent for both variants of the free energy descent \[Eq. (37) and (38)\], we found it more convenient and faster to perform a discrete steepest descent.

All programs have been implemented using the tensor-algebra environment called Yorick , an interpreted environment which is capable of fast numerical computations and which can be equipped with the FFTW library . Execution times are given for a Digital Personal workstation 500au, but the programs run equally well on modest personal computers with Pentium-II processors and less than 64 Mb of memory.

As test case we have considered three situations. Case A is the simplest one and corresponds to a stationary trap $`\mathrm{\Omega }=0`$, an intense nonlinearity $`g=100`$ and starting the minimization with a radially symmetric Gaussian of unit width, that is, $`\psi _0e^{|𝐱|^2}`$, a shape which is similar to that of the true ground state.

Cases B and C involve a rotating condensate with $`\mathrm{\Omega }=0.6`$ and $`g=100`$. Under these conditions the functional achieves both an absolute minimum, which is the same that found in case A, and a local minimum. The local minimum is a solution of vortex type, a topological defect, whose behavior near zero is $`\psi (x_1+ix_2)/|𝐱|^2`$.

For this reason we designed two test cases with different initial conditions. Case B starts with a Gaussian profile with a centered vortex or $`\psi _0|𝐱|e^{|𝐱|^2}(x_1+ix_2)/|𝐱|^2.`$ In this case all methods are trapped on the local minimum with the centered vortex.

Finally, a third set of simulations, case C, uses the same parameters $`\{\mathrm{\Omega }=0.6,g=100\}`$ but starts from a nonsymmetric initial state, $`\psi _0|𝐱|e^{|𝐱|^2}((x_1y_1)+i(x_2y_1))/|𝐱𝐲|^2`$ which is close to the vortex solution but belongs to the basin of attraction of the ground state. Here all minimization methods require more computational work since they must find a path out of the local minimum, and it is precisely in this case where the differences between methods are best shown.

In Table 1 we summarize the results of the simulations. In case A it is apparent that the Sobolev preconditioning has a positive influence over convergence, with an astonishing result of 55 steps for the steepest descent with free energy. A similar behavior is found in case B.

In case C, the Sobolev preconditioning enhances convergence speed by two orders of magnitude. An intuitive explanation of why the steepest descent with a Sobolev gradient takes less steps to converge will be discussed in detail in Appendix B.

## 5 Applications to Nonlinear Optics

### 5.1 The model

In this section we will consider a model for a pair of incoherently interacting light beams. To be precise we will study the light field of each beam, $`u(𝐱,t)`$ and $`w(𝐱,t)`$, propagating through a weakly nonlinear saturable optical medium. This system may be modeled by the Cauchy problem

$`i_tu`$ $`=`$ $`\mathrm{}u+{\displaystyle \frac{u}{1+\kappa (|u|^2+|w|^2)}},`$ (39)
$`i_tw`$ $`=`$ $`\mathrm{}w+{\displaystyle \frac{w}{1+\kappa (|u|^2+|w|^2)}}.`$ (40)
for the complex functions $`u,w:^2\times ^+`$, which vanish at infinity and satisfy the initial data $`u(𝐱,0)=u_0(𝐱)`$ and $`w(𝐱,0)=w_0(𝐱)`$. Here $`\kappa ^+`$, $`\mathrm{}=\frac{^2}{x^2}+\frac{^2}{y^2}`$ is the Laplacian operator which accounts for the diffraction of light and the nonlinear term $`(1+|u|^2+|w|^2)^1`$ models the saturable interaction among the beams.

Let us define the two component vector

$$\stackrel{~}{U}(𝐱,t)=\left(\begin{array}{c}u(𝐱,t)\\ w(𝐱,t)\end{array}\right),$$
(41)
then, the energy functional for the system reads

$$E(\stackrel{~}{U})=\left[\stackrel{~}{U}^{}\mathrm{}\stackrel{~}{U}+G(|\stackrel{~}{U}|^2)\right]d^nx.$$
(42)
where $`G(\rho )=\frac{1}{\kappa ^2}\left(\mathrm{ln}\left(1+\kappa \rho \right)\kappa \rho \right)`$. The analysis of this section may be generalized to more general nonlinearities with only minor changes to $`G(\rho )`$.

### 5.2 Stationary solutions

We are interested on stationary solutions, which are of the form

$$\stackrel{~}{U}(𝐱,t)=\left(\begin{array}{cc}e^{i\mu _ut}& 0\\ 0& e^{i\mu _wt}\end{array}\right)U(𝐱)=e^{iMt}U(𝐱).$$
(43)
The equations for the stationary solutions now are of elliptic type

$$MU=\mathrm{}U+G^{}(|U|^2)U.$$
(44)
with zero Dirichlet boundary conditions at infinity. Again this formulation poses a nonlinear eigenvalue problem for the pair $`\{M,U\}`$.

The stationary solutions are critical points of the energy functional subject to a constraint on the $`L^2`$-norm of each component. That is, defining

$`N_u`$ $`=`$ $`{\displaystyle |u|^2d^nx},`$ (45)
$`N_w`$ $`=`$ $`{\displaystyle |w|^2d^nx},`$ (46)
the first order variation of the energy around $`U_0`$ for fixed $`N_u`$ and $`N_w`$ must be zero

$$\frac{\delta E}{\delta U}|_{N_u,N_w}=0.$$
(47)
The particular fixed values of $`\{N_u,N_w\}`$ represent the total intensity of each beam of light.

### 5.3 Ground states

In principle there are many different stationary solutions of Eq. (44). However, if we focus on the ground states (stationary solutions of minimal energy) we can apply some of the methods that we mentioned in §4 with minor modifications to account for the higher dimensionality of the problem. For instance, one may define a free energy functional

$$F(U)=E(U)+\frac{1}{2}\left(N_u\lambda _u\right)^2+\frac{1}{2}\left(N_w\lambda _w\right)^2,$$
(48)
and minimize it using a discrete steepest descent with Sobolev preconditioning.

By performing this minimization using different parameters $`\{\lambda _u,\lambda _w\}`$ one obtains nodeless localized solutions for $`u`$ and $`w`$ as shown in Fig. 2. The precise values of $`\lambda _u`$ and $`\lambda _w`$ determine the norm of each component.

Let us remark that, up to a global factor, the ground state has the same shape in both envelopes. That is,

$`u_0(𝐱)`$ $`=`$ $`N\sqrt{\xi }\rho (|𝐱|,N),`$ (49)
$`w_0(𝐱)`$ $`=`$ $`N\sqrt{1\xi }\rho (|𝐱|,N),\xi [0,1].`$
The common shape $`\rho `$ depends on the total intensity $`N=N_u+N_w`$, and it is the one that fixes the values of $`\mu _u`$ and $`\mu _w`$. For this reason if we had chosen the traditional definition of the free energy $`F_{QM}(U)=E(U)+\mu _uN_u+\mu _wN_w`$ we would have found an infinite degeneracy which prevents convergence to the desired values of $`\{N_u,N_w\}`$. This problem is removed by the use of our nonlinear Lagrange multipliers and the $`\{\lambda _u,\lambda _w\}`$ parameters.

### 5.4 Excited states

In the field of guided light waves there has been a great interest on the properties of solutions of Eq. (44) which are not ground states, the so called excited states. Minimization methods based on the energy functional \[See §4\] cannot be applied to this task since excited states are not necessarily local minima of the energy. Instead, a variational principle somehow equivalent to Eq. (44) must be defined.

Let us rewrite Eq. (44) as the application of a nonlinear operator

$$f(U)\left[\mathrm{}M+I(|U|^2)\right]U=0,$$
(50)
define the error functional

$$F(U)f(U)^{}f(U)0$$
(51)
and take our variational principle to be *“find $`U_0`$ such that $`F`$ reaches an absolute minimum $`F(U_0)=0`$.”*

With this principle and a Sobolev preconditioning our descent technique becomes

$$\frac{\nu (𝐱,\tau )}{\tau }=_SF(\nu )=\left(1\mathrm{}\right)^1F,$$
(52)
and it may be proven that $`U_0(𝐱)=lim_\tau \mathrm{}\nu (𝐱,\tau )`$. The ordinary gradient $`F`$ reads

$`G`$ $``$ $`(\mathrm{}+I(|U|^2))U,`$ (53)
$`F`$ $`=`$ $`(\mathrm{}+I(|U|^2))G+I^{}(|U|^2)\left(U^{}G+G^{}U\right)U,`$ (54)
Its calculation using a discrete Fourier basis is straightforward.

There are several advantages of this approach. The first one is that $`F(U)`$ makes no distinction between ground and excited states: any stationary state with the right eigenvalues $`\mu _i=M_{ii}`$ is a local minimum of this new functional. The second one is that we do not need to add any Lagrange multipliers to $`F(U)`$ since they are already present in the $`M`$ operators. Finally we expect that the number of minima of $`F(U)`$ be discrete and separated so as to avoid problems with descent methods being trapped on critical points that are not minima. Indeed, it is very easy to know when this accidental trapping happens since any absolute minimum of $`F`$ must also be a zero of it $`F(U_0)=0`$.

A remarkable feature of the method that we have outlined above is that to distinguish among the different excited states we have to change both the ad-hoc eigenvalues, $`\{\mu _u,\mu _w\}`$, and the initial data of the minimization method.

We have concentrated on two types of excited states of particular interest for applications: the first one is called *vortex vector soliton* , its features being summarized in Fig. 3(a-b). Choosing different initial data for the minimization process we obtain an second type of excited states which are called *dipole-mode vector solitons* . An example of these asymmetric solutions is depicted in Fig. 3(d). Depending on the parameters $`\mu _{u,w}`$ we find different norms of the solution. The fact that our method allows the finding of these nonsymmetric stationary states is very interesting since conventional approaches to the problem have severe difficulties. In fact, the application of the Sobolev preconditioning not only greatly enhances the convergence but it is necessary to get convergence.

Nevertheless the method works equally well for more complicated stationary solutions. In Fig. 4 we show several exotic solutions that are obtained by solving Eq. (52) with different initial conditions. All of those states are dynamically unstable and could have not been found with traditional minimization methods.

### 5.5 Performance and grid refinement

In this subsection we want to analyze the performance of the method from §5.4 and to introduce a multigrid–like technique to improve convergence rates while looking for more accurate solutions.

The idea of the multigrid technique is to use solutions from coarse-grain grids to calculate better approximations on finer grids . Roughly the algorithm consists in setting a coarse-grain initial data, solving the equation or the variational principle with this initial data until the error is small enough, interpolating the result over a finer grid and iterating using this new grid until both the error and the spatial discretization of the solution are the ones we desire.

Since we are already working with Fourier modes over discrete grids the logical choice for our algorithm is indeed Fourier interpolation. The idea is to use the expansion from Eq. (19) *outside of the original grid*, that is,

$$\psi ^{(new)}(x)\underset{n}{}c_ne^{i𝐤_n𝐱},𝐱\mathrm{\Omega }.$$
(55)
This approximation is then discretized over a finer grid, resulting an expansion with a larger number of modes

$$\psi _m^{(new)}=\underset{n}{}c^{(new)}e^{i𝐤_𝐧^{(\mathrm{𝐧𝐞𝐰})}𝐱_m},$$
(56)
which may be used as the starting point for further iteration.

In our case we have used only two different grids, one with 32$`\times `$32 points (*coarse grid*) and a other with 64$`\times `$64 points (the *fine grid*). We have used two different initial data for our minimization method: either a pair of Hermite modes which resemble the desired shape, or the solution of the coarse grid interpolated over the finer grid. As both the evolution of the error in Fig. 5 and the computation times in Table 2 show, interpolation saves a significant amount of time. Indeed, the interpolated solution with only 32$`\times `$32 Fourier modes is already a good approximation for the fine grid, as it shows the small change of the error in Fig. 5

## 6 Conclusions

In this paper we have shown a numerically efficient way to improve the convergence of several minimization methods using the so called Sobolev gradients and applied it to different problems which involve Nonlinear Schrödinger Equations. We have also proven that these gradients represent a preconditioning over the traditional definition of gradients on $`L^2`$. In Appendix B we suggest a generalization of this method to different vector spaces and scalar products.

We have presented two different methods for solving our minimization problems: a traditional one, the *imaginary time evolution*, and a new one, *the minimization of a nonlinear free energy*. Both methods have been shown to be suitable for a Sobolev preconditioning, giving us two new methods that we call *preconditioned imaginary time evolution* and a *preconditioned free energy*.

We have implemented all four methods using a discrete Fourier basis and a Fast Fourier Transform. With these tools we have shown that the Sobolev preconditioning becomes an inexpensive additional step over existing methods. The four resulting solvers have been applied to several realistic problems and in all tests the nonlinear free energy with the Sobolev preconditioning showed the best convergence rates. Indeed the Sobolev preconditioning has an important effect on convergence which can be as good as gaining two orders of magnitude over the traditional techniques. Furthermore, opposite to what happens with finite differences  the preconditioning may be applied without significant computational cost.

We have derived a new method for finding excited states of coupled nonlinear Schrödinger equations. This method introduces a new variational principle which is not based on an energy functional but on finding the zeros of a nonlinear operator which corresponds to the equation to be solved. We have also shown how to improve convergence using Fourier interpolation and a two-grid method as a source for better initial approximations of the iterative method. This approach may be extended to more sophisticated multigrid methods.

## Appendix A Existence of minimizers for the nonlinear free energy functional.

In this appendix we want to prove the existence of minimizers for the free energy (36). Let us write the free energy functional in the following form

$$F(\psi )=\overline{\psi }A_\mathrm{\Omega }\psi d^nx+\frac{g}{2}|\psi |^4d^nx+\frac{1}{2}\left(N\lambda \right)^2,$$
(57)
where $`N=|\psi |^2d^nr`$ and

$$A_\mathrm{\Omega }=\frac{1}{2}\mathrm{}+V(𝐱)\mathrm{\Omega }L_z,\mathrm{\Omega }[0,1)$$
(58)
is a positive hermitian operator $`A_\mathrm{\Omega }\mu _{min}>0.`$

Let us also define the following spaces of functions defined over $`^n`$. We will work in $`L^p`$ spaces

$$L^p=\{\psi :^2/\psi _p:=\left(|\psi (𝐱)|^pd^nx\right)^{1/p}<+\mathrm{}\}.$$
(59)
Also of interest will be the space of functions over which $`A_\mathrm{\Omega }`$ is well defined

$$H_\mathrm{\Omega }=\{\psi :^2/\psi _\mathrm{\Omega }:=\left(\overline{\psi }(𝐱)A_\mathrm{\Omega }\psi (𝐱)\right)^{1/2}<+\mathrm{}\}.$$
(60)
Let us remark that the dimensionality of the space, $`^2`$ is important, since for $`n2`$ it is easy to show that $`H_\mathrm{\Omega }L^4`$.

With the preceding notation we will state the following relevant theorem:

###### Theorem 1.

The free energy functional $`F`$ given by Eq. (57) has at least one absolute minimum in the set given by the inequalities

$$0<\psi _\mathrm{\Omega }^2\lambda \psi _2^2\lambda (\lambda \mu _{min}).$$

###### Proof.

The first step of the proof will be to show that the domain of $`F`$ is indeed the whole space $`H_\mathrm{\Omega }`$. Using the positivity of the $`A_\mathrm{\Omega }`$ operator we show that

$$\psi _2\sqrt{\mu _{min}}\psi _\mathrm{\Omega },$$
(61)
which means that $`H_\mathrm{\Omega }L^2`$. We need Sobolev’s inequality

$$\psi _k\psi _2^{1d}\psi _2^d$$
(62)
where $`d=n/2n/k`$, $`n`$ is the dimensionality of the space and for us $`d=1/2`$. Applying this inequality we obtain the following bound

$$\psi _4\sqrt{\psi _2\psi _2}\mu _{min}^{\frac{1}{4}}\psi _\mathrm{\Omega },$$
(63)
which means that $`H_\mathrm{\Omega }L^4`$. Since $`_2,_4,_\mathrm{\Omega }<+\mathrm{}`$ inside $`H_\mathrm{\Omega }`$ we conclude that $`F`$ as given by Eq. (57) is well defined over the whole space.

$`F`$ is also continuous. To prove it let us rewrite the free energy functional as

$$F(\psi )=\psi _\mathrm{\Omega }^2+\frac{g}{2}\psi _4^4+\frac{1}{2}(\psi _2^2\mu )^2,$$
(64)
and using the bounds (61) and (63) it is easy to show that

$$|F(\psi )F(\xi )|\delta (\epsilon ),\xi \psi \xi _\mathrm{\Omega }ϵ$$
(65)
where the constant $`\delta (\epsilon )`$ is given by $`\epsilon `$ and $`\psi _4`$.

We will also need to show that $`F`$ is coercive

$$\underset{\psi \mathrm{}}{lim}\frac{F(\psi )}{\psi }\alpha >0.$$
(66)
Using Eq. (64) one may show that indeed

$$\frac{F(\psi )}{\psi _\mathrm{\Omega }}\psi _\mathrm{\Omega },$$
(67)
and thus the quotient tends to infinity as the norm grows. The continuity and coercitivity of $`F`$ allow to use theorem 1 (section 1.2) of , which states the existence of at least one minimizer $`\{infF(u):uX\}`$ of $`F`$ provided it is a weakly lower semicontinuous and coercive application $`F:X`$.

So we know that there is at least one minimum, but we do not know where to look for it. Let us now show how the $`\lambda `$ parameter allows us to select different targets for the minimization problem. To do so we define a real one-dimensional function for any given direction $`\psi H_\mathrm{\Omega }`$

$$f(k)F(k\psi ).$$
(68)
This function is nothing but a polynomial over $`k`$

$$f(k)=kNa_\psi +k^2N^2\frac{u_\psi }{2}+\frac{1}{2}\left(kN\lambda \right)^2,$$
(69)
where $`a_\psi =\overline{\psi }A\psi /N`$ and $`u_\psi =\frac{g}{2}|\psi |^4/N^2`$ are constants that depend only on the precise direction $`\psi `$.

By differentiating the polynomial and imposing $`k=0`$ we find that $`f^{}(0)<0`$ at the origin for all possible directions. This means that, as we mentioned above, our Lagrange “multiplier” $`\lambda `$ allows us to avoid the useless solution, $`\psi =0`$.

Furthermore, we can restrict the location of the minimizer to a surface of a certain norm. By differentiating $`f(k)`$ we reach

$$(a_\psi \lambda )+Nk(u_\psi +1)=0.$$
(70)
This equation has a single solution which gives us the norm of minimal energy along the $`\psi `$-direction

$$N_{min}(\psi )=\mathrm{max}\{0,\frac{\lambda a_\psi }{u_\psi +1}\}\lambda \mu _{min},$$
(71)
a value which is bounded above by our Lagrange multiplier $`\lambda `$.

From Eq. (70) it also follows that for any absolute minimum of the functional, $`\psi _{min}`$, the expected value of the $`A_\mathrm{\Omega }`$ operator must be bounded by the $`L^2`$ norm

$$\overline{\psi }_{min}A\psi _{min}=\psi _{min}_\mathrm{\Omega }^2\lambda N(\psi _{min}).$$
(72)
Otherwise the trivial solution $`\psi =0`$ would have less energy than $`\psi _{min}`$. We can thus, instead of working with an unknown surface, delimit the location of the minimum to a set which is given by two inequalities, (71) and (72)

$$W=\{\psi H_\mathrm{\Omega }:\mathrm{\hspace{0.33em}0}<N(\psi )\lambda \mu _{min},\psi _\mathrm{\Omega }^2\lambda N(\psi )\}.$$
(73)
There are several practical consequences of this theorem. First, it states that the problem of finding the minima of $`E(\psi )`$ subject to fixed norm has one solution, i.e., there exists at least one ground state. Second, but equally important, it proves our Lagrange penalizer $`\frac{1}{2}(N\mu )^2`$ to be specially well suited for this problem since it avoids both the useless solution $`\psi =0`$ and those with too large norm. And finally it gives us a bounded in which the minimizer must be. Indeed we can extend this result by proving that we can restrict our search space to a compact superset of $`W`$.

###### Theorem 2.

Any absolute minimum, $`\psi `$, of the functional $`F`$ given by Eq. (57) lays inside the compact set of $`L^2`$

$$\overline{U}=\{\psi L^2:\psi _\mathrm{\Omega }^2\lambda \psi _2^2\lambda (\lambda 1)\}.$$

###### Proof.

For the type of spaces that we work with, a set is compact iff it is closed and we can build an $`\epsilon `$-net for any positive number $`\epsilon `$. The $`\epsilon `$-net is a finite set $`K_\epsilon =\{\nu _1,\mathrm{},\nu _k\}`$ such that for each $`\psi \overline{U}`$ there is an element $`\nu K_ϵ`$ verifying $`\psi \nu _\mathrm{\Omega }<\epsilon `$. Thus compactness is equivalent to the possibility of building an arbitrarily good approximation of our minimizer using a finite but sufficiently large basis of functions.

It is evident that the set $`\overline{U}`$ is closed in the subspace $`H_\mathrm{\Omega }`$ of $`L^2`$. Let $`\{u_n\}\overline{U}`$ be a convergent sequence and let $`u`$ be their limit. Since for each element of the sequence

$$u_n_\mathrm{\Omega }\sqrt{\lambda }u_n_2<\lambda (\lambda 1),$$
(74)
it is also obvious that

$$u_\mathrm{\Omega }\sqrt{\lambda }u_2\lambda (\lambda 1),$$
(75)
also thus the limit belongs to $`\overline{U}`$.

The compactness of the closed set $`\overline{U}`$ essentially follows from the fact that the eigenstates of $`A_\mathrm{\Omega }`$ form a complete basis of $`H_\mathrm{\Omega }`$, and that the eigenvalues of $`A_\mathrm{\Omega }`$ form a monotonously growing unbounded set of positive numbers. These eigenstates are of the form

$$\varphi _{n,l}=P_n^l(|𝐱|)\frac{x_1+ix_2}{|𝐱|^2}e^{|𝐱|^2/2},$$
(76)
where $`P_n^l`$ are Laguerre’s generalized polynomials. And the corresponding eigenvalues are

$$A_\mathrm{\Omega }\varphi _{n,l}=\mu _{n,l}\varphi _{n,l}=(2n+(1\mathrm{\Omega })l+1)\varphi _{n,l},n,l\{0\}.$$
(77)
Let us choose any natural number $`k`$ such that $`k+1>\lambda `$. We can split the whole space as a direct sum $`H_\mathrm{\Omega }=H_0^{k+1}H_{k+1}^{\mathrm{}}`$ where

$$H_j^k=\mathrm{lin}\{\varphi _{n,l}:j2n+l<k\}.$$
(78)
The important point is that since $`H_0^{k+1}`$ is isomorph to $`^m`$ for some natural number $`m`$, and $`\overline{U}_k=\overline{U}H_0^{k+1}`$ lays inside a compact $`m`$-dimensional ball of radius $`\sqrt{\lambda 1}`$, then we can find any $`\epsilon `$-net for $`\overline{U}_k`$. Furthermore, by separating

$$\psi =\psi _a+\psi _b,\psi _a\overline{U}_k,\psi _bH_{k+1}^{\mathrm{}},$$
(79)
and using the definition of $`\overline{U}_k`$ we show that the projection of $`\psi `$ outside of $`\overline{U}_k`$ can be made arbitrarily small

$$\psi _b_2^2\frac{\lambda }{k+1}N.$$
(80)
A direct consequence of this is that for $`k>\lambda N/\epsilon `$, a $`\epsilon `$-net of $`\overline{U}_k`$ is also an $`\epsilon `$-net of $`\overline{U}`$, which proves the compactness.

Finally, by inspecting the eigenvalues of $`A_\mathrm{\Omega }`$ and using Eq. (71-72) it is not difficult to see that the absolute minimum of $`F`$ must lay in $`\overline{U}`$.

## Appendix B Extending the Sobolev gradients

We have shown that redefining the gradient turns out to be a kind of preconditioning over the original choice of the direction of descent. Let us assume that our functional has the following form

$$E(\psi )=\overline{\psi }A\psi +f(|\psi |^2,𝐱),$$
(81)
where $`A`$ is a non–negative hermitian operator that may involve some derivatives. Let us also assume that $`H`$ is a suitable space equipped with the following scalar product

$$\psi ,\varphi =\overline{\psi }\left(1+A\right)\varphi ,$$
(82)
which is indeed a scalar product because $`A|_H0.`$

Let us rewrite the energy functional in the following way

$$E(\psi )=\psi ,\psi +f(|\psi |^2,𝐱)|\psi |^2.$$
(83)
The Sobolev gradient in $`H`$ is

$$_AE=\psi +(1+A)^1\left[_1f\psi \right],$$
(84)
while the so called ordinary gradient is $`E=A\psi +_1f.`$ Hence the preconditioning nature of the method is recovered

$$_AE=\left(1+A\right)^1E.$$
(85)
We can thus think that the new choice of the scalar product aims at making the linear part of the energy functional close to some quadratic form, $`\psi ,B\psi `$, such that the new operator $`B`$ is almost the unity. In our example, indeed, $`B\psi =\psi .`$ We that this preconditioning will both enhance the directions of decreasing energy and also have a smoothing effect on the nonlinear part. However, the problem is complicated and we know of no proof that these arguments be of general applicability.


text:
# Emergence of Bulk CsCl Structure in (CsCl)nCs+ Cluster Ions

## I Introduction

A general goal of cluster physics is to study the emergence of bulk behavior right from the molecular limit, by building clusters of increasing size and following the size evolution of selected properties. From the theoretical point of view, this ambitious plan has been largely impeded because of the slow and nonmonotonic size evolution observed in many properties. The predicted cluster structures are not simply related to the corresponding bulk structures in many cases, which precludes the possibility of a meaningful extrapolation to the bulk limit. Moreover, cluster structure is difficult to determine theoretically due to the huge increase in the number of isomers with size, and experimentally due to the small number of scatterers compared with the bulk case. Nevertheless, recent advances involving ion mobility measurements,  electron diffraction from trapped cluster ions, or the use of photoelectron spectra as a fingerprint of structure have been successful in elucidating the structures of several ionic and covalent clusters.

Abundance patterns obtained from the mass spectra of binary ionic clusters like the alkali halides and alkaline-earth oxides point towards a prompt establishment of bulk rock-salt symmetry. Theoretical calculations have shown, however, that small sodium iodide and lithium halide clusters adopt ground state structures based on the stacking of hexagonal rings.  In the case of alkaline-earth oxide clusters, the large and coordination-dependent values of the oxide polarizabilities favor the formation of structures with a large proportion of ions in surface sites, inducing a delay in the emergence of bulk structural properties. Turning to the alkali halides, bulk CsCl, CsBr, and CsI crystallize in the CsCl-type structure, while both experimental mass spectra and theoretical calculations indicate that small clusters of those materials adopt ground state structures which are fragments of a rock-salt lattice. This implies that there has to be a structural phase transition as the cluster size is increased. Ion mobility measurements performed by Löffler suggest that (CsCl)nCs+ cluster ions with n=32 are specially compact, which might be explained by the presence of isomers with the shape of a perfect three-shell rhombic dodecahedron (that is with bulk CsCl symmetry) in the cluster beam. The electron diffraction experiments performed recently in the group of Parks show that there is a substantial proportion of isomers with bulk CsCl symmetry for the same size.

In this theoretical work we analyce the above mentioned size-induced phase transition in (CsCl)nCs+ cluster ions. We consider only those sizes that correspond to geometrical shell closings for the CsCl-type (perfect rhombic dodecahedra with n=32,87,184,335,552) and rock-salt (perfect cubes with n=13,62,171,364,665) structural series. In doing so, we try to avoid any nonmonotonic size evolution in the calculated properties. In the upper part of Fig. 1 we display the relative number of atoms with a given coordination as a function of N-1/3, where N=2n+1 is the total number of atoms in the cluster. In the lower part we show the number of atoms with nonbulk coordination relative to the total number of surface atoms. For the largest sizes considered the proportion of bulklike atoms is dominant, and within the surface the proportion of face-like atoms is already much larger than those of edge and vertex-like atoms. From those sizes to the bulk, the only meaningful size evolution of these proportions will be a slow approach to zero of the face-like atoms. We thus expect to capture all the physical information relevant to the phase transition by studying this set of clusters and the corresponding bulk phases, which have been studied both with the same theoretical model. In this way inaccuracies related to the use of different methodologies are avoided and a meaningful extrapolation to the bulk limit can be done. In a second part of the work, we explicitely analyze the structures adopted by (CsCl)nCs+ cluster ions in the size range n=31–33, in order to explain the experimental findings of Refs.  and .

The rest of the paper is organized as follows: Section II includes just a brief description of the theoretical model employed in the calculations, as a full account of it has been given in previous publications. In Section III we present and discuss the results of the calculations, and Section IV summarizes the main conclusions.

## II Theoretical model

Cluster energies have been obtained by performing Perturbed Ion (PI) plus polarization calculations. This is a well tested method that describes accurately both bulk and cluster limits. Its theoretical foundation lies in the theory of electronic separability. Very briefly, the cluster wave function is broken into local group functions (ionic in nature in our case) that are optimised in a stepwise procedure. In each iteration, the total energy is minimized with respect to variations of the electron density localized in a given ion, with the electron densities of the other ions kept frozen. In the subsequent iterations each frozen ion assumes the role of nonfrozen ion. When the self-consistent process finishes, the outputs are the total cluster energy and a set of localized wave functions, one for each geometrically nonequivalent ion of the cluster. This procedure leads to a linear scaling of the computational effort with cluster size, which allows the investigation of large clusters with an explicit inclusion of the electronic structure. The cluster binding energy can be decomposed into ionic additive contributions

$$E_{clus}^{bind}=\underset{Rclus}{}(E_{add}^RE_0^R),$$
(1)
being E$`{}_{}{}^{R}{}_{add}{}^{}`$ the contribution of the ion R to the total cluster energy and E$`{}_{}{}^{R}{}_{0}{}^{}`$ the energy of the ion R in vacuo. In this way the contribution of ions with different coordinations to the binding energy can be separately analyzed, which is particularly convenient for our study. Each additive energy can be decomposed in turn as a sum of deformation and interaction terms

$$E_{clus}^{bind}=\underset{Rclus}{}(E_{def}^R+\frac{1}{2}E_{int}^R),$$
(2)
where E$`{}_{}{}^{R}{}_{def}{}^{}`$ is the self-energy of the ion R, measured relative to the vacuum state, and E$`{}_{}{}^{R}{}_{int}{}^{}`$ contains electrostatic, exchange and repulsive overlap energy terms. The polarization contribution to the cluster binding energy is not computed in the actual version of the PI code, as it assumes (for computational simplicity) that the electronic charge distribution of each ion in the cluster is spherically symmetric. Thus, a polarization correction to the PI energy is computed semiempirically as described in Refs.. Bulk polarizabilities are used for both Cs+ and Cl- ions. This is a good approximation for the Cs+ cations. The main effect on the anion polarizabibities when passing from the bulk to a cluster environment is an increase of the polarizabilities of those ions located on the cluster surface, due to the lower average coordination compared to the bulk. However, we have checked that our main conclusions are not affected by an increase in the surface chloride polarizabilities as large as 10–20 %, which are typical values for halides. The short-range induction damping parameters have been obtained through the scaling procedure validated in Ref.. The reliability for cluster calculations of the mixed quantum-mechanical/semiempirical energy model thus obtained has been checked and shown to be high in previous publications.

## III Results

### A The rock-salt to CsCl-type structural transition

Fig. 2 shows the size evolution of the binding energy per ion. First of all, we note that the PI model properly reproduces the stability trend in the bulk, predicting the CsCl structure as the most stable one. This is a tough problem for semiempirical methods, as Pyper has shown that a full account of the coordination number dependence of the self-energy and overlap contributions is necessary to obtain the correct ground state structure. The values of the binding energy, plotted as a function of N-1/3, lie neatly on a straight line. The regression coefficients obtained from a fit are 0.9998 in all cases if we exclude from the fitting the NaCl-type cluster with n=13, which is the smallest one. We have calculated after the fitting procedure the energy of the 5$`\times `$5$`\times `$7 cuboid (also included in Fig. 2), and checked that it lies on the fitted NaCl-type energy curve. This shows that a consideration of perfect cubes (or cuboids) on one hand, and rhombic dodecahedra on the other hand removes the nonmonotonic behavior from the size evolution of the binding energies. Our results predict that the rhombic dodecahedra become definitely more stable after the completion of the fifth shell of atoms, that is for n=184. The four-shell rhombic dodecahedron and the 5$`\times `$5$`\times `$7 cuboid are essentially degenerate, so both of them will contribute to the enhanced abundance observed experimentally for n=87.  We have not found any experimental mass spectrum for values of n as high as 184, but we predict that a new set of magic numbers, corresponding to the closing of rhombic dodecahedral atomic shells, should emerge from this size on. The magic numbers corresponding to the closing of perfect cubic shells will probably not disappear still at that specific size from the mass spectra, because they do not coincide with the CsCl shell closings, and complete cubes can remain more stable than incomplete rhombic dodecahedra until larger values of n are reached. Polarization has little influence on these general results, and only affects the energetic ordering of the two essentially degenerate isomers mentioned above.

Now we turn to an analysis of the physical factors responsible for this transition. In Fig. 3 we show the binding energy per ion, averaged over subsets of ions with a fixed coordination. The contribution of bulklike ions to the binding energy favors always the formation of CsCl-type structures. However, the contribution of facelike ions favors the formation of rock-salt fragments. As soon as the proportion of bulk ions is larger than that of surface atoms, which occurs after the completion of the fifth rhombic dodecahedral atomic shell, fragments of the CsCl-type lattice become more stable. The energy contribution of those ions in edge positions is approximately the same for both structural families except for the smallest clusters; finally, corner atoms favor the CsCl-type structures, but their small relative number results in a very small contribution to the total energy for those sizes where the transition occurs.

Fig. 3 has reduced the structural phase transition in (CsCl)nCs+ cluster ions to an essentially bulk effect. By this we mean that CsCl-type structures become more stable as soon as the proportion of bulklike atoms is dominant. To complete our discussion we have then to address the stability question in the bulk. This is more easily understood by analyzing the reasons why other alkali halides like NaCl or CsF do not crystallize in the CsCl-type lattice. The largest contribution to the binding energy of an ionic crystal is the Madelung energy term E$`{}_{M}{}^{}=\frac{A_M}{R_e}`$, with $`A_M`$ the Madelung constant and $`R_e`$ the equilibrium interionic distance. The Madelung constant of the CsCl-type lattice (1.762675) is larger than that of the rock-salt lattice (1.747565), so were the value of Re the same for both structures, the CsCl-type would always be more stable. We have solved for the electronic structure of NaCl and CsF crystals in the CsCl-type structure at a nonequilibrium value of the interionic distance, chosen in such a way that the Madelung energy term is exactly the same as in the corresponding rock-salt lattice at equilibrium. In the case of NaCl, Eadd(Na+) favors the CsCl-type structure, but Eadd(Cl-) largely favors the rock-salt phase. The main reason is the large anion-anion overlap at that artificial distance, that is the Na+ cation is so small compared to the Cl- anion that eight anions can not be packed efficiently around a cation. In CsF the situation is reversed, and it is the cation-cation overlap that is too large. This demonstrates that the stability situation in the bulk is a purely packing effect: in CsCl, CsBr and CsI, the large value of the cation-anion size ratio allows for an equilibrium interionic approach in the CsCl-type structure close enough as to obtain a Madelung energy term more negative than in the rock-salt phase, without a large overlap interaction between like ions. The same is true for the bulklike ions in the clusters studied, and so when those ions begin to dominate the energetics, the bulklike fragments become more stable.

We have made a prediction above that can be tested experimentally, namely the emergence of a new set of magic numbers from n=184 on. Here we propose another, perhaps more indirect, experimental test. In Fig. 4, the eigenvalues of the 3p orbitals of Cl- (with opposite sign) are plotted as a function of N-1/3. We have a band of eigenvalues for each size because the anions occupy nonequivalent positions in the clusters. As the clusters under study are formed by closed shell ions whose wave functions are strongly localized, it can be assumed that an electron is extracted from a specific localized orbital when the cluster is ionized. This is the lowest bound 3p orbital, which corresponds always to a chloride anion with a low coordination. Thus the dashed lines represent the size evolution of the vertical ionization potential IP (in the Koopmans’ approximation) for both structural families. For the rock-salt series, that size evolution is approximately linear in N-1/3, but for the CsCl-type series it shows a more or less oscillating behavior, which should be detected in experimental measurements of the vertical IP if rhombic dodecahedra actually are the ground state structures from a given size on. We can explain these different electronic behaviors in a very simple way: in the rock-salt clusters the eight corner sites are always occupied by Cs+ cations. The weakest bound electron corresponds always to a Cl- anion with fourfold coordination, namely anyone of those closer to the corner cation sites. On the other hand, rhombic dodecahedra have fourteen corner sites. When the number of atomic shells is even, all these sites are occupied by Cs+ cations, but when that number is odd, eight of them are cationic sites and the other six anionic sites. Thus the nonmonotonic behavior of the vertical IP is due to the different local coordination of the Cl- anion to which the weakest bound electron is attached as the number of atomic shells increases.

### B Structures of (CsCl)nCs+ (n=31–33) and comparison to experiment

We finish our study with an explicit consideration of (CsCl)nCs+ clusters in the size range n=31–33, the range covered in the experiments of Löffler and Parks. Specifically, we have considered the most compact 7$`\times `$3$`\times `$3, 4$`\times `$4$`\times `$4 and 5$`\times `$4$`\times `$3 rock-salt structures, and the three-shell rhombic dodecahedron, with some atoms added or removed from different positions. The binding energies are shown in Table I. The ground state (GS) structure of (CsCl)31Cs+ is a complete 7$`\times `$3$`\times `$3 cuboid. The 4$`\times `$4$`\times `$4 fragment with an anion removed from a corner position is slightly less stable, and the lowest energy rhombic dodecahedral isomer we have obtained has a still lower stability. For n=32, the complete three-shell rhombic dodecahedron becomes the GS isomer. All the different incomplete rock-salt fragments have a smaller binding energy. For n=33, the different rock-salt isomers are essentially degenerate, but the CsCl-type structure is found again at a higher energy. This sequence of GS structures for (CsCl)nCs+ clusters is consistent with the experimental findings. The relative mobility is a local maximum for n=32, as the perfect three-shell rhombic dodecahedron is evidently more compact than the complete 7$`\times `$3$`\times `$3 (CsCl)31Cs+ cuboid or any of the incomplete rock-salt structures obtained for n=33. Also, the energetical ordering of the isomers is consistent with the large proportion of CsCl-type isomers found for n=32 in the electron diffraction experiments.

## IV Summary

We have reported a computational study of the size-induced rock-salt to CsCl-type structural phase transition in (CsCl)nCs+ cluster ions. For this purpose, the Perturbed Ion (PI) method, supplemented with a semiempirical account of polarization effects, has been employed. Only cluster ions with an atomic closed-shell configuration have been considered in order to avoid nonmonotonic behavior in the calculated properties. Moreover, we have employed the same theoretical model to study both cluster and bulk limits, which allows for a meaningful extrapolation strategy. The main result is that rhombic dodecahedral isomers become definetely more stable than rock-salt structures after the completion of the fifth rhombic-dodecahedral atomic shell, that is for a size n=184. Thus, it is predicted that a new set of magic numbers, reflecting the establishment of the new structural symmetry, should emerge from that size on. The size evolution of the vertical ionization potential of the cluster ions should also be a good experimental fingerprint of the transition. In order to explain the nature of the transition, an analysis of the binding energy into ionic components has been performed. The result is quite simple: bulklike ions always prefer to have a CsCl-type environment, even for the smallest cluster sizes (this has been shown to be a purely packing effect), while surface-like atoms prefer to adopt rock-salt structures. The transition occurs as soon as the proportion of bulklike atoms is large enough to dominate the energetics of the whole cluster. One of the possibilities advanced by Parks consistent with his experimental results  is the existence of isomers with mixed symmetry. Our results indicate that the formation of isomers with a CsCl-type core and a rock-salt-type surface could be energetically favored, if the strain accumulated in the bonds at the interface region separating both phases can be conveniently relaxed. This point deserves further investigation.

The structures adopted by (CsCl)nCs+ cluster ions have been more carefully studied in the size range n=31–33, which has been covered in the experimental investigations. Our results are consistent with the experimental findings, and show that the three-shell rhombic dodecahedron is the lowest energy isomer for n=32.

ACKNOWLEDGMENTS
This work has been supported by DGES (Grant PB98-0368) and Junta de Castilla y León (VA70/99). The author is indebted to J. M. López for a careful reading of the manuscript.

Captions of Tables.

Table I Binding energy, in eV/ion, of the different rock-salt and CsCl-type structures for the size range n=31–33.

Captions of Figures.

Figure 1 Size evolution of the number of atoms with a given coordination, relative to the total number of atoms (upper half) or to the total number of surface atoms (lower half). The left half refers to CsCl-type symmetry and the right half to rock-salt symmetry.

Figure 2 Size evolution of the binding energy per ion for both CsCl-type and rock-salt structural families, with (lower half) and without (upper half) the inclusion of polarization corrections. The value of N-1/3 at the transition point has been indicated with an arrow.

Figure 3 Size evolution of the binding energy contributions from ions with different coordinations. Full circles represent ions in the CsCl-type structures and squares represent ions in the rock-salt structures.

Figure 4 Size evolution of the 3p orbital eigenvalues of chloride anions. The dashed line represents the variation of the vertical ionization potential in the Koopmans’ approximation with size.

\[