text:
# Multiresolution analysis in statistical mechanics. I. Using wavelets to calculate thermodynamic properties

## I Introduction

Spin models are popular tools for theoretical calculations and for numerical simulations, as their universality classes allow a huge range of different systems—as varied as binary metal alloys, surface adsorption, and neural networks—to be modeled simultaneously. For example, even the “trivial” one-dimensional Ising model can be used to model the helix-coil transition in biopolymers; the deep connection between magnetic models and polymer chains allow us to predict scaling behavior and other properties across an even wider range of materials. Lattice models are still widely used in modeling the thermodynamics of complex systems, because their regular structure simplifies the type and nature of interactions among components of the system. Moreover, the difficulty in obtaining analytical solutions of lattice systems, and the relative ease of computational simulations thereof, make them ideal test cases for new simulation algorithms.

Although simulations of lattice models are relatively straightforward to implement, they share the same drawbacks as off-lattice models. The chief drawback is that as the number of particles grows large, the time required to sample the system accurately increases rapidly. A popular approach for addressing this problem is to coarse-grain the system: that is, we “rescale” the problem by increasing the basic size of a simulation element. For example, we might coarse-grain an atomic representation of a polymer chain into a “united atom” model, where a chain molecule is treated as if it consisted only of the backbone. More creative approaches redefine the problem to be addressed: for example, Mattice and coworkers have produced a method which maps a polymer chain atomistically onto a high-coordination lattice; this lattice is then used as the basis for a Monte Carlo simulation; the resulting configuration is then used to map back to continuous space to provide “evolution in time.”

This paper illustrates the use of the wavelet transform as a mathematical basis for performing thermodynamic computations of lattice models. The wavelet transform is an important tool in multiresolution analysis, which analyzes a system simultaneously at several length or frequency scales selected to reflect the actual physical processes underlying the observed behavior as closely as possible. The wavelet transform possesses a number of convenient properties, including orthogonality, compactness, and reconstruction; we will make these concepts more precise in Section II B below. The orthogonal nature of most wavelet constructions makes them a logical choice for use in ab initio density-matrix quantum chemistry computations, in which the selection of an accurate basis set is crucial to the convergence and efficiency of the calculations. Wavelet decompositions have been applied principally in electrical engineering, particularly in the field of signal processing. In this context, white noise and Markov processes have been studied using multiscale methods. To date, however, wavelet analysis does not seem to have been extensively applied to models in statistical mechanics. Huang uses wavelet analysis to observe the statistical distribution of multiplicity fluctuations in a lattice gas, while Gamero et al. employ wavelets to introduce their notion of multiresolution entropy, although their primary goal is dynamic signal analysis rather than statistical mechanics simulations, while O’Carroll attempts to establish a theoretical foundation connecting wavelets to the block renormalization group. A more in-depth review of the connection between wavelets and renormalization theory is provided in a recent monograph by Battle.

## II Wavelet transform fundamentals

### A The conceptual picture

The wavelet transform is a hierarchical method for decomposing a data set into averages and differences. Like the Fourier transform, it can be used to provide a decomposition in both real space and reciprocal space ($`𝐤`$-space), or time space and frequency space. Unlike the Fourier transform, however, it is capable of providing simultaneously localized transformations in both real and reciprocal space. A function localized in position space, such as a finite impulse function, cannot be represented by a few terms of its Fourier series: many terms are required before good convergence is achieved. By contrast, in wavelet space, this same function can be almost completely described by just a handful of wavelet coefficients. Although the first wavelet was discovered almost a century ago by Haar, they have become an important computational technique only in the last decade, following the work of Mallat, Daubechies, and others.

The wavelet transform, like any other transform, takes a mathematical object and transforms it into another: we can represent its action by writing

$$\stackrel{~}{u}=𝒲\left[u\right];$$
(1)
the specific form of $`𝒲`$ depends both on the type of wavelet we have selected, and the object $`u`$ which we wish to transform. All versions of the wavelet transform $`𝒲`$, however, are derived from the same source: a set of coefficients which define the transform. If $`u`$ is a discrete data set, such as a signal sampled at regular intervals, then $`𝒲`$ is usually represented as a matrix; while if $`u`$ represents a continuous data set, such as the same signal measured at all times, then $`𝒲`$ acts as an integral operator. While the matrix form of $`𝒲`$ is often called a “filter bank” and the integral form a “wavelet transform,” we will not distinguish between them in what follows, as the theory developed here for discrete lattices and filter banks should carry over to continuous systems and wavelet transforms essentially unchanged.

Similar to the Fourier transform, the wavelet transform decomposes the object $`x`$ into two separate components, as two different functions, a scaling function $`\varphi `$ and a wavelet function $`\psi `$, both operate on $`x`$. However, the two functions separate its components not into cosines and sines, but into averages and differences, with a “wavelength” equal to the “window” over which the scaling and wavelet functions are nonzero. In another important distinction, the wavelet transform is recursive, so that it can be applied in succession to any set of averages which is produced using that wavelet transform, to produce another level of averages and another level of details.

We shall now make the above concepts more mathematically precise. Let us define $`u`$ to be a discrete set of samples $`u=(u\left(1\right),u\left(2\right),\mathrm{},u\left(n\right))`$. Then applying the scaling and wavelet functions $`\varphi `$ and $`\psi `$ to $`u`$ create a set of averages $`s\left(i\right)`$ and a set of differences $`\delta \left(i\right)`$:

$`s\left(i\right)`$ $`={\displaystyle \underset{k=0}{\overset{r1}{}}}\varphi \left(k\right)u\left(i+k\right),`$ (2)
$`\delta \left(i\right)`$ $`={\displaystyle \underset{k=0}{\overset{r1}{}}}\psi \left(k\right)u\left(i+k\right),`$ (3)
where $`r`$ is a finite integer which defines the length scale, often referred to as the “size of the support,” over which $`\varphi `$ and $`\psi `$ are nonzero. The index $`i`$ runs from $`1`$ to $`n;`$ generally the data set is padded with zeros to ensure that all sums in (2) and (3) are well-defined, although periodicity is sometimes used instead.  The coefficients $`\varphi \left(k\right)`$ and $`\psi \left(k\right)`$ in (2) and (3) are related, and are central in controlling the features of the wavelet transform. Note the wavelet transform is inherently redundant: for every sample $`u\left(i\right)`$ in the original set $`u`$, we now have two values, a local average $`s\left(i\right)`$ and a local difference $`\delta \left(i\right)`$. Since the new data are simply linear combinations of the original values, it is superfluous to retain both sets; at the same time, it is obvious that we cannot simply discard one set of data and recover all the original information using only the other data set. Instead, we choose to keep only the odd-numbered $`s\left(i\right)`$’s and $`\delta \left(i\right)`$’s, eliminating the even-numbered samples; this process is called downsampling. Downsampling removes half of the $`s\left(i\right)`$ and half of the $`\delta \left(i\right)`$, regardless of the length $`r`$ of the wavelet. Now we are left with $`n`$ data: $`s\left(1\right),s\left(3\right),\mathrm{},s\left(n1\right)`$ and $`\delta \left(1\right),\delta \left(3\right),\mathrm{},\delta \left(n1\right)`$. These $`n`$ data points can be stored as the level-one wavelet transform $`\stackrel{~}{u}`$ of $`u`$, by assigning $`s\left(1\right),s\left(3\right),\mathrm{},s\left(n1\right)`$ to $`\stackrel{~}{u}^{\left(1\right)}\left(1\right),\stackrel{~}{u}^{\left(1\right)}\left(2\right),\mathrm{},\stackrel{~}{u}^{\left(1\right)}\left(n/2\right)`$, and the corresponding $`\delta \left(i\right)`$’s as $`\stackrel{~}{u}^{\left(1\right)}\left(n/2+1\right),\mathrm{},\stackrel{~}{u}^{\left(1\right)}\left(n\right)`$. \[The superscript $`\left(1\right)`$ denotes that the wavelet transform has been applied once to this data set.\] We can either stop at this level of description, or continue by further decomposing the averages: then the new object $`u^{\left(1\right)}`$ to be transformed is $`u^{\left(1\right)}\left(1\right)=\stackrel{~}{u}\left(1\right),\mathrm{},u^{\left(1\right)}\left(n/2\right)=\stackrel{~}{u}\left(n/2\right)`$, and so on. Note that although $`\stackrel{~}{u}^{\left(1\right)}`$ contains the averages $`s\left(i\right)`$’s and the differences $`\delta \left(i\right)`$’s obtained in the previous step, successive transforms only apply to averages obtained in the previous step. This process can be repeated until we have reduced our set of averages to a single point; no further averaging is possible. We assume henceforth that the data set $`u^{\left(k\right)}`$ has been sufficiently downsampled to retain only the minimum data set required.

### B Properties and examples of wavelet and scaling functions

Until this point, we have not introduced any specific wavelet or scaling functions. Before we do so, we note that the choice of a wavelet transform to apply to a given system usually hinges on the desired properties which one wishes to include in the transformed data. Three principal properties are almost universally required for filter banks and wavelet families:

1. Perfect reconstruction: No data is distorted by performing analysis followed by synthesis, so that the only permissible change is a delay in recovery of the original sample.
2. Orthogonality: Wavelets computed at different length scales or at different spatial locations are mutually orthogonal; thus fluctuations in the system are localized at the scales where they are most relevant.
3. Compact support: Properly designed wavelets are identically zero except for a finite interval, which means that exact results can be obtained using only a finite number of terms.

Other properties, such as orthonormality, symmetry in the functional form of the wavelet or a certain number of vanishing moments, can be taken into account when constructing the wavelet transform.

The two most commonly encountered selections are the Haar and Daubechies wavelets, named after their respective discoverers. The Haar pair is the oldest and simplest set of wavelets: the coefficients of the scaling function are $`\varphi =(\varphi \left(0\right),\varphi \left(1\right))=\frac{1}{\sqrt{2}}(1,1)`$, while the coefficients of the wavelet function are $`\psi =(\psi \left(0\right),\psi \left(1\right))=\frac{1}{\sqrt{2}}(1,1)`$. No other wavelet can be described with two points, and therefore no other wavelet has a support as compact as the Haar wavelet. The scaling function $`\varphi `$ simply averages the values stored at neighboring points, while $`\psi `$ finds the difference between those values; the extra factor of $`\sqrt{2}`$ is incorporated to ensure orthonormality between overlapping $`\varphi \left(k\right)`$ and $`\psi \left(k\right)`$. A simple example of the action of the Haar wavelet is shown in Fig. 1.

The Daubechies wavelets are a family of orthonormal functions whose construction was explicitly designed to have orthogonality as well as vanishing higher-order moments. Daubechies was able to show that the Haar wavelet is in fact the “first” member of the Daubechies family; that is, the Haar wavelet is the Daubechies wavelet with the shortest support. The second such member has four terms in its definition: the scaling function is defined by $`\varphi =(\varphi \left(0\right),\varphi \left(1\right),\varphi \left(2\right),\varphi \left(3\right))=\frac{1}{4\sqrt{2}}(1+\sqrt{3},3+\sqrt{3},3\sqrt{3},1\sqrt{3}).`$ The wavelet function reverses the order of the coefficients and inverts the sign of every other component, which allows the orthonormality properties to be satisfied: $`\psi =(\psi \left(0\right),\psi \left(1\right),\psi \left(2\right),\psi \left(3\right))=(\varphi \left(3\right),\varphi \left(2\right),\varphi \left(1\right),\varphi \left(0\right))`$. We can see that the Haar wavelet obeys the same pattern as the Daubechies wavelet: $`\psi _H=(\varphi \left(1\right),\varphi \left(0\right))`$. This pattern can be extended, using different coefficients but the same general sign rules, for wavelets with $`6`$, $`8`$, $`10`$, $`\mathrm{}`$ coefficients. The resulting wavelet and scaling functions become increasingly smooth, and therefore are better suited for data sets in which there is only a gradual change in the data set with position—in thermodynamic systems, this would be more useful for, say, a spin-$`N`$ Ising model than a spin-$`\frac{1}{2}`$ Ising model (presuming that $`N\frac{1}{2}`$).

### C Matrix formulation of the wavelet transform

For discrete systems, a conceptually simple method of implementing the wavelet transform is to set up the transform as a matrix equation. The input $`u`$ is converted into a column vector $`𝐮,`$ so that the coefficients $`s\left(i\right)`$ and $`\delta \left(i\right)`$ are obtained via the dot product of $`𝐮`$ with vectors $`𝐡\left(i\right)`$ and $`𝐥\left(i\right)`$, where the vectors are padded so that the first nonzero element is located at position $`i`$:

$`𝐡\left(i\right)`$ $`=(0,\mathrm{},0,\varphi \left(0\right),\varphi \left(1\right),\mathrm{},\varphi \left(r1\right),0,\mathrm{},0),`$
$`𝐥\left(i\right)`$ $`=(0,\mathrm{},0,\psi \left(0\right),\psi \left(1\right),\mathrm{},\psi \left(r1\right),0,\mathrm{},0).`$
The vectors $`𝐬`$ and $`\delta `$, which contain the wavelet-transformed coefficients of the decomposition, can be obtained by forming the matrices $`𝐇`$ and $`𝐋`$ and right-multiplying by the vector $`𝐮`$:

$`𝐬=\mathrm{𝐇𝐮}\text{ and }\delta =\mathrm{𝐋𝐮}\text{,}`$

where the rows of the matrices $`𝐇`$ and $`𝐋`$ are the vectors $`(𝐡\left(1\right),\mathrm{},𝐡\left(n\right))`$ and $`(𝐥\left(1\right),\mathrm{},𝐥\left(n\right))`$, respectively. As mentioned in the previous section, we do not need to keep all of the $`s\left(i\right)`$’s and $`\delta \left(i\right)`$’s in order to obtain perfect reconstruction of our signal; thus, we can obtain all the necessary coefficients in a single matrix multiplication by combining the relevant rows of $`𝐇`$ and $`𝐋`$ into a single matrix $`𝐖^{\left(1\right)}`$, whose rows are $`(𝐡\left(1\right),𝐡\left(3\right),\mathrm{},𝐡\left(n1\right),𝐥\left(1\right),𝐥\left(3\right),\mathrm{},𝐥\left(n1\right))`$. Thus, the wavelet transformation (1) can be written as

$$\left[\begin{array}{c}s\left(1\right)\\ \mathrm{}\\ s\left(n1\right)\\ \delta \left(1\right)\\ \mathrm{}\\ \delta \left(n1\right)\end{array}\right]=\left[\begin{array}{c}𝐡\left(1\right)\\ \mathrm{}\\ 𝐡\left(n1\right)\\ 𝐥\left(1\right)\\ \mathrm{}\\ 𝐥\left(n1\right)\end{array}\right]\left[\begin{array}{c}u\left(1\right)\\ u\left(2\right)\\ \\ \mathrm{}\\ \\ u\left(n\right)\end{array}\right].$$
(4)
We will denote the product on the left-hand side of (4) as $`\stackrel{~}{𝐮}`$.

As stated above, the wavelet process can be applied recursively: the set of averages $`(s\left(1\right),s\left(3\right),\mathrm{},s\left(n1\right))`$ can be treated as a new data sample $`𝐮^{\left(1\right)}=(u^{\left(1\right)}\left(1\right),\mathrm{},u^{\left(1\right)}\left(n/2\right))`$, and operated on by an $`\frac{N}{2}\times \frac{N}{2}`$ reduction of $`𝐖^{\left(1\right)},`$ which we denote $`𝐖^{\left(2\right)}`$, to produce a new set of $`n/4`$ averages $`(s^{\left(2\right)}\left(1\right),s^{\left(2\right)}\left(3\right),\mathrm{},s^{\left(2\right)}\left(n/21\right))`$ and corresponding new set of $`n/4`$ differences $`(\delta ^{\left(2\right)}\left(1\right),\delta ^{\left(2\right)}\left(3\right),\mathrm{},\delta ^{\left(2\right)}\left(n/21\right))`$. To reconstruct the original data set, we combine these $`n/2`$ values along with the $`n/2`$ differences $`(\delta ^{\left(1\right)}\left(1\right),\delta ^{\left(1\right)}\left(3\right),\mathrm{},\delta ^{\left(1\right)}\left(n1\right))`$ obtained from applying $`𝐖^{\left(1\right)}`$. This process can be repeated as many times as desired, dividing the $`m`$ non-downsampled averages $`𝐬^{\left(k\right)}`$ into $`m/2`$ averages $`𝐬^{\left(k+1\right)}`$ and $`m/2`$ differences $`\delta ^{\left(k+1\right)}`$. However, since at each iteration the matrix $`𝐖^{\left(k\right)}`$ only operates on selected elements of the vector $`\stackrel{~}{𝐮}^{\left(k\right)}=𝐖^{\left(k1\right)}𝐮^{\left(k1\right)}`$, computations for multiple levels can be performed at the same time. Thus, if we wish to apply the wavelet transform $`K`$ times, we can write this operation as an extended matrix product:

$$\stackrel{~}{𝐮}^{\left(K\right)}=\mathrm{𝐖𝐮}=\underset{k=1}{\overset{K}{}}𝐐^{\left(k\right)}𝐮\text{,}$$
(5)
where the $`𝐐^{\left(k\right)}`$ are a family of matrices of the form

$$𝐐^{\left(k\right)}=\left[\begin{array}{cc}𝐖^{\left(k\right)}& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝐈\end{array}\right].$$
(6)
In (6), $`𝐐^{\left(k\right)}`$ is always an $`N\times N`$ matrix, while the matrix $`𝐖^{\left(k\right)}`$ has size $`\left(N/2^{k1}\right)\times \left(N/2^{k1}\right)`$.

To recover the original data sample $`𝐮`$ following a wavelet transform $`𝐖`$, we can simply multiply $`\stackrel{~}{𝐮}`$ by the inverse of the wavelet matrix $`𝐖`$. The matrix $`𝐖`$ is unitary: that is, its inverse $`𝐖^1`$ is equal to its transpose $`𝐖^T`$. Consequently, if $`𝐖`$ is known, all that is necessary to reverse the transformation is to left-multiply $`\stackrel{~}{𝐮}^{\left(k\right)}`$ by the transpose $`𝐖^T`$. Moreover, the computation (4) of the wavelet transform can usually be carried out “in place” by manipulating local coordinates; in this manner, the computation is carried out even more rapidly than a standard multiplication, and without the increased storage costs associated with matrix multiplications.

### D Multidimensional wavelet transforms

Since virtually all problems in lattice thermodynamics are in multiple dimensions, it is necessary to take the wavelet transform of a multivariate function or data set. Several methods have been developed to carry out such transformations; among them are Cohen and Daubechies’s separable wavelets, which form the multidimensional scaling and wavelet functions $`\varphi (x,y)`$, $`\psi _{xx}(x,y)`$, $`\psi _{xy}(x,y)`$, and $`\psi _{yy}(x,y)`$ from products of the one-dimensional scaling and wavelet functions $`\varphi \left(x\right)`$ and $`\psi \left(x\right)`$ . A more general algorithm, the lifting algorithm, has been developed by Sweldens. It divides the wavelet transform into two steps: the first computes the wavelet coefficients $`\delta \left(i\right);`$ the second step uses the wavelet coefficients to speed up the calculation of the scaling coefficients. “Initialization” of the lifting algorithm requires the use of an appropriately selected basis function.

A particularly convenient basis function for the multidimensional lifting transform is the generalized orthogonal Haar wavelets outlined by Sweldens.  An extension of the one-dimensional wavelet transform, they can be created in any number of dimensions, and have the same basic orthonormality properties as the one-dimensional Haar functions, although the orthonormality constant becomes $`2^{d/2},`$ where $`d`$ is the dimensionality of the system. \[The two-dimensional version is shown in Figure 2.\] Moreover, the use of the Haar wavelets as a starting point for further iterations of the lifting algorithm allow the development of additional, “better” lifted wavelets with more desirable properties, such as smoothness.

It can further be seen that the “oversampling” problem which exists in one dimension will be magnified in multiple dimensions: since the number of wavelet functions which are produced from each data point increases by a factor of two with each additional dimension added, we must reduce the number of points maintained for each wavelet function by that same factor. Thus, in two dimensions, we keep only every fourth point; in three dimensions, every eighth point, and so on. The wavelet transform for multidimensional systems can thus still be written in the form of (5) and (6), after we have written the multidimensional data set in terms of a column vector. This can be accomplished by wrapping around the edges of the matrix in creating $`𝐮`$: for example, after inserting element $`(1,N)`$ of a two-dimensional data set into $`𝐮`$, we next store element $`(2,1)`$, and so forth. The other significant difference in the structure of these equations is that the size of the submatrix $`𝐖^{\left(k\right)}`$ in (6) is now $`N/2^{\left(k1\right)d}\times `$ $`N/2^{\left(k1\right)d}`$ instead of $`N/2^{k1}\times N/2^{k1}`$.

## III Wavelet analysis of lattice thermodynamics

### A Applying the wavelet transform to Hamiltonians

The standard model for studying the thermodynamic behavior of lattice systems is the spin-$`\frac{1}{2}`$ Ising model, which contains both nearest-neighbor pairwise interactions as well as interactions between lattice sites and an external field. The Hamiltonian for this system is normally written in the form

$$\beta =h\underset{i}{}\sigma _i+J\underset{ij}{}\sigma _i\sigma _j,$$
(7)
where $`h`$ is the strength of the external field in the direction of the spins $`\sigma _i`$, and $`J`$ is the strength of the interaction between nearest-neighbor pairs of spins on the lattice; these pairs are indicated by the subscript $`ij`$ in the second summation in (7). The inverse temperature $`\beta =\left(k_BT\right)^1`$; for convenience we let $`k_B=1`$, so that temperature, external field, and nearest-neighbor interactions are all dimensionless quantities. The model can be further extended by the inclusion of a position-dependent external field, or by the inclusion of pairwise interactions beyond nearest neighbors; in this case, (7) can be written in the more general form

$$\beta =\underset{i}{}h_i\sigma _i+\underset{i}{}\underset{j}{}J_{ij}\sigma _i\sigma _j,$$
(8)
where in (8) $`h_i`$ is the strength of the external field at lattice site $`i`$, and $`J_{ij}`$ is the strength of the interaction between sites $`i`$ and $`j`$. Analytical solutions of (7) are well-known for one-dimensional systems and two-dimensional systems when $`h=0`$; it has recently been suggested that analytic solutions do not exist for more complicated systems.

While (7) and (8) are compact representations of the Hamiltonian of the system, the expansion of the lattice variables $`s_i`$ and $`s_j`$ as a sum of wavelet coefficients makes these equations impractical for applying the wavelet transformation. Since the system is described discretely, we want to use discrete wavelets, and therefore a matrix formulation of the Hamiltonian would be convenient. Using graph theory, this is readily accomplished: let the vectors $`𝐮=(\sigma _1,\sigma _2,\mathrm{},\sigma _N)`$ and $`𝐡=(h_1,\mathrm{},h_N)`$ denote the values of each of the $`N`$ lattice variables in the system and the set of external-field interaction strengths, respectively, constructed in the row-wise manner described in the previous section. Furthermore, define the matrix $`𝐉`$ such that element $`J_{ij}`$ is the strength of the interaction between site $`i`$ and site $`j`$. If these sites do not interact, then $`J_{ij}=0`$. Then, the Hamiltonian (8) can be written in the form of a matrix equation:

$$\beta =𝐡^T𝐮+𝐮^T\mathrm{𝐉𝐮},$$
(9)
where the superscript $`T`$ denotes the transpose of the vector (or matrix) which precedes it.

The matrix $`𝐖`$ which defines the wavelet transform satisfies by construction $`𝐖^T𝐖=𝐈`$, where $`𝐈`$ is the identity matrix. Therefore, to apply the wavelet transform, we simply insert $`𝐖^T𝐖`$ between each pair of terms in (9), thereby obtaining

$$\beta =\left(𝐡^T𝐖^T\right)\left(\mathrm{𝐖𝐮}\right)+\left(𝐮^T𝐖^T\right)\left(\mathrm{𝐖𝐉𝐖}^T\right)\left(\mathrm{𝐖𝐮}\right).$$
(10)
Using the general matrix property that $`𝐁^T𝐀^T=\left(\mathrm{𝐀𝐁}\right)^T`$, (10) can be rewritten in terms of the wavelet-transformed vectors $`\stackrel{~}{𝐡}^{\left(K\right)}=\mathrm{𝐖𝐡}`$, $`\stackrel{~}{𝐮}^{\left(K\right)}=\mathrm{𝐖𝐮}`$, and the wavelet-transformed matrix $`\stackrel{~}{𝐉}^{\left(K\right)}=\mathrm{𝐖𝐉𝐖}^T`$:

$$\beta \stackrel{~}{}=\left(\stackrel{~}{𝐡}^{\left(K\right)}\right)^T\stackrel{~}{𝐮}^{\left(K\right)}+\left(\stackrel{~}{𝐮}^{\left(K\right)}\right)^T\stackrel{~}{𝐉}^{\left(K\right)}\stackrel{~}{𝐮}^{\left(𝐊\right)}.$$
(11)
It is important to note that in writing equation (11), we have not made any explicit assumptions about the form of the matrix $`𝐖`$, other than to require it to be a matrix describing a wavelet transform. As a result, we can perform several levels of multiresolution simultaneously through appropriate preparation of the matrix $`𝐖`$ in the manner outlined above. Using the inverse wavelet transform $`𝐖^1=𝐖^T`$ to recover an original configuration $`𝐮`$ is a one-to-one mapping only if we are provided with $`n`$ distinct data points, as contained in $`\stackrel{~}{𝐮}^{\left(K\right)}`$.

Examining (11), we see that the result of applying the wavelet transform to a set of spins $`𝐮`$ is to create a new representation $`\stackrel{~}{𝐮}^{\left(K\right)}`$, which contains the same information about the spins as does the original state vector $`𝐮`$. However, the vector $`\stackrel{~}{𝐮}^{\left(K\right)}`$ contains $`n/2^{dK}`$ averages, where $`d`$ is the dimensionality of the lattice; these averages can be viewed as “block spins” in a sense similar to that of Kadanoff. The remaining elements of $`\stackrel{~}{𝐮}^{\left(K\right)}`$ contain the local differences in the spins; that is, they can be used to describe the specific set of spins which give rise to a particular block spin $`\stackrel{~}{s}_i^{\left(K\right)}`$.

### B Computing thermodynamic functions

The canonical partition function

$$Z=\underset{𝐮𝕊}{}\mathrm{exp}\left(\beta \left(𝐮\right)\right),$$
(12)
where $`𝕊`$ is the configuration space of the system, can be used to derive all the thermodynamic properties of a lattice system. Applying the wavelet transform $`𝐖`$ to the state vector $`𝐮`$ results in a new state vector $`\stackrel{~}{𝐮}^{\left(K\right)}`$ belonging to the configuration space $`\stackrel{~}{}S^{\left(K\right)}`$. Provided that $`𝐖`$ satisfies the perfect reconstruction property, if the summation over $`𝐮\stackrel{~}{}S`$ in (12) is replaced with the summation over $`\stackrel{~}{𝐮}^{\left(K\right)}\stackrel{~}{}S^{\left(K\right)}`$, the results will be identical. This result follows naturally, since perfect reconstruction necessarily implies that there is a unique state vector $`\stackrel{~}{𝐮}^{\left(K\right)}\stackrel{~}{}S^{\left(K\right)}`$ for each state vector $`𝐮𝕊`$, and by construction of the wavelet-transformed Hamiltonian (11), $`\beta =\beta \stackrel{~}{}^{\left(K\right)}.`$

Ensemble averages involving wavelet-transformed variables are in general no more complicated than computations involving the original variables. In general, the transformation of a function (or functional) of the original variables $`f\left(𝐮\right)`$ or $`f\left[u\right]`$ will generally have the same characteristics after applying the wavelet transform to obtain $`\stackrel{~}{f}^{\left(K\right)}\left(\stackrel{~}{𝐮}^{\left(K\right)}\right)`$ or $`\stackrel{~}{f}^{\left(K\right)}\left[\stackrel{~}{u}^{\left(K\right)}\right]`$. Moreover, the standard properties of ensemble averaging, such as linearity, also apply, which makes calculations of moments of the distribution particularly simple.

The above formulation can be applied to any system whose Hamiltonian can be written in the form of equation (8), or as the sum of contributions, each of which is of that form. While the wavelet transform can be applied to any lattice system whose Hamiltonian is a function of the components of the state vector $`𝐮`$—or, indeed, to any Hamiltonian which is a functional of $`𝐮\left(𝐫\right)`$ for continuous systems—in most of these cases, it is necessary to rely on the more cumbersome series expansions, or integral transforms in the case of continuous systems.

In addition, for Ising spin variables, the use of the wavelet transform poses an additional challenge. While it is entirely straightforward to describe the possible values an individual lattice site can take—for a spin-$`q`$ Ising model, the allowed values of an individual spin $`\sigma `$ are $`q`$, $`q+1,\mathrm{},q1,q`$—the rules which determine whether an arbitrarily selected “transformed” state vector $`\stackrel{~}{𝐮}^{\left(K\right)}`$ represents a real state vector $`𝐮`$ are cumbersome to manipulate and to implement, and have proven a formidable challenge in prior research as well.

## IV Analysis

In this section, we use two variants of the two-dimensional Ising model as a basis for our calculations: we look at $`4\times 4`$ and $`32\times 32`$ Ising lattices; the former is used for calculations when it is desired to perform calculations over all states explicitly, while the latter is illustrative of larger systems, for which exact calculations are intractable. The emphasis in this paper will be on the use of wavelets to yield approximate answers in significantly faster time than is possible with a Monte Carlo simulation incorporating all degrees of freedom explicitly. The methodology by which the wavelet transform can be extended to Monte Carlo simulations of lattice systems is the focus of the companion paper.

Before proceeding to the results of the calculations, we make note of the time required to execute the simulations. Each of Figures 3 through 6 plots the observed variables in the temperature range $`T=0.50`$ to $`T=5.00`$, with a step size of $`\mathrm{\Delta }T=0.01`$. The computations for the original problem, with $`2^{16}=65536`$ states to consider explicitly, required more than 6.8 seconds per point to execute the required calculations on a 733 MHz Pentium III; the same calculations using one and two applications of the wavelet transform required just 0.061 and 0.026 seconds per point to consider $`5^4=625`$ and $`17`$ states, respectively.

### A Weighting functions for wavelet-transformed statistics

Until now, we have worked with wavelet transforms which preserve the number of degrees of freedom between the original and transformed problems. This approach yields results for the wavelet-transformed system in exact agreement with those for the original system. However, as mentioned above, in such cases the transformed equations are usually harder to model than the original ones. Hence it is desirable to use the wavelet transform not only as a means of describing a lattice system, but also to derive an approximation scheme whereby estimates of thermodynamic properties can be made efficiently, while still offering error estimates that bound the true results.

The approach we adopt in this paper is to ignore all local differences: that is, we assume that $`\delta _i^{\left(k\right)}=0`$ for all values of $`i`$ and $`k`$. As a secondary assumption, we assume that correlation functions which include wavelet coefficients are also equal to zero: that is, we assume $`\delta _i^{\left(k\right)}A()=0`$ for any choice of property $`A()`$. These extreme assumptions represent a “worst-case scenario” for the use of the wavelet transform method; any more accurate representation of the behavior of the wavelet coefficients will lead to similarly more accurate results in the calculations we present below. Note that under certain circumstances, these approximations are accurate: for Hamiltonians of the form (8) which do not exhibit interactions with an external field, there exist configurations with equal energies but opposite signs for $`\delta _i^{\left(k\right)}`$; consequently, when we take the average over all configurations, the ensemble average of $`\delta _i^{\left(k\right)}`$ will vanish.

Let us assume that we have applied the wavelet transform to our original Hamiltonian, $`\left(𝐬\right),`$ and have obtained a new function $`\stackrel{~}{}^{\left(K\right)}\left(\stackrel{~}{𝐮}^{\left(K\right)}\right)`$. If we then make our approximation that all the wavelet coefficients are negligible, then we can reduce $`\stackrel{~}{}^{\left(K\right)}\left(\stackrel{~}{𝐮}^{\left(K\right)}\right)`$ to a new function $`\stackrel{~}{}^{\left(K\right)}\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)`$, where $`\stackrel{~}{𝐬}^{\left(K\right)}`$ represents the set of all averages $`(s_1^{\left(K\right)},\mathrm{},s_m^{\left(K\right)})`$ that were preserved by the wavelet transform. We know that the partition function of the system before we performed the wavelet transform is given by (12); after using the wavelet transform, we hypothesize that the partition function of the new system is given by

$$\stackrel{~}{Z}=\underset{\stackrel{~}{𝐬}^{\left(K\right)}\stackrel{~}{}S^{\left(K\right)}}{}w\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)\mathrm{exp}\left(\beta \stackrel{~}{}^{\left(K\right)}\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)\right),$$
(13)
where the $`w\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)`$ is the weight of configuration $`\stackrel{~}{𝐬}^{\left(K\right)}`$ in the configuration space $`\stackrel{~}{}S^{\left(K\right)}`$; since multiple configurations $`𝐮𝕊`$ can correspond to the same $`\stackrel{~}{𝐬}^{\left(K\right)}`$, we cannot set $`w\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)=1`$ as we did in the untransformed Ising model. Clearly, the new partition function (13) will be identical to the original partition function (12) if we define

$$w\left(\stackrel{~}{𝐬}\right)=\underset{𝐮:𝐮\stackrel{~}{𝐬}^{\left(K\right)}}{}\mathrm{exp}\left(\beta \left(\left(𝐮\right)\stackrel{~}{}^{\left(K\right)}\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)\right)\right),$$
(14)
where the notation $`𝐮:𝐮\stackrel{~}{𝐬}^{\left(K\right)}`$ denotes that the sum is to be performed over all configurations $`𝐮`$ which have the same set of wavelet-transformed averages $`\stackrel{~}{𝐬}^{\left(K\right)}`$. However, evaluating (14) to obtain the weighting functions is no more tractable than the computation of the original partition function. Thus, any computational efficiency to be gained is by finding an economical approximation for (14). One such approach is to take $`w\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)`$ to be equal to the number of states $`𝐬`$ whose wavelet transform yields $`\stackrel{~}{𝐬}^{\left(K\right)}`$, so that (13) is the standard form of the canonical partition function for a system with degenerate energy levels. Let us call the number of such states $`𝐬`$ with equivalent averages $`\stackrel{~}{𝐬}^{\left(K\right)}`$ the degeneracy of state $`\stackrel{~}{𝐬}^{\left(K\right)}`$, and denote this degeneracy as $`g\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)`$. The restriction of $`\stackrel{~}{𝐮}^{\left(K\right)}`$ to the averages $`\stackrel{~}{𝐬}^{\left(K\right)}`$ prevents a unique reconstruction of $`𝐮,`$ unless $`g\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)=1.`$

### B Order parameter

A natural variable to compute is the order parameter, generically denoted $`\eta `$, which for lattice spin systems is the average magnetization,

$$m=\frac{1}{N}\underset{i}{}\sigma _i;$$
(15)
for other members of the Ising universality class, the order parameter can represent the overall density $`\rho `$ or the difference between the densities of two phases. While the computation of the order parameter is straightforward in simulations, its calculation can be made difficult in the case of zero external field because of symmetries in the configuration space: for every configuration with magnetization $`m_i`$, there exists another configuration with the same energy and magnetization $`m_i`$. As a result, when all the states are combined using (15), we do not observe spontaneous magnetization, but instead find $`m=0`$ at all temperatures. Thus, we consider the absolute value of the magnetization, $`\left|m\right|`$ in place of the magnetization $`m`$. The results of this calculation for the $`4\times 4`$ Ising model with the wavelet transform, setting all $`\delta \left(i\right)`$’s to zero, and without the wavelet transform are shown in Figure 3. We note that the error is essentially negligible for temperatures below $`T=1`$, and decreases again for large values of $`T`$, where differences in energy levels become negligible and the average magnetization of a state is the only contributing factor to $`\left|m\right|`$.

### C Free-energy considerations

The Helmholtz free energy is just the logarithm of the partition function: $`A=k_BT\mathrm{ln}Z`$; if we estimate the partition function using an expression such as (13), we naturally expect the approximated value $`A`$ to differ from its true value. When we examine the behavior of the $`4\times 4`$ Ising model under the wavelet transform, we find as expected that the two free energy surfaces are similar, although they are clearly not identical. In particular, the exact numerical values obtained from the two equations are not the same. However, since the assignment $`A=0`$ is arbitrary in any system, we can choose to define $`A=0`$ as either the maximum or the minimum free energy obtained in each system. Under these conditions, the energy scales for the exact calculation and the wavelet transform calculation are essentially identical, particularly at the so-called “fixed points” of the system—that is, for points which are not affected by renormalization transformations. For the two-dimensional Ising model, these points are at $`T=0`$, and at infinite (positive or negative) values of the external field interaction $`h`$.

This agreement at the fixed points should not be surprising: the fixed points correspond to unique configurations of the system, such that for the configurations $`\stackrel{~}{𝐬}^{\left(K\right)}`$ which result from the wavelet transform of the fixed points $`𝐮^{}`$, the degeneracy of the states is unity. Consequently, the approximation $`w\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)=g\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)`$ is correct for the dominant configuration in the system, and therefore the behavior of the approximate partition function $`Z^{}`$ is almost identical to that of the true partition function $`Z`$ in the vicinity of the fixed points.

As can be seen in Figure 4, the free energy converges to the same values in the low-temperature limit, where only a few states which are essentially unaffected by the wavelet transform contribute to the partition function of the system. At high temperatures, the agreement is less exact, because of the approximations made in evaluating the Boltzmann factors of the block spin configurations.

### D Entropy of a wavelet-transformed lattice system

It is relatively straightforward to show that coarse-graining a system with a countable phase space—that is, a phase space whose states can be completely enumerated—requires a correction factor to ensure satisfactory agreement with results from the fine-scale system. As an example, consider the entropy of a wavelet-transformed spin-$`\frac{1}{2}`$ Ising model. The original lattice has $`2^{N_t}`$ total configurations, and therefore in the high-temperature limit, the entropy approaches

$$S_{\mathrm{max}}=k_BN_t\mathrm{ln}2.$$
(16)
If we ignore the weighting factor of the wavelet-transformed system by setting $`w\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)=1`$ for all states $`\stackrel{~}{𝐬}^{\left(K\right)}`$, the resulting system will have $`\left(1+_{i=1}^{K1}N^{\left(i\right)}\right)^{N^{\left(K\right)}}`$, and the maximum possible entropy for this system is

$$S_{\mathrm{max}}^{\left(K\right)}=k_BN^{\left(K\right)}\mathrm{ln}\left(\underset{i=1}{\overset{K1}{}}N^{\left(k\right)}+1\right).$$
(17)
Comparing (17) to the original limit of $`S_{\mathrm{max}}=k_BN_t\mathrm{ln}2`$, we must have $`S_{\mathrm{max}}>S_{\mathrm{max}}^{\left(K\right)}`$, since the configuration space is smaller, which indicates that our coarse-graining has effectively reduced the available entropy of the system. We can conclude from this that unless the volume of phase space is conserved by the coarse-graining procedure, some entropy will be lost at high temperatures, regardless of the accuracy of the coarse-graining. Our choice of the degeneracy $`g\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)`$ for the weighting function $`w\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)`$ preserves the volume of phase space in the limit of $`T\mathrm{}`$ for the systems under study here.

Looking at the specific case of the $`4\times 4`$ Ising model, we compute the entropy as a function of temperature for zero external field for the exact problem, and for the wavelet-transformed problem using one and two iterations of the two-dimensional Haar wavelet. The results are shown in Figure 5. As before, we find that errors vanish in the low-temperature limit, where only the lowest-energy states make a contribution to the partition function; since there are two such states, namely those with all spins up and all spins down, we find that the $`T0`$ limit of the entropy is

$`S\left(T0\right){\displaystyle \underset{i=1}{\overset{2}{}}}{\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{1}{2}}=\mathrm{ln}2`$

In the high-temperature limit, the entropies agree because, by construction, we have $`_{\stackrel{~}{𝐬}\stackrel{~}{}S}w\left(\stackrel{~}{𝐬}\right)=2^{16},`$and therefore, since the Boltzmann factor goes to unity for all states as $`T\mathrm{}`$, $`Z_{\stackrel{~}{𝐬}\stackrel{~}{}S}w\left(\stackrel{~}{𝐬}\right)`$, and the entropy tends toward the value $`S_{\mathrm{max}}`$ given by (16). Disagreement in the intermediate temperature regime is largely the result of grouping together states with different energies into a single transformed state with a single energy. In particular, the grouping of higher-energy states together at a lower energy increases the probability of being in those states, and therefore increases the total entropy of the system. This explains the seemingly anomalous result of coarse-graining increasing the entropy at intermediate temperatures observed in Figure 5.

For more general systems, we can consider the formal definition of the entropy for a continuous distribution,

$$S=k_B_{𝐮𝕊}𝑑𝐮p\left(𝐮\right)\mathrm{ln}p\left(𝐮\right),$$
(18)
where $`p\left(𝐮\right)`$ is the probability of observing the configuration $`𝐮`$. For the original Ising model, each configuration $`𝐮`$ is unique, and therefore the probability $`p\left(𝐮\right)`$ is just the Boltzmann weight $`Z^1\mathrm{exp}\left(\beta \left(𝐮\right)\right),`$ and thus the entropy is given by

$$S=k_B\mathrm{ln}Z+\frac{k_B}{Z}_{𝐮𝕊}𝑑𝐮\beta \left(𝐮\right)\mathrm{exp}\left(\beta \left(𝐮\right)\right).$$
(19)
After performing the wavelet transform, and discarding finer-scale details, we must account for the weighting factor $`w\left(\stackrel{~}{𝐬}\right)`$ in our expression for the entropy. Consequently, the entropy for such a system is given by

$`S^{\left(K\right)}`$ $`=`$ $`k_B{\displaystyle _{\stackrel{~}{𝐬}\stackrel{~}{}S}}𝑑\stackrel{~}{𝐬}w\left(\stackrel{~}{𝐬}\right)p\left(\stackrel{~}{𝐬}\right)\mathrm{ln}p\left(\stackrel{~}{𝐬}\right)`$ (20)
$`=`$ $`k_B\mathrm{ln}\stackrel{~}{Z}+{\displaystyle \frac{k_B}{\stackrel{~}{Z}}}{\displaystyle _{\stackrel{~}{𝐬}\stackrel{~}{}S}}𝑑\stackrel{~}{𝐬}w\left(\stackrel{~}{𝐬}\right)\beta \stackrel{~}{}^{\left(K\right)}\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)e^{\beta \stackrel{~}{}^{\left(K\right)}\left(\stackrel{~}{𝐬}^{\left(K\right)}\right)},`$ (21)
where the coarse-grained partition function $`\stackrel{~}{Z}`$ is given by (13).

### E Fluctuation properties

For a given thermodynamic system, the constant-volume heat capacity $`C`$ is defined as the fluctuation in the internal energy of the system:

$$C=\frac{E^2E^2}{k_BT^2}.$$
(22)
It is well-known that in the vicinity of a continuous phase transition, the heat capacity $`C`$ diverges. At the same time, it is also known that no system of finite size can have $`C\mathrm{}`$, since the energy, and hence the variance $`E^2E^2`$ in the energy of the system will also remain finite. However, we can still observe evidence of power-law divergence near the critical temperature. For small finite systems, we do not observe evidence of a divergence in the heat capacity; instead, the heat capacity is a smooth function of temperature $`T`$. The wavelet transform largely preserves the behavior of the original model; we observe the same general functional form in both the original model as well as the transformed systems, as seen in Figure 6. However, for large finite systems, we can still detect the characteristic power-law divergence in the vicinity of the critical point; such a divergence is shown for a $`32\times 32`$ Ising model as Figure 7.

Applying the wavelet transform to a given thermodynamic system, we expect that we will still find evidence of a critical point; however, because of the mean-field like behavior of the wavelet transform method, the critical point should be found at a higher temperature in the wavelet-transformed system than in the original system. Computing the heat capacity for large lattices cannot be done by exhaustive enumeration; therefore, we save discussion of our numerical results for a subsequent paper, and present a brief heuristic argument. As the size of a block spin increases, corresponding to a greater number of iterations of the wavelet transform, we expect that the suppression of fluctuations in the system will lead to an increase in the temperature of the phase transition, while the actual numerical value of the heat capacity itself decreases. The decrease in the number of degrees of freedom in the system is accompanied by an increasing trend towards homogeneity of the system—the remaining configurations begin to look more and more similar to one another. As a result, the variance measured by $`E^2E^2`$ decreases with increasing amounts of coarse-graining. This trend is much less pronounced for small systems, as shown in Figure 6. When the number of degrees of freedom is small, the ability to perform an exact enumeration ensures that the only errors observed are those introduced by the wavelet approximation itself. Thus, since the free energy and entropy of small systems are accurately reproduced by the transform, only minor deviations in the phase transition temperature for small systems are expected.

## V Similarities to Renormalization Group theories

It should already be apparent that close affinities exist between the method described here and position-space renormalization methods. The construction of wavelet-based averages naturally corresponds to Kadanoff’s concept of a “block spin” transformation: both the present method and Kadanoff’s approach rely on combining a region of contiguous spins into a single new “block spin.” However, with the approach outlined above, we do not seek to impose the requirement that the block spins must be restricted to have the same spin values as the original spins. For example, in a wavelet-transformed spin-$`\frac{1}{2}`$ Ising model, the block spins can take on values other than $`+1`$ and $`1`$. In addition, we do not impose “majority rules” or other “tiebreakers” in the case where there are equal numbers of “up” and “down” spins contained in a single block.

A less apparent connection can also be established with Migdal’s “bond-moving” approximation, or similarly with Wilson’s recursion method. In particular, Migdal’s method can be compared to applying a separable wavelet transformation, in that the recursion is performed in only one spatial direction at a time. However, using the methods of either Migdal or Wilson, after performing the bond-moving or decimation transformation, the new interaction strengths $`J_{ij}^{}`$ are determined by manipulating the resulting Hamiltonian and recasting it in the form of the original, leading at times to very complicated, nonlinear formulas for the $`J_{ij}^{}`$ as a function of $`J_{ij}`$. By contrast, using wavelets, the values for the transformed coupling coefficients $`\stackrel{~}{J}_{ij}`$ are obtained directly by transforming the interaction matrix $`𝐉`$, since $`\stackrel{~}{𝐉}=\mathrm{𝐖𝐉𝐖}^T`$. The trade-off that must be made for the algorithmic transparency of obtaining the new coupling constants via a matrix multiplication is that we cannot solve for the fixed points of the recursion relation, which makes the production of a renormalization “flow diagram” using the wavelet method a difficult problem.

The behavior of the wavelet transform and of the renormalization group differs in another, important manner. By the nature of its construction, the group of renormalization transformations is only a semi-group, inasmuch as reversing the mapping to move from a coarsened description to a more detailed description is impossible. Using the wavelet transform, a reverse mapping is theoretically possible: reversibility fails because of the approximations invoked to reduce the number of wavelet coefficients which are kept, rather than as an inherent limitation of the wavelet transform itself.

## VI Conclusions

Using the wavelet transform as a basis for calculations of lattice systems yields an impressive reduction in the computational time required to obtain estimates of thermodynamic properties, at the price of modest errors in accuracy, given the relatively small size of the systems considered here, and the approximations introduced to simplify our calculations. The wavelet transform provides a systematic approach to coarse-grain systems to any level of resolution. The method produces correct results in the vicinity of fixed attractors, with decreasing accuracy as one approaches the critical point in parameter space.

In addition, the use of the wavelet transform on such systems is simple to implement: after selecting a particular set of wavelets, all of the operations can be reduced to matrix multiplications, as shown in Section III A. The computational requirements needed to implement the transformation are small: if the transformation cannot be accomplished in the explicit form of a matrix multiplication, it is possible to organize the calculation to be performed in-place. Moreover, the amount of coarse-graining that can be achieved using the wavelet transform can be adjusted dynamically. These properties of the wavelet transform lead naturally to an extension of the method to systems for which simulations are required; the resulting algorithm is the focus of the companion paper.

## VII Acknowledgments

Funding for this research was provided in part by a Computational Sciences Graduate Fellowship (AEI) sponsored by the Krell Institute and the Department of Energy.


format: text

text:
# Optimization of evaporative cooling towards a large number of Bose-Einstein condensed atoms

## I Introduction

Evaporative cooling is an essential experimental technique in the recent demonstrations of Bose-Einstein condensation (BEC) with magnetically trapped atomic vapor . This cooling method is based on both the selective removal of energetic atoms through evaporation and collisional rethermalizations among remaining atoms . The successful utilization of interatomic elastic collisions makes evaporative cooling highly efficient. A real system, however, involves other atomic collisions, such as elastic collisions with background gas, inelastic collisions due to dipolar relaxation and three-body recombination processes. These undesirable collisions result in a loss of the trapped atoms and seriously reduce the effectiveness of evaporative cooling.

For alkali-metal atoms, such as 87Rb and 23Na, three-body recombination becomes a dominant loss mechanism during evaporative cooling , and limits the final number of condensed atoms, $`N_0`$, in the experiments. It is known that 87Rb atoms suffer much severer losses than 23Na atoms due to their much large rate coefficient of three-body recombination process . BEC experiments have clearly reflected this fact: the reported $`N_0`$ value in 87Rb is $`10^6`$ at most , while that in 23Na is up to $`5\times 10^7`$ . The small BEC of 87Rb is a great disadvantage with respect to its future application as an “atom laser” and has to be overcome.

In this paper, we report the theoretical optimization of evaporative cooling, focusing on BEC experiments with 87Rb atoms. Several theoretical works have considered the optimization of evaporative cooling to increase the cooling efficiency . Especially, Ref.  optimized the evaporative cooling in atomic hydrogen not only theoretically but also experimentally. However, these studies are based on kinetic theory of a classical gas which fails around the BEC transition point, and therefore don’t give us any information on the number of condensed atoms produced by evaporative cooling. Our analysis is based on quantum kinetic theory of a Bose gas  and can provide quantitative results on the increase in the number of condensed atoms by optimizing the evaporative cooling. We optimize the temporal frequency sweep of the radio-frequency (rf) magnetic field under a fixed trapping potential by following the variational method developed by Sackett et al. for evaporative cooling in a classical gas system . We show that an accelerated sweep of the rf-field frequency around the BEC transition point can largely enhance the final number of condensed atoms compared with conventional cooling with an exponential sweep of the frequency. It is also demonstrated that highly efficient evaporative cooling can be achieved by an initial exponential and then a rapid linear sweep of frequency. Our present results provide a useful guideline to experimentalists aiming to achieve large condensates in 87Rb atoms.

## II Kinetic theory of evaporative cooling

### A Thermodynamics of trapped bosonic atoms

We begin with a review of our theoretical formulation . During evaporative cooling, a magnetic potential, $`U(𝐫)`$, is truncated depending on the frequency of the applied rf-magnetic field. The thermalized distribution of noncondensed atoms in such a truncated potential is well approximated by the truncated Bose-Einstein distribution function such that

$$\stackrel{~}{f}(𝐫,p)=\frac{1}{\mathrm{exp}(ϵ_p/k_BT)/\stackrel{~}{\xi }(𝐫)1}\mathrm{\Theta }(\stackrel{~}{A}(𝐫)ϵ_p),$$
(1)
where $`T`$ is the temperature, $`ϵ_p=p^2/2m`$ is the kinetic energy of atoms with momentum $`p`$ and mass $`m`$, and $`\mathrm{\Theta }(x)`$ is a step function. In this paper, we add tilde to the notations of the quantities evaluated through this truncated distribution function. $`\stackrel{~}{\xi }(𝐫)`$ denotes the local fugacity including mean-field potential energy; $`\stackrel{~}{\xi }(𝐫)=\mathrm{exp}\{[\mu U(𝐫)2v\stackrel{~}{n}(𝐫)]/k_BT\}`$, where $`\mu `$ is the chemical potential, $`v=4\pi a\mathrm{}^2/m`$ is the interaction strength of atoms in proportional to the s-wave scattering length $`a`$, and $`\stackrel{~}{n}(𝐫)`$ is the density profile of atoms. The step function eliminates the momentum states whose kinetic energy exceeds the effective potential height $`\stackrel{~}{A}(𝐫)=ϵ_tU(𝐫)2v\stackrel{~}{n}(𝐫)`$, where $`ϵ_t`$ is the truncation energy of a magnetic potential. We calculate the density profile of atoms using this truncated distribution function in a self-consistent way such that

$$\stackrel{~}{n}(𝐫)=4\pi h^3\stackrel{~}{f}(𝐫,p)p^2𝑑p,$$
(2)
and the internal energy density similarly as

$$\stackrel{~}{e}(𝐫)=4\pi h^3ϵ_p\stackrel{~}{f}(𝐫,p)p^2𝑑p+v\stackrel{~}{n}^2(𝐫)+U(𝐫)\stackrel{~}{n}(𝐫).$$
(3)
The spatial integrations of these density functions give the total number of atoms, $`\stackrel{~}{N}=\stackrel{~}{n}(𝐫)𝑑𝐫`$, and the total internal energy, $`\stackrel{~}{E}=\stackrel{~}{e}(𝐫)𝑑𝐫`$ respectively.

After the BEC transition, the atom density in the condensed region is obtained by the sum of a condensate $`n_0(𝐫)`$ and saturated noncondensed atoms $`\stackrel{~}{n}_n`$ such that $`\stackrel{~}{n}(𝐫)=n_0(𝐫)+\stackrel{~}{n}_n`$. The condensate is described by the Thomas-Fermi distribution $`n_0(𝐫)=n_pU(𝐫)/v`$ with the peak density of the condensate at the center of potential $`n_p`$, while $`\stackrel{~}{n}_n`$ is evaluated through the truncated Bose-Einstein distribution function with the condition $`\stackrel{~}{\xi }(𝐫)=1`$.

### B Dynamics of evaporative cooling

The dynamics of trapped atoms during evaporative cooling is investigated on the basis of quantum kinetic theory of a Bose gas, developed previously by the present authors . Trapped atoms are removed from the potential by not only evaporation but also undesirable collisions, as mentioned in the introduction. The change rates (i.e., loss rates) of the total number of atoms and the total internal energy are calculated respectively as the sum of the contributions of all processes such that,

$`{\displaystyle \frac{d\stackrel{~}{N}}{dt}}`$ $`=`$ $`{\displaystyle \dot{n}_{\mathrm{ev}}(𝐫)d𝐫}`$ (5)
$`{\displaystyle \underset{s=1}{\overset{3}{}}}G_s{\displaystyle K_s(𝐫)[\stackrel{~}{n}(𝐫)]^sd𝐫}+\left({\displaystyle \frac{\stackrel{~}{N}}{ϵ_t}}\right)_{T,\mu }\dot{ϵ}_t,`$
$`{\displaystyle \frac{d\stackrel{~}{E}}{dt}}`$ $`=`$ $`{\displaystyle \dot{e}_{\mathrm{ev}}(𝐫)d𝐫}`$ (7)
$`{\displaystyle \underset{s=1}{\overset{3}{}}}G_s{\displaystyle K_s(𝐫)\stackrel{~}{e}(𝐫)[\stackrel{~}{n}(𝐫)]^{s1}d𝐫}+\left({\displaystyle \frac{\stackrel{~}{E}}{ϵ_t}}\right)_{T,\mu }\dot{ϵ}_t.`$
The rates $`\dot{n}_{\mathrm{ev}}`$ and $`\dot{e}_{\mathrm{ev}}`$ denote the evaporation rates of density functions derived from a general collision integral of a Bose gas system . The parameters $`G_1`$, $`G_2`$, and $`G_3`$ are the decay rate constants of trapped atoms due to the background gas collisions, dipolar relaxation, and three-body recombination, respectively. $`K_s`$ represents the correlation function which describes the $`s`$-th order coherence of trapped atoms, and we assume the expressions of $`K_s`$ for an ideal Bose gas system . The terms proportional to the change rate of truncation energy, $`\dot{ϵ}_t=dϵ_t/dt`$, give the contribution of extra atoms “spilled over” when $`ϵ_t`$ changes continuously in the forced evaporative cooling . The dynamics of the evaporative cooling including the loss of trapped atoms is described by the coupled differential equations in Eq. (5). In our theoretical framework based on truncated Bose-Einstein distribution function, all thermodynamic quantities are given as the complicated functions of three independent variables: temperature $`T`$, chemical potential $`\mu `$, and truncation energy $`ϵ_t`$. We derive here the useful thermodynamic relations in terms of the change rates of these variables:

$`{\displaystyle \frac{d\stackrel{~}{N}}{dt}}`$ $`=`$ $`\left({\displaystyle \frac{\stackrel{~}{N}}{T}}\right)_{\mu ,ϵ_t}\dot{T}+\left({\displaystyle \frac{\stackrel{~}{N}}{\mu }}\right)_{T,ϵ_t}\dot{\mu }+\left({\displaystyle \frac{\stackrel{~}{N}}{ϵ_t}}\right)_{T,\mu }\dot{ϵ}_t,`$ (8)
$`{\displaystyle \frac{d\stackrel{~}{E}}{dt}}`$ $`=`$ $`\left({\displaystyle \frac{\stackrel{~}{E}}{T}}\right)_{\mu ,ϵ_t}\dot{T}+\left({\displaystyle \frac{\stackrel{~}{E}}{\mu }}\right)_{T,ϵ_t}\dot{\mu }+\left({\displaystyle \frac{\stackrel{~}{E}}{ϵ_t}}\right)_{T,\mu }\dot{ϵ}_t,`$ (9)
where $`\dot{T}`$ and $`\dot{\mu }`$ are the change rate of temperature and that of chemical potential respectively. Substituting these relations into Eq. (5), we clearly see that the terms proportional to $`\dot{ϵ}_t`$ are canceled out and both $`\dot{T}`$ and $`\dot{\mu }`$ are determined by three variables $`T`$, $`\mu `$, and $`ϵ_t`$.

### C Optimization of cooling trajectory

The optimum cooling trajectory in the evaporative cooling process is obtained approximately by following the instantaneous optimization procedure investigated in Ref. . We determine the optimized cooling trajectory so as to achieve the maximum increase in phase-space density $`\stackrel{~}{\rho }`$ with the smallest decrease in the number of trapped atoms $`\stackrel{~}{N}`$. In the nonuniform gas system, the phase-space density is evaluated at the peak position and given by $`\stackrel{~}{\rho }=\stackrel{~}{n}(\mathrm{𝟎})\lambda ^3`$ with the thermal de Broglie wavelength $`\lambda =h/\sqrt{2\pi mk_BT}`$. At each time, we calculate the optimum truncation energy $`ϵ_t^{}`$ by maximizing the derivative, $`d(\mathrm{ln}\stackrel{~}{\rho })/d(\mathrm{ln}\stackrel{~}{N}`$), which describes the efficiency of evaporative cooling . We use a more useful form of this derivative rewritten in terms of change rates of temperature, chemical potential, and truncation energy as,

$`{\displaystyle \frac{d(\mathrm{ln}\stackrel{~}{\rho })}{d(\mathrm{ln}\stackrel{~}{N})}}`$ $`=`$ $`{\displaystyle \frac{\left(\frac{\stackrel{~}{\rho }}{T}\right)_{\mu ,ϵ_t}\dot{T}+\left(\frac{\stackrel{~}{\rho }}{\mu }\right)_{T,ϵ_t}\dot{\mu }+\left(\frac{\stackrel{~}{\rho }}{ϵ_t}\right)_{T,\mu }\dot{ϵ}_t}{\left(\frac{\stackrel{~}{N}}{T}\right)_{\mu ,ϵ_t}\dot{T}+\left(\frac{\stackrel{~}{N}}{\mu }\right)_{T,ϵ_t}\dot{\mu }+\left(\frac{\stackrel{~}{N}}{ϵ_t}\right)_{T,\mu }\dot{ϵ}_t}}`$ (11)
$`\times \left({\displaystyle \frac{\stackrel{~}{N}}{\stackrel{~}{\rho }}}\right).`$
In the optimized evaporative cooling, the truncation parameter $`\eta =ϵ_t/k_BT`$ changes more slowly than the variables $`T`$, $`\mu `$, and $`ϵ_t`$. We neglect the change rate of this truncation parameter, and consequently $`\dot{ϵ}_t`$ is related to $`\dot{T}`$ as $`\dot{ϵ}_t/ϵ_t=\dot{T}/T`$ in Eq. (11). The value of optimum truncation energy $`ϵ_t^{}`$ is given as a function of both temperature $`T`$ and chemical potential $`\mu `$. Substituting this $`ϵ_t^{}`$ into Eq. (5), we obtain the coupled differential equations describing the dynamics of the optimized evaporative cooling. The time-evolution calculations are carried out numerically by repeating a series of steps consisting of instantaneous optimization and dynamical short-time evolution of Eq. (5) with the temporal value of $`ϵ_t^{}`$. Furthermore we assume quick rethermalization, where the system is always described by the truncated distribution function . The obtained function $`ϵ_t^{}(t)`$ corresponds to the optimized sweeping function of the rf-field frequency.

In the next section, we will show the numerical results calculated by using relevant experimental parameters and discuss the experimental possibility of making large condensates in 87Rb atoms.

## III Numerical Results and Discussions

### A Optimized evaporative cooling

In our numerical calculations, we employed parameters corresponding to the experiments performed at the University of Tokyo . The cloverleaf magnetic trap with bias field $`B_0`$, radial gradient $`B^{^{}}`$, and axial curvature $`B^{^{\prime \prime }}`$ was modeled by the magnetic field of Ioffe-Pritchard type such that $`B(r,z)=\sqrt{(\alpha r)^2+(\beta z^2+B_0)^2}`$, where $`\alpha =\sqrt{B^{}_{}{}^{}2B^{^{\prime \prime }}B_0/2}`$ and $`\beta =B^{^{\prime \prime }}/2`$. We used the typical values adopted experimentally, such as $`B_0=0.9\mathrm{G}`$, $`B^{^{}}=180\mathrm{G}/\mathrm{cm}`$, and $`B^{^{\prime \prime }}=170\mathrm{G}/\mathrm{cm}^2`$. The corresponding trapping frequencies are $`\omega _r=2\pi \times 171`$ Hz in the radial direction, and $`\omega _z=2\pi \times 11.8`$ Hz in the axial direction. The collisional parameters for the trapped $`\mathrm{F}=1,\mathrm{m}_\mathrm{F}=1`$ state were set as follows : s-wave scattering length $`a=5.8`$ nm, loss rate constants due to background gas collisions $`G_1=0.01\mathrm{s}^1`$, dipolar relaxation $`G_2=1.5\times 10^{16}\mathrm{cm}^3\mathrm{s}^1`$, and three-body recombination $`G_3=4.3\times 10^{29}\mathrm{cm}^6\mathrm{s}^1`$. We chose the temperature $`T=150\mu `$K, the peak density $`\stackrel{~}{n}(\mathrm{𝟎})=4.0\times 10^{11}`$ cm-3, and the number of trapped atoms $`\stackrel{~}{N}=1.0\times 10^9`$ as the initial evaporative cooling conditions.

The exponential sweep of the rf-field frequency is usually adopted in BEC experiments . We also performed calculations assuming that the frequency, $`\nu `$, is swept exponentially from 35 to 0.67 MHz in 60 s as in Ref.. The comparison between the different sweeps of rf-field frequency, optimized or exponential, will clarify the merit of the present optimization procedure.

Figure 1 shows each sweep of the rf-field frequency. We see that the sweeping time (i.e., the cooling time) in the optimized case is about half that of exponential-sweep cooling for 60 s. The shorter sweeping time in the optimized cooling is advantageous for BEC experiments since the long time evaporative cooling easily becomes inefficient due to the instability of the magnetic trap. Furthermore, it is found that the optimized frequency sweep becomes rapid at the final stage of the cooling process, i.e., after about 27 s, and such a sweep is quite opposite to the exponential one which slows down over time. As shown in the inset of Fig. 1, the optimization procedure accelerates the evaporation around the BEC transition point. Similar acceleration of evaporation at the final stage of optimized cooling has been also demonstrated by the previous studies based on kinetic theory of a classical gas .

In Figs. 2(a) and 2(b), we show the time evolution of the number of trapped atoms $`\stackrel{~}{N}`$ and that of temperature $`T`$, respectively. In the optimized cooling, both $`\stackrel{~}{N}`$ and $`T`$ rapidly decrease at the final stage of cooling due to the rapid evaporation mentioned above. On the other hand, $`\stackrel{~}{N}`$ and $`T`$ decrease monotonically in the exponential cooling case. At the BEC transition point indicated by the dots in the figure, both the number of trapped atoms and temperature are calculated to be $`\stackrel{~}{N}=7.2\times 10^7`$ and $`T=1.2\mu `$K at 29.01 s for the optimized cooling, and $`\stackrel{~}{N}=7.1\times 10^6`$ and $`T=0.56\mu `$K at 48.9 s for the exponential cooling.

Next, we show how the optimization affects the cooling efficiency in the whole evaporative cooling process. In Fig. 3, the cooling trajectories of both cases are plotted in the $`\stackrel{~}{N}`$-$`\stackrel{~}{\rho }`$ plane and several corresponding time values are also added to the cooling trajectories to demonstrate the time evolution. At the early stage of cooling where the loss due to background gas collisions is dominant, both trajectories are quite similar, indicating that the exponential sweep of rf-field frequency is sufficiently good. However, the difference becomes larger as the trajectories approach the quantum degenerate region where the phase-space density exceeds the critical value of $`\rho _c=2.612`$ and Bose statistics of atoms becomes dominant. Clearly, more atoms are lost in conventional exponential-sweep cooling than in the optimized case. For example, we have ten times larger number of atoms in the optimized cooling at the BEC transition point when the trajectory reaches the horizontal line of $`\stackrel{~}{\rho }/\rho _c=1`$. Note that the point where the two trajectories split up in Fig. 3 (i.e., $`\stackrel{~}{\rho }/\rho _c0.01`$) corresponds to the time of around 28 s where the rf-field frequency curves cross in Fig. 1. The sudden increase in phase-space density just after the BEC transition in both trajectories reflects the rapid growth of condensate due to bosonic stimulation.

The number of condensed atoms, $`N_0`$, is largely increased by the optimization. In the optimized cooling, $`N_0`$ increased monotonically and we obtained $`N_0=2.1\times 10^7`$ after 29.04-s cooling. In the exponential-sweep cooling, on the other hand, $`N_0`$ had the maximum value of $`1.74\times 10^5`$ at 57.0 s and then gradually decreased. We finally obtained $`N_0=1.65\times 10^5`$ for 60-s exponential sweep cooling. More than a 100-fold enhancement of the number of condensed atoms was achieved in our present calculations. The lifetime of such a large condensate is evaluated to be less than 100 ms due to the enormous three-body recombination loss, as we will see in Fig. 4. However, it is possible to resolve this problem by adiabatically expanding the trapping potential. We also add that the value of $`N_0=1.7\times 10^5`$ quantitatively agrees with the experimental results reported in Ref. and confirms the validity of our numerical calculations.

The undesirable loss of trapped atoms during the cooling process gives us the physical origin of the rapid frequency sweep in the optimized cooling. Figure 4 shows the time evolution of loss rates of total number of trapped atoms at the final stage of the optimized cooling. It is seen that the loss due to the three-body recombination greatly increases and selectively removes atoms with high density (i.e., atoms with lower energy). This means that the serious heating of the system is caused by three-body recombination loss. In order to get a large number of condensed atoms, evaporative cooling should therefore be accelerated to compensate for such undesirable heating.

Here, we briefly discuss the validity of the quick rethermalization assumption in the rapid evaporation necessary for optimized cooling. The evolution of the system in the evaporative cooling process is well scaled by the characteristic time $`\tau `$ given by the inverse of normalized loss rate as $`\tau =(\dot{N}_{\mathrm{loss}}/\stackrel{~}{N})^1`$. The thermalization of the system, on the other hand, is characterized by the elastic collision time $`\tau _{\mathrm{coll}}=(n\sigma v)^1`$; $`n`$ is the mean atom density, $`\sigma =8\pi a^2`$ is the elastic collision cross section, and $`v=4\sqrt{k_BT/\pi m}`$ is the mean velocity of relative motion of the atoms . At the BEC transition point in the optimized cooling, these characteristic times are calculated to be $`\tau =140`$ ms and $`\tau _{\mathrm{coll}}=0.4`$ ms, respectively. The relation $`\tau \tau _{\mathrm{coll}}`$ confirms that quick rethermalization is still valid during the rapid evaporation that occurs for optimized cooling.

### B Efficient evaporative cooling based on general sweeping functions of rf-frequency

The optimized frequency sweep in Fig. 1 is the complicated function of the experimental parameters such as the temperature, the density of gas, the atomic collisional parameters, the magnetic trapping potential, and so on. This means that we need to calculate the optimized cooling for each BEC experiment depending on its experimental conditions. Here we try to demonstrate the efficient evaporative cooling based on general sweeping functions of rf-field frequency, which is useful for many experimentalists who want to achieve large condensates.

From the results in Figs. 3 and 4, the exponential sweep is highly efficient until the three-body recombination loss becomes dominant. The rapid frequency sweep at the final stage of the optimized cooling in Fig. 1 might be approximated by a linear function. Thus it is expected that a combination of an exponential sweep of rf frequency and a rapid linear one realizes the high cooling efficiency. To test this cooling method, we carried out the time-evolution calculation by applying the same numerical parameters used in the optimization calculation of the proceeding subsection. As shown in Fig. 5(a), the frequency is first swept exponentially from 35 MHz to 2.26 MHz for 27 s and then swept linearly to 0.67 MHz for three different times of 1 s, 5 s, and 10 s respectively. The exponential curve in Fig. 5(a) is in conformity with the corresponding part of the dashed line in Fig. 1.

Figure 5(b) shows the cooling trajectories in the $`\stackrel{~}{N}`$-$`\stackrel{~}{\rho }`$ plane for the evaporative cooling using the sweeping functions of Fig. 5(a), together with those for the optimized and exponential cooling plotted in Fig. 3. We can see that the linear frequency sweep at the final stage of the evaporative cooling strongly improves the cooling efficiency as expected. Especially, the cooling trajectory in the case of 27-s exponential and 1-s linear frequency sweep (dashed line) is very close to the optimized one (solid line). The final number of condensed atoms in this cooling method is calculated to be $`N_0=1.8\times 10^7`$ and comparable to the value in the optimized cooling such as $`N_0=2.1\times 10^7`$. One can achieve the evaporative cooling with sufficiently high efficiency by an initial exponential and a sequential rapid linear sweep of rf-field frequency.

### C Experimental realization of efficient evaporative cooling

Recently, a group at Gakushuin University performed the efficient evaporative cooling for 87Rb atoms with F=2, $`\mathrm{m}_\mathrm{F}=2`$ trapped state . They applied a frequency sweep similar to the optimized one which we calculated according to their experimental situations. The time of evaporative cooling was divided into three segments and in each segment the optimized sweeping function was approximated by the linear one: the frequency was swept linearly like 32 $``$13 $``$ 2 $``$1.1 MHz for 10, 10, and 5 s respectively. A large number of condensed atoms, $`N_0=1.4\times 10^6`$, were finally obtained after this 25-s evaporative cooling starting from the experimental conditions such as initial temperature $`T=380\mu K`$, initial total number of atoms $`N=4\times 10^8`$, and initial phase-space density $`\rho =5.8\times 10^8`$.

On the other hand, we theoretically investigated this evaporative cooling. It was found that the cooling trajectory through phase space is close to that of optimized cooling and the final 5-s linear sweep of rf-frequency enhances the phase space density more than four orders of magnitude. The number of condensed atoms for this 25-s evaporative cooling was calculated to be $`1.7\times 10^6`$. Our calculated results agree well with the experiment at Gakushuin University.

## IV Conclusion

We have theoretically studied the optimization of evaporative cooling on the basis of quantum kinetic theory of a Bose gas, focusing on BEC experiments with rubidium atoms. The calculated results indicate that the acceleration of evaporative cooling around the BEC transition point is very effective against serious three-body recombination loss. The number of condensed atoms is expected to be largely increased by this optimization procedure. Furthermore, it is shown that highly efficient evaporative cooling can be achieved by an initial exponential and then a rapid linear sweep of frequency. These results are consistent with the experiments and provide a useful guideline to experimentalists aiming to achieve large condensates in 87Rb atoms.

In the present work, we optimized the sweeping function of the rf-field frequency only. A more complete optimization together with the adiabatic change of trapping potential curvature might produce a large and long-lived condensate. On the other hand, our theory neglects the influence of gravity which reduces the cooling efficiency at low temperatures . An interesting future work would be to include the gravity in our present kinetic analysis and study how the optimized cooling trajectory is modified by it.

## Acknowledgments

We are grateful to Y. Yoshikawa, K. Araki, T. Kuwamoto, T. Hirano of Gakushuin University, M. Kozuma of the Tokyo Institute of Technology, T. Sugiura, T. Kuga of the University of Tokyo, M. W. Jack and T. Hong of NTT Basic Research Laboratories for valuable discussions.


text:
# Knotted Solitons

## 1. Introduction

The term “soliton” entered applied mathematics in 1965. It was coined by M. Kruskal and N. Zabusky for a special solution of nonlinear Korteweg-de Vries (KdV) equation, depicting solitary wave . Use of convention of particle physics language shows that the author envisioned the particle-like interpretation for the object which they called soliton.

The attention of mathematical physicists to solitons was attracted after the inverse scattering method was devised by G. Gardner, J. Green, M. Kruskal and R. Miura for solving the KdV equation  and its extension to Nonlinear Schroedinger Equation was found by V. Zakharov and A. Shabat . In the 1970’s this method and its generalizations got a lot of attention and involved quite a few active participants. Rather complete review can be found in . In the end of that decade the quantum variant of the method was constructed and particle-like interpretations of solitons got natural confirmation in terms of quantum field theory, see review in . The mathematical structure of the quantum method was deciphered in pure algebraic way leading in the 1980’s to notion of quantum groups with new applications in pure mathematics and mathematical physics.

The value of solitons for the particle physics consists in the possibility of going beyond the paradigm of the perturbation theory. Indeed, soliton solutions correspond to full nonlinear equations and disappear in their linearized form. Characteristic for solitons is that they interact strongly if the excitations of the linearized fields interact weakly. Another attractive feature is the appearance of elementary topological characteristics for solitons topological charges.

This was understood already in the middle of 1970’s by several groups as I underlined in my lectures, when I was touring USA in 1975 (see e. g. ). However, all these tantalizing features of solitons had one very important drawback: the developed methods applied only in 1 + 1 dimensional space-time.

Naturally the search for 3 + 1 dimensional generalizations became eminent. General considerations showed that many features of 1 + 1 dimensional systems, such as complete integrability and existence of exact many-particle solutions could not be generalized to 3 + 1 dimensions. However, the mere existence of “one-particle” soliton solutions was not excluded. One particular example was introduced by Skyrme in a pioneer paper  long before the soliton rush. Another example was proposed by G. ’t Hooft and A. Polyakov in 1975 . In the following years their solutions got real applications in nuclear and high energy physics.

In both examples the solitons are “point-like”, namely their deviation from the vacuum is concentrated around central point in space. Moreover they have spherical symmetry, allowing the separation of variables in the corresponding equations, reducing them to ODE, which one can treat on a usual PC.

In my lectures , already mentioned, I proposed one more possibility for 3 + 1 dimensional system, allowing solitons. The model, which superficially looks as a slight modification of Skyrme model, has quite distinct features. The center of the would-be soliton is not a point, but a closed contour, possibly linked or knotted. However my proposal remained unnoticed. The reason was evident: the maximal symmetry for such a soliton is axial, reducing 3-dimensional nonlinear PDE to 2-dimensional one. Existing computers were not able to treat such a problem. Thus my proposal was in slumber for 20 years until my colleague Antti Niemi became interested and agreed to sacrifice a year to learn computing and devising the programm. The preliminary results published in  attracted the attention of professionals in computational physics and now we have an ample evidence, confirming my proposal , .

The development which followed showed unexpected universality of my model. The variables, used in it, were shown to enter the list of degrees of freedom for several systems, having realistic physical applications , .

In this talk I shall describe all these developments in detail. First I shall introduce the model, then briefly discuss its numerical treatment and finish with the description of the applications.

## 2. The field configurations and Hopf invariant

The space time is 4-dimensional Minkowski space $`M`$ with linear coordinates $`x^\mu `$, $`\mu =0,1,2,3`$, $`x^0`$ being time and $`x^k`$, $`k=1,2,3`$ space variables. The field $`\stackrel{}{n}(x)`$ is defined on $`M`$ and has values on 2-dimensional sphere $`𝕊^2`$:

$$\stackrel{}{n}:M𝕊^2.$$

The boundary condition on spatial infinity is introduced

$$\stackrel{}{n}|_{r=\mathrm{}}=\stackrel{}{n}_0,$$
(2.1)
where $`r=((x^1)^2+(x^2)^2+(x^3)^2)^{1/2}`$ and $`\stackrel{}{n}_0`$ is a fixed vector, e.g. corresponding to the north pole

$$\stackrel{}{n}_0=(0,0,1).$$

We shall consider mostly the time independent configurations, corresponding to a soliton at rest. The boundary condition (2.1) effectively compactifies the space $`^3`$, turning it into sphere $`𝕊^3`$, thus the stationary configurations realize the map

$$\stackrel{}{n}:𝕊^3𝕊^2,$$
(2.2)
which are known to be classified by Hopf invariant, sort of topological charge.

In general, the density of topological charge is the zero component $`J_0`$ of the current $`J_\mu `$, which is conserved

$$_\mu J_\mu =0$$

independently of the equations of motion. Mathematically it is more natural to use the 3-form $`J`$ dual to 1-form $`J^{}=J_\mu \mathrm{d}x^\mu `$ and define the topological charge as an integral of $`J`$ over space section

$$Q=_^3J.$$

In our case the 3-form $`J`$ is constructed as follows. The pull-back of the volume 2-form on $`𝕊^2`$ via map (2.2) defines the closed 2-form on the space time

$$H=H_{\mu \nu }\mathrm{d}x^\mu \mathrm{d}x^\nu ,$$

where antisymmetric tensor $`H_{\mu \nu }`$ is expressed via field configuration $`\stackrel{}{n}(x)`$ as follows

$$H_{\mu \nu }=(_\mu \stackrel{}{n}\times _\nu \stackrel{}{n},\stackrel{}{n}).$$
(2.3)
Here I use usual notations of vector analysis in 3-space. In fact $`H`$ is exact

$$H=\mathrm{d}C$$

and current 3-form is given by

$$J=\frac{1}{4\pi }HC.$$

In more detail, we have the relations

$$H_{ik}=_iC_k_kC_i$$

and

$$Q=\frac{1}{4\pi }\epsilon _{ikj}H_{ik}C_jd^3x.$$

For regular configurations $`Q`$ gets integer values. This integer has a nice interpretation in the description of which I shall use the terminology of magnetostatic.

Tensor $`H_{ik}`$ can be interpreted as a field strength of the stationary magnetic field in Maxwell theory. The corresponding lines of force are defined via equations

$$\frac{d}{ds}x_i=\frac{1}{2}\epsilon _{ikj}H_{kj},$$

where $`s`$ is a local parameter along the line. It is easy to see that components of $`\stackrel{}{n}(x)`$ along these lines are constant

$$\frac{d}{ds}\stackrel{}{n}(x)=0,$$

giving two “integrals of the motion”. In other words, the Maxwell lines of force are the preimages of points on $`𝕊^2`$ under the map (2.2). Hopf invariant is the intersection number of any pair of such lines.

All these facts are well known and can be found in textbooks (see e.g. ). However I decided to include them into my text to make it more selfcontained.

## 3. The dynamical model

I introduce the dynamical model by giving the relativistic action functional

$$𝒜=a(_\mu \stackrel{}{n})^2d^4x+b(H_{\mu \nu })^2d^4x.$$

In the usual convention of high-energy physics $`𝒜`$ is dimensionless, so the parameter $`a`$ has dimension $`[\text{length}]^2`$ and parameter $`b`$ is dimensionless. Corresponding static energy $`E`$ has the same form as $`𝒜`$ with space-time coordinates substituted by space variables only

$$E=a(_k\stackrel{}{n})^2d^3x+b(H_{ik})^2d^3x.$$
(3.4)
and has proper dimension $`[\text{length}]^1`$. The structure of $`E`$ is similar to that of Skyrme model, where the field variable having values in $`𝕊^3`$ is used and corresponding topological charge is just a degree of map.

Usual check based on the scale transformation is favorable for (3.4) in the same way as in Skyrme model. Indeed

$$E=E_2+E_4,$$

where $`E_2`$ and $`E_4`$ are quadratic and quartic in derivatives of $`\stackrel{}{n}`$ correspondingly. Thus under scaling $`x\lambda x`$ we have

$$E_2\lambda E_2,E_4\frac{1}{\lambda }E_4$$

and the virial theorem states that on the minimal configuration (if any)

$$E_2=E_4.$$

In terms of quantum theory $`E_2`$ has a standard interpretation of the energy of nonlinear sigma-model whereas $`E_4`$ is rather exotic. On the contrary in the magnetic interpretation, mentioned above, $`E_4`$ is a natural term — it is just the Maxwell magnetic energy, whereas the nature of $`E_2`$ is not that clear. However in what follows the presence of both $`E_2`$ and $`E_4`$ is crucial for the existence of solitons as the scaling argument already showed.

This is confirmed also by a beautiful estimate, obtained in

$$Ec|Q|^{3/4},$$

which shows that in the sectors with nonzero $`Q`$ the minimum of energy is strongly positive. Thus the soliton solutions should be obtained by the minimizing of $`E`$ with $`Q0`$ fixed.

Unfortunately until now there exists no proof of the compactness of the minimizing sequence in general case. For the case of axial symmetry uncouraging result are obtained in . So the main argument for the evidence of solitons in my model is based on the numerical work.

## 4. Numerical work

To find the numerical evidence of the existence of localized solitons it is not necessary to solve the nonlinear elliptic equation, obtained by the variational principle

$$\frac{\delta E}{\delta \stackrel{}{n}}=0.$$
(4.5)
Instead one can introduce an auxiliary time $`s`$ and consider the parabolic equation

$$\frac{d\stackrel{}{n}}{ds}=\frac{\delta E}{\delta \stackrel{}{n}}$$
(4.6)
with initial value $`\stackrel{}{n}_{\text{init}}`$

$$\stackrel{}{n}|_{s=0}=\stackrel{}{n}_{\text{init}}$$

being a configuration with the prescribed Hopf invariant. Of course to simulate (4.6) on the computer one is to use some difference scheme. If for large $`s`$ solution of (4.6) stabilizes it gives the solution of (4.5). In other words the soliton appears as an attractor for the evolution equation.

There are of course many important practical details how to discretize equation, how to take into account the normalization condition $`\stackrel{}{n}^2=1`$ and how to choose the initial configuration $`\stackrel{}{n}_{\text{init}}`$. The main papers  and  use different prescription for all this, however quite satisfactorily the final results coincide. I refer to these papers for the details of calculations and proceed to describe the results.

The iterative process was performed for the configuration with $`Q=1,2,\mathrm{}7`$. The results are as follows: for $`Q=1`$ and $`Q=2`$ the solutions are axial symmetric. The center line — the preimage of the point $`n=(0,0,1)`$ — is a circle. The surfaces $`n_3=\alpha `$, $`1<\alpha <1`$ are toroidal and they are spanned by the lines of force wrapping the torus once for $`Q=1`$ and twice for $`Q=2`$. In other words the soliton can be viewed as a filament of lines of force, closed and twisted once or twice.

The solution for $`Q=3`$ is similar but not axial symmetric any more, the corresponding “cable” is warped. For $`Q=4`$ the soliton is a link of two twisted filaments. Especially beautiful case is $`Q=7`$, the central line of the corresponding soliton is a trefoil knot.

The file  contains impressive moving pictures illustrating the convergence of the iterations. I plan to show these movies in my talk, but unfortunately can not do it in a written text.

Thus the numerical work gives the compelling evidence of the existence of string-like solitons in my model. There remains an important mathematical challenge to provide the rigorous existence theorem. Another interesting direction is to find some realistic applications of the model. Some progress in this direction is already obtained and I proceed to the description of it.

## 5. The applications

Nonlinear fields such as $`\stackrel{}{n}(x)`$ rarely enter the dynamical models directly. However they can appear as a part of degrees of freedom in a suitable parameterization of the original fields. For example in condensed matter theory one uses the complex valued functions $`\psi _\alpha (x)`$, $`\alpha =1,\mathrm{},N`$ to describe the density amplitudes of Bose gas or the gap function of superconductor. The interaction supports the configurations, for which

$$\rho ^2=\underset{\alpha =1}{\overset{N}{}}|\psi _\alpha |^2$$
(5.7)
is nonvanishing. In this case it is natural to use $`\rho `$ as one of the independent variables and introduce new variables

$$\chi _\alpha =\psi _\alpha /\rho $$

such that

$$\underset{\alpha =1}{\overset{N}{}}|\chi _\alpha |^2=1.$$
(5.8)
In this way the compact target (I use the slang of the string theory) $`𝕊^{2N1}`$ appears.

When magnetic interaction is introduced the invariance with respect to the phase transformation

$$\psi _\alpha (x)e^{i\lambda (x)}\psi _\alpha (x)$$

is invoked. This means that the target $`𝕊^{2N1}`$ changes

$$𝕊^{2N1}𝕊^{2N1}/U(1).$$

In particular for $`N=2`$ we have

$$𝕊^3/U(1)𝕊^2$$

and the field $`\stackrel{}{n}(x)`$ naturally appears. Quite satisfactorily the tensor $`H_{ik}`$ also emerges as a contribution to the magnetic field strength.

Let us illustrate it in more detail. From the beginning we shall treat the stationary system, so no electric field will be used.

The magnetic field is described in a usual way by means of the vector potential $`A_k(x)`$ and its interaction with $`\psi `$-fields is introduced via covariant derivatives

$$_k\psi =_k\psi +iA_k\psi .$$

The energy density (of Landau-Ginsburg-Gross-Pitaevsky type) looks as follows

$$E=\underset{\alpha =1}{\overset{2}{}}|_k\psi _\alpha |^2+\frac{1}{2}F_{ik}^2+V(|\psi _\alpha |),$$
(5.9)
where

$$F_{ik}=_iA_k_kA_i$$

is the field strength of the magnetic field. The energy is invariant with respect to the gauge transformations

$$\psi _\alpha e^{i\lambda }\psi _\alpha ,A_kA_k_k\lambda .$$

We shall make the change of the field variables so that only gauge invariant ones will remain. For that observe that the first term in the RHS of (5.9) is a quadratic form in $`A`$

$$\underset{\alpha =1}{\overset{2}{}}|_k\psi _\alpha |^2=\underset{\alpha =1}{\overset{2}{}}|_k\psi _\alpha |^2+A_kJ_k+\rho ^2A_k^2,$$

where we use variable $`\rho `$ from (5.7) and introduce current

$$J_k=i\underset{\alpha =1}{\overset{2}{}}(\overline{\psi }_\alpha _k\psi _\alpha _k\overline{\psi }_\alpha \psi _\alpha ).$$

It is easy to check, that under the gauge transformations the current $`J_k`$ changes as follows

$$J_kJ_k+2\rho ^2_k\lambda ,$$

so that the sum

$$C_k=A_k+\frac{1}{2\rho ^2}J_k$$

is gauge invariant. We shall use this variable instead of $`A_k`$. Another gauge invariant combination is given by the quadratic form

$$\stackrel{}{n}=(\overline{\chi }_1,\overline{\chi }_1)\stackrel{}{\tau }\left(\begin{array}{c}\chi _1\\ \chi _2\end{array}\right),$$
(5.10)
where $`\stackrel{}{\tau }=(\tau _1,\tau _2,\tau _3)`$ are Pauli matrices

$$\tau _1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\tau _2=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),\tau _3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$

Normalization (5.8) for $`\chi _\alpha `$ implies that $`\stackrel{}{n}`$ is a real unit vector. In fact the map

$$(\chi _1,\chi _2)\stackrel{}{n}$$

defined in (5.10) is a standard Hopf map. Variable $`\stackrel{}{n}`$ is manifestly gauge invariant and the set of variables $`(\rho ,\stackrel{}{n},C_k)`$ is our gauge invariant choice, substituting for the initial set $`(\psi _\alpha ,A_k)`$. The energy density can be explicitly expressed via $`\rho ,\stackrel{}{n}`$ and $`C_k`$ as follows

$$E=(_k\rho )^2+\rho ^2((_k\stackrel{}{n})^2+C_k^2)+\frac{1}{2}(_kC_i_iC_k+H_{ik})^2+v(\rho ,n_3).$$

The most notable feature is the appearance of the tensor $`H_{ik}`$, defined in (2.3). The model, described in the main text, emerges if we put $`\rho =\text{const}`$ and $`C=0`$. Hopefully nontrivial $`\rho `$ and $`C`$, at least confined to some range, do not spoil the soliton picture. This problem is under discussion now, see , .

Let us stress, that the use of two fields $`\psi _\alpha `$, $`\alpha =1,2`$ is most essential in this example. If $`N=1`$ only variables $`\rho `$ and $`C`$ remain after the reduction, similar to just described. If $`N>2`$ the $`CP(N1)`$ field generalizing $`\stackrel{}{n}`$ has no topological characteristics.

Another application, considered recently , deals with the parameterization for the $`SU(2)`$ Yang-Mills field $`A_\mu ^a(x)`$, $`\mu =0,1,2,3`$, $`a=1,2,3`$. The Yang-Mills Lagrangian is invariant with respect to the nonabelian gauge transformations

$$\delta A_\mu ^a=_\mu \epsilon ^a+f^{abc}A_\mu ^b\epsilon ^c.$$

However in some treatments one reduces this invariance by the partial gauge fixing to the abelian one

$$\delta B_\mu =i\epsilon B_\mu ,\delta A_\mu ^3=_\mu \epsilon ,$$

where $`B_\mu =A_\mu ^1+iA_\mu ^2`$ is a complex vector field. I shall not discuss the reason for this reduction here and proceed assuming that it is done. Observe, that two vector fields $`A_\mu ^1`$, $`A_\mu ^2`$ in generic situation define a plane in Minkowski space and introduce an orthonormalized basis in this plane $`e_\mu ^\alpha `$, $`\alpha =1,2`$.

$$e_\mu ^\alpha e_\mu ^\beta =\delta _{\alpha \beta }.$$

Let $`e_\mu =e_\mu ^1+ie_\mu ^2`$. The basis is defined up to rotation

$$e_\mu e^{i\omega }e_\mu .$$

The fields $`B_\mu `$ can be written in terms of this basis as

$$B_\mu =\psi _1e_\mu +\psi _2\overline{e}_\mu $$

and thus two complex valued fields $`\psi _1`$ and $`\psi _2`$ appear. The situation becomes quite similar to the previous example and indeed in  the complete parameterization of the Yang-Mills variables is introduced with appearance of $`\stackrel{}{n}`$-field and corresponding $`H`$-tensor. This is an indication that the Yang-Mills theory can have string-like excitations. However the situation is not that simple. The classical Yang-Mills theory is conformally invariant and has no dimensional parameters. Thus no hope for the localized regular classical solution exists. Nevertheless this complication could be lifted by quantum corrections. The famous “dimensional transmutation”, which leads to the appearance of dimensional parameter in quantum effective action, could favor the nonvanishing value of the corresponding $`\rho `$-variable. All these considerations at the moment are rather speculative and need much more work to become reasonable. Personally I am quite impressed by this possibility and continue to work on it.

## 6. Conclusions

I think that the topic of my talk is quite instructive. It connects different domains in mathematics and mathematical physics: nonlinear PDE, elementary topology, quantum field theory, numerical methods. It illustrates the essential unity of mathematics, theoretical and applied. Finally it could lead to the realistic physical applications. For all these reasons I decided to present it to the ICM2002.