text:
# Untitled Document

Efficient collinear third-harmonic generation in a single two-dimensional nonlinear photonic crystal

Todor Karaulanov, Solomon M. Saltiel
Faculty of Physics, University of Sofia, 5 J. Bourchier Boulevard, Sofia 1164, Bulgaria,
Tel +359-2 62 15 30 e-mails: tod\_kulan@hotmail.com, saltiel@phys.uni-sofia.bg

## Abstract

We propose novel multi-phase-matched process that starts with generation of a pair of symmetric second-harmonic waves. Each of them interacts again with the fundamental wave to produce two constructively interfering third harmonic waves collinear to the fundamental input wave.

Efficient collinear third-harmonic generation in a single two-dimensional nonlinear photonic crystal

Todor Karaulanov, Solomon M. Saltiel
Faculty of Physics, University of Sofia, 5 J. Bourchier Boulevard, Sofia 1164, Bulgaria,

Nonlinear photonic crystals have attracted much interest lately because of their unique properties and many potential applications in nonlinear optics and optics communications. Two dimensional nonlinear photonic crystals (2DNPC)  allow realization of simultaneous phase matching of several processes , however in most of the cases generated wave is noncollinear to the input fundamental wave that was considered as a disadvantage for 2DNPC. Here we propose for the first time phase-matching geometry for third harmonic (TH) generation in 2DNPC, for which generated TH wave is collinear to the input wave. In addition to that the proposed interaction at low input intensities is four times more efficient than the conventional schemes.

Reciprocal lattice of the 2DNPC with the main reciprocal wave-vectors $`G_a,G_b,G_a`$ are shown on the Fig. 1. Multiple phase matching geometry is presented on Fig. 2. The process starts with generation of pair second harmonic waves with wave vectors $`k_{2u},k_{2d}`$. Phase matching is assured by the reciprocal fundamental vectors $`G_a,G_a`$ . Each second harmonic wave interacts again with the fundamental wave via pair phase matched interaction with participation of the reciprocal vectors generating this way pair collinear TH waves with one and the same wave vector $`k_3`$. The reciprocal vectors that phase match the processes of sum frequency mixing $`\omega +2\omega =3\omega `$ are combination of the two fundamental reciprocal vectors $`G_c=NG_bG_a`$ ; $`G_c=NG_bG_a`$. The two TH waves interfere constructively in the direction of the input fundamental wave leading this way to higher overall TH efficiency. The process is described by system of 5 differential equations, that have been solved in several approximation. Design and optimization of TH frequency converters made from 2DNPC implemented in LiNbO3 and LiTaO3 is presented. The frequency converter can be also useful for testing the homogeneity of the sample. The better is the homogeneity the higher will be TH efficiency.
1. V. Berger, Phys. Rev. Lett. 81, 4136 (1998).
2. S. Saltiel and Yu. S. Kivshar, Opt. Lett. 25, 1204 (2000).


format: text

text:
# Statistical Dynamics of On-line Independent Component Analysis

## 1 Introduction

On-line learning algorithms are often used for dealing with very large data sets or in dynamic situations in which data is changing according to a non-stationary process. Independent component analysis (ICA) is often applied under one or both of these conditions and a number of on-line ICA algorithms have been developed (see, e.g. Hyvärinen, 1999). In on-line learning the model parameters are updated after the presentation of each training example. Although there is good understanding of the asymptotic properties of on-line learning, much of the learning takes place far from the asymptotic regime and here a theoretical understanding of the process is often poor. Training methods are typically based on experimental observations in the absence of a successful theoretical analysis of on-line learning in this regime. The efficiency of on–line training is very sensitive to the choice of training parameters such as the learning rate and this dependence can slow down learning and influence the ability of learning to converge successfully to desired states. A deeper theoretical understanding of the on-line learning process is needed to provide reliable methods for setting the parameters and optimising the training process.

Most classical theoretical results on on-line learning are from stochastic approximation theory (Kushner and Clark, 1978). For a review of recent advances and modern theoretical approaches, see e.g. Saad (1998). Theories describing on-line learning have mainly been developed in two different asymptotic regimes. Most work has been done in the limit of the long times, in which case the model parameters are close to a stable fixed point of the learning dynamics. Here one can work out the asymptotic distribution of the parameters for constant learning rate or study their convergence to a fixed point as the learning rate is reduced according to some annealing schedule. The conditions under which the parameters converge at an optimal rate are quite well understood (White, 1989). Work has also been carried out in the limit of small learning rates, in which case the dynamics can be shown to follow the mean gradient or flow globally (Kushner and Clark, 1978). Unifying these two strands to some extent one can study the long-time global behaviour under an annealed learning schedule (Kushner, 1987).

It is unclear how much practical relevance the classical asymptotic limits have, since often in practical applications learning is not asymptotic in the learning times or in the learning rate. In this work we pursue a different type of asymptotic analysis by considering the limit of large system size. This limit has been studied extensively by researchers applying statistical mechanics methods to the study of learning systems (Engel and Van den Broeck, 2001). For example, the dynamics of gradient decent in multilayer perceptrons (MLPs) (Biehl and Schwarze, 1995, Saad and Solla, 1995) and the dynamics of Sanger’s principal component analysis (PCA) algorithm (Biehl, 1994, Biehl and Schlösser, 1998) have been studied in this limit. In these examples the dynamics displays interesting and non-trivial transient behaviour with unstable, sub-optimal fixed points appearing due to a symmetry in the parameter space of the models. The dynamics of natural gradient learning in MLPs has been studied using these techniques, showing that the transient fixed point becomes less serious and transient learning performance is improved compared to standard online gradient descent (Rattray et al., 1998, Rattray and Saad, 1999). Since a natural gradient algorithm has been shown to provide an asymptotically efficient on-line learning algorithm for ICA (Amari et al., 1996, Amari, 1998) we are interested in whether the statistical mechanics approach developed here can help understand its transient performance.

We have recently developed a novel theoretical framework for studying the dynamics of on-line ICA in the limit of large data dimension (Rattray, 2002, Rattray and Basalyga, 2002). We study a Hebbian learning rule, which is a simple and popular algorithm for on-line ICA with nice stability conditions (Hyvärinen and Oja, 1998) and is closely related to a popular fixed-point batch algorithm (Hyvärinen and Oja, 1997). The algorithm is particularly amenable to analysis in the limit of large input dimension because it can be used to extract a small number of independent components from high-dimensional data. We will see later that this allows the dynamics to be represented by a relatively small number of variables when the data dimension becomes large. As well as studying the standard version of the algorithm we also develop a natural gradient variant and study its dynamics. Because the parameter space is constrained to the Stiefel manifold of orthogonal rectangular matrices, the standard equivariant or natural gradient ICA algorithms due to Cardoso and Laheld (1996) and Amari et al. (1996) are not appropriate here. Instead it is possible to use the ideas developed by Edelman et al. (1999) in order to construct an algorithm which follows the gradient on a Stiefel manifold. Our variant uses the gradient defined by Amari (1999) but includes the orthogonalisation term from Hyvärinen and Oja (1998) in order to keep the model parameters on the Stiefel manifold.

In order to obtain a tractable ICA model, we consider an idealised data set in which a small number of non-Gaussian sources are mixed into a large number of Gaussian sources (Rattray, 2002). By using the methods of statistical mechanics, we provide a solution to the dynamics of Hebbian ICA in the limit of large input dimension. This generalises on previous results (Rattray, 2002, Rattray and Basalyga, 2002), which were limited to the simplest single source case (except for the late time asymptotics, which were solved for the general case). We find that the transient dynamics of Hebbian ICA can be described as a stochastic process in which the system moves through a sequence of metastable fixed points, each of which can be described as a multi-dimensional diffusion in polynomial potential. By solving the dynamics of the multi-source case we can characterise the symmetry breaking process which is critical to performance of the learning process.

It is interesting to observe that the dynamics of ICA is very different from the dynamics of learning in MLPs and in PCA studied previously (Saad and Solla, 1995, Biehl and Schlösser, 1998). In these cases the dynamics was observed to be “self-averaging” so that it followed a smooth trajectory in the limit of large data dimension $`N`$ and the learning happened on an $`O(N)`$ time-scale. The ICA dynamics displays significant fluctuations even in this limit and learning occurs on a much slower time-scale, typically requiring of the order of $`N^2`$ or $`N^3`$ iterations depending on the details. We find that the natural gradient variant of the Hebbian algorithm has very nice asymptotic convergence properties with uniform convergence in the case of equal source statistics. However, the algorithm is shown to have equivalent transient performance to the standard algorithm.

This paper is organised as follows. The data model is described in Section 2. The on-line Hebbian ICA algorithm is introduced in Section 3. In Section 4 we introduce the macroscopic variables which provide a compact description of the dynamics for large $`N`$. In Section 5 we show that the learning dynamics near a sub-optimal fixed point close to the initial conditions can be considered as a diffusion process. The transient dynamics through a sequence of metastable states is analysed in Section 6. In Section 7 we introduce the natural gradient version of Hebbian ICA algorithm and study its dynamics. General conclusions are made in Section 8.

## 2 Data Model

We consider the following idealised linear data model introduced in Rattray (2002). The $`N`$-dimensional data $`𝒙`$ is generated from a noiseless linear mixture of a small number $`M`$ of non-Gaussian sources $`𝒔`$ and a large number $`NM`$ of uncorrelated Gaussian components, $`𝒏𝒩(\mathrm{𝟎},𝑰_{NM})`$,

$$𝒙=𝑨\left[\begin{array}{c}𝒔\\ 𝒏\end{array}\right]=𝑨_s𝒔+𝑨_n𝒏,$$
(1)
where $`𝑨=[𝑨_s𝑨_n]`$ is the $`N\times N`$ mixing matrix, $`𝑰_N`$ denotes an $`N\times N`$ identity matrix and $`𝒩(𝒂,𝚺)`$ denotes a Gaussian distribution with mean $`𝒂`$ and covariance matrix $`𝚺`$. Without loss of generality it can be assumed that the sources each have unit variance. In order to apply the Hebbian ICA algorithm we assume also that the data is already sphered, i.e. the data has zero mean and an identity covariance matrix. This can be achieved for on-line learning by an adaptive sphering algorithm, such as the one introduced by Cardoso and Laheld (1996). These model assumptions lead to the constraint that $`𝑨`$ must be an orthogonal matrix, e.g.

$`\left[𝑨_s𝑨_n\right]\left[\begin{array}{c}𝑨_s^\mathrm{T}\\ 𝑨_n^\mathrm{T}\end{array}\right]`$ $`=`$ $`𝑨_s𝑨_s^\mathrm{T}+𝑨_n𝑨_n^\mathrm{T}=𝑰,`$ (4)
$`\left[\begin{array}{c}𝑨_s^\mathrm{T}\\ 𝑨_n^\mathrm{T}\end{array}\right]\left[𝑨_s𝑨_n\right]`$ $`=`$ $`\left[\begin{array}{cc}𝑨_s^\mathrm{T}𝑨_s& 𝑨_s^\mathrm{T}𝑨_n\\ 𝑨_n^\mathrm{T}𝑨_s& 𝑨_n^\mathrm{T}𝑨_n\end{array}\right]=\left[\begin{array}{cc}𝑰& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝑰\end{array}\right].`$ (11)
## 3 On-line Hebbian ICA Learning Rule

The goal of ICA is to find the de-mixing matrix $`𝑾`$ such that the projections,

$$𝒚𝑾^\mathrm{T}𝒙,$$
(12)
will coincide with the non-Gaussian sources $`𝒔`$ up to scaling and permutations. The best possible solution is one in which the $`K`$ projections will learn as many as possible of the $`M`$ non-Gaussian sources. Note that specialisation of each projection to a particular source mostly depends on the details of the initial conditions.

We consider a simple Hebbian learning rule (Hyvärinen and Oja, 1998), which extracts non-Gaussian sources from the data by maximising some measure of non-Gaussianity of the projections. The change of the $`N\times K`$ de-mixing matrix $`𝑾`$ at the each learning step is given by,

$$𝚫𝑾=\eta 𝒙\mathit{\varphi }(𝒚)^\mathrm{T}𝝈+\alpha 𝑾(𝑰𝑾^\mathrm{T}𝑾).$$
(13)
Here, $`\eta `$ is the learning rate and $`𝝈`$ is a $`K\times K`$ diagonal matrix with elements

$$\sigma _{ii}=\text{sign}\left(\text{E}_{s_i}[s_i\varphi (s_i)\varphi ^{}(s_i)]\right),$$
(14)
which ensures stability of the correct solution as $`y_is_i`$ (Hyvärinen and Oja, 1998). The first term on the right of (13) maximises some measure of non-Gaussianity of the projections. The second term provides orthogonalisation of the de-mixing matrix. The choice of learning rate $`\eta `$ greatly influences the performance of this algorithm. Choosing too large a learning rate results in slow and inefficient learning but choosing too high a value may result in the learning dynamics becoming trapped in a poor solution, as we will see later. The learning dynamics is less sensitive to the choice of the parameter $`\alpha `$ and we set $`\alpha =0.5`$ in simulations. The function $`\mathit{\varphi }(𝒚)=[\varphi (y_i)]`$ is some smooth non-linear function which is applied to every component of the vector $`𝒚`$. An even non-linearity, e.g. $`\varphi (y)=y^2`$, is usually used to detect asymmetric non-Gaussian signals, while an odd non-linearity, e.g. $`\varphi (y)=y^3`$ or $`\varphi (y)=\mathrm{tanh}(y)`$, is used to extract symmetric non-Gaussian signals.

## 4 Learning Dynamics for Large Data Dimension

In the case of high-dimensional data (when $`N`$ becomes very big) it is difficult to analyse the dynamics of the $`N\times K`$ de-mixing matrix. In order to provide a compact description of the system dynamics in the limit $`N\mathrm{}`$, we introduce new macroscopic variables

$`𝑹𝑾^\mathrm{T}𝑨_s\text{ and }𝑸𝑾^\mathrm{T}𝑾,`$ (15)
the dimension of which are $`K\times M`$ and $`K\times K`$ respectively. These overlap matrices contain all necessary information about the relationship between the projections $`𝒚`$ and the sources $`𝒔`$. Therefore the system can be described by a small number of macroscopic quantities in the limit of large $`N`$ as long as $`K`$ and $`M`$ remain small. We will usually order the indices retrospectively so that the dynamics approaches the optimal solution with $`R_{ij}=\delta _{ij}`$ (assuming orthogonal $`𝑾`$ such that $`𝑸=𝑰`$). However, it should be remembered that there are equivalent optima related to this solution by a permutation of indices or changes in sign of the components of $`𝑹`$. In order to learn successfully the algorithm has to break symmetry and specialise to a particular solution.

Using Equation (4) one can show that,

$`𝒚=𝑾^\mathrm{T}(𝑨_s𝒔+𝑨_n𝒏)=𝑹𝒔+𝒛,`$ (16)
where $`𝒛𝒩(\mathrm{𝟎},𝑪)`$ and $`𝑪=𝑸𝑹𝑹^\mathrm{T}`$. In order to analyse the dynamics we will have to compute expectations with respect to the distribution of $`𝒚`$. The covariance matrix $`𝑪`$ is symmetric and positive definite so that we can always write $`𝒛`$ in the form $`𝒛=𝑳𝝁`$ where $`𝝁𝒩(\mathrm{𝟎},𝑰)`$ are Gaussian variables and the matrix $`𝑳`$ can be found by special decomposition. A particularly useful decomposition for our purposes is the Cholesky decomposition (see, e.g. Press et al., 1992) in which case $`𝑪=𝑳𝑳^\mathrm{T}`$ where $`𝑳`$ is lower-triangular with non-zero components satisfying the following recurrence relations

$$L_{kk}=\sqrt{C_{kk}\underset{j=1}{\overset{k1}{}}L_{kj}^2}$$
(17)
and for the $`k`$-th row of $`𝑳`$ we have

$$L_{ki}=\frac{C_{ki}_{j=1}^{i1}L_{ij}L_{kj}}{L_{ii}},$$
(18)
for $`i=1,2,\mathrm{},k1`$.

From Equation (13) we can calculate the changes in $`𝑹`$ and $`𝑸`$ after a single learning step,

$`\mathrm{\Delta }𝑹`$ $`=`$ $`\eta 𝝈\mathit{\varphi }(𝒚)𝒔^\mathrm{T}+\alpha (𝑰𝑸)𝑹,`$ (19)
$`\mathrm{\Delta }𝑸`$ $`=`$ $`\eta 𝝈(𝑰+\alpha (𝑰𝑸))\mathit{\varphi }(𝒚)𝒚^\mathrm{T}+\alpha ^2(𝑰𝑸)^2𝑸`$ (20)
$`+`$ $`\eta 𝝈𝒚\mathit{\varphi }(𝒚)^\mathrm{T}(𝑰+\alpha (𝑰𝑸))+2\alpha (𝑰𝑸)𝑸`$
$`+`$ $`\eta ^2\mathit{\varphi }(𝒚)𝒙^\mathrm{T}𝒙\mathit{\varphi }(𝒚)^\mathrm{T},`$
where we used the constraint in Equation (11) to set $`𝒙^\mathrm{T}𝑨_s=𝒔^\mathrm{T}`$.

The dynamics is not very sensitive to the exact value of $`\alpha `$ as long as $`\alpha \eta `$. We will see later that the learning rate must be chosen very small for large $`N`$ so that typically we will have this situation in practice. As $`\alpha `$ increases relative to $`\eta `$, $`𝑸`$ approaches $`𝑰`$ since the orthogonalisation term in Equation (13) dominates. If one defines $`𝑸𝑰𝒒/\alpha `$ and sets $`\alpha `$ large relative to $`\eta `$ then the fixed point of Equation (20) to leading order is,

$$𝒒=\frac{1}{2}\left[\eta 𝝈\left(\mathit{\varphi }(𝒚)𝒚^\mathrm{T}+𝒚\mathit{\varphi }(𝒚)^\mathrm{T}\right)+\eta ^2N\mathit{\varphi }(𝒚)\mathit{\varphi }(𝒚)^\mathrm{T}\right],$$
(21)
where we have dropped terms lower than $`O(\eta ^2N)`$ and $`O(\eta )`$. Substituting Equation (21) into Equation (19) leads to the following update equation for $`𝑹`$,

$`\mathrm{\Delta }𝑹=\eta 𝝈\left[\mathit{\varphi }(𝒚)𝒔^\mathrm{T}\frac{1}{2}\left(\mathit{\varphi }(𝒚)𝒚^\mathrm{T}+𝒚\mathit{\varphi }(𝒚)^\mathrm{T}\right)𝑹\right]\frac{1}{2}\eta ^2N\mathit{\varphi }(𝒚)\mathit{\varphi }(𝒚)^\mathrm{T}𝑹.`$ (22)
This simplification procedure is an example of adiabatic elimination of fast variables (see, for example, Gardiner, 1985). In a more rigorous treatment one should consider the mean and covariance of $`\mathrm{\Delta }𝑹`$ and $`\mathrm{\Delta }𝑸`$ using the appropriate large $`N`$ scalings which are described in the next sections. One finds that fluctuations in $`𝑸`$ are negligible in the limit of large $`N`$ and therefore the dynamics of $`𝑸`$ can be described by a differential equation in this limit with stable fixed point given by Equation (21).

In Figure 1 we show some typical dynamical trajectories. We observe the following types of fixed points in the $`𝑹`$-dynamics:

* Optimal fixed points (see Figure 1 (a) and (b))

Asymptotically, when $`y_is_i`$, the optimal solution is given by $`R_{ij}^{}=\delta _{ij}`$, which is a fixed point of Equation (22) as $`\eta 0`$. Note that all other possible solutions can be obtained by a trivial permutation of indices and/or changes in sign. Asymptotically the learning rate should be annealed in order to approach this fixed point at an optimal rate. For a detailed account of the asymptotic dynamics of the Hebbian ICA algorithm under an annealed learning rate, see Rattray (2002).
* Sub-optimal fixed points causing trapping near the initial conditions (see Figure 1 (c))

Due to the $`O(\eta ^2)`$ fluctuation term in Equation (22), the algorithm has a special class of sub-optimal fixed points near $`𝑹=\mathrm{𝟎}`$ which causes the presence of a stochastic trapping state near the initial conditions. We will discuss this situation in detail in the next Section.
* Transient fixed points (see Figure 1 (a) and (b))

When $`T<M`$ non-Gaussian sources have been learned, the dynamics can become localised in metastable states somewhere between the initial conditions and the final, optimal fixed point. In this case the fixed points of Equation (22) are

$$R_{ij}^{}=\delta _{ij}\text{I}\left[iT\right],$$
(23)
where

$$\text{I}\left[\text{predicate}\right]=\{\begin{array}{c}1\text{ if predicate is true},\hfill \\ 0\text{if predicate is false}.\hfill \end{array}$$
(24)
Again, similar fixed points can be obtained from (23) by a simple permutation of indices. In this case the system has to leave such metastable states in order to complete learning. We will analyse the dynamics near these metastable states in Section 6.

## 5 Stochastic Trapping State near the Initial Conditions

In similar studies of on-line learning, macroscopic quantities like the overlap $`𝑹`$ usually have a “self-averaging” property such that the variance of these macroscopic quantities tends to zero in the limit $`N\mathrm{}`$ (see e.g. Saad and Solla, 1995, Biehl and Schlösser, 1998). A random and uncorrelated choice for $`𝑨`$ and the initial entries of de-mixing matrix $`𝑾`$ leads us to expect $`𝑹=O(N^{1/2})`$ initially. Larger initial values of $`𝑹`$ could only be obtained with some prior knowledge of the mixing matrix which we will not suppose. In this case one can no longer assume that fluctuations are negligible as $`N\mathrm{}`$. Moreover, as we will see below, the mean and variance of the change in $`𝑹`$ at each iteration are of the same order. That means that the overlap $`𝑹`$ does not self-average and the fluctuations have to be considered even in the limit. Therefore, it is more natural to model the on-line learning dynamics near the initial conditions as a diffusion process (see, for example, Gardiner, 1985, van Kampen, 1992). In order to establish a clear picture of the dynamics we have to choose an appropriate scaling for macroscopic quantities and learning parameters. In the following discussion we set $`𝒓𝑹\sqrt{N}`$ where $`𝒓`$ is assumed to be an $`O(1)`$ quantity.

### 5.1 Diffusion in a Potential

For a diffusion process the probability density $`p(𝒓,t)`$ of a random variable $`𝒓=\left[r_{ij}\right]`$ at time $`t`$ obeys the Fokker-Planck equation

$`{\displaystyle \frac{p(𝒓,t)}{t}}={\displaystyle \underset{ij}{}}{\displaystyle \frac{}{r_{ij}}}\left(A_{ij}p(𝒓,t)\right)+{\displaystyle \frac{1}{2}}{\displaystyle \underset{ijkl}{}}{\displaystyle \frac{^2}{r_{ij}r_{kl}}}\left(B_{ijkl}p(𝒓,t)\right),`$ (25)
where the coefficients $`A_{ij}`$, which are called “drift” coefficients, represent the expectation of the change of variable $`r_{ij}`$ with respect to the stochastic process in question. They can be written as

$$A_{ij}=\text{E}[\mathrm{\Delta }r_{ij}]=\frac{U(𝒓)}{r_{ij}}N^p,$$
(26)
where $`N`$ and $`p`$ are the input dimension and the scaling order for our system and $`U(𝒓)`$ is some differentiable function of $`𝒓`$. This function is analogous to the potential function for the case of a particle undergoing a diffusion in a potential. The coefficients $`B_{ijkl}`$ are the covariance of the change of variable $`r_{ij}`$,

$$B_{ijkl}=\text{Cov}[\mathrm{\Delta }r_{ij},\mathrm{\Delta }r_{kl}]=\text{E}[\mathrm{\Delta }r_{ij}\mathrm{\Delta }r_{kl}]\text{E}[\mathrm{\Delta }r_{ij}]\text{E}[\mathrm{\Delta }r_{kl}]=D_{ijkl}N^p,$$
(27)
where $`D_{ijkl}`$ is called the “diffusion term”. Usually this would be a matrix but notice that in our case the dynamical variables are in a matrix and therefore the diffusion term has four indices. However, we can think of each pair as a single index in a vectorised system. If we do so then for our case the diffusion term can be considered a diagonal matrix of magnitude $`D`$,

$$D_{ijkl}=D\delta _{ik}\delta _{jl},$$
(28)
where $`D`$ is called the “diffusion coefficient”. It is typical for the diffusion process that the mean and covariance of the random variable are of the same order $`O(N^p)`$. This means that, for large $`N`$, the system is equivalent to diffusion in the potential $`U(𝒓)`$ with a characteristic time-scale $`(\delta t)^1=N^p`$.

In our case the potential $`U(𝒓)`$ has a minimum at $`r_{ij}=0`$ surrounded by potential barriers of differing heights. Examples are shown in Figure 2 for an odd and even non-linearity on the left and right respectively. The escape time (the mean first passage time) from the minimum of the potential at $`𝒓=\mathrm{𝟎}`$ is mainly determined by the effective size of the minimal potential barrier $`\mathrm{\Delta }U`$ (see, for example, Gardiner, 1985, van Kampen, 1992),

$$T_{\text{escape}}=A\mathrm{exp}\left(\frac{\mathrm{\Delta }U}{D}\right),$$
(29)
where prefactor $`A`$ is proportional to the characteristic time-scale and depends on the curvature of the potential. Because diffusion through a higher potential barrier is exponentially less likely as the difference between barrier heights increases, the escape time will effectively be inversely proportional to the number of escape points (the number of barriers with the same minimal height $`\mathrm{\Delta }U`$) when the Arrhenius factor $`\mathrm{\Delta }U/D`$ is large.

In the next two sections we will find the form of potential barrier and estimate the escape time from the trapping state for the two main classes of non-linear function $`\varphi (y)`$.

### 5.2 Odd Non-linearity

If we need to extract symmetrical non-Gaussian signals from the data then we have to use an odd non-linearity, e.g. $`\varphi (y)=y^3`$ or $`\varphi (y)=\mathrm{tanh}(y)`$ are common choices. In this case the appropriate scaling for the learning rate will be $`\eta =\nu N^2`$, where $`\nu `$ is an $`O(1)`$ scaled learning rate parameter. After expanding Equation (22) near $`𝒓=\mathrm{𝟎}`$ we obtain the following expressions for the mean and covariance of the change in $`𝒓`$ at each iteration,

$`\text{E}[\mathrm{\Delta }r_{ij}]=\left(\frac{1}{2}\varphi ^2(\mu )\nu ^2r_{ij}+\frac{1}{6}\kappa _4^j\varphi ^{\prime \prime \prime }(\mu )\sigma _{ii}\nu r_{ij}^3\right)N^3+O(N^4),`$ (30)
$`\text{Cov}[\mathrm{\Delta }r_{ij},\mathrm{\Delta }r_{kl}]=\varphi ^2(\mu )\nu ^2\delta _{ik}\delta _{jl}N^3+O(N^4),`$ (31)
where $`\kappa _4^j`$ is the fourth cumulant of the $`j`$-th source distribution (measuring kurtosis) and brackets denote averages over a Gaussian variable $`𝝁𝒩(\mathrm{𝟎},𝑰)`$. In this case the system can be described by a Fokker-Planck equation for large $`N`$ with a characteristic time-scale $`(\delta t)^1=N^3`$. To compute the expectations we have made use of the Cholesky decomposition of $`𝒚`$ as described in Equations (17) and (18). This allows us to remove the dependence of the parameters in the Gaussian averages by writing $`𝒚=𝑳𝝁`$ where the lower diagonal (i.e. non-zero) elements of $`𝑳`$ are found to be,

$$L_{ij}=\delta _{ij}\left(\frac{1}{2}\delta _{ij}\underset{m}{\overset{M}{}}r_{im}^2+\underset{m}{\overset{M}{}}r_{im}r_{jm}\right)N^1+O(N^2)\text{for}ij.$$

The dynamics is equivalent to diffusion in the following potential

$$U(𝒓)=\underset{i=1}{\overset{K}{}}\underset{j=1}{\overset{M}{}}\left(\frac{1}{4}\varphi ^2(\mu )\nu ^2r_{ij}^2\frac{1}{24}\kappa _4^j\varphi ^{\prime \prime \prime }(\mu )\sigma _{ii}\nu r_{ij}^4\right)$$
(32)
with a diagonal diffusion matrix of magnitude $`D=\varphi ^2(\mu )\nu ^2`$. Notice that the potential is a sum of contributions $`U(𝒓)=_{ij}U(r_{ij})`$ which depend only on a single element in $`𝒓`$. Since the diffusion matrix is diagonal, this means that each element of $`𝒓`$ undergoes an independent diffusion process equivalent to the one-dimensional case described by Rattray (2002). The effective size of the potential barriers for each element $`r_{ij}`$ is given by,

$$\frac{\mathrm{\Delta }U(r_{ij})}{D}=\frac{3\varphi ^2(\mu )\nu }{8|\kappa _4^j\varphi ^{\prime \prime \prime }(\mu )|}.$$
(33)
Escape over one such potential barrier results in the corresponding source $`j`$ being learned by projection $`i`$. This breaks the symmetry of the system as one projection specialises to a particular source. Once this happens the system can again become trapped in a metastable state with other sources remaining unlearned. The dynamics in the neighbourhood of this more general class of fixed point is described in Section 6.

The shape of the potential for an example with two sources (with equal kurtosis) and one projection ($`K=1`$, $`M=2`$) is shown on the left of Figure 2. In this case we have four minimal potential barriers which the system has to overcome. For the special case when all non-Gaussian sources have the same kurtosis $`\kappa _4^i=\kappa _4`$ ($`i=1,2,\mathrm{},M`$), the escape time, i.e. the mean first passage time for a single source to be learned, is given by,

$$T_{\text{escape}}^{\text{odd}}\frac{N^3}{2MK}\mathrm{exp}\left[\frac{3\varphi ^2(\mu )\nu }{8|\kappa _4\varphi ^{\prime \prime \prime }(\mu )|}\right],$$
(34)
where $`\sigma =\text{sign}(\varphi ^{\prime \prime \prime }(\mu )\kappa _4)`$ is a necessary condition for successful learning. Notice that the time-scale for escape diverges with the learning rate. On the left of Figure 3 we show numerical simulations for this situation. We have generated data of dimension $`N=50`$ with $`M`$ non-Gaussian sources each with a uniform distribution and we extract a single projection $`K=1`$. The figure shows the dependence of the escape time on the number of non-Gaussian sources, on a log-log plot. The solid line shows the inverse scaling with $`M`$ predicted by the theory and this is consistent with the simulation results. A source was considered learned when the associated overlap $`R_{ij}`$ was observed to be greater than a threshold magnitude (0.1 greater than the position of the potential barrier).

If the sources have different kurtosis then the algorithm will be most likely to learn the source with highest kurtosis first since this will correspond to the escape point with the lowest potential barrier. The mean first passage time will be dominated by the contribution from this barrier for large learning rates and the escape time will not depend on the number of sources.

### 5.3 Even Non-linearity

If the non-Gaussian signals are asymmetrical, then we can use an even non-linearity, for example $`\varphi (y)=y^2`$. In this case the appropriate scaling for the learning rate is $`\eta =\nu N^{3/2}`$. After expanding Equation (22) near $`𝒓=\mathrm{𝟎}`$ we find that the mean and covariance of the change in $`𝒓`$ at each iteration are given by (to leading order in $`N^1`$),

$`\text{E}[\mathrm{\Delta }r_{ij}]=(\frac{1}{2}\varphi ^2(\mu )\nu ^2r_{ij}+\frac{1}{2}\kappa _3^j\varphi ^{\prime \prime }(\mu )\sigma _{ii}\nu r_{ij}^2\frac{1}{2}\varphi (\mu )^2\nu ^2{\displaystyle \underset{li}{\overset{K}{}}}r_{lj})N^2,`$ (35)
$`\text{Cov}[\mathrm{\Delta }r_{ij},\mathrm{\Delta }r_{kl}]=\varphi ^2(\mu )\nu ^2\delta _{ik}\delta _{jl}N^2,`$ (36)
where $`\kappa _3^j`$ is the third cumulant of the $`j`$-th source distribution (third central moment), which measures skewness, and brackets denote averages over Gaussian variables $`𝝁𝒩(\mathrm{𝟎},𝑰)`$. Again the system can be described by a Fokker-Planck equation for large $`N`$ but now with shorter characteristic time-scale $`(\delta t)^1=N^2`$. The system is locally equivalent to a diffusion process in the cubic potential

$`U(𝒓)={\displaystyle \underset{i=1}{\overset{K}{}}}{\displaystyle \underset{j=1}{\overset{M}{}}}(\frac{1}{4}\varphi ^2(\mu )\nu ^2r_{ij}^2\frac{1}{6}\kappa _3^j\varphi ^{\prime \prime }(\mu )\sigma _{ii}\nu r_{ij}^3+\frac{1}{2}\varphi (\mu )^2\nu ^2{\displaystyle \underset{li}{\overset{K}{}}}r_{lj}r_{ij}),`$ (37)
with a diagonal diffusion matrix of magnitude $`D=\varphi ^2(\mu )\nu ^2`$ as before. In this case the sum in the potential contains “cross-terms” which depend on more than one element in $`𝒓`$. The dynamics is therefore not equivalent to the one-dimensional case and features of the potential will depend on the particular value of $`K`$ and $`M`$ considered, making analysis less straightforward than for the odd non-linearity.

The shape of the potential for an example with two sources of equal skewness and one projection ($`K=1`$, $`M=2`$) is shown on the right of Figure 2. In this case we have a ledge in the potential with two points of minimum height $`\mathrm{\Delta }U`$ (escape points). In the general case we find that the effective size of the minimal barriers is given by,

$$\frac{\mathrm{\Delta }U}{D}=\frac{f(K)\nu ^2}{\varphi ^2(\mu )(\kappa _3^i)^2},$$
(38)
where the function $`f(K)`$ has a complex form which depends on the choice of non-linear function $`\varphi (\mu )`$. For example, the shape of this function for $`\varphi (\mu )=\mu ^2`$ is shown on the left in Figure 4. We see that the size of potential barrier decreases with increasing number of projections $`K`$ and this appears to be a general feature of the function. This suggests that parallel algorithms for extracting asymmetrical signals may prove more efficient than deflationary ones which separate one signal at a time.

For the case of a single projection ($`K=1`$) the potential does decompose into a sum of independent terms and each component of $`𝒓`$ evolves independently. For the case of sources with equal skewness $`\kappa _3^i=\kappa _3`$ ($`i=1,2,\mathrm{},M`$), the mean first passage time, i.e. the time until one of the source signals is learned, is then given by,

$$T_{\text{escape}}^{\text{even}}\frac{N^2}{M}\mathrm{exp}\left[\frac{1}{12}\left(\frac{\varphi ^2(\mu )\nu }{\kappa _3\varphi ^{\prime \prime }(\mu )}\right)^2\right].$$
(39)
Numerical results from simulations of this scenario with $`N=50`$ are shown on the right of Figures 3 and 4. The asymmetrical sources used in the simulations are binary and each have the same skewness with $`\kappa _3=1.5`$. In Figure 3 we show the escape time as a function of the number of sources $`M`$ on a log-log plot. The solid line shows the inverse scaling predicted by the theory and is consistent with the experimental results. On the right of Figure 4 we show how the escape time (on a log scale) depends on the learning rate parameter $`\nu `$. The slope of the solid line is the theoretical prediction from Equation (39) and is consistent with the simulation results.

## 6 Transient Dynamics

Consider the more general situation when $`T<M`$ non-Gaussian sources have already been learned by the system. The corresponding fixed points of Equation (22) are $`R_{ij}^{}=\delta _{ij}\text{I}\left[iT\right]`$, where $`\text{I}\left[iT\right]`$ is defined by (24). We have already considered the special case when $`T=0`$ which is appropriate close to the initial conditions when no sources have yet been learned. A typical learning dynamics proceeds by passing through these states one by one until all the sources are learned. We show such a dynamical trajectory on the left of Figure 5. The indices have been labeled retrospectively so that the sources are learned in order although this labeling is clearly arbitrary and only chosen for notational convenience. It appears that the typical time to learn each source increases as more sources are learned, a phenomenon which is explained below.

We introduce new $`O(1)`$ scaled variables,

$$\nu =\eta N^d,𝒗=(𝑹𝑹^{})\sqrt{N},$$
(40)
where $`N`$ is the input dimension and $`d`$ is the scaling order for the learning rate. For the case of an even non-linearity we set $`d=\frac{3}{2}`$ while for the case of an odd non-linearity we choose $`d=2`$. In our new variables the fixed point is $`𝒗=\mathrm{𝟎}`$. We can compute the mean and covariance of these variables as we did for the $`𝒓`$ variables close to the initial conditions in the previous Section. In the present case it is convenient to consider four categories of variables separately.

1. $`iT,jT`$.

$`\text{E}[\mathrm{\Delta }v_{ij}]`$ $`=`$ $`\left[\left(\xi _i+\frac{1}{2}\xi _j\right)v_{ij}\frac{1}{2}\xi _iv_{ji}\right]\nu N^{1p}+O(N^p),`$
$`\text{Cov}[\mathrm{\Delta }v_{ij},\mathrm{\Delta }v_{kl}]`$ $`=`$ $`O(N^p).`$ (41)
2. $`i>T,jT`$

$$\text{E}[\mathrm{\Delta }v_{ij}]=\frac{1}{2}\xi _jv_{ij}\nu N^{1p}+O(N^p),\text{Cov}[\mathrm{\Delta }v_{ij},\mathrm{\Delta }v_{kl}]=O(N^p).$$
(42)
3. $`iT,j>T`$

$$\text{E}[\mathrm{\Delta }v_{ij}]=\xi _iv_{ij}\nu N^{1p}+O(N^p),\text{Cov}[\mathrm{\Delta }v_{ij},\mathrm{\Delta }v_{kl}]=O(N^p).$$
(43)
4. $`i>T,j>T`$

$$\text{E}[\mathrm{\Delta }v_{ij}]=\frac{U(𝒗)}{v_{ij}}N^p+O(N^{p1}),$$
(44a)
$$\text{Cov}[\mathrm{\Delta }v_{ij},\mathrm{\Delta }v_{kl}]=D\delta _{ik}\delta _{jl}N^p+O(N^{p1}).$$
(44b)

In these equations, $`p`$ is the scaling order of our system with $`p=2`$ for the even non-linearity and $`p=3`$ for the odd non-linearity. We define,

$$\xi _i=\sigma _{ii}\left(\text{E}_{s_i}[s_i\varphi (s_i)\varphi ^{}(s_i)]\right).$$

In Equation (44b) the exact expressions for the potential $`U(𝒗)`$ and for the diffusion coefficient $`D`$ have the same form as those which were found in Section 5 given by Equations (30), (31), (35) and (36).

In the first three groups of variable we observe that the fluctuations are of negligible order, so that the dynamics can be described by linear differential equations in the large $`N`$ limit with a relatively fast time-scale of $`\delta t^1=N^{p1}`$. Equations (41), (42) and (43) therefore converge exponentially to the fixed point as long as the condition in Equation (14) is met. However, the variables in the fourth group display a similar diffusive dynamics to that considered in the previous Section. The dynamics for these variables is completely equivalent to the $`𝒓`$-dynamics close to the initial conditions. Therefore we observe the same behaviour in these variables, with localisation at the fixed point until one component escapes over the potential barrier resulting in another source being learned. Once all the sources are learned we effectively have $`T=M`$ and only the first three groups of variables above remain. In this case we can increase the learning rate to an $`O(N^1)`$ quantity in principle without the stochastic effects dominating, but then the learning rate should be annealed as described by Rattray (2002) in order to converge asymptotically to the optimal solution.

The picture on the left in Figure 5 is a good illustration of the transient dynamics. We show numerical simulations for the typical dynamics of a Hebbian ICA algorithm extracting four ($`K=M=4`$) uniformly distributed non-Gaussian sources from an $`N=50`$ dimensional data set. On the right we show a log-log plot of the time $`\mathrm{}t_i^{\text{escape}}`$ required to learn the next source signal in the case when $`i`$ non-Gaussian sources have already been learned by the system and we see that it is consistent with the expected trend shown by the solid line. The total learning time for extracting all the non-Gaussian sources in this case will be

$$T_{\text{escape}}^{\text{ total}}=\underset{i=0}{\overset{M1}{}}\mathrm{}t_i^{\text{escape}}\mathrm{exp}\left[\frac{\mathrm{\Delta }U}{D}\right]\underset{i=0}{\overset{M1}{}}\frac{1}{2(Mi)^2},$$
(45)
where $`\mathrm{\Delta }U/D`$ is given by Equation (33).

## 7 Natural Gradient ICA

Natural gradient algorithms have been developed for ICA which use the structure of the parameter space to define a Riemannian gradient descent direction (Amari et al., 1996). Along with closely related relative gradient algorithms (Cardoso and Laheld, 1996) these methods provide some advantages over standard gradient descent methods, such as greater simplicity, robustness and asymptotic efficiency (Amari, 1998). However, these algorithms have mainly been defined for the special case where the mixing matrix is square and invertible.

The algorithm we use here searches the space of tall thin orthogonal matrices. This allows it to extract a relatively small number of independent components from a high-dimensional data set possibly containing Gaussian components. Standard natural gradient ICA algorithms are not appropriate in this case and we therefore need a different approach. One possibility would be to use a parameterisation of the set of orthogonal matrices. This approach is considered by Moon and Gunther (2002) who provide an interesting reinterpretation of natural gradient as a pullback. This allows them to define natural gradient algorithms for various structured matrices. Although they restrict their attention to square, invertible matrices their ideas could be extended to tall thin matrices. However, the available parameterisations appear to be quite complex in this case and computing the gradient even more so.

The approach of Moon and Gunther (2002) is to use a set of coordinates which are intrinsic to the manifold. An alternative approach is to use the original variables subject to constraints, i.e. work in the space of tall thin matrices $`𝑾`$ but impose the orthogonality constraint,

$$𝑾^\mathrm{T}𝑾=𝑰.$$
(46)
This is the approach taken by Edelman et al. (1999) and it leads to a much more straightforward gradient definition for ICA which is described by Amari (1999). The constraint surface is known as a Stiefel manifold and for a function $`F(𝑾)`$ defined on the Stiefel manifold, the “natural” gradient of $`F`$ at the point $`𝑾`$ of the manifold is defined by

$$\stackrel{~}{\mathbf{}}_WF=\mathbf{}_WF𝑾\mathbf{}_WF^\mathrm{T}𝑾,$$
(47)
where the standard gradient $`\mathbf{}_WF`$ is the $`K`$-by-$`M`$ matrix of partial derivatives of $`F`$ with respect to the elements of $`𝑾`$. The loss function used in Hebbian ICA is some non-quadratic function of the projections $`F(𝒚)`$ and the standard gradient of this function is given by

$$\mathbf{}_WF=𝒙\varphi (𝒚)^\mathrm{T}𝝈.$$
(48)
Then, according to Equation (47), the natural gradient of this function on the Stiefel manifold will be (Amari, 1999),

$$\stackrel{~}{\mathbf{}}_WF=𝒙\varphi (𝒚)^\mathrm{T}𝝈𝑾𝝈\varphi (𝒚)𝒚^\mathrm{T}.$$
(49)
A disadvantage of using non-intrinsic variables is that the algorithm is not guaranteed to stay on the manifold. This is especially problematic for stochastic gradient algorithms which only approximately follow the gradient direction. We therefore add the same orthogonalisation term used in the standard Hebbian algorithm. The natural gradient algorithm then has the following form,

$$𝚫𝑾=\eta [𝒙\mathit{\varphi }(𝒚)^\mathrm{T}𝝈𝑾𝝈\varphi (𝒚)𝒚^\mathrm{T}]+\alpha 𝑾(𝑰𝑾^\mathrm{T}𝑾).$$
(50)
The update increment for the overlap matrices $`𝑹𝑾^\mathrm{T}𝑨_s`$ and $`𝑸𝑾^\mathrm{T}𝑾`$ at every learning iteration is found to be,

$`\mathrm{\Delta }𝑹`$ $`=`$ $`\eta \left(𝝈\mathit{\varphi }(𝒚)𝒔^\mathrm{T}𝒚\varphi (𝒚)^\mathrm{T}𝝈𝑹\right)+\alpha (𝑰𝑸)𝑹,`$ (51)
$`\mathrm{\Delta }𝑸`$ $`=`$ $`\eta 𝝈(𝑰𝑸+\alpha (𝑰𝑸)\alpha (𝑰𝑸)𝑸)\mathit{\varphi }(𝒚)𝒚^\mathrm{T}+\alpha ^2(𝑰𝑸)^2𝑸`$ (52)
$`+`$ $`\eta 𝝈𝒚\mathit{\varphi }(𝒚)^\mathrm{T}(𝑰𝑸+\alpha (𝑰𝑸)\alpha 𝑸(𝑰𝑸))+2\alpha (𝑰𝑸)𝑸`$
$`+`$ $`\eta ^2\mathit{\varphi }(𝒚)𝒙^\mathrm{T}𝒙\mathit{\varphi }(𝒚)^\mathrm{T}.`$
After adiabatic elimination of the $`𝑸`$ variables by a similar procedure as we carried out for the case of Hebbian ICA (see Section 3) we have the following dynamical equations for the overlaps,

$`\mathrm{\Delta }𝑹=\eta 𝝈\left(\mathit{\varphi }(𝒚)𝒔^\mathrm{T}𝒚\mathit{\varphi }(𝒚)^\mathrm{T}𝑹\right)\frac{1}{2}\eta ^2N\mathit{\varphi }(𝒚)\mathit{\varphi }(𝒚)^\mathrm{T}𝑹.`$ (53)
The typical learning dynamics of this natural gradient version of Hebbian algorithm is shown on the left of Figure 6, where we used the odd non-linearity $`\varphi (y)=y^3`$ to extract four symmetrical sources with kurtosis $`\kappa _4=4/9`$ from 50-dimensional data ($`K=M=4`$, $`N=50`$).

Following the procedure outlined in Section 6 we can expand near the general fixed points with $`TM`$ sources learned. Using the same variables we expand around $`𝒗=\mathrm{𝟎}`$ and find the following results

1. $`iT,jT`$.

$$\text{E}[\mathrm{\Delta }v_{ij}]=\left[\left(\xi _i+\xi _j\right)\right]v_{ij}\nu N^{1p}+O(N^p),\text{Cov}[\mathrm{\Delta }v_{ij},\mathrm{\Delta }v_{kl}]=O(N^p).$$
(54)
2. $`i>T,jT`$

$$\text{E}[\mathrm{\Delta }v_{ij}]=\xi _jv_{ij}\nu N^{1p}+O(N^p),\text{Cov}[\mathrm{\Delta }v_{ij},\mathrm{\Delta }v_{kl}]=O(N^p).$$
(55)
3. $`iT,j>T`$

$$\text{E}[\mathrm{\Delta }v_{ij}]=\xi _iv_{ij}\nu N^{1p}+O(N^p),\text{Cov}[\mathrm{\Delta }v_{ij},\mathrm{\Delta }v_{kl}]=O(N^p).$$
(56)
4. $`i>T,j>T`$

$$\text{E}[\mathrm{\Delta }v_{ij}]=\frac{U(𝒗)}{v_{ij}}N^p+O(N^{p1}),$$
(57a)
$$\text{Cov}[\mathrm{\Delta }v_{ij},\mathrm{\Delta }v_{kl}]=D\delta _{ik}\delta _{jl}N^p+O(N^{p1}).$$
(57b)

As before, $`p=2`$ for the even non-linearity and $`p=3`$ for the odd non-linearity and the potential in the last case is the same as for standard Hebbian ICA. We see that the equations for the first set of variables ($`iT,jT`$) in Equation (54) are simplified in comparison to Equation (41) and no longer contain cross-terms. This means, for example, that the algorithm will enjoy uniform asymptotic convergence if all sources have identical statistics and $`M=K`$. Generally speaking the eigenvalues determining the convergence of the variables in the first three groups have lower variance and the asymptotic convergence of the natural gradient algorithm will be faster than that of the Hebbian algorithm. However, Equations (57b) and (44b) are identical and therefore the transient dynamics of the algorithms will be very similar. These are the variables which provide the rate limiting factor and learning time-scale during the transient.

The plot on the right of Figure 6 shows the escape time from each of the transient fixed points encountered during the dynamics. These simulation results confirm that the transient dynamics is very similar to the standard Hebbian algorithm results (see Figure 5) as predicted by our theory.

## 8 Conclusion

The dynamics of on-line ICA learning has been studied in the limit of large data dimension. We have analysed a Hebbian learning algorithm which is appropriate for extracting a prescribed number of components from high dimensional data possibly containing Gaussian components. We also studied a natural gradient variant of the algorithm which uses the gradient defined on the Stiefel manifold of orthogonal matrices.

We find that the learning time-scale of both algorithms is mainly determined by the transient dynamics. Learning takes place by a sequence of symmetry breaking steps in which a new source is learned and these steps can be described as a diffusion and escape process. The learning time-scale is found to be longer than expected from the analysis of related algorithms such as on-line back-propagation and Sanger’s PCA algorithm (e.g. Saad and Solla, 1995, Biehl and Schlösser, 1998). To learn each symmetric source typically requires of the order of $`N^3`$ learning iterations while to learn an asymmetrical source using an even non-linearity typically requires of the order of $`N^2`$ learning iterations. Both algorithms exhibit equivalent transient dynamics and we only find an advantage in using the natural gradient variant asymptotically.

### Acknowledgments.

This work was supported by an EPSRC award (ref. GR/M48123).


text:
# Granular packings with moving side walls

## I Introduction

There has been a recent resurgence in interest in the formation and structure of granular packings in the physics community Jaeger et al. (1996); Duran (2000). One particular facet of granular packings to receive attention recently is their stress profiles, both because of the industrial applications of silos Nedderman (1992), and because of new experimental techniques to measure effective stress in packings Mueth et al. (1998); Vanel and Clément (1999); Vanel et al. (2000). The first theoretical attempt to understand stress in a silo geometry dates back over 100 years to Janssen Janssen (1895), who obtained a one-parameter form for the vertical stress in a silo. Several assumptions were made to arrive at this result. One was to treat the granular material as a continuous medium where a fraction $`\kappa `$ of vertical stress is converted to horizontal stress. Another assumption was that the forces of friction between particles and walls are at the Coulomb failure criterion: $`F_t=\mu _wF_n`$, where $`F_t`$ is the magnitude of the tangential friction force, $`F_n`$ is the normal force at the wall, and $`\mu _w`$ is the coefficient of friction for particle-wall contacts. This assumption is also known as incipient failure. For a cylindrical container of radius $`R`$ with static wall friction $`\mu _w`$ and granular pack of total height $`z_0`$, the Janssen analysis predicts the vertical stress $`\sigma _{zz}(z)`$ at a height $`z`$ is

$$\sigma _{zz}(z)=\rho gl\left[1\mathrm{exp}\left(\frac{z_0z}{l}\right)\right]$$
(1)
where the decay length $`l=\frac{R}{2\kappa \mu _w}`$. $`\kappa `$ represents the fraction of the weight carried by the side walls, $`\rho `$ is the volumetric density, $`g`$ is gravity, and $`z_0`$ is height of the top of the packing.

Numerous experiments have been carried out to verify this theory, but precise experiments are difficult. Recently, extremely well-controlled experiments on granular packings have been done Vanel and Clément (1999); Vanel et al. (2000) that show a deviation from the ideal Janssen form, which we shall hereafter refer to as the Vanel-Clément form. This phenomenological form with one additional free parameter includes a hydrostatic-like region of linear dependence of pressure with depth (defined by the length $`a`$), followed by a region that conforms to the Janssen theory:

$`z_0z<a`$ $`:`$ $`\sigma _{zz}(z)=\rho g(z_0z)`$
$`z_0z>a`$ $`:`$ $`\sigma _{zz}(z)=\rho g\left(a+l\left[1\mathrm{exp}\left({\displaystyle \frac{z_0za}{l}}\right)\right]\right)`$ (2)
This form was also found in extensive molecular dynamics simulations of granular packings in both two and three dimensions Landry et al. (2003a, b). These packings were created both through pouring and sedimentation and then allowed to settle under the influence of gravity.

Many questions about stress in granular packings still remain unanswered. Even after a packing has been formed, it may be perturbed in many ways that radically change its stress profile and physical structure. Many studies have focused on tapping as a means to compress the packing and its logarithmic response time Knight et al. (2003); Nowak et al. (1998); Josserand et al. (2000); Philippe and Bideau (2002). Another method to perturb a packing is to move the side walls Bertho et al. (2002); Ovarlez et al. (2003); Ovarlez and Clément (2003). The effect of this movement is not well understood and is the focus of this study.

Recently, experiments have been conducted on granular packings in cylindrical containers with movable side walls Bertho et al. (2002); Ovarlez et al. (2003); Ovarlez and Clément (2003). The experiments make use of a movable cylinder enclosing a granular packing supported by an independent base. These experiments find over a wide range of wall velocities good agreement with the Janssen form for the vertical stress after the system has relaxed following cessation of wall movement. This is in contrast to the earlier experiments on packings with fixed side walls Vanel et al. (2000); Vanel and Clément (1999). Here we present large-scale three-dimensional ($`3D`$) discrete element, molecular dynamics simulations of granular packings in cylindrical containers (silos) with movable side walls. Our aim is to understand how the motion of the side wall modifies the stress profiles in granular packings and to compare our results with the recent experimental findings Bertho et al. (2002); Ovarlez et al. (2003); Ovarlez and Clément (2003). We analyze how wall movement and its cessation affect the stress profile of the packings. We also investigate in depth the effects of wall movement on the internal structure and particle positions of these packings, which cannot be easily measured experimentally. The behavior of the system under wall movement is similar across a wide range of wall velocities. Finally we show that wall movement in the direction opposite to gravity drives tangential forces to the Coulomb criterion everywhere, leading directly to the Janssen form for the stress profiles.

The simulation technique and model are presented in Sec. II. The stress profiles and their features are discussed in Sec. III. In Sec. IV we examine the particle motion during and after wall movement, while Sec. V we discuss the force distribution of the resultant packings. A brief summary and conclusions are presented in Sec. VI.

## II Simulation Method

We present discrete element, molecular dynamics (MD) simulations in $`3D`$ of model systems of $`N=50000`$ mono-dispersed spheres of fixed mass $`m`$ and diameter $`d`$. The system is constrained by a cylinder of radius $`R=10d`$, centered on $`x=y=0`$, with its axis along the vertical $`z`$ direction. The cylinder is bounded below by a layer of randomly-arranged immobilized particles approximately $`2d`$ high to provide a rough base. This work builds on previous MD simulations of packings in silos, where further details of the model can be found Landry et al. (2003a).

The spheres interact only on contact through a spring-dashpot interaction in the normal and tangential directions to their lines of centers. Contacting spheres $`i`$ and $`j`$ positioned at $`𝐫_i`$ and $`𝐫_j`$ experience a relative normal compression $`\delta =|𝐫_{ij}d|`$, where $`𝐫_{ij}=𝐫_i𝐫_j`$, which results in a force

$$𝐅_{ij}=𝐅_n+𝐅_t.$$
(3)
The normal and tangential contact forces are given by

$$𝐅_n=f(\delta /d)(k_n\delta 𝐧_{ij}\frac{m}{2}\gamma _n𝐯_n)$$
(4)
$$𝐅_t=f(\delta /d)(k_t𝚫𝐬_t\frac{m}{2}\gamma _t𝐯_t)$$
(5)
where $`𝐧_{ij}=𝐫_{ij}/r_{ij}`$, with $`r_{ij}=|𝐫_{ij}|`$. $`𝐯_n`$ and $`𝐯_t`$ are the normal and tangential components of the relative surface velocity, and $`k_{n,t}`$ and $`\gamma _{n,t}`$ are elastic and viscoelastic constants, respectively. $`f(x)=1`$ for Hookean (linear) contacts while for Hertzian contacts $`f(x)=\sqrt{x}`$. $`𝚫𝐬_t`$ is the elastic tangential displacement between spheres, obtained by integrating tangential relative velocities during elastic deformation for the lifetime of the contact. The magnitude of $`𝚫𝐬_t`$ is truncated as necessary to satisfy a local Coulomb yield criterion $`F_t\mu F_n`$, where $`F_t|𝐅_t|`$ and $`F_n|𝐅_n|`$ and $`\mu `$ is the particle-particle friction coefficient. Frictionless spheres correspond to $`\mu =0`$. Particle-wall interactions are treated the same, but the particle-wall friction coefficient $`\mu _w`$ is set independently. The side wall of the container is smooth, and thus the particle-wall normal force is always perpendicular to the $`xy`$ plane. A more detailed description of the model is available elsewhere Silbert et al. (2001); Landry et al. (2003a).

These simulations are run with a fixed set of parameters: $`k_n=2\times 10^5mg/d`$, $`k_t=\frac{2}{7}k_n`$, and $`\gamma _n=50\sqrt{g/d}`$. For Hookean springs we set $`\gamma _t=0`$. In these simulations, it takes far longer to drain the energy out of granular packs using the Hertzian force law, since the coefficient of restitution $`ϵ`$ is velocity-dependent Schäfer et al. (1996) and goes to zero as the velocity goes to zero. We thus use Hookean contacts Her , which for the above parameters give $`ϵ=0.88`$. The convenient time unit is $`\tau =\sqrt{d/g}`$, the time it takes a particle to fall its radius from rest under gravity. For this set of parameters, the timestep $`\delta t=10^4\tau `$. The particle-particle friction and particle-wall friction are the same: $`\mu =\mu _w=0.5`$. The simulations are run using a parallel distributed memory code on 20 DEC Alpha processors. One million timesteps takes approximately 7 hours. Our longest simulation for slow wall velocities, $`v=5\times 10^4d/\tau `$, ran for $`t=1.14\times 10^4\tau `$, which corresponds to approximately 800 hours.

The simulations all begin with the same initial packing to minimize sample to sample fluctuations. This packing was generated by pouring particles into a container using method P1 as described in ref. Landry et al. (2003a). The beginning packing is quiescent. Over the course of a simulation, the cylindrical side wall of the packing moves for a time $`t_s`$, which usually is $`10^3\tau `$, or over a fixed distance $`\mathrm{\Delta }z`$. After this period, the walls cease to move and the packing settles. We consider a packing quiescent when the kinetic energy per particle $`E_k10^8mgd`$. The time scale for this relaxation is very short, usually less than $`10\tau `$. The cylindrical side wall is moved either up ($`+z`$) or down ($`z`$) with a constant velocity $`v_s`$ varying from $`10^1`$ to $`10^5d/\tau `$. As in the experiments Bertho et al. (2002); Ovarlez et al. (2003); Ovarlez and Clément (2003), only the side wall moves — the rough base is immobile throughout the course of the simulation.

The wall velocities used here for upward velocities are in the range used in the experiments by Bertho et al. Bertho et al. (2002). In that study, glass beads with $`d=2\mathrm{mm}`$ were moved with a velocity $`v_s`$ ranging from $`2\times 10^2`$ to $`35\mathrm{mm}/\mathrm{s}`$ or $`1.4\times 10^4`$ to $`0.25d/\tau `$ over distances up to $`\mathrm{\Delta }z=70\mathrm{mm}=35d`$. Ovarlez et al. Ovarlez et al. (2003) use glass beads with $`d=1.5`$ mm and a fixed $`v_s=1.5\times 10^3`$ mm/s or $`1.2\times 10^5d/\tau `$ over a distance $`\mathrm{\Delta }z1.5\times 10^2\mathrm{mm}=10^2d`$. This is a very low velocity over a very short distance. In an earlier study Vanel et al. (2000); Vanel and Clément (1999), Vanel et al. used a higher velocity, in conjunction with tapping, of $`v_s=2\times 10^2\mathrm{mm}/\mathrm{s}`$ for particles of the same diameter. Another distinction between the experiments by Ovarlez et al. Ovarlez et al. (2003) and Vanel et al. Vanel et al. (2000); Vanel and Clément (1999) was the time at which measurement occurs. In the Ovarlez et al. experiments, measurement occurs right at the end of wall movement, while in the Vanel et al. experiments, the packing is allowed to settle before measurements are made and tapping is sometimes also applied. Ovarlez and Clément also studied Ovarlez and Clément (2003) moving the wall downwards with a $`v_s`$ ranging from $`5\mathrm{nm}/\mathrm{s}`$ to $`100\mu \mathrm{m}/\mathrm{s}`$, which corresponds to wall velocities of $`4\times 10^8d/\tau `$ to $`8\times 10^3d/\tau `$. This range of velocities overlaps with our velocity range for downward wall movement, but also extends to much slower velocities.

Packings are examined both during wall movement and after cessation of wall movement and settling. Figure 1 shows the structure of a packing for a smaller system of $`20000`$ particles (used for illustration) with wall velocity $`v_s=10^1d/\tau `$ in the upward ($`+z`$) direction. This is a very high velocity, so during wall movement, there is significant particle rearrangement, and a number of particles originally in contact with the wall move upward. Particles between $`z=5`$ and $`z=15`$ have been colored light gray (green online) to provide a visual picture of particle movement over time. The height of the pile changes significantly during wall movement (4.7 %), but does not change after the wall movement has ceased and the packing reaches equilibrium. For much slower wall velocities applied for the same time, very little particle movement is observed, and the height of the packing does not change.

## III Stress Profiles

In this section we discuss the effect of wall movement on stress profiles of packings. In general, wall movement increases the tangential force between particles and the wall. The wall moves in relation to the particles, so the integrated displacement $`\mathrm{\Delta }s`$ increases as the wall moves. Since the tangential force is proportional to $`\mathrm{\Delta }s`$, $`F_t`$ increases over time. Wall movement thus drives $`F_t`$ towards the Coulomb limit, $`F_t=\mu _wF_n`$, for particles in contact with the wall. We know from previous work Landry et al. (2003a) that the prime factor determining the form of the stress profile in these packings is the tangential force between particle and wall. If $`F_t`$ is close to the Coulomb criterion, then we expect the stress to follow the Janssen form throughout the pile. In addition, the more particle-wall interactions that are at the Coulomb criterion, the stronger this effect and the lower the value of the saturation stress in the packing, i.e. the value of the stress in the depths of the packing, where the stress becomes depth-independent. We believe this effect is the origin of the experimental measurements observing remarkable agreement with the Janssen theory after the side walls of a packing are moved Bertho et al. (2002); Ovarlez et al. (2003). Below, we investigate the specific effects of wall movement on the stress profile.

Figure 2 shows the change in the stress profile of a packing with wall velocity $`v_s=0.01d/\tau `$, a relatively large velocity. The stress profile with the largest saturation stress corresponds to $`t=0`$, before wall movement has begun. This stress profile obeys the Vanel-Clément form Vanel and Clément (1999); Landry et al. (2003a), with a linear stress profile at the top of the pile crossing over to a Janssen form in the depth of the pile. After moving the wall for $`t=100\tau `$, the saturation stress has decreased by more than a factor of three to $`\sigma _{zz}11mg/d^2`$ and the height of the pile has increased. As the wall movement continues, the saturation stress slowly increases. After a long time $`t700\tau `$, the system reaches steady state and the stress profile ceases to change markedly. Over the course of the wall movement, the height of the pile has increased by approximately 3%, and there has been considerable particle rearrangement. In addition, the linear stress disappears almost immediately after the initiation of wall movement. This change in the saturation stress is consistent with that observed in granular experiments performed by Bertho et al. Bertho et al. (2002), where the apparent mass (effectively the saturation stress) dropped quickly after the initiation of wall movement and then slowly increased with time until it reached a saturation value.

The behavior of the stress profile is similar for low wall velocities, although not identical. Unfortunately, computational limitations are much more severe for these packings than for high wall velocity packings, because the computational time required for the side wall to travel the same distance is so much larger. Figure 3a shows the change in the stress profile for a relatively slow wall velocity $`v_s=10^4d/\tau `$. Each profile is $`2\times 10^2\tau `$ after the previous one, starting at $`t=0`$. For this wall velocity, the stress drops quickly at the base of the pile, but initially remains unchanged at the top of the pile. As wall movement continues, the new reduced stress profile propagates up the pile. Eventually, after moving the walls for $`t=2.2\times 10^3\tau `$, which corresponds to a very small vertical distance $`\mathrm{\Delta }z=0.22d`$, the entire stress profile follows the Janssen form. At this point in the simulation, the height of the pile has not changed. This minimum saturation stress is the same as that observed in Figure 2, $`\sigma _{zz}11mg/d^2`$, and is stable for a significant period of subsequent wall movement. This suggests that this is the limit for low saturation stress in this packing, and its value is controlled purely by geometric factors.

We explored the full behavior of the stress profile during wall movement using a low wall velocity of $`v_s=5\times 10^4d/\tau `$ for $`1.14\times 10^4\tau `$, so that the side walls moved a total distance of $`\mathrm{\Delta }z=5.7d`$, which is comparable to the distance walls were moved for higher velocities. As shown in Figure 3b, prolonged wall movement at low velocities does eventually move the stress profile away from the minimum saturation stress. Prolonged wall movement also increases the height of the pile by 3.3 %, which is similar to the height change for packings with high wall velocities. As the wall moves upward, changes in the stress profile propagate upwards from the bottom of the pack to the top. This occurs not only at early times, when the minimum saturation stress is propagated up the pile, but also at later times, when enhanced stress is also propagated up the pile as further wall movement increases the eventual saturation stress. This drop to a minimum saturation stress followed by an increase to a larger, stable saturation stress is the same for all observed wall velocities. The drop to a minimum saturation stress occurs at wall movement distances $`\mathrm{\Delta }z0.3d`$ for all wall velocities measured. In addition, the increase in saturation stress from the minimum occurs when the wall has moved a distance $`\mathrm{\Delta }z1.6d`$. At this distance, particles near the wall have completely moved past nearby particles, which means that whatever contacts these particles had initially have been destroyed. The minimum saturation stress is thus the optimized stress network for a given initial condition, and wall movement of greater than $`d`$ destroys the initial stress network, and forms another, which is no longer optimal. This process is not steady, however. The stress changes relatively abruptly at certain times, then remains essentially unchanged for long periods (on order of $`3\times 10^3\tau `$ for $`v_s=5\times 10^4d/\tau `$), then abruptly changes again. One exception to this is that the height of the packing increases smoothly with time until it reaches a stable height at $`9\times 10^3\tau `$ of $`140.5d`$, after which it remains unchanged under further wall movement.

Figure 4 shows the final stress profiles for the same packings after moving the wall for $`t=10^3\tau `$ for high and low wall velocity and then allowing them to settle. The wall movement has very different effects on the final stress profile depending on its strength. As seen in Figure 4a, wall velocities of $`v_s10^3d/\tau `$ increase the height of the pile relative to the original packing, even after relaxation. The change in height is substantial and similar for many different wall velocities, 3.7%, meaning there has been a large change in the density of the packing. The final stress after relaxation is related to the magnitude of the velocity $`v_s`$, for $`v_s10^3d/\tau `$. Those packings with larger $`v_s`$ had larger saturation stress in the final packing. As we shall see in Section IV, the larger $`v_s`$, the more the particles rearrange. In addition, the larger $`v_s`$, the larger the particle rearrangement after cessation of wall movement is. This suggests that particle rearrangement strongly influences the final saturation stress. Those packings with large $`v_s`$ have enough particle rearrangement to completely disrupt the force network. The higher the velocity, the more contacts both at the wall and in the bulk are broken and the larger the eventual saturation stress. This behavior was also observed by Bertho et al. Bertho et al. (2002): after the wall stopped moving, the packing settled and the apparent mass increased. In addition, they observed the same trend after relaxation. The larger the $`v_s`$, the larger the final saturation stress.

By contrast, small wall velocities ($`v_s<10^3d/\tau `$) do not change the height of the pile over the same time period, largely because the wall does not travel far enough to substantially disrupt the force network. Large-scale rearrangements do not occur for this time duration. In this case, the particles against the wall are fully mobilized and the resultant saturation stress is small, because most of the pressure is supported by the walls. When the packing is allowed to relax, there are no large particle rearrangements, because the wall movement is not large enough to move particle positions significantly, and the wall movement has not put very much energy into the packing. The packing stays in the minimum saturation stress configuration. In addition, if the wall does not move for long enough to force the entire packing into the minimum saturation stress configuration, the intermediate stress configuration is stable under relaxation, as shown in Fig. 4b for the packing with wall velocity $`v_s=10^4d/\tau `$.

If the walls are moved long enough at a low velocity to approach the stable final saturation stress (higher than the minimum saturation stress) the height of the pile does increase, as shown in Figure 3b for $`v_s=5\times 10^4d/\tau `$. After relaxation, however, the stress remains unchanged. The slow wall movement does not introduce enough energy into the packing for the stress to change during relaxation, even though for these long times, there is significant particle rearrangement before relaxation.

Figure 5 shows a comparison of the stress profile of two distinct packings. The first is the original quiescent packing. The second packing is this same packing after the wall has moved for $`t=10^3\tau `$ at $`v_s=0.1d/\tau `$ and then the packing has settled. The nature of the stress profile has been radically changed. While the original packing fits very well the Vanel-Clément form Vanel and Clément (1999), the same packing after wall movement fits the one parameter Janssen form Janssen (1895) well with $`\kappa =0.44`$. As we shall see in Section V, the wall movement has forced the tangential force of the particles at the wall to the Coulomb criteria everywhere, eliminating the linear stress region. Another example is given for a slow velocity $`v_s=5\times 10^4d/\tau `$, where $`\kappa =0.63`$. In this case, however, the Janssen fit overshoots the actual stress near the top of the pack, the opposite of the case where it fails for the original stress region. For all velocities observed, the linear stress region is destroyed, but only the high velocity cases fit the Janssen stress form well.

We have compiled $`\kappa `$’s for the various wall movement rates and show them in Figure 6. These values can be compared with Figure 4 in Bertho et al. Bertho et al. (2002). In both cases, as the wall velocity increases, the resultant Janssen length $`l`$ increases. Since $`\kappa =R/2\mu _wl`$, $`\kappa `$ decreases with increasing wall velocity $`v_s`$. The $`\kappa `$’s observed are much lower than those obtained for the original packings with the Vanel-Clément form, which were slightly greater than 1 Landry et al. (2003a). In our case, the intermediate stresses between the minimal stress and final stress do not follow the Janssen form and $`\kappa `$ for those stress profiles has no meaning. Using the observed $`l`$ and $`R=7.5d`$ given in ref. Bertho et al. (2002), we find $`\kappa `$’s ranging from $`\kappa =0.349`$ for $`v_s=20mm/s=0.14d/\tau `$ to $`\kappa =0.405`$ for $`v_s=0.2mm/s=1.4\times 10^3d/\tau `$. These values are slightly lower but close to our observed values, shown in Fig. 6.

Downward motion of the wall is a very different phenomenon. In this case, the wall movement merely increases the stress at the base of the pile, without changing the stress profile elsewhere, as shown in Figure 7, which shows the stress profile for a particle packing with $`v_s=10^4d/\tau `$ applied for $`t_s=10^3\tau `$. This stress buildup increases with duration of wall movement and occurs for all the velocities probed, from $`v_s=10^1d/\tau `$ to $`10^5d/\tau `$. In addition, this stress buildup is robust and does not disappear when the packing settles after cessation of wall movement. The extra stress is locked in. There is a slight difference between high and low wall velocities. Low wall velocities do not change the height of the packing, while high wall velocities ($`v_s10^3d/\tau `$) increase the density at the bottom of the pile and lower the overall height of the packing by a slight amount.

## IV Particle Rearrangement

We have studied the motion of particles during and after wall movement to understand the effects of wall movement on the position and density of the packing. As the walls are moved upward with velocity $`v_s`$, they drag particles in contact with the walls upwards by means of the frictional force between them. As particles on the edge of the cylindrical container move upward, particles in the middle of the packing move downward and outward toward the walls to fill the voids created by the upward-moving particles.

We present two sets of data for both high and low velocities. The first data set is for a high wall velocity of $`v_s=10^2d/\tau `$, and is presented in two ways. The first is a histogram of the motion of the particles in $`z`$ from their initial starting position to the final position after wall movement for $`t=10^3\tau `$, $`\mathrm{\Delta }z=10d`$. The measured movement of the particles in the $`z`$-direction $`\delta z`$ is recorded and a histogram is generated. Only those particles with an initial position $`z10d`$ are included. In the second method we average $`\delta z`$ over $`z`$ and present an image map of the motion of these particles in $`z`$. Both the top ($`z>120d`$) and bottom ($`z<10d`$) regions of the initial packing are excluded to avoid edge effects. These data sets are presented for the packing after wall movement for $`t=10^3\tau `$ and then after a relaxation to equilibrium in Figure 8.

For high wall velocities, there is significant particle movement. In this case, particles travel upwards as much as $`7d`$. However, all of this movement is restricted to a small ring around the cylinder walls. In the center of the pack, particle movement is less than $`d`$. Particles slide past one another and the initial stress networks are quickly destroyed. New force networks are constantly created and destroyed over the course of the wall movement.

The initial drop in the saturation stress when the side wall has moved approximately $`0.3d`$ is not enough to move the particles past each other. This minimum saturation stress thus reflects a particle configuration where the initial contacts have been strengthened. Once particles at the wall move past their original contacts, the original stress network is completely disrupted, and later networks as they form and dissolve support less stress, causing the saturation stress to increase.

When the wall movement stops, the packing settles, and the new $`\delta z`$ is much smaller in magnitude than the $`\delta z`$ during wall movement. Also, the particles at the wall move much less than the particles in the center, which suggests that the particles at the wall are supported by friction at the walls.

Figure 9 shows $`\delta z`$ for the low wall velocity $`v_s=5\times 10^4d/\tau `$ applied for the much longer time $`t=1.14\times 10^4\tau `$ for a total wall movement of $`\mathrm{\Delta }z=5.7d`$. The change in particle heights for $`\delta z`$ is similar to that seen in high velocity runs. Particles along the walls move upward much more than those in the center. However, unlike in the high-velocity case, there is no second ring of particles one or two $`d`$ in from the wall that also shows a large $`\delta z`$ relative to the center. In this case, the particles not in contact with the wall show much less difference in $`\delta z`$ as a function of radius.

Relaxation is very interesting in this case. The difference in $`\delta z`$ between particles at the wall and particles in the bulk is much less than in the high velocity case. In addition, relaxation occurs in discrete jumps with many particles moving together. This is a general feature of $`\delta z`$ relaxation histograms for low wall velocity packings. In the low velocity case, energy is imparted into the system very slowly and is dissipated more quickly. There is thus much less energy available to rearrange the packing after the wall movement ceases, and the relaxation of the particles is much less. In addition, because individual particles have very little kinetic energy when the wall movement ceases, motion happens coherently as many particles rearrange at once.

## V Coulomb Criterion

We also examined the force distributions in the packings during and after wall movement. These provide further evidence that the final packing state after relaxation is very different from the initial state, and also demonstrate the different behavior of the two regimes after relaxation. A useful quantity is the ratio of normal forces $`F_n`$ to tangential forces $`F_t`$. A key assumption of the Janssen analysis is that the ratio of these forces $`\zeta =F_t/\mu _wF_n=1`$ for particle-wall forces in the packing. These distributions are shown in Figure 10 for several different packings. We consider three different distributions: forces in the bulk of the sample far from the wall, particle-wall forces in center of the packing, and particle-wall forces at the top of the packing. In this case, $`\zeta =F_t/\mu F_n`$ in the bulk of the sample, and $`\zeta =F_t/\mu _wF_n`$ for particle-wall forces.

The original packing has a strong peak near $`\zeta =0.2`$ in the bulk of the material, showing that the vast majority of forces are far from the Coulomb criterion. At the wall, the peak in the distribution is much closer to the Coulomb criterion at $`\zeta 0.9`$. At the top of the packing, in the hydrostatic region Landry et al. (2003a), the distribution of forces is far from the Coulomb criterion.

After wall movement, the distribution of forces is radically changed. For both the high velocity ($`v_s=10^2d/\tau `$) and low velocity ($`v_s=5\times 10^4d/\tau `$) cases, all three distributions are driven toward the Coulomb criterion, with the peak in the distributions at the Coulomb criterion. In both cases, the original distribution is not completely destroyed. For the bulk distribution, a subsidiary peak appears in the new bulk distribution where the original peak appeared. This change in the distribution of forces is the main cause for the more Janssen-like stress distributions observed.

We also present results for the $`v_s=10^2d/\tau `$ case after cessation of wall movement and relaxation in Fig. 10c. In this case, the main change in the distribution after relaxation is the increase in particles forces at very low $`\zeta `$. This occurs because some contacts between particles disappear during relaxation as particles move relative to their neighbors. When new contacts are made, these are far from the Coulomb criterion by definition, since particles are not, in general, rotating relative to each other during relaxation. By contrast, the low velocity case $`v_s=5\times 10^4d/\tau `$, Fig. 10 exhibits no change in the distribution of forces after cessation and relaxation. In this case, particles move together, and contacts are not destroyed or created, so the force distribution does not change. This lack of change in the force distribution may explain why low velocity packings do not exhibit perfect agreement with the Janssen form after relaxation.

## VI Conclusions

We have explored the effects of side wall movement on granular packings. Much of the resultant structure of the pack does not depend strongly on the magnitude of the wall movement, as long as the packing is moved for an equivalent distance. Small differences emerge in the behavior after cessation, where high wall velocity packings substantially rearrange and increase their saturation stress while low wall velocity packings remain essentially unchanged. The main effect of wall movement is to drive both the particle-wall and particle-particle contacts to the Coulomb criterion, so that the ratio of tangential forces to normal forces is maximized. This condition forces the packing in the high wall velocity case to obey the Janssen form, which takes the Coulomb criterion as one of its main assumptions. For low wall velocities, the final form is very different from the initial form, though it does not perfectly match the Janssen form.

This work was supported by the Division of Materials Science and Engineering, Basic Energy Sciences, Office of Science, U.S. Department of Energy. This collaboration was performed under the auspices of the DOE Center of Excellence for the Synthesis and Processing of Advanced Materials. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.