text:
# Role of surface disorder on the magnetic properties and hysteresis of nanoparticles

## 1 Introduction

The progressive reduction of particle sizes with technological application down to the nanometric range results in magnetic behaviours that are predominantly dominated by the particle surface spins. The magnetic behaviour of the particle surface differs from that corresponding to the core because of the lower coordination number and the existence of broken exchange bonds. Among other effects, this results in increased anisotropy compared to bulk values due to the surface anisotropy with uniaxial normal direction that comes from symmetry breaking at the boundaries of the particle . Moreover, depending on the character of the exchange interaction between spins in the particle and the geometry of their underlying lattice, surface spins may suffer random canting and form a spin-glass state as has been argued in studies of fine particles of different ferrimagnetic oxides . These systems display anomalous magnetic properties at low temperatures. In particular, they display high closure fields with high values of the differential susceptibility and shifts in the hysteresis loops after field cooling  whose origin is still a matter of controversy. While some of these features have been attributed the existence of dipolar interactions among the particles , there are experimental evidences that finite-size and surface effects are crucial in order to understand this phenomenology. In this study, we present the results of Monte Carlo (MC) simulations of the magnetic properties of individual nanoparticles which aim at clarifying what is the specific role played by increased radial anisotropy at the surface on the magnetization processes.

## 2 Model

The model considered is an extension to Heisenberg spins of our previous MC simulations for an Ising spin model of a $`\gamma `$-Fe2O3 particle . The particles considered are spheres of diameter $`D`$ (in units of the linear unit cell size $`a`$). The simulations have been run for ferrimagnetic particles with maghemite lattice structure and compared with those corresponding to ferromagnetic (FM) particles with the same underlying lattice. In order to characterize the specific contribution of the surface spins to the magnetic properties of the particle, we have distinguished two different regions in the particle: the surface formed by the outermost unit cells and an internal core of diameter $`D_{Core}=D2`$.

The magnetic ions are modeled by Heisenberg classical spins with the following interaction Hamiltonian:

$`H/k_B={\displaystyle \underset{i,j}{}}J_{ij}\stackrel{}{S}_i\stackrel{}{S}_j{\displaystyle \underset{i=1}{\overset{N}{}}}\stackrel{}{H}\stackrel{}{S}_i{\displaystyle \underset{i=1}{\overset{N}{}}}\left[k_C(S_i^z)^2+k_S(\stackrel{}{S}_i\widehat{n}_i)^2\right],`$ (1)
which is an extension of our previous study with Ising spins . Here, $`\stackrel{}{S}_i`$ are three dimensional classical vector spins placed on a discrete lattice and $`\stackrel{}{H}`$ is the magnetic field, that in the following will be given in temperature units as $`\stackrel{}{h}=\mu \stackrel{}{H}/k_B`$, with $`\mu `$ the moment of the magnetic ion. The $`J_{ij}`$’s are the nearest neigbour exchange constants, which might vary from site to site depending on the specific material considered, for the case of maghemite, which is a two sublattice (O,T) ferrimagnet, their values are $`J_{TT}=21`$ K, $`J_{OO}=8.6`$ K, $`J_{TO}=28.1`$ K . The last term in Eq. 1 accounts for the magnetic anisotropy energy, which has been taken as uniaxial along the z-axis for the core spins and along the radial direction for the spins at the surface. $`k_C`$ and $`k_S`$ are the core and surface anisotropy constants respectively in temperature units. Although the values of bulk anisotropy constants $`K_C`$ can be obtained indirectly by magnetic measurements, the surface contribution $`K_S`$ to the anisotropy energy is more difficult to evaluate. However, for maghemite $`K_C4.7\times 10^4`$ erg/cm3  and $`K_S`$ has been estimated as $`K_S0.06`$ erg/cm2 from Mössbauer experiments . Therefore, in the simulations for the maghemite particle, we will vary the core anisotropy values in the range $`k_C=0.011`$ K and those of the surface anisotropy in the range $`k_S=1100`$ K.

The simulations have been performed using the standard Metropolis algorithm  for continuous spins. As for the details, let us only mention that, in this case, special care has to be taken in the update algorithm with respect to the case of Ising spins. In order to sample efficiently the configuration space both in the low and high anisotropy limits, we have chosen a sequential combination of three kinds of trial angles in a similar way as Hinzke and Nowak . The update algorithm consists of a series of three uniform trial steps at random, one small trial step inside a cone around the initial spin direction and a reflection of the spin direction at the x-y plane to guarantee the correct performance of the simulation at high values of the anisotropy.

## 3 Results

In order to gain insight on the equilibrium properties of particles with surface anisotropy , we have performed a study of the equilibrium configurations in zero magnetic field at low temperatures. For this purpose, the system is initially prepared in a high temperature state with spins randomly oriented and then it is subsequently cooled in constant temperature steps $`\delta T=1`$ down to zero temperature. At each temperature we measure the thermal averages of the quantities of interest such as the total energy, sublattice magnetizations ($`M_O`$, $`M_T`$), surface $`M_S`$ and core $`M_C`$ contributions to the total magnetization $`M_{Total}`$, and also $`M_n=_{i=1}^N\left|\stackrel{}{S}_i\widehat{n}_i\right|`$, the sum of projections of the spins onto the local anisotropy axis.

For a better understanding of the role of surface anisotropy on the zero temperature configurations, we will first show the results of simulations for a spherical particle with $`D=3`$ and $`k_C=1`$, and the same structure than maghemite but with equal ferromagnetic (FM) interactions ($`J_{ij}=1`$) among the spins instead of the real values. In Fig. 1 we show the surface and core contributions to thermal dependence of the magnetization during the above-mentioned simulated annealing procedure for a range of values of $`k_S`$ spanning two orders of magnitude. The increase of $`M_z`$ as $`T`$ is reduced indicates that ordered states with FM order start to appear below $`T10`$ K which finally converge towards the equilibrium state at zero temperature . As reflected by the $`M_z`$ and $`M_n`$ values attained at $`T=0`$, the increase in $`k_S`$ results in a change in the kind of low $`T`$ order. For $`k_S10`$, quasi-uniform magnetization configurations are stable with the core spins pointing along the z axis ($`M_n1`$ and $`|M_z|1`$ as $`T`$ approaches zero in the right column curves of Fig. 1) and groups of surface spins slightly deviated towards the radial direction ($`M_n0.5`$ and $`|M_z|1`$ as $`T`$ approaches zero in the left column curves of Fig. 1). However, as $`k_S`$ is further increased from $`10`$, hedgehog-like structures of surface spins start to form (note that $`M_n`$ at $`T=0`$ tends to $`1`$ at the surface) inducing the deviation of core spins towards the radial direction (indicated by the approach of $`M_n`$ to $`0.5`$ and of $`M_z`$ to $`0`$ as $`T`$ decreases for core spins). Snapshots of configurations representative of low and high $`k_S`$ values are depicted in Fig. 3 (upper rows) in qualitative agreement with simulations performed by Labay et al. for a FM particle with sc lattice .

A similar phenomenon occurs in the case of a ferrimagnetic maghemite spherical particle with the real antiferromagnetic (AF) exchange constants given in the preceding section. However now, due to the dominant AF intersublattice coupling, there is a tendency to achieve antiparallel alignment among spins in different sublattices that causes a greater degree of disorder at the surface with respect to the FM case. Notice also that, even for small $`k_S`$, $`M_z`$ does not approach $`1`$ as $`T`$ tends to zero as in the FM case, but instead tends to the value of magnetization of the noncompensated spins ($`M_{Unc}=(N_ON_T)/N_{Total}`$), that for a particle of diameter $`D=3`$ is $`M_{Unc}=0.285(0.412)`$ for surface (core) spins . Moreover, for the ferrimagnetic particle, there is a change in the magnetic order at low $`T`$ as the $`k_S/k_C`$ ratio increases. When surface anisotropy is small, the particle orders into a quasi-AF state in which spins in each sublattice are almost aligned along the core easy-axis. When increasing $`k_S`$, surface spins start to form throttled structures as those seen for $`k_S`$ in Fig. 3 that depart the spins from the z-axis. Finally, for $`k_S>>k_C`$, hedgehog-like configurations are favoured by the dominant radial anisotropy contribution (see the configuration for $`k_S=100`$).

In this case, the direct visualization of equilibrium configurations presented in Fig. 3 shows that the reduction of the saturation magnetization with particle size observed experimentally in different fine particles of ferrimagnetic oxides, can be attributed to the random canting of surface spins caused by the competing antiferromagnetic interaction between sublattices . Moreover, as the results of our simulations confirm, the degree of disorder at the surface is larger than for the FM particle due to the complex interplay between the AF intralattice interactions and the local anisotropy easy-axis.

In order to study the influence of surface anisortropy on the magnetization processes, we have also simulated hysteresis loops at zero temperature. The loops have been computed by starting from a saturated state achieved after application of a high enough field of $`h=200`$ K along the z axis and decreasing the field in constant steps $`\delta h=1`$K, during which the magnetization was averaged over $`1000`$ MC steps after thermalization. The results for a spherical particle with $`D=3`$ and FM interactions are shown in Fig. 4. For a FM particle, the hysteresis loops are dominated by the surface contribution for all the values of $`k_S`$ studied as indicated by the non-squaredness of the loops around the coercive field. For high values of the surface anisotropy ($`k_S=50,100`$), a magnetic field as high as $`h=20`$ K is able to saturate the core but the surface spins instead point along the radial direction during the magnetization process. This is more clearly reflected in the right panels of Fig. 4, where we see that for high $`k_S`$, $`M_n`$ remains close to $`1`$ at the surface during all the reversal process, while the core spins depart from their easy directions ($`M_n1`$) dragged by the surface towards the radial direction close to the coercive field, where $`M_n0.5`$. The hysteresis loops for the FM particle are to be compared with those for the ferrimagnetic maghemite particle presented in Fig. 5 for $`k_S=10,100`$. First, although the recersal process is still dominated by the surface spins, notice that now the loops for the ferrimagnetic particle become more elongated and have higher closure fields, ressembling those observed experimentally . The increase of $`k_S`$ at the surface results in an increase in the coercive field $`h_C`$ and a reduction in the high field susceptibility. These features are due to the dominance of surface anisotropy over exchange interactions that create surface disordered states which become more difficult to reverse by the magnetic field. More importantly, there is a change in the magnetization reversal mechanism when increasing $`k_S`$ from $`10`$ to $`100`$. In the first case (left panels in Fig. 5), the core and surface have similar $`h_C`$ and closure fields. The particle core reverses in a quasi-uniform form with the spins pointing mostly along the z-axis ($`M_n1`$) except near $`h_c`$, and with the surface spins following the core reversal (with $`M_n<1`$ indicating alignment close to the z direction) . However, at higher $`k_S`$ (right panels in Fig. 5), surface spins remain close to the local radial direction ($`M_n1`$) during all the reversal process, driving the core spins out of their local easy axis and making their reversal non-uniform ($`0.5<M_n<1`$) due to the appearance of the hedgehog-like structures during the reversal .

In conclusion, we have studied the differences in the equilibrium states and reversal processes of ferrimagnetic with respect to FM particles, showing in both cases that finite-size effects, together with the existence of surface disorder due to increased anisotropy at the surface, are necessary ingredients for the low temperature anomalous magnetic properties observed experimentally in nanoparticle systems.

## Acknowledgements

We acknowledge CESCA and CEPBA under coordination of C4 for the computer facilities. This work has been supported by SEEUID through project MAT2000-0858 and CIRIT under project 2001SGR00066.

## FIGURE CAPTIONS

FIG. 1: Spherical ($`D=3`$) particle with FM interactions ($`J_{ij}=1`$): thermal dependence of the surface (left column) and core (right column) contributions to the projection of the magnetization along the z axis $`M_z`$ and to $`M_n`$ as defined in the text during a progressive cooling from a high $`T`$ at a constant rate $`\delta T=1`$ K.

FIG. 2: Spherical ($`D=3`$) particle with the real maghemite AF interactions quoted in the text: thermal dependence of the surface (left column) and core (right column) contributions to the projection of the magnetization along the z axis $`M_z`$ and to $`M_n`$ as defined in the text during a progressive cooling from a high $`T`$ at a constant rate $`\delta T=1`$ K.

FIG. 3: Snapshots of the equilibrium configurations attained after cooling from a high temperature disordered state for FM (upper row) and ferrimagnetic maghemite (lower row) spherical particles with $`D=3`$ and maghemite lattice structure.

FIG. 4: Hysteresis loops for a spherical particle of diameter $`D=3`$ with FM interactions. The core anisotropy constant is $`k_C=1`$, the results for several values of the surface anistropy constant are displayed $`k_S=1,5,10,20,50,100`$. The magnetization component along the applied field direction is shown in lower right panel, the core contribution is displayed in the upper right panel. Right panels show the suface and core contributions to $`M_n`$ as defined in the text.

FIG. 5: Hysteresis loops for a spherical particle of diameter $`D=3`$ with the real values of exchange constants of maghemite. The main panels show the total magnetization along the field direction and the corresponding surface (dashed lines) and core (circles) contributions. The upper subpanels show the surface (circles) and core (continuous lines) contributions to $`M_n`$. Left panel: $`k_S=10`$, right panel: $`k_S=100`$.


format: text

text:
# Random Walks on Complex Networks

## Abstract

We investigate random walks on complex networks and derive an exact expression for the mean first passage time (MFPT) between two nodes. We introduce for each node the random walk centrality $`C`$, which is the ratio between its coordination number and a characteristic relaxation time, and show that it determines essentially the MFPT. The centrality of a node determines the relative speed by which a node can receive and spread information over the network in a random process. Numerical simulations of an ensemble of random walkers moving on paradigmatic network models confirm this analytical prediction.

From biology over computer science to sociology the world is abundant in networks. Watts and Strogatz Watts.Strogatz98  demonstrated that many of these real-world networks pose characteristic features like the small world and clustering property. Since classical network models do not display these features, new models have been developed to understand the structure underlying real-world networks Albert.Barabasi02 ; Dorogovtsev.Mendes02 . The important finding is that there exists a scale-free (SF) network characterized by a power-law degree distribution Albert.Jeong.Barabasi99 , that is also a characteristic for the World-Wide-Web (WWW) Albert.Jeong.Barabasi99 , and for many other networks in various disciplines Albert.Barabasi02 .

The SF network has a heterogeneous structure. For instance, the WWW analyzed in Ref. Albert.Jeong.Barabasi99  includes nodes with degrees ranging from $`𝒪(1)`$ to $`𝒪(10^3)`$. The heterogeneity leads to intriguing properties of SF networks. In the context of percolation, SF networks are stable against random removal of nodes while they are fragile under intentional attacks targeting on nodes with high degree Cohen.E.A.H00 ; Albert.Jeong.Barabasi00 . Statistical mechanical systems on SF networks also display interesting phase transitions Pastor-Satorras.Vespignani01 ; Dorogovtsev.Goltsev.Mendes02 ; Igloi.Turban02 . In a transport process, each node in SF networks does not contribute equally likely. The importance of each node in such a process is measured with the betweenness centrality Newman01 , which has a broad power-law distribution Goh01.02 .

In this paper, we study a random walk on general networks with a particular attention to SF networks. The random walk is a fundamental dynamic process Hughes95 . It is theoretically interesting to study how the structural heterogeneity affects the nature of the diffusive and relaxation dynamics of the random walk RWref . Those issues will be studied further elsewhere Noh.Rieger03 . The random walk is also interesting since it could be a mechanism of transport and search on networks Adamic.L.P.H01 ; Guimera.etal02 ; Holme03 . Those processes would be optimal if one follows the shortest path between two nodes under considerations. Among all paths connecting two nodes, the shortest path is given by the one with the smallest number of links comment1 . However the shortest path can be found only after global connectivity is known at each node, which is improbable in practice. The random walk becomes important in the extreme opposite case where only local connectivity is known at each node. We also suggest that the random walk is a useful tool in studying the structure of networks.

In the context of transport and search, the mean first passage time (MFPT) is an important characteristic of the random walk. We will derive an exact formula for the MFPT of a random walker from one node $`i`$ to another node $`j`$, which will be denoted by $`T_{ij}`$, in arbitrary networks. In the optimal process it is just given by the number of links in the shortest path between two nodes, and both motions to one direction and to the other direction are symmetric. However, a random walk motion from $`i`$ to $`j`$ is not symmetric with the motion in the opposite direction. The asymmetry is characterized with the difference in the MFPT’s. It is revealed that the difference is determined by a potential-like quantity which will be called the random walk centrality (RWC). The RWC links the structural heterogeneity to the asymmetry in dynamics. It also describes centralization of information wandering over networks.

We consider an arbitrary finite network (or graph) which consists of nodes $`i=1,\mathrm{},N`$ and links connecting them. We assume that the network is connected (i.e. there is a path between each pair of nodes $`(i,j)`$), otherwise we simply consider each component separately. The connectivity is represented by the adjacency matrix $`𝐀`$ whose element $`A_{ij}=1(0)`$ if there is a link from $`i`$ to $`j`$ (we set $`A_{ii}=0`$ conventionally). In the present work, we restrict ourselves to an undirected network, namely $`A_{ij}=A_{ji}`$. The degree, the number of connected neighbors, of a node $`i`$ is denoted by $`K_i`$ and given by $`K_i=_jA_{ij}`$.

The stochastic process in discrete time that we study is a random walk on this network described by a master equation. The transition probabilities are defined by the following rule: A walker at node $`i`$ and time $`t`$ selects one of its $`K_i`$ neighbors with equal probability to which it hops at time $`t+1`$, thus the transition probability from node $`i`$ to node $`j`$ is $`A_{ij}/K_i`$ comment3 . Suppose the walker starts at node $`i`$ at time $`t=0`$, then the master equation for the probability $`P_{ij}`$ to find the walker at node $`j`$ at time $`t`$ is

$$P_{ij}(t+1)=\underset{k}{}\frac{A_{kj}}{K_k}P_{ik}(t).$$
(1)
The largest eigenvalue of the corresponding time evolution operator is $`1`$ corresponding to the stationary distribution $`P_j^{\mathrm{}}=lim_t\mathrm{}P_{ij}(t)`$, i.e. the infinite time limit comment2 . An explicit expression for the transition probability $`P_{ij}(t)`$ to go from node $`i`$ to node $`j`$ in $`t`$ steps follows by iterating Eq. (1)

$$P_{ij}(t)=\underset{j_1,\mathrm{},j_{t1}}{}\frac{A_{ij_1}}{K_i}\frac{A_{j_1j_2}}{K_{j_1}}\mathrm{}\frac{A_{j_{t1}j}}{K_{j_{t1}}}.$$
(2)
Comparing the expressions for $`P_{ij}`$ and $`P_{ji}`$ one sees immediately that

$$K_iP_{ij}(t)=K_jP_{ji}(t).$$
(3)
This is a direct consequence of the undirectedness of the network. For the stationary solution, Eq. (3) implies that $`K_iP_j^{\mathrm{}}=K_jP_i^{\mathrm{}}`$, and therefore one obtains

$$P_i^{\mathrm{}}=\frac{K_i}{𝒩}$$
(4)
with $`𝒩=_iK_i`$. Note that the stationary distribution is, up to normalization, equal to the degree of the node $`i`$ — the more links a node has to other nodes in the network, the more often it will be visited by a random walker.

How fast is the random walk motion? To answer to this question, we study the MFPT. The first-passage probability $`F_{ij}(t)`$ from $`i`$ to $`j`$ after $`t`$ steps satisfies the relation

$$P_{ij}(t)=\delta _{t0}\delta _{ij}+\underset{t^{}=0}{\overset{t}{}}P_{jj}(tt^{})F_{ij}(t^{}).$$
(5)
The Kronecker delta symbol insures the initial condition $`P_{ij}(0)=\delta _{ij}`$ ($`F_{ij}(0)`$ is set to zero). Introducing the Laplace transform $`\stackrel{~}{f}(s)_{t=0}^{\mathrm{}}e^{st}f(t)`$, Eq. (5) becomes $`\stackrel{~}{P}_{ij}(s)=\delta _{ij}+\stackrel{~}{F}_{ij}(s)\stackrel{~}{P}_{jj}(s)`$, and one has

$$\stackrel{~}{F}_{ij}(s)=(\stackrel{~}{P}_{ij}(s)\delta _{ij})/\stackrel{~}{P}_{jj}(s).$$
(6)
In finite networks the random walk is recurrent Hughes95 , so the MFPT is given by $`T_{ij}=_{t=0}^{\mathrm{}}tF_{ij}(t)=\stackrel{~}{F}_{ij}^{}(0)`$.

Since all moments $`R_{ij}^{(n)}_{t=0}^{\mathrm{}}t^n\{P_{ij}(t)P_j^{\mathrm{}}\}`$ of the exponentially decaying relaxation part of $`P_{ij}(t)`$ are finite, one can expand $`\stackrel{~}{P}_{ij}`$ as a series in $`s`$ as

$$\stackrel{~}{P}_{ij}(s)=\frac{K_j}{𝒩(1e^s)}+\underset{n=0}{\overset{\mathrm{}}{}}(1)^nR_{ij}^{(n)}\frac{s^n}{n!}.$$
(7)
Inserting this series into Eq. (6) and expanding it as a power series in $`s`$, we obtain that

$$T_{ij}=\{\begin{array}{cc}\frac{𝒩}{K_j},& \text{ for }j=i\hfill \\ \frac{𝒩}{K_j}\left[R_{jj}^{(0)}R_{ij}^{(0)}\right],& \text{ for }ji\hfill \end{array}.$$
(8)
A similar expression is derived in Ref. Hughes95  for the MFPT of the random walk in periodic lattices.

It is very interesting to note that the average return time $`T_{ii}`$ does not depend on the details of the global structure of the network. It is determined only by the total number of links and the degree of the node. Since it is inversely proportional to the degree, the heterogeneity in connectivity is well reflected in this quantity. In a SF network with degree distribution $`P(K)K^\gamma `$, the MFPT to the origin $`T_o`$ also follows a power-law distribution $`P(T_o)T_o^{(2\gamma )}`$. The MFPT to the origin distributes uniformly in the special case with $`\gamma =2`$.

Random walk motions between two nodes are asymmetric. The difference between $`T_{ij}`$ and $`T_{ji}`$ for $`ij`$ can be written as (using Eq. (8))

$$T_{ij}T_{ji}=𝒩\left(\frac{R_{jj}^{(0)}}{K_j}\frac{R_{ii}^{(0)}}{K_i}\right)𝒩\left(\frac{R_{ij}^{(0)}}{K_j}\frac{R_{ji}^{(0)}}{K_i}\right),$$

where the last term vanishes due to Eq. (3). Therefore we obtain

$$T_{ij}T_{ji}=C_j^1C_i^1,$$
(9)
where $`C_i`$ is defined as

$$C_i\frac{P_i^{\mathrm{}}}{\tau _i},$$
(10)
where $`P_i^{\mathrm{}}=K_i/𝒩`$ and the characteristic relaxation time $`\tau _i`$ of the node $`i`$ is given by

$$\tau _i=R_{ii}^{(0)}=\underset{t=0}{\overset{\mathrm{}}{}}\{P_{ii}(t)P_i^{\mathrm{}}\}.$$
(11)
We call $`C_i`$ the random walk centrality since it quantifies how central a node $`i`$ is located regarding its potential to receive informations randomly diffusing over the network. To be more precise: Consider two nodes $`i`$ and $`j`$ with $`C_i>C_j`$. Assume that each of them launches a signal simultaneously, which is wandering over the network. Based on Eq. (9), one expects that the node with larger RWC will receive the signal emitted by its partner earlier. Hence, the RWC can be regarded as a measure for effectiveness in communication between nodes. In a homogeneous network with translational symmetry, all nodes have the same value of the RWC. On the other hand, in a heterogeneous network the RWC has a distribution, which leads to the asymmetry in the random dynamic process.

The RWC is determined by the degree $`K`$ and $`\tau `$. The order of magnitude of the characteristic relaxation time $`\tau `$ is related to the second largest eigenvalue (nota bene comment2 ) of the time evolution operator in (1): $`P_{ii}(t)=P_i^{\mathrm{}}+_{\lambda =2}^Na_i^{(\lambda )}b_i^{(\lambda )}\mathrm{\Lambda }_\lambda ^t`$, where $`a^{(\lambda )}`$ and $`b^{(\lambda )}`$ are the left and right eigenvectors, respectively, of the time evolution operator belonging to the eigenvalue $`\mathrm{\Lambda }_\lambda `$. If we order the eigenvalues according to the modulus ($`|\mathrm{\Lambda }_2||\mathrm{\Lambda }_3|\mathrm{}|\mathrm{\Lambda }_N|`$) the asymptotic behavior is $`P_{ii}(t)P_i^{\mathrm{}}a_i^{(2)}b_i^{(2)}\mathrm{\Lambda }_2^t`$ and $`\tau _ia_i^{(2)}b_i^{(2)}/|\mathrm{ln}|\mathrm{\Lambda }_2||`$. Thus the relaxation time $`\tau _i`$ has a node dependence only through the weight factor, which is presumably weak. On the other hand, the degree dependence is explicit.

We examined the distribution of the RWC in the Barabási-Albert (BA) network Barabasi.Albert99 . This is a model for a growing SF network; at each time step, a new node is added creating $`m`$ links with other nodes which are selected with the probability proportional to their degree. We grew the network, solved the master equation numerically with the initial condition $`P_i(t=0)=\delta _{ik}`$, and calculated the relaxation time $`\tau _k`$ for each $`k`$. Figure 1 (a) shows the plot of $`\tau `$ vs. $`K`$ in the BA network of $`N=10^4`$ nodes grown with the parameter $`m=2`$. The degree is distributed broadly over the range $`2K400`$. On the other hand, the relaxation time turns out to be distributed very narrowly within the range $`1\tau 2`$. We also studied BA networks of different sizes, but did not find any significant broadening of the distribution of $`\tau `$. So the RWC distribution is mainly determined by the degree distribution. In Fig. 1 (b) we show the plot of $`C`$ vs. $`K`$ in the same BA network. It shows that the RWC is roughly proportional to the degree. Note however that the RWC is not increasing monotonically with the degree due to the fluctuation of $`\tau `$ as seen in Fig. 1 (a).

The RWC is useful when one compares the random walk motions between two nodes, e.g., $`i`$ and $`j`$ with $`C_i>C_j`$. On average a random walker starting at $`j`$ arrives at $`i`$ before another walker starting at $`i`$ arrives at $`j`$. Now consider an intermediate node $`k`$, which may be visited by both random walkers. Since $`T_{ij}>T_{ji}`$, it is likely that a random walker starting at node $`k`$ will arrive at node $`i`$ earlier than at node $`j`$. Although this argument is not exact since we neglected the time spent on the journey to the intermediate node, it indicates that nodes with larger RWC may be typically visited earlier than nodes with smaller RWC by the random walker. If we interpret the random walker as an information messenger, nodes with larger RWC are more efficient in receiving information than nodes with smaller RWC.

We performed numerical simulations to study the relation between the RWC and this efficiency. To quantify it, we consider a situation where initially all nodes in a network are occupied by different random walkers. They start to move at time $`t=0`$, and we measure $`n_i`$, the fraction of walkers which have passed through the node $`i`$, as a function of time $`t`$. It is assumed that the walkers do not interact with each other. They may be regarded as a messenger delivering an information to each node it visits. Then, with the information distribution uniformly initially, $`n_i`$ is proportional to the amount of information acquired by each node. The argument in the previous paragraph suggests that typically nodes with larger values of RWC have larger value of $`n_i`$ at any given time.

The BA network Barabasi.Albert99  and the hierarchical network of Ravasz and Barabási cite:hnet  were considered in the simulations. The hierarchical network is a deterministic network growing via iteration; at each iteration the network is multiplied by a factor $`M`$. The emergent network is scale-free when $`M3`$. Since it is a deterministic network, several structural properties are known exactly Noh03 . We measured $`n_i`$ in the BA network with $`m=2`$ and $`N=512`$ nodes and in the hierarchical network with $`M=5`$ and $`N=5^4`$ nodes for $`0t2048`$, which are presented in the left and the right column of Fig. 2, respectively. The value of $`n_i`$ is color-coded according to the reference shown in Fig. 2.

The time evolution of $`n_i`$ is presented in three different ways. In the first row, the nodes are arranged in ascending order of the node index $`i`$. In the BA network, the node index corresponds to the time step at which the node is added to the network. The indexing scheme for the hierarchical network is explained in Ref. Noh03 . In the second row, the nodes are arranged in descending order of the degree $`K`$ and in the third row they are arranged in descending order of the RWC $`C`$. At a given time $`t`$, the plot in the first row shows that $`n`$ is non-monotonous and very irregular as a function of the node index. As a function of the degree it becomes smooth, but still non-monotonic tendencies remain. However, as a function of the RWC, it becomes much smoother and almost monotonous. We calculated for each node $`i`$ the time $`\tau _i^{}`$ at which $`n_i`$ becomes greater than $`1/2`$. In the BA network, among all node pairs $`(i,j)`$ satisfying $`\tau _i^{}<\tau _j^{}`$, only $`3\%`$ violate the relation $`C_i>C_j`$, whereas the number of pairs that violate the relation $`K_i>K_j`$ is five times larger.

In summary, we studied the random walk processes in complex networks. We derive an exact expression for the mean first passage time (see Eq. (8)). The MFPT’s between two nodes differ for the two directions in general heterogeneous networks. We have shown that this difference is determined by the random walk centrality $`C`$ defined in Eq. (10). Among random walk motions between two nodes, the walk to the node with larger value of $`C`$ is faster than the other. Furthermore, it is argued that in a given time interval nodes with larger values of $`C`$ are visited by more random walkers which were distributed uniformly initially. We confirmed this by numerical simulations on the BA and the hierarchical network. One may regard the random walkers as informations diffusing through the network. Our results imply that information does not distribute uniformly in heterogeneous networks; the information is centralized to nodes with larger values of $`C`$. The nodes with high values of $`C`$ have the advantage of being aware of new information earlier than other nodes. On the other hand, it also implies that such nodes are heavily loaded within an information distribution or transport process. If the network has a finite capacity, the heavily loaded nodes may cause congestions Holme03 . Therefore much care should be taken of the nodes with high $`C`$ values in network management. In the current work, we consider the random walks on undirected networks. The generalization to directed networks would be interesting. And in order to study congestion, the random walk motions with many interacting random walkers would also be interesting. We leave such generalizations to a future work.

Acknowledgement: This work was supported by the Deutsche Forschungsgemeinschaft (DFG) and by the European Community’s Human Potential Programme under contract HPRN-CT-2002-00307, DYGLAGEMEM.


text:
# Are current-induced forces conservative?

## Abstract

The expression for the force on an ion in the presence of current can be derived from first principles without any assumption about its conservative character. However, energy functionals have been constructed that indicate that this force can be written as the derivative of a potential function. On the other hand, there exist compelling specific arguments that strongly suggest the contrary. We propose physical mechanisms that invalidate such arguments and demonstrate their existence with first-principles calculations. While our results do not constitute a formal resolution to the fundamental question of whether current-induced forces are conservative, they represent a substantial step forward in this direction.

Electrons flowing in a conductor can transfer momentum to the ions. From a static point of view, the steady-state electron density in a current-carrying conductor is microscopically different from that in the absence of current. In either picture, the passage of current generates additional forces on ions. This may result in atomistic rearrangements and diffusion of atoms, a phenomenon known as electromigration. sorbe  The advent of microelectronics has sparked off a considerable interest in electromigration since it has been found to be a major failure mechanism in metallic interconnects. blech  In recent years, a new type of electronics is being considered which envisions the use of nanoscale components. Since nanoscale conductors can carry current densities much larger than their macroscopic counterparts, the problem of electromigration is again of key importance. Recent work has studied the effect of current-induced forces in atomic lango ; agr ; ita ; yas ; tod1 ; mioyang ; mads  and molecular wires aya ; md1 . Combined with early investigations on the subject, locdipole  this work has shed new light on the microscopic mechanism of current-induced mechanical effects.

While a consensus on the physical origin of current-induced forces seems to have been reached, a fundamental question remains unanswered: are current-induced forces conservative? The difficulty in answering this question is intimately related to the difficulty in defining a force in a non-equilibrium problem. md2 ; tod00

In this paper we address this question in the following way: we consider representative arguments against the conservative nature of current-induced forces and suggest physical mechanisms that invalidate these arguments. We then corroborate the existence of such mechanisms with first-principles calculations. While this by no means proves that these forces are conservative, we believe it is a substantial step forward in this direction and hope that our suggestions will generate much-needed discussions. prec

Let us start from the expression for forces in the presence of current. This expression can be derived in different ways. One approach is to start from a fully dynamical quantum description of the many-body problem comprising electrons and ions, and then take the limit as the ionic motion becomes classical. md2  In this derivation, the force expression comes about as the classical limit of Ehrenfest’s theorem applied to the rate of change of the momentum of an ion. md2  Formally, this limit corresponds to high ionic masses or high kinetic energies tod01 . A different approach is to start from an appropriate mixed quantum-classical Lagrangian. The force expression then results from the Euler-Lagrange equation of motion for the classical ions tod01 . Both approaches result in the same expression for the force on a given ion, with position $`𝐑`$, due to the self-consistent electronic density $`\rho (𝐫)`$ under current flow tod01 ,

$$𝐅=𝑑𝐫\frac{v}{𝐑}\rho (𝐫),$$
(1)
which can be equivalently written as md2

$$𝐅=\underset{i}{}\psi _i\left|\frac{v}{𝐑}\right|\psi _i\underset{\mathrm{\Delta }0}{lim}_\sigma 𝑑E\psi _\mathrm{\Delta }\left|\frac{v}{𝐑}\right|\psi _\mathrm{\Delta }.$$
(2)
where $`v`$ is the electron-ion interaction. comment  The first term on the RHS of Eq. (2) is similar to the usual Hellmann-Feynman contribution to the force due to localized electronic states $`|\psi _i`$. The second term is the contribution to the force due to continuum states. md2  The wavefunctions $`|\psi _\mathrm{\Delta }`$ are eigendifferentials in the continuum $`\sigma `$md2  Due to the algebraic equivalence of the above two equations, in the following, when discussing current-induced forces, we will refer to either one without distinction.

Let us now turn to the question of whether or not these forces are conservative. It is essential to emphasize that the force obtained either from the Ehrenfest’s theorem md2  or a Lagrangian formulation tod01 , is not necessarily conservative because the above first-principles dynamical derivations do not contain any assumption about the existence, or otherwise, of a free energy functional, of which this force might be the derivative. However, as a separate exercise, it was suggested both in Ref. md2,  and in Ref. tod00,  that forces in the presence of current can also be written as derivatives of appropriate energy functionals. For instance, in Ref. md2,  the stationarity of an action was used to map the force problem into a variational one in which the force can be written as the derivative of an energy function. In Ref. tod00,  an entropy function was used, whose maximization (constrained to the energy of the current-carrying system and the numbers of electrons in it that originate from the right and left reservoirs) again yields the force as the derivative of a potential function. These independent derivations suggest that current-induced forces are conservative.

On the other hand, there exist compelling specific arguments that strongly suggest that current-induced forces cannot be conservative. These arguments consist in gedanken experiments, aimed at showing that the work done by current-induced forces along specific closed paths is not zero and that, therefore, these forces are not conservative. Below we refer to one notable such construction due to Sorbello. sorbe  This example is of charming simplicity and is particularly convincing. The argument goes like this. Picture a single defect atom in a current-carrying wire prec1  between bulk electrodes under an applied bias $`V`$. First, let the atom move, quasi statically, from a point $`A`$ to a point $`B`$ inside the wire (see schematic in Fig. 1). The current-induced force on the defect in the wire does some non-zero work on the way. Now imagine taking the atom from $`A`$ to $`B`$ along a path that goes outside the wire (see Fig. 1). Since outside the wire there apparently is no current, no current-induced force should be present. Hence, the work done by current-induced forces on the atom along the path from $`A`$ to $`B`$ inside the wire cannot be the same as along the path outside the wire. This would be sufficient to prove that current-induced forces are not conservative.

However, there is something missing in this argument: the quantum-mechanical nature of current outside the wire. To see this, picture the atom at a point in space outside the wire. Let $`|\psi `$ be the global current-carrying electron states for the system. If there always is some (however small) quantum-mechanical coupling between the atom and the wire, there will always be some non-zero tunneling probability for electrons to get from the wire onto the atom, and back. The atom then acts as a resonant center and there will always exist an energy window where the states $`|\psi `$ contain resonant paths that emerge from the wire, pass through the atom, and go back into the wire. Due to the scattering of electrons at the atom, these resonant paths may be expected to create a charge dipole, centered at the atom, much like the local resistivity dipole in electromigration. sorbe  Increasing the distance $`d`$ between the atom and the wire will make the resonant energy window narrower, while the resonance will remain pinned against the energy window $`eV`$, due to self-consistency. Now, the total current carried by the resonance is given by the energy integral, over the energy window $`eV`$, of a transmission function, $`T(E)`$. As the resonance gets narrower in energy, so does $`T(E)`$. $`T(E)`$ is always bounded and hence, as $`d`$ increases, the total current carried by the resonance will decrease. On the other hand, the charge dipole across the resonant center originates from the density of states, $`\mathrm{DOS}(E)`$, due to the populated current-carrying electron states in the energy window $`eV`$. By contrast with $`T(E)`$, $`\mathrm{DOS}(E)`$ is unbounded: if $`d`$ increases, the resonant peak in $`\mathrm{DOS}(E)`$ gets narrower but taller in such a way that its energy integral remains constant. Since current-induced forces are related to the integral of the density of states \[see e.g. Eq. (2)\], we expect the force on the atom due to the resonant dipole to be constant with $`d`$. In other words, there is always a current-induced force on the atom, even outside the wire, and this force does not vanish with increasing $`d`$.

A key assumption of our reasoning is that there always is some coupling between the atom in the vacuum and the metal. If this coupling is truncated artificially, then the current and the “resonant-dipole force” will both be discontinuously extinguished. However, such a truncation would at the same time create a bound state (at the atom) in the continuum, and would thus modify the nature of the electronic spectrum and of the corresponding Born-Oppenheimer (BO) manifold. Such a discontinuous modification corresponds to a breakdown of the BO approximation, in which case the notion of conservative forces is lost by definition.

Our analysis above relies only on generic properties of resonant transmission. We therefore now seek to demonstrate the existence of a resonant dipole in the simplest possible resonant configuration: that of a single atom in the vacuum between two bulk electrodes. Our goal is to show that this resonant dipole gives rise to a finite current-induced force that does not vanish with increasing electrode separation, even when the current dies out. To this end, we have performed first-principles calculations of the current-induced force on a single Si atom between two jellium electrodes. prec12  The interior electron density of the electrodes is taken equal to the value for metallic gold ($`r_s3`$). prec13  The electrodes are kept at a certain bias with the left electrode positively biased so that electrons flow from right to left. In Fig. 2 we plot the density of states as a function of energy for a bias of 0.1 V and different distances between the jellium edges of the two electrodes. It is clear from this figure that a resonance develops in the energy window between the left and right chemical potentials. The resonance becomes narrower and taller with increasing distance. It is also evident from Fig. 2 that this resonance is pinned towards the right chemical potential. Such pinning can be intuitively understood as follows. By increasing the distance between the electrodes the resonance becomes narrower and narrower. In order to preserve charge neutrality, the amount of charge accommodated by the resonance needs to be the same at all distances. The only way to achieve this is for the “center” of the resonance to move closer to the right electrochemical potential.

Now that we have shown the existence of the resonance, let us look at the spatial distribution of the charge dipole it creates across the atom. In order to do this we follow Lang’s approach lango  and calculate first the difference between the electronic distribution at a given bias $`\rho (𝐫,V)`$ and the corresponding distribution at zero bias $`\rho (𝐫,0)`$. Call this quantity $`\stackrel{~}{\rho }(𝐫,V)`$. It contains three contributions: (i) the polarization of the biased electrode-electrode system without the atom, (ii) the polarization of the electron cloud at the atom due to the field generated by (i), and (iii) the current-induced resonant dipole. The first two terms generate a force of equal magnitude but opposite sign so they do not contribute to the total force. (An isolated atom in a uniform electric field does not have a net force on its nucleus). To eliminate (i) and (ii), we subtract (i) from $`\stackrel{~}{\rho }(𝐫,V)`$ and take the difference between the result of that subtraction at a given electrode-electrode distance and the same quantity at a larger separation scaled to the field of the previous distance. lango  This gives a negative fraction of the current-induced dipole, so a change of sign is necessary. Call this quantity $`\rho _f(𝐫,V)`$. Finally, to isolate the antisymmetric (dipole) component of this distribution we evaluate $`[\rho _f(𝐫,V)\rho _f(𝐫,V)]/2`$. The result is plotted in Fig. 3 for a separation of 10 a.u. (with respect to 16 a.u.) and clearly shows the current-induced resonant dipole across the atom. The dipole has positive charge on the left of the atom and negative charge on its right so that the current-induced force pushes the atom to the right, i.e. against the electron flow.

Finally, we have to show that the current decreases with electrode separation while the force does not. This is plotted in Fig. 4 for different electrode separations. It is clear from Fig. 4 that the current decreases substantially with electrode separation while the force saturates at a value of about 2 mRy/a.u. Notice again that this force is positive, i.e. it pushes the atom against the electron flow, in agreement with the sign of the resonant dipole. These numerical results confirm our conjecture: as long as the atom is coupled to the electrodes, a resonant dipole forms across it, such that a non-zero current-induced force is exerted on the atom that persists with separation even as the current dies out.

We conclude by stressing again that our goal here is to show that existing arguments against the conservativeness of current-induced forces ignore specific physical mechanisms. These mechanisms predict the existence of current-induced forces in unexpected places, such as an atom in the vacuum next to a current-carrying metal. We have demonstrated the existence of such forces using first-principles calculations. These results do not prove that current-induced forces are conservative but we believe that they give extra support in favor of such a conclusion. More theoretical work is clearly necessary to fully resolve this fundamental problem.

We thank J. Lagerqvist for help in generating the figures in this paper. We are grateful to Adrian Sutton for invaluable discussions and for the inspiration he has provided through his deep and long-standing interest in this problem. Two of us (MD and YCC) acknowledge support from the National Science Foundation Grants No. DMR-01-02277 and DMR-01-33075, Carilion Biomedical Institute, and the Petroleum Research Fund. The calculations reported in this paper were performed on the beowulf cluster of the Laboratory for Advanced Scientific Computing and Applications at Virginia Tech.


text:
# Charge stripes due to electron correlations in the two-dimensional spinless Falicov-Kimball model

## I Introduction

We find it fitting to write a paper on the spinless Falicov-Kimball (FK) model Falicov and Kimball (1969) to celebrate Elliott Lieb’s seventieth birthday. Elliott, and his collaborators, provided two seminal results on this model: (i) the first, with Tom Kennedy, proved that there was a finite temperature phase transition to a checkerboard charge-density-wave (CDW) phase in two or more dimensions for the symmetric half-filled case Kennedy and Lieb (1986); Lieb (1986), and (ii) the second, with Daniel Ueltschi and Jim Freericks, proved that the segregation principle holds for all dimensions Freericks et al. (2002a, b) (which states that if the total particle density is less than one, then the ground state is phase separated if the interaction strength is large enough Freericks and Falicov (1990)). The Kennedy-Lieb result (along with an independent Brandt-Schmidt paper Brandt and Schmidt (1986, 1987)) inspired dozens of follow-up papers by researchers across the world. The Freericks-Lieb-Ueltschi paper generalized Lemberger’s proof Lemberger (1992) from one dimension to all dimensions, which finally proved the decade old Freericks-Falicov conjecture Freericks and Falicov (1990). Both papers are important, because they are the only examples where long-range order and phase separation can be proved to occur in a correlated electronic system.

The Falicov-Kimball model has an interesting history too. Leo Falicov and John Kimball invented the spin-one-half version of the model in 1969 to describe metal-insulator transitions of rare-earth compounds Falicov and Kimball (1969). It turns out that John Hubbard actually “discovered” the spinless version of the FK model four years earlier in 1965 Hubbard (1965), when he developed the alloy-analogy solution to the Hubbard model Hubbard (1963) (the so-called Hubbard III solution). This latter version was rediscovered by Kennedy and Lieb in 1986 Kennedy and Lieb (1986) when they formulated it as a simple model for how crystallization can be driven by the Pauli exclusion principle.

In this contribution, we focus on another problem that can be analyzed in the FK model—the problem of stripe formation in two dimensions. The question of the relation between charge stripes, correlated electrons, and high-temperature superconductivity has been asked ever since static stripes were first seen in the nickelate Chen et al. (1993); Tranquada et al. (1994); Sachan et al. (1995); Tranquada et al. (1995a) and cuprate Tranquada et al. (1995b); Tranquada et al. (1996, 1997); Mook et al. (2002) materials starting in 1993. Two schools of thought emerged to describe the theoretical basis for stripe formation in the Hubbard model. The Kivelson-Emery scenario Emery et al. (1990); Carlson et al. (1998); Pryadko et al. (1998, 1999) says that at large $`U`$ the Hubbard model is close to a phase separation instability but the long-range Coulomb force restricts the phase separation on the nanoscale; a compromise results in static stripe-like order. The Scalapino-White scenario White and Scalapino (1998a, b); White and Scalapino (1999, 2000) says that stripes can form due to a subtle balance between kinetic-energy effects and potential-energy effects, mediated by spin fluctuations. No long range Coulomb interaction or phase separation is needed to form these stripes. There are numerous numerical studies that have tried to shed light onto this problem. Unfortunately, they have conflicting results. High temperature series expansions on the related $`tJ`$ model Putikka et al. (1992); Putikka and Luchini (2000) show that phase separation exists, but only when $`J`$ is large enough, so it is not present in the large-correlation-strength limit of the Hubbard model (where $`J0`$). Monte Carlo calculations Cosentini et al. (1998); Calandra et al. (1998); Hellberg and Manousakis (2000) and exact diagonalization studies Dagotto (1994); Hellberg and Manousakis (1999) give different results: some calculations predict the stripe formation to occur, others show a linkage between the stripe formation and the boundary conditions selected for the problem. Mean-field-theory analyses Zaanen and A. M. Oleś (1996); D. Góra et al. (1999) seem to predict stripe formation without any phase separation. One way to make sense of these disparate results is that both the energy of the intrinsic stripe phases and the energy for phase separation are quite close to one another, so any small change (induced by finite-size effects, statistical errors, effects of correlations not included in the perturbative expansions, or due to terms dropped or added to the Hamiltonian) can have a large effect on the phase diagram by producing a small relative change in the energetics of the different many-body states (because of their near degeneracy).

We take an alternate point of view here. We choose to examine a model that can be analyzed rigorously, and can be continuously connected to the Hubbard model. We choose the regime that connects directly with the $`S_z=0`$ states of the Hubbard model. The model we analyze is the spinless Falicov-Kimball model on a square lattice

$$=t\underset{ij}{}(c_i^{}c_j+c_j^{}c_i)+U\underset{i=1}{\overset{|\mathrm{\Lambda }|}{}}w_ic_i^{}c_i,$$
(1)
where $`c_i^{}`$ ($`c_i`$) creates (destroys) a spinless conduction electron at site $`i`$, $`t`$ is the hopping matrix element ($`ij`$ denotes a summation over nearest-neighbor pairs on a square lattice), $`w_i=0`$ or 1 is a classical variable denoting the localized electron number at site $`i`$, and $`U`$ is the on-site Coulomb interaction energy. The Fermionic operators satisfy anticommutation relations $`(c_i^{},c_j^{})_+=0`$, $`(c_i,c_j)_+=0`$, and $`(c_i^{},c_j)_+=\delta _{ij}`$. The symbol $`|\mathrm{\Lambda }|`$ denotes the total number of lattice sites in the square lattice $`\mathrm{\Lambda }`$. We will always be dealing with periodic configurations of localized electrons, which means we can always consider our lattice to have a large but finite number of lattice sites and periodic boundary conditions. A short presentation of these results has already appeared Lemański et al. (2002).

The Falicov-Kimball model can be viewed as a Fermionic quantum analogue of the Ising model, while the Hubbard model can be viewed as the Fermionic quantum analogue of the Heisenberg model (indeed in the large-$`U`$ limit at half filling, the Falicov-Kimball model maps onto an effective nearest-neighbor Ising model, while the Hubbard model maps onto an effective nearest-neighbor Heisenberg model). The way to link the Falicov-Kimball model to the Hubbard model is to imagine a generalization of the Hubbard model where the down-spin hopping matrix element differs from the up-spin hopping matrix element. Then as $`t_{}0`$, the down spins become heavy and are localized on the lattice; the quantum-mechanical ground state is determined by the configuration of down-spin electrons that minimizes the energy of the up-spin electrons. This is precisely the Falicov-Kimball model!

In order to maintain the connection to the Hubbard model in zero magnetic field, we must choose the conduction electron density $`\rho _e=_{i=1}^{|\mathrm{\Lambda }|}c_i^{}c_i/|\mathrm{\Lambda }|`$ to be equal to the localized electron density $`\rho _f=_{i=1}^{|\mathrm{\Lambda }|}w_i/|\mathrm{\Lambda }|`$, which we do here. We study the evolution from the checkerboard phase at half filling ($`\rho _e=\rho _f=1/2`$) to the segregated phase, which appears when $`\rho _e=\rho _f`$ is small enough. Since these two phases are drastically different from each other, the transition is likely to include many different intermediate phases. Indeed, the ground state phase diagram of the Falicov-Kimball model can be quite complex. There are many different periodic phases that can be stabilized for different values of $`U`$ or $`\rho _e=\rho _f`$. As $`U`$ becomes large though, the phase diagram simplifies, as the segregated phase becomes the ground state for wider and wider ranges of the electron densities.

## II Formalism

Our strategy to examine the FK model is a brute-force approach which is straightforward to describe, but tedious to carry out. We employ the so-called restricted phase diagram approach, where we consider the grand-canonical thermodynamic potential of the system for all possible periodic phases of the localized electrons, selected from a finite set of candidate phases. In this work, we consider 23,755 phases, which corresponds to the set of all inequivalent phases with a unit cell that includes 16 or fewer lattice sites. In order to calculate the thermodynamic potential, we first must determine the electronic band structure for the conduction electrons for each candidate periodic phase. We employ a Brillouin-zone grid of $`110\times 110`$ momentum points for each bandstructure. This requires us to diagonalize up to $`16\times 16`$ matrices at each discrete momentum point in the Brillouin zone and results in at most 16 different energy bands. Hence, our calculations can be viewed as finite-size cluster calculations with cluster sizes ranging from $`110\times 110\times 1`$ up to $`110\times 110\times 16`$ depending on the number of atoms in the unit cell. An example of such a bandstructure is shown in Fig. 1. The eigenvalues of the band structure are summed to determine the ground-state energy for each number of conduction electrons. The Gibbs thermodynamic potential is then calculated for all possible values of the chemical potentials of the conduction and localized electrons through the formula

$$G(\{w_i\})=\frac{1}{12100N_0}\underset{ϵ_j<\mu _e}{}ϵ_j(\{w_i\})\mu _e\rho _e\mu _f\rho _f,$$
(2)
with $`\mu _e`$ and $`\mu _f`$ denoting the chemical potentials for the conduction and localized electrons, respectively, and $`N_0`$ denoting the number of atomic sites in the unit cell for the given configuration of localized electrons. The symbol $`ϵ_j(\{w_i\})`$ denotes the energy eigenvalues of the band structure for the given configuration of localized electrons. Since the thermodynamic potential is concave, the phase diagram can be directly determined in the chemical potentials plane Watson and Lemański (1995); Gajek et al. (1996a, b). Next, we convert the grand canonical ensemble into a canonical ensemble to determine the ground-state phase diagram as functions of $`\rho _e`$ and $`\rho _f`$. We find the ground state is often a phase separated mixture of two or three different phases, which can be periodic phases, or the segregated phase. This step of the analysis is quite complicated, because small areas of stability in the grand canonical phase diagram can correspond to large regions in the canonical phase diagram, and vice versa. Finally, we restrict the analysis to the case $`\rho _e=\rho _f`$ and plot the phase diagram as a function of the total filling for each chosen value of $`U`$. This computational algorithm is illustrated schematically in Fig. 2.

We find that of the initial 23,755 candidate phases, only 111 can be found in the ground-state phase diagram for the values of $`U`$ that we considered. Any phase energetically excluded from appearing in the restricted phase diagram must also be excluded from the complete phase diagram. What we do not know is how our computed phase diagram will change as more candidate phases are introduced (although the majority of these additional phases also won’t appear in the phase diagram).

## III Results

The different phases that are stabilized in our restricted phase diagram can be grouped into different families that represent different types of geometric arrangements of the localized electrons. Unfortunately, there is no way to rigorously categorize these phases, so the grouping we have chosen arises in part from our personal taste in determining which phases appear most similar. Nevertheless, the groupings we have made are in some sense “obvious”, and we believe the analysis presented here is a useful way to categorize and summarize the data. We will concentrate on describing different kinds of striped phases that are present in the phase diagram and we will motivate some of the physical principles behind their appearance in the phase diagram.

We separate the different stable phases into 10 different groups. Every stable phase is labeled by a number and depicted in Fig. 3. The small dots indicate the absence of a localized electron, while the large dots indicate the positions of the localized electrons. In the lower left corner, we shade in the unit cell of the configuration and we show with the two solid lines the translation vectors of the unit cell that allow the square lattice to be tiled by the unit cell. The different families of configurations are as follows: (i) the empty lattice ($`\rho _e0`$ and $`\rho _f=0`$) denoted E which contains no localized electrons (configuration 1); (ii) the full lattice ($`\rho _e=0`$ and $`\rho _f=1`$) denoted F which contains a localized electron at each site (configuration 2); (iii) the checkerboard phase ($`\rho _e=\rho _f=1/2`$) denoted Ch which has the localized electrons occupying the A sublattice only of the square lattice in a checkerboard arrangement (configuration 3); (iv) diagonal non-neutral stripe phases ($`\rho _e1\rho _f`$) denoted DS which consist of diagonal checkerboard phases separated by empty diagonal stripes of slope 1 (configuration 4); (v) axial non-neutral checkerboard stripes ($`\rho _e1\rho _f`$) denoted AChS which consist of checkerboard regions arranged in stripes oriented parallel to the $`x`$-axis and separated by empty stripes with slope 0 (configurations 5–10); (vi) diagonal neutral stripe phases ($`\rho _e=1\rho _f`$) denoted DNS which consist of localized electrons arranged in the checkerboard phase separated by fully occupied striped regions of slope 1, or equivalently, checkerboard phases with diagonal antiphase boundaries (configurations 11–19); (vii) axial non-neutral stripe phases ($`\rho _e1\rho _f`$) denoted AS which consist of fully occupied vertical (or horizontal) stripes separated by empty stripes, which are translationally invariant in the vertical (or horizontal) direction (configurations 20–54); (viii) neutral phases ($`\rho _e=1\rho _f`$) denoted N which consist of neutral phases in an arrangement that does not look like any simple stripe phase (some neutral phases can be described in a stripe picture, such as configuration 61 which has a slope 1/3 empty lattice stripe, but we prefer to refer to them as non-stripe phases) (configurations 55–70); (ix) four-molecule phases ($`\rho _e1\rho _f`$) denoted 4M which can be described as a “bound” four-molecule square of empty sites tiled inside an occupied lattice framework (configurations 71–74); (x) two-dimensional non-neutral phases ($`\rho _e1\rho _f`$) denoted 2D which consist of phases with the localized electrons arranged in a fashion that is not stripe-like and requires a two-dimensional unit cell to describe them (once again, some phases like configuration 75 could be described as a slope 3/2 stripe, but appears to us more like a 2D phase) (configurations 75–111).

Generically, we find the canonical phase diagram does not contain pure phases from one of the 111 stable phases, but rather forms mixtures of two or three periodic phases, or one or two periodic phases and the empty lattice (which is often needed to get the conduction-electron filling correct in the mixture). When we are doped sufficiently far from half filling, we are in the segregated phase, which is a mixture of the E and F phases.

We consider 5 different values of $`U`$ in our computations: $`U=1`$, 2, 4, 6, and 8. The phase diagram is quite complex, with many of the different 111 phases appearing for different values of $`U`$. We summarize which phases appear in Table 1.

We begin our discussion with the weak-coupling value $`U=1`$ where 50 phases appear. The phase diagram is summarized in Fig. 4. We use a solid line to indicate the region of the particle density where a particular phase appears in the ground state (either as a pure phase or as a mixture). The phases that appear in a mixture at a given density are found by determining the solid vertical lines that intersect a horizontal line drawn to pass through the given particle density. The phase diagram has shading included to separate the regions of the different categories of phases. The numeric labels are shown to make it easier to determine the actual phases present in the diagram. We plot similar phase diagrams for $`U=2`$ (38 phases), 4 (42 phases), 6 (30 phases), and 8 (25 phases) in Figs. 58, respectively. A schematic phase diagram that illustrates the generic features of the phase diagram in the electron density, interaction-strength plane appears in Fig. 9.

As can be seen from these figures, the generic phase diagram is quite complex, and by looking at the different phases in Fig. 3, many of the phases have stripe-like structures to them. To begin our discussion of these results, we must first recall the rigorous results known for this model. When $`\rho _e=\rho _f=1/2`$, the ground state is the checkerboard phase (configuration 3) for all $`U`$. This can be seen in all of the phase diagrams plotted. When $`\rho _e=\rho _f1/2`$, the ground state becomes the phase separated segregated phase when $`U`$ is large enough. So there is a simplification in the phase diagram as we increase $`U`$, and the most complex phase diagram appears in the $`U0`$ limit. That limit is also the most difficult computationally, because the differences in the energies between different configurations also becomes small for small $`U`$, and the numerical accuracy must be huge in order to achieve trustworthy results. This is why we do not report any phase diagrams with $`U<1`$ here.

Looking at the $`U=8`$ case shown in Fig. 8, we see that as we move away from half filling, we initially find mixtures between the checkerboard phase, other diagonal stripe phases, and the empty lattice. When we examine the structure factors associated with the diagonal stripe phases, we find that they tend to have more weight along the Brillouin zone diagonal than elsewhere. Hence, these diagonal stripe phases are being stabilized by a “near-nesting” instability of the noninteracting Fermi surface, and the overall mixtures are required to maintain the average fillings of the conduction and localized electrons. As we move farther from half filling, the checkerboard phase disappears from the mixtures, and then a series of neutral phases enter the mix which retain some appearance of diagonal stripes, but with more and more “defects” to the stripes that make them look more two-dimensional. We find the localized electron density of these phases increases as we reduce the total filling, which is what we expect as we move toward the segregated phase which involves a mixture of the E and F configurations. Note that the formation of many different stripe phases, occurs without needing the long-range Coulomb interaction to oppose the tendency towards phase separation, when we are close to half filling. Indeed, the ground state is often a phase separated mixture, but it is a mixture of stripe-like phases, which occur automatically, without the need to add any other physics to the system. This regime, is the closest to the Kivelson-Emery picture, but we see it has more complex behavior than what they envisioned when they examined the Hubbard model.

Moving on to the $`U=6`$ case in Fig. 7, we find a significant change in the phase diagram. The grouping of diagonal stripes near half filling disappears and we instead find the ground state to initially be a mixture between the checkerboard phase, a truly two-dimensional screen-like phase (configuration 108) and the empty lattice. Here, if we include a long-range Coulomb interaction, we would likely form diagonal stripes, but the mixture would be more complicated because it would include this screen-like structure as well. As we dope further away, we see a smaller number of the neutral phases, which look somewhat like diagonal stripes with a large number of defects in them, and then we go to a very different class of mixtures, dominated by the presence of the axial stripe phase in configuration 33. As that phase becomes destabilized, we find a cascade of many other axial stripes entering, before the segregated phase takes over. This transition from diagonal stripes to axial stripes as a function of the electron filling, also occurs because of a “near-nesting” effect. The structure factors of the axial stripe phases are peaked predominantly along the zone edge, and as we dope further from half filling, this is where nesting is more likely to occur. The cascade of stable phases that enter after configuration 33 is destabilized, have a progression of the peaks in their structure factor moving towards the zone center, which is also expected, since they are progressively heading towards the segregated phase. A similar kind of transition from diagonal stripes to axial stripes is seen in the Hubbard model studies, with the critical density lying near 0.375, as we see here too.

By the time we decrease to $`U=4`$ shown in Fig. 6, we find even more interesting behavior. Now, when we are near half filling, we find two more configurations, a nonneutral phase (configuration 59) and a two-dimensional phase (configuration 109) joining with the checkerboard phase and configuration 108 in the initial mixtures. Each of these phases looks like a “square-lattice screen” with differing size “holes” in the screen. These two-dimensional structures are not stripe-like and it would be interesting to see if they could appear in the Hubbard model. As we dope further away, we enter the axial stripe region, now dominated by configuration 20 first, then there is a cascade to configuration 33, then a cascade to the segregated phase. This value of $`U`$ is a truly intermediate value, where many different mechanisms for ordering are present and the system can change very rapidly in response to a modification in the density.

As $`U=2`$ (Fig. 5), we see more modifications in the phase diagram. Now we see other diagonal stripe phases mixing with the checkerboard phase near half filling. This region would correspond to the Scalapino-White regime, where the stripe formation is driven more by interplays between the kinetic and potential energies and nesting effects (driven by charge fluctuations in the FK model and spin fluctuations in the Hubbard model). In addition, a much larger number of the 2D phases enter also close to half filling, illustrating the prevalence of these “screen-like” phases as well. The axial stripes also enter as we dope further away from half filling, but the configurations 20 and 33 are not nearly as stable as they are for slightly larger $`U`$. Here, we see the four-molecule phases being stabilized just before the system phase separates into the segregated phase.

Finally, for $`U=1`$, shown in Fig. 4, the predominance of the diagonal stripes, near half filling increases now supplemented by the axial checkerboard stripes, but then there is a plethora of different 2D phases that also enter as the system is doped somewhat farther from half filling, then we see a similar evolution, first to AS and then to 4M phases before the segregated phase. Here there is a tremendous complexity to the phase diagram, with many different mixtures being present due to the competition between kinetic energy and potential energy minimization brought about by the many-body aspects of the problem.

The general picture, illustrated schematically in Fig. 9, now emerges: near half filling, we often find diagonal stripes and screen-like two-dimensional phases, then a rapid transition to the segregated phase for large $`U`$. As $`U`$ is reduced, we can dope farther away from half filling before segregating, which allows many other phases to enter. In particular, there is a large region of stability for axial stripes, and as $`U`$ is reduced further, we see the emergence of axial checkerboard stripes close to half filling, near the diagonal stripes, and four-molecule phases appearing near the segregation boundary.

## IV Conclusions

In this manuscript we have numerically studied how the FK model makes the transition from the checkerboard phase at half filling to the segregated phase as the density is lowered. Since these two phases are very different from one another, there are many different pathways that one might imagine the system to take in making this crossover. Indeed we find that the pathway varies dramatically as a function of $`U`$. For large $`U`$, we have a relatively simple transition between diagonal stripe-like phases which become more two-dimensional as the localized electron density increases, until the system gives way to the segregated phase. As $`U`$ is lowered, we first see two-dimensional-“screen”-like phases enter, then we see axial stripes emerge, followed by four-molecule phases and axial checkerboard stripes. The complexity of the phase diagram greatly increases as the interaction strength decreases.

It is interesting to ask how we might expect these results to change if we allowed more configurations into our restricted phase diagram. We don’t know this answer in particular, but we do know, that of the 23,755 candidate phases only a small fraction (111 or 0.5%) appear in the phase diagram, so we don’t expect too many additional phases to appear.

Another interesting question to ask is how do these results for the FK model shed light on the stripe-formation problem in the Hubbard model. By continuity, we expect these results not to change too dramatically as we turn on a small hopping for the localized electrons (although now we must summarize our results in terms of correlation functions for the two kinds of electrons, since both are now mobile). But we also know for many fillings, there will be a “phase transition” as a function of the hopping, since the ground state of the Hubbard model is not ferromagnetic for all fillings and large $`U`$ (which is what the segregated phase maps to in the Hubbard model). The results are likely to be closer to what happens in the Hubbard model close to half filling, because the analogue of the antiferromagnetic phase is the checkerboard phase, and that is present for all $`U`$ in the Hubbard model at $`T=0`$. In general, we also feel that the FK model phase diagram must be more complicated than the Hubbard model phase diagram because of the mobility of both electrons in the latter. We feel one of the most important results of this work is that there may be two-dimensional phases that are not stripe like that form ground-state configurations for some values of the filling in the Hubbard model, and such configurations will be worthwhile to investigate with the numerical techniques that currently exist.

In conclusion, we are delighted to be able to shed some light on the interesting question for the FK model of how one makes a transition from the checkerboard phase at half filling to the segregated phase away from half filling. Since Elliott Lieb has had an important impact in proving the stabilization of these two phases, we find it fitting to ask the questions about how the two phases inter-relate. Perhaps these numerical calculations can further inspire new rigorous work that helps to identify the pathway between these two phases.

## V Acknowledgments

R.L. and G.B. acknowledge support from the Polish State Committee for Scientific Research (KBN) under Grant No. 2P03B 131 19 and J.K.F. acknowledges support from the NSF under grant No. DMR-0210717. We also acknowledge support from Georgetown University for a travel grant in the fall of 2001.

Tables

Figure Captions

Figure 1.

Figure 2.

Figure 3.

Figure 4.

Figure 5.

Figure 6.

Figure 7.

Figure 8.

Figure 9.


text:
# Deformed Harry Dym and Hunter-Zheng Equations

## 1 Introduction:

The search for exactly solvable equations has acquired enormous importance since the nonlinear Kortweg-de Vries equation was shown to be integrable . Exactly solvable equations, linear or nonlinear, constitute a very special class of dynamical systems with many interesting properties. The significance of determining new integrable systems, as well as studying their properties and solutions, can not be overestimated. The structure of integrable systems (or partial differential equations) is highly restrictive and, in general does not allow for deformations. In particular, the presence of arbitrary constant parameters (which can not be transformed away by some symmetry) is quite rare. One notable exception is the two boson equation , which has an arbitrary constant parameter present, is known to be integrable. Different values of this parameter reduce the model to other known soluble models. In this paper, we propose a new system of equations, with two arbitrary constant parameters, that is exactly soluble and reduces to various other known physical models in different limits of these parameters.

Let us recall that the Harry Dym (HD) equation

$$u_t=(u_{xx}^{1/2})_x,$$
(1)
and the Hunter-Zheng (HZ) equation

$$\left(u_t+uu_x\right)_x=\frac{1}{2}u_x^2,$$
(2)
are known to be integrable and, in fact, belong to the same hierarchy corresponding to negative and positive order flows, respectively . The HD equation has interest in the study of the Saffman-Taylor problem which describes the motion of a two-dimensional interface between a viscous and a non-viscous fluid . The HZ equation, on the other hand, arises in the study of massive nematic liquid crystals and in the study of shallow water waves . In this paper we will show that the following system of two equations,

$$u_t=\sqrt{2}\left(\frac{1\lambda u_{xx}}{\sqrt{2u_{xx}\alpha \lambda u_{xx}^2}}\right)_x$$
(3)
and

$$u_{xxt}=\alpha u_x2u_xu_{xx}uu_{xxx}+\frac{\lambda }{6}\left(u_x^3\right)_{xx},$$
(4)
are integrable for arbitrary values of the constant parameters $`\alpha ,\lambda `$. Furthermore, they belong to the same hierarchy corresponding to the negative and positive flows respectively, much like the HD and the HZ equations. It is interesting to note from (4) that when $`\alpha 0`$, it can be scaled to unity through the scaling $`xx/\alpha `$, $`uu/\alpha `$. Therefore, looking at (4) alone, it would appear that there are only two meaningful values for $`\alpha `$, namely, $`\alpha =0,1`$. Similarly, under a scaling $`tt/\alpha ^{3/2}`$, $`xx/\alpha `$, $`uu/\alpha `$, the parameter $`\alpha `$ in the dHD equation (3) can be scaled to unity when it is not zero. However, there is no scaling which will scale $`\alpha `$ to unity (when it is nonzero) simultaneously in both the equations and, consequently, if the two systems belong to the same hierarchy, $`\alpha `$ has to be thought of as an arbitrary parameter.

We note that, for $`\alpha =\lambda =0`$, eqs. (3) and (4) reduce respectively to the HD equation (1) and the HZ equation (2). For lack of a better name, we will refer to (3) as the deformed Harry Dym equation (dHD) and (4) as the deformed Hunter-Zheng equation (dHZ). (The name “generalized Harry Dym” equation has already been used earlier in the context of a multi-component Harry Dym equation .) For $`\alpha 0`$ and $`\lambda =0`$ equation (4) was studied by Alber et al  and its solutions contain solitons of the type known as umbilic solitons. Equation (4) with $`\alpha =1`$ and $`\lambda 0`$ has appeared more recently in the literature  as describing short capillary-gravity waves and its solutions are known to become multi-valued in a finite amount of time. While a matrix Lax pair for this system was provided in  and solubility of the model was argued based on a map to the sine-Gordon equation, interesting properties, such as the bi-Hamiltonian structure, the infinite number of conserved charges were not discussed at all. Our results clarify these aspects and provide direct support for the integrability of this new system. The deformed equation (3) for arbitrary $`\alpha `$ and $`\lambda 0`$ is truly a new integrable system (which to the best of our knowledge has not been studied in the literature) and leads to the above mentioned models in different limits. This can, therefore, be considered as a very rich system, much like the two boson hierarchy.

Our paper is organized as follows. In section 2, we obtain the first Hamiltonian structure of the dHD and dHZ equations. In section 3, a second Hamiltonian structure, compatible with the first one, is obtained. We show that both the dHD and dHZ equations are bi-Hamiltonian and, therefore, integrable and belong to the same hierarchy corresponding to negative and positive flows. In section 4, we present a Lax pair description of both the system of equations and construct the conserved charges of negative order from the Lax operator. We obtain a non-standard Lax description for the dHD equation for the special values of the deformation parameters $`\lambda =0,\lambda =3/\alpha `$. We argue that, for general values of the parameters, a non-standard Lax description will involve a pseudo-differential operator of infinite order. In section 5 we present our conclusions.

## 2 dHZ and dHD as Hamiltonian Systems:

To describe the dHD and the dHZ equations in a compact manner, let us introduce the following notation. Let us define

$$F_{(\alpha ,\lambda )}^22u_{xx}\alpha \lambda u_{xx}^2,A\frac{1\lambda u_{xx}}{F}.$$
(5)
Then, it follows from the definitions in (5) that

$`\kappa `$ $``$ $`(1\alpha \lambda )=\lambda F^2+(1\lambda u_{xx})^2=(\lambda +A^2)F^2,`$
$`{\displaystyle \frac{F_x}{F^3}}`$ $`=`$ $`{\displaystyle \frac{1}{\kappa }}AA_x,`$
$`{\displaystyle \frac{u_{xxx}}{F^3}}`$ $`=`$ $`{\displaystyle \frac{1}{\kappa }}A_x.`$ (6)
These and other relations following from these prove very useful in the analysis of the two systems.

Given the dHZ equation (4), we obtain from the definition in (5) that

$$F_t=\left[F\left(u\frac{\lambda }{2}u_x^2\right)\right]_x.$$
(7)
Similarly, in terms of the new variables in (5), the dHD equation (3) can be written as

$$u_t=\sqrt{2}A_x,$$
(8)
and it follows that, under the evolution of dHD,

$`F_t`$ $`=`$ $`\sqrt{2}AA_{xxx},`$
$`A_t`$ $`=`$ $`\sqrt{2}\kappa {\displaystyle \frac{A_{xxx}}{F^3}}.`$ (9)
It is now straightforward to note from eqs. (7) and (9) that

$$H_1=\sqrt{2}𝑑xF$$
(10)
is conserved under both the dHD and dHZ flows.

The dHZ equation can be obtained from a variational principle, $`\delta 𝑑t𝑑x`$, with the Lagrangian density

$$=\frac{1}{2}u_xu_t+\frac{\alpha }{2}u^2+\frac{1}{2}uu_x^2\frac{\lambda }{24}u_x^4.$$
(11)
This is a first order Lagrangian density and, consequently, the Hamiltonian structure can be readily read out, or we can use, for example, Dirac’s theory of constraints  to obtain the Hamiltonian and the Hamiltonian operator associated with (11). The Lagrangian is degenerate and the primary constraint is obtained to be

$$\mathrm{\Phi }=\pi \frac{1}{2}u_x,$$
(12)
where $`\pi =/u_t`$ is the canonical momentum. The total Hamiltonian can be written as

$`H_T`$ $`=`$ $`{\displaystyle 𝑑x\left(\pi u_t+\beta \mathrm{\Phi }\right)}`$ (13)
$`=`$ $`{\displaystyle 𝑑x\left[\frac{\alpha }{2}u^2\frac{1}{2}uu_x^2+\frac{\lambda }{24}u_x^4+\beta \left(\pi \frac{1}{2}u_x\right)\right]},`$
where $`\beta `$ is a Lagrange multiplier field. Using the canonical Poisson bracket relation

$$\{u(x),\pi (y)\}=\delta (xy),$$
(14)
with all others vanishing, it follows that the requirement of the primary constraint to be stationary under time evolution,

$$\{\mathrm{\Phi }(x),H_T\}=0,$$

determines the Lagrange multiplier field $`\beta `$ in (13) and the system has no further constraints.

Using the canonical Poisson bracket relations (14), we can now calculate

$$K(x,y)\{\mathrm{\Phi }(x),\mathrm{\Phi }(y)\}=\frac{1}{2}_y\delta (xy)\frac{1}{2}_x\delta (xy).$$
(15)
This shows that the constraint (12) is second class and that the Dirac bracket between the basic variables has the form

$$\{u(x),u(y)\}_D=\{u(x),u(y)\}𝑑z𝑑z^{}\{u(x),\mathrm{\Phi }(z)\}J(z,z^{})\{\mathrm{\Phi }(z^{}),u(y)\}=J(x,y),$$

where $`J`$ is the inverse of the Poisson bracket of the constraint (15),

$$𝑑zK(x,z)J(z,y)=\delta (xy).$$

This last relation determines

$$_xJ(x,y)=\delta (xy),$$

or

$$J(x,y)=𝒟_1\delta (xy),$$

where

$$𝒟_1=^1,$$
(16)
and can be thought of as the alternating step function in the coordinate space. We can now set the constraint (12) strongly to zero in (13) to obtain

$$H_2H_T=𝑑x\left(\frac{\alpha }{2}u^2+\frac{1}{2}uu_x^2\frac{\lambda }{24}u_x^4\right).$$
(17)
Therefore, the dHZ equation can be written in the Hamiltonian form

$$u_t=𝒟_1\frac{\delta H_2}{\delta u},$$

with $`𝒟_1`$ and $`H_2`$ given by (16) and (17), respectively.

From the results in  we know that the HD and HZ equations belong to the same hierarchy of equations. Here, too, we will see that both the dHD and the dHZ equations belong to the same hierarchy. In particular, we note that

$$u_t=𝒟_1\frac{\delta H_1}{\delta u},$$

with $`H_1`$, given by (10), yields the deformed Harry Dym equation (8). As a result, the dHD equation also is Hamiltonian with the same Hamiltonian structure of the dHZ equation in (16) and, consequently, has a Lagrangian description given by

$$=\frac{1}{2}u_xu_t\sqrt{2}F.$$

The reader can easily check that $`H_2`$ is also conserved by both the dHD and dHZ equations, much like $`H_1`$. Note that for $`\alpha =\lambda =0`$, these two charges reduce to the corresponding ones in the HD and HZ systems.

## 3 dHD and dHZ as bi-Hamiltonian Systems:

It is well known that a system can be shown to be integrable if it is bi-Hamiltonian . This corresponds to the system having a Hamiltonian description with two distinct Hamiltonian structures that are compatible. Therefore, we try to see if the dHD and the dHZ equations can be described as bi-Hamiltonian systems. For this, we have to find a second Hamiltonian description for the two systems.

In order to determine a second Hamiltonian structure for the two systems, we recall  that, when $`\alpha =\lambda =0`$, the corresponding $`H_1`$ is a Casimir of the second Hamiltonian structure for the HD and HZ systems, namely,

$$𝒟_2^{\left(\alpha =\lambda =0\right)}=^2u_{xx}^1+^1u_{xx}^2$$
(18)
satisfies

$$𝒟_2^{\left(\alpha =\lambda =0\right)}\frac{\delta H_1^{\left(\alpha =\lambda =0\right)}}{\delta u}=0.$$

Since we know $`H_1`$ for the deformed systems which is a generalization of $`H_1^{\left(\alpha =\lambda =0\right)}`$, we look for a Hamiltonian structure for which it is a Casimir. With some work, it can be determined that $`H_1`$ given in (10) is a Casimir of

$$𝒟_2𝒟_2^{(\alpha ,\lambda )}=\frac{1}{2}\left(^2F^2^1+^1F^2^2\right)+\lambda ^2u_{xxx}^1u_{xxx}^2.$$
(19)
Note that this structure reduces to (18) when $`\alpha =\lambda =0`$. The skew symmetry of this Hamiltonian structure is manifest. The proof of the Jacobi identity for this structure as well as its compatibility with $`𝒟_1`$ in (16) can be determined through the standard method of prolongation described in ref. , which we discuss briefly.

Performing the change of variables

$$w=u_{xx},$$

the Hamiltonian structures (16) and (19) assume the forms

$`𝒟_1`$ $`=`$ $`^3,`$
$`𝒟_2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(F^2+F^2\right)+\lambda w_x^1w_x.`$
We can construct the two bivectors associated with the two structures as

$`\mathrm{\Theta }_{𝒟_1}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle 𝑑x\left\{\theta 𝒟_1\theta \right\}}={\displaystyle \frac{1}{2}}{\displaystyle 𝑑x\theta }\theta _{xxx},`$
$`\mathrm{\Theta }_{𝒟_2}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle 𝑑x\left\{\theta 𝒟_2\theta \right\}}`$
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle 𝑑x\left\{\alpha \theta \theta _x+2w\theta \theta _x\lambda w^2\theta \theta _x+\lambda w_x\theta (^1w_x\theta )\right\}}.`$
Using the prolongation relations,

$`\text{pr}\stackrel{}{v}_{𝒟_1\theta }(w)`$ $`=`$ $`\theta _{xxx},`$
$`\text{pr}\stackrel{}{v}_{𝒟_1\theta }(w^2)`$ $`=`$ $`2w\text{pr}\stackrel{}{v}_{𝒟_1\theta }(w),`$
$`\text{pr}\stackrel{}{v}_{𝒟_1\theta }(w_x)`$ $`=`$ $`\left(\text{pr}\stackrel{}{v}_{𝒟_1\theta }(w)\right)_x,`$
$`\text{pr}\stackrel{}{v}_{𝒟_2\theta }(w)`$ $`=`$ $`\alpha \theta _x+2w\theta _x\lambda w^2\theta _x+w_x\theta \lambda ww_x\theta +\lambda w_x(^1w_x\theta ),`$
$`\text{pr}\stackrel{}{v}_{𝒟_2\theta }(w^2)`$ $`=`$ $`2w\text{pr}\stackrel{}{v}_{𝒟_2\theta }(w),`$
$`\text{pr}\stackrel{}{v}_{𝒟_2\theta }(w_x)`$ $`=`$ $`\left(\text{pr}\stackrel{}{v}_{𝒟_2\theta }(w)\right)_x,`$ (20)
it is straightforward to show that the prolongation of the bivector $`\mathrm{\Theta }_{𝒟_2}`$ vanishes,

$$\text{pr}\stackrel{}{v}_{𝒟_2\theta }\left(\mathrm{\Theta }_{𝒟_2}\right)=0,$$

implying that $`𝒟_2`$ satisfies Jacobi identity. Using (20), it also follows that

$$\text{pr}\stackrel{}{v}_{𝒟_1\theta }\left(\mathrm{\Theta }_{𝒟_2}\right)+\text{pr}\stackrel{}{v}_{𝒟_2\theta }\left(\mathrm{\Theta }_{𝒟_1}\right)=0,$$

showing that $`𝒟_1`$ and $`𝒟_2`$ are compatible. Namely, not only are $`𝒟_1,𝒟_2`$ genuine Hamiltonian structures, any arbitrary linear combination of them is as well. Any physical system that is Hamiltonian with respect to these two structures, therefore, defines a pencil system and is integrable. It is worth pointing out here that when $`\alpha =\lambda =0`$, the second Hamiltonian structure (18) represents the centerless Virasoro algebra  (with dimension zero operators). The structure in (19) appears to be a highly nonlocal generalization of this algebra, but we are not familiar with any study of such an algebra in the literature.

To show that the dHD and dHZ are bi-Hamiltonian systems, we note that the charges

$`H_1`$ $`=`$ $`{\displaystyle 𝑑xu_x^2},`$ (21)
$`H_2`$ $`=`$ $`{\displaystyle \frac{1}{2\kappa }}{\displaystyle 𝑑xFA_x^2}`$ (22)
are also conserved by both the dHZ and dHD equations. While the charge in (21) is the unmodified charge of the HD and HZ systems (this seems to be a simple coincidence), the charge in (22) is a true generalization of the corresponding charge of the HD and HZ systems. With these charges, it is easy to see that the dHZ equation can be written in a truly bi-Hamiltonian form

$$u_t=𝒟_1\frac{\delta H_2}{\delta u}=𝒟_2\frac{\delta H_1}{\delta u}.$$

Similarly, the dHD equation can also be written in the bi-Hamiltonian form

$$u_t=𝒟_1\frac{\delta H_1}{\delta u}=𝒟_2\frac{\delta H_2}{\delta u}.$$

Thus, we see that both the dHD as well as dHZ systems are bi-Hamiltonian with the same two compatible Hamiltonian structures and are, therefore, integrable.

## 4 The Lax Representation:

When a system is bi-Hamiltonian, we can naturally define a hierarchy of commuting flows through the relation

$$u_t=K_n[u]=𝒟_1\frac{\delta H_{n+1}}{\delta u}=𝒟_2\frac{\delta H_n}{\delta u},n=0,1,2,\mathrm{}.$$
(23)
In the present case, both the Hamiltonian structures have Casimirs. We have already seen that $`H_1`$ is a Casimir of $`𝒟_2`$ and it can be checked that the trivial Hamiltonian

$$H_0=𝑑xu_{xx}$$
(24)
formally defines the Casimir of $`𝒟_1`$ (namely, if we write formally $`\frac{\delta H_0}{\delta u}=^2`$, it is annihilated by $`𝒟_1`$). As a result, the system of flows can be extended to both positive and negative integer values for $`n`$. In this way, we see that much like in the HD and HZ systems , the dHD and dHZ systems also belong to the same hierarchy corresponding to the negative and positive flows.

Let us introduce the recursion operator following from the two Hamiltonian structures as

$$R=𝒟_2𝒟_1^1.$$

Then, it follows from (23) that

$$K_{n+1}=RK_n,$$

and

$$\frac{\delta H_{n+1}}{\delta u}=R^{}\frac{\delta H_n}{\delta u},$$
(25)
where

$$R^{}=^1u_{xx}^1+(\alpha +u_{xx})^2\lambda u_{xx}^1u_{xx}^1$$
(26)
is the adjoint of $`R`$. The conserved charges for the hierarchy can, of course, be determined in principle recursively from (25). However, in practice, integrating the recursion relation is highly nontrivial. Therefore, we look for a Lax representation for the system of dHD and dHZ equations which will allow us to construct the conserved charges directly.

It is well known  that for a bi-Hamiltonian system of evolution equations, $`u_t=K_n[u]`$, a natural Lax description

$$\frac{M}{t}=[M,B],$$

is easily obtained where, we can identify

$`M`$ $``$ $`R^{},`$
$`B`$ $``$ $`K_n^{}.`$
Here $`K_n^{}`$ represents the Fréchet derivative of $`K_n`$, defined by

$$K_n^{}[u]v=\frac{d}{dϵ}K_n[u+ϵv]|_{ϵ=0}..$$

For the dHD and dHZ system of equations in (3) and (4) respectively, we have

$`B^{\text{ dHD}}`$ $``$ $`K_2=\sqrt{2}\kappa ^2F^3,`$
$`B^{\text{ dHZ}}`$ $``$ $`K_1=(\alpha +u_{xx})^1+\left(u{\displaystyle \frac{\lambda }{2}}u_x^2\right).`$
The two systems have the same $`M=R^{}`$ given in (26). It can now be checked that

$`{\displaystyle \frac{M}{t}}`$ $`=`$ $`[M,B^{\text{ dHD}}],`$
$`{\displaystyle \frac{M}{t}}`$ $`=`$ $`[M,B^{\text{ dHZ}}],`$ (27)
do indeed generate the dHD and the dHZ equations and, thereby, provide a Lax pair for the system.

One of the advantages of a Lax representation is that they directly give the conserved charges of the system. From the structure of (27), it follows that $`\mathrm{Tr}M^{\frac{2n+1}{2}}`$ are conserved, where “Tr” represents Adler’s trace . We note that

$`\text{Tr}M^{\frac{2n+1}{2}}`$ $`=`$ $`0,n1`$
$`\text{Tr}M^{\frac{1}{2}}`$ $`=`$ $`{\displaystyle 𝑑xF},`$
$`\text{Tr}M^{\frac{1}{2}}`$ $`=`$ $`{\displaystyle \frac{1}{2\kappa }}{\displaystyle 𝑑xFA_x^2},`$
$`\text{Tr}M^{\frac{3}{2}}`$ $`=`$ $`{\displaystyle \frac{3}{\kappa }}{\displaystyle 𝑑x\left(4\frac{A_{xx}^2}{F}+\frac{1}{\kappa }FA_x^4\right)},`$ (28)
$`\mathrm{}`$ $`.`$
The first two nontrivial charges correspond respectively to $`H_1,H_2`$ given in eqs. (10) and (22), constructed earlier by brute force. In fact, all $`H_{n1}`$ with positive $`n0`$ can be constructed from $`\mathrm{Tr}M^{\frac{2n1}{2}}`$ and by construction (namely, because of the nature of (27)), they are conserved under both the dHD and dHZ flows. Unfortunately, as is clear from (28), this procedure does not yield the charges $`H_n`$ with positive integer values. These are, in general, non-local and even in the HD and HZ case, construction of these charges relies primarily on the recursion relation (25). It remains an interesting question to construct these charges in a more direct manner.

The Harry Dym equation has a Gelfand-Dikii representation for the Lax pair, while the HZ equation does not. We will now discuss the existence of such a Lax representation for the dHD equation. A spectral problem associated with the dHD equations can be obtained from the recursion relation (25) (see  and references therein) as follows. Introducing a spectral parameter $`\mu `$ and defining

$$\psi ^2(x,t,\mu )=\underset{n=0}{\overset{\mathrm{}}{}}\mu ^n\frac{\delta H_n}{\delta u},$$

we note that the recursion relation (25) can be written compactly as (recall that $`H_0`$ is a Casimir of $`𝒟_1`$)

$$(𝒟_1\mu 𝒟_2)\psi ^2=0,$$

or

$$(1\mu R^{})\psi ^2=0,$$
(29)
which defines an eigenvalue problem for the eigenfunction $`\psi ^2`$ with eigenvalue $`1/\mu `$. A linear eigenvalue problem can be derived from this if we can factorize the operator $`(1\mu R^{})`$.

Let us note that the operator $`R^{}`$ in (26) can be rewritten in the form

$$R^{}=\frac{1}{2}\left[^1\left(F_{(\alpha ,\lambda )}^2+X\right)^1+\left(F_{(\alpha ,\lambda )}^2+X\right)^2\right],$$
(30)
where

$$X=\underset{n=1}{}X_n^n,$$

and the coefficients $`X_n`$ can be determined recursively to be

$`X_1`$ $`=`$ $`0,`$
$`X_2`$ $`=`$ $`{\displaystyle \frac{\lambda }{4}}u_{xxx}^2,`$
$`X_3`$ $`=`$ $`{\displaystyle \frac{1}{2}}X_{2,x},`$
$`\mathrm{}`$ $`.`$
When $`\lambda =0`$, it follows that $`X=0`$ and that $`F_{\left(\alpha ,\lambda =0\right)}^2`$ is a simple function. In this case, the eigenvalue problem (29) can be factorized as

$$(1R^{})\psi ^2=2^1\varphi ^2\varphi ^3\left(^2\frac{\mu }{4}F_{\left(\alpha ,\lambda =0\right)}^2\right)\varphi =0,$$

where we have identified

$$\varphi ^2=(^2\psi ^2).$$

This shows that if the linear equation

$$\left(^2\frac{\mu }{4}F_{\left(\alpha ,\lambda =0\right)}^2\right)\varphi =0$$

is satisfied, then (29) will hold and this identifies the Lax operator for the system to be

$$L=\frac{1}{F_{\left(\alpha ,\lambda =0\right)}^2}^2.$$

In fact, it can be readily checked that when $`\lambda =0`$, the hierarchy of dHD equations can be obtained from the non-standard Lax equation

$$\frac{L}{t_n}=4\sqrt{2}[L,(L^{(2n1)/2})_2].$$

The conserved quantities for this system can be obtained from $`\mathrm{Tr}L^{(2n1)/2}`$, $`n=0,1,2,\mathrm{}`$.

On the other hand, when $`\lambda 0`$, the coefficients $`X_n`$ are nontrivial and $`X`$ represents a pseudo-differential operator. The factorization, in such a case, is not so simple and, in principle would involve an infinite series of terms. For arbitrary values of $`\kappa `$, the terms in the series can possibly be determined recursively. However, this is not very interesting. We simply note here that for the special value $`\kappa =4`$, the infinite series of terms seems to have a simpler compact form and the Lax operator, in such a case, has the form

$$L_{\kappa =4}=\frac{1}{F}\frac{1}{F}\frac{1}{4}A_x^1A_x,$$

and, in this case, the hierarchy of dHD equations can be obtained from the non-standard Lax representation

$$\frac{L_{\kappa =4}}{t_n}=4\sqrt{2}[L_{\kappa =4},\left(L_{\kappa =4}^{(2n1)/2}\right)_2].$$

The conserved quantities, in this case, also follow from $`\mathrm{Tr}L_{\kappa =4}^{(2n1)/2}`$ and up to multiplicative factors, they have the forms given in (28) with $`\kappa =4`$. A simple Lax description for arbitrary $`\lambda `$, however, remains an open question.

## 5 Conclusion:

In this paper, we have studied the general system of dHD and dHZ equations. We have shown that both these systems are bi-Hamiltonian and, therefore, integrable and belong to the same hierarchy corresponding to negative and positive flows. The Lax pair for the two system of equations have been derived and conserved charges corresponding to negative integer values follow from the Lax operator. A simple construction of the charges for positive integer values remains an open question. For $`\lambda =0`$, we have constructed a non-standard Lax representation for the dHD equation, which involves a purely differential Lax operator. For arbitrary values of $`\lambda `$, we have argued that a non-standard Lax representation will necessarily involve a Lax operator which is a pseudo-differential operator of infinite order. For the particular case of $`\kappa =4`$, however, this takes a simpler compact form.

## Acknowledgments

One of us (AD) would like to thank the members of the physics Departments at UFSC (Brazil) and Wrocław University (Poland) for hospitality, where parts of this work was carried out. This work was supported in part by US DOE grant no. DE-FG-02-91ER40685 as well as by NSF-INT-0089589.


text:
# The effective mass of two–dimensional 3He

## ACKNOWLEDGMENTS

This work was supported, in part, by the Austrian Science Fund under project P12832-TPH and from DGI (Spain) Grant No. BFM2002-00466 and Generalitat de Catalunya Grant No. 2001SGR-00222.