text:
# Electromechanical Limits of Polymersomes

## Abstract

Self-assembled membranes of amphiphilic diblock copolymers enable comparisons of cohesiveness with lipid membranes over the range of hydrophobic thicknesses $`d=315`$ $`nm`$. At zero mechanical tension the breakdown potential $`V_c`$ for polymersomes with $`d=15`$ $`nm`$ is $`9`$ $`V`$ compared to $`1`$ $`V`$ for liposomes with $`d=3`$ $`nm`$. Nonetheless, electromechanical stresses at breakdown universally exhibit a $`V_c^2`$ dependence, and membrane capacitance shows the expected strong $`d`$-dependence, conforming to simple thermodynamic models. The viscous nature of the diblock membranes is apparent in the protracted post-poration dynamics.

A primary task for any biological membrane is to separate inside from out, partitioning ions and other solutes that generate a transmembrane potential $`V_m`$. Electrically excitable cells, particularly neurons, control and propagate this potential for purposes of signaling and inter-communication . Efforts to exploit and better understand such phenomena have most recently motivated the semi-patterned growth of nerve cells on circuits  as well as the generation of artificial networks with soft nanotubes pulled from vesicles . Further opportunities in such directions now seem possible with wholly synthetic block copolymers that, analogous to lipids in water, self-assemble into membranes, minimizing interfacial exposure of hydrophobic segments (Fig. 1) and thereby generating vesicles termed polymersomes .

Physical limits of self-assembled lipid membranes impose important constraints on electrochemical signals. Indeed, the operating range of biomembrane excitations, such as the action potential of a neuron, is generally less than $`0.1`$ $`V`$ over time scales of milliseconds. At the same time, a finite resting tension $`\tau `$ is exerted on most cell membranes, including neurons , through substrate adhesion and additional cellular mechanisms. From optical trap techniques and other methods, $`\tau `$ is now generally believed to be in the range of $`0.01`$ to $`0.1`$ $`pN/nm`$ , though additional transient stresses on cells can readily magnify such tensions . For example, electrocompressive stresses arise with a voltage drop across a membrane and have long been recognized as coupled, most simply, to lateral tensions through an integrated form of the Lippmann equation :

$$\tau +\frac{1}{2}C_mV_m^2=\tau _{net}$$
(1)
where $`C_m`$ is the capacitance per unit area of the membrane, and $`\tau _{net}`$ is an equivalent tension. Alternative descriptions of electromechanical coupling have been put forward and include a thermodynamic study of block copolymer lamellae perforation that was specialized in its comparisons to lipid membranes . Direct tests of lipid membrane cohesion by Needham and Hochmuth  have clearly demonstrated the electromechanical limits of lipid membranes via the combination of micropipette aspiration and electroporation techniques. From these and related tests it has become clear that transmembrane potentials higher than $`1`$ $`V`$ and mechanical tensions higher than $`10`$ $`pN/nm`$ are not sustainable in lipid membranes . Clearly, a hydrophobic thickness of only $`d=34`$ $`nm`$ for such membranes is a determinant of electromechanical stability.

Formation of polymersomes in dilute aqueous solutions has been described elsewhere . Generally, the phase behavior of amphiphiles composed of strongly segregating segments is determined by the hydrophilic fraction $`f`$ and mean molecular weight $`\overline{M}_n`$ . Lipids provide initial guidance in the synthesis  of biomimetic super-amphiphiles; Table I shows $`f2939\%`$ for vesicle-forming phospholipids and diblock copolymers, although the latter are much higher in $`\overline{M}_n`$ and somewhat more polydisperse. For the copolymers used here, the hydrophilic segment of polyethylene oxide (PEO) forms a highly hydrated brush  while aggregation is driven by a hydrophobic block of either polyethylethylene (PEE) or its unsaturated homologue polybutadiene (PBD). As will be reported elsewhere , scaling of $`d`$ with $`\overline{M}_n`$ (as well as elastic properties) is consistent with the hydrophobic core behaving as a fluid-like melt.

When small pieces of our copolymers are added to water, vesicular tubes spontaneously sprout and grow (Fig. 1(b)). Vesicle formation processes - kinetics, yields, and morphologies - exhibit temperature dependences that likely stem from hydration-facilitated melting of PEO . The vesicles formed exhibit a wide variety of shapes that range from multi-armed stars to hundreds of micron-long axon-like tubes (Fig. 1) as well as flaccid spheroids. Transformations between morphologies are easy to achieve through osmotic adjustments of the external solution  and thus reveal the semi-permeable nature of the copolymer membranes. Our first generation polymersomes composed of a PEO-PEE diblock copolymer (designated OE7 in Table I) possess a core thickness $`d8\pm 1`$ $`nm`$ and have already been shown to be mechanically tough . In this Letter, we elaborate the more general electromechanical responses of polymersomes including two new covalently-crosslinkable PEO-PBD copolymers of similar block fraction but higher molecular weight (OB16 and OB18 in Table I). Specific comparisons are made with our own measurements of the highly representative lipid stearoyl-oleoyl-phosphatidylcholine (SOPC), which are comparable to the results reported elsewhere . Thus, for the first time we can thoroughly study the cohesiveness of membranes as a function of their hydrophobic core thickness over the range $`315`$ $`nm`$.

Flaccid polymersomes were progressively aspirated into a micropipette (Fig. 2 inset) at stress rates of $`0.1`$ $`pNnm^1sec^1`$ up to the point of rupture. Values of the lysis tension $`\tau _c`$ for SOPC as well as the three polymersomes studied are shown in Table II. The results generally show that increasing $`\overline{M}_n`$ leads to an increase in stability, consistent with general ideas of meso-phase stability for strongly segregated copolymers . As such, the hydrophobic core is expected to behave as a dense fluid. Assuming incompressibility, the thickness at rupture is given by $`d_c=d/(1+\alpha _c)`$, where $`\alpha _c`$ is the area strain at rupture as measured directly from images of aspirated vesicles. For lipids, $`\alpha _c`$ is very small, typically $`0.030.05`$ , so that $`d_cd`$. For polymersomes, however, $`\alpha _c`$ ranges from $`0.2\pm 0.06`$ for OE7 to $`0.44\pm 0.08`$ for OB18. Thus, polymersomes should thin considerably. A simple linear fit of $`\tau _c`$ vs. $`d_c`$ is found to intersect the $`y`$-axis at essentially zero tension (Fig. 2). Moreover, the slope of this fit defines a lysis stress $`\mathrm{\Sigma }_c=\tau _c/d_c27`$ $`atm`$. Such a stress is typical of cavitation pressures for homogeneous fluids  suggesting that lysis of membranes can occur through nucleation of water-filled cavities inside the hydrophobic core.

Electromechanical experiments were performed to determine the breakdown potential $`V_c`$ at different membrane tensions. Once again a flaccid vesicle was aspirated to a prescribed, subcritical tension while being manipulated to a position between two platinum wire electrodes separated by $`1`$ $`mm`$. A $`60`$ $`\mu sec`$ square-shaped pulse was applied across the electrodes at intervals that varied from seconds to minutes depending on the post-poration dynamics of the vesicle (see below). The applied potential was increased at discrete intervals until membrane rupture was observed. A typical electroporation event at low applied tension is shown in Fig. 3. The electrodes were arranged parallel to the pipette, generating an electric field $`𝐄`$ as illustrated, and the suction pressure was sufficient to hold the vesicle in the micropipette. A first image taken at zero applied field ($`0`$ $`msec`$) demonstrates the integrity of the membrane by showing a phase-dense vesicle interior (sucrose solution) against a phase-light exterior (equi-osmolar glucose). Dramatic rupture within 1-2 video frames of the applied pulse was invariably indicated by a ’jet’ of released sucrose; such jets always occurred at focal point(s) on the membrane where the surface normal is parallel/antiparallel to $`𝐄`$. However, the dynamics of pore formation differed considerably from one membrane system to another. Membranes composed of OB16 (Fig. 3) typically showed two very large, antipodal pores that grew to many microns in size over seconds. Continued growth often completely rendered the OB16 polymersomes, despite the low tensions. Phospholipid vesicles, in contrast, showed more rapid and even reversible pore dynamics at vanishing tension, as described by others . Polymersomes of the smallest copolymer, OE7, exhibited liposome-like dynamics whereas polymersomes composed of the largest copolymer, OB18, were far more protracted in their poration dynamics. With OB18, encapsulated sucrose was invariably lost over periods of tens of seconds, with weak, barely visible jets implying relatively small but sustained pores.

Kinetic diversity in membrane poration is understood at a first level in terms of the interplay between edge energy, characterized by a line tension $`\lambda `$ that tends to close the pore, and the work done by any dynamic membrane tensions that tend to force the pore open. A balance between these two energies gives the metastable pore size $`r_0=\lambda /\tau `$ above which the pore will grow - provided sufficient tension is maintained. Sandre et al.  recognized the important role of viscosities in vesicle pore dynamics. In the limit of vanishing tension, the rate of pore closure was postulated to scale as $`\lambda /\eta _s`$ where $`\eta _s`$ is the surface viscosity. Although $`\lambda `$ is expected - from simple bending energy considerations of a hydrophilic pore - to increase in linear proportion to $`d`$, separate measurements of lateral diffusion of copolymers (Lee et al. (unpublished)) suggest that $`\eta _s`$ increases strongly with $`\overline{M}_n`$ and is at least ten-fold greater in the present polymersomes compared to lipid membranes. Moreover, the $`\overline{M}_n`$ of both OB16 and OB18 exceed the relevant entanglement molecular weights ($`2`$ to $`3`$ $`kDa`$), so that the polymeric nature of the present systems provides a clear basis for the slowed pore dynamics.

The transmembrane potential $`V_m`$ has been calculated previously for a spherical shell  as $`V_m=1.5R_vE\mathrm{cos}\theta [1\mathrm{exp}(t/t_c)]`$ where $`R_v`$ is the radius of the vesicle and $`\theta `$ is the angle between the membrane surface normal and $`𝐄`$. The charging time $`t_c`$ is given by $`t_c=R_vC_m(\rho _i+0.5\rho _e)`$, where $`\rho _i`$ and $`\rho _e`$ are internal and external solution resistivities, respectively. In this work, $`t_c`$ is orders of magnitude smaller than the pulse duration used. To assess coupling of the mechanical and electrical stresses as suggested by Eq. (1), polymersomes are aspirated to a set tension and subsequently subjected to electrical stresses as explained above. At each imposed tension, $`V_c`$ is determined and the resulting points are plotted in the $`\tau V_m`$ plane (Fig. 4). The data show that the higher the membrane tension, the lower the breakdown potential. For OB18, the results are remarkable: this $`15`$ $`nm`$ thick membrane can transiently withstand the potential of a $`9`$ $`V`$ battery.

Fitting the rupture data ($`\tau ,V_c`$) for each self-assembled membrane system by Eq. (1) readily generates an amphiphile-specific estimate of $`C_m`$ at rupture (Table II). Within this phase space spanned by $`\tau `$ and $`V_m`$, the area under each parabolic segment provides a phenomenological measure of the range of electromechanical function or robustness. For SOPC, this robustness is very small (see Table II); for polymersomes, in contrast, the robustness can be orders of magnitude larger. Fitted values of polymersome $`C_m`$ are nonetheless very low compared to those reported for lipid membranes, $`1`$ $`\mu F/cm^2`$ . This is qualitatively consistent with the enhanced thickness of polymersome membranes since $`C_mϵ/d_c`$, where $`ϵ`$ is the hydrophobic dielectric constant. On the other hand, $`V_c`$ at $`\tau =0`$ increases considerably with $`d_c`$ (see Eq. (1) and Fig. 4). By definition, the surface charge at rupture $`\sigma _c=C_mV_c`$ should then be independent of $`d_c`$. Indeed, Table II indicates that $`\sigma _c`$ varies by only a factor of about two within this structurally diverse set of membrane systems studied.

Polymersomes made from OB16 and OB18 offer the further possibility of crosslinking the membrane core. Free radical polymerization in solution proves highly effective  with dramatic increases in the robustness of the membranes. Simple osmotic rupture tests - where large vesicles are observed to burst after suspension in a sufficiently hypotonic solution - indicate that $`\tau _c`$ of cross-linked membranes are of order $`10^3`$ $`pN/nm`$. Thus $`\mathrm{\Sigma }_c1000`$ $`atm`$. Using such values in Eq. (1) and assuming values for $`C_m`$ listed in Table II, $`V_c`$ for the cross-linked OB’s are estimated to exceed $`125`$ $`V`$ for $`60`$ $`\mu sec`$ pulses.

Electrically excitable cells, particularly neurons, control and propagate electrochemical potentials for signaling purposes , but the transmembrane voltages employed are always constrained by physical limits of a self-assembled lipid membrane system. The results here with a series of self-assembling biomimetic diblock copolymers considerably expand these limits, and clearly indicate the overall electromechanical robustness of hyperthick membranes. As might be expected, $`\tau _c`$ and $`V_c`$ both increase with $`d`$ while $`C_m`$ decreases and $`\sigma _c`$ hardly changes. Though a microscopic theoretical treatment of membrane electroporation is lacking, the capacitance and insulation properties of such self-assemblies provide a more physically reliable basis for microelectronics applications such as sensor platforms and those that might exploit conducting copolymer segments. In addition, drug delivery applications could benefit through formulations of robust copolymer vesicles that reversibly reseal or not when challenged by sufficient membrane tensions and/or voltages.

ACKNOWLEDGMENTS HAE and HB thank R. Parthasarathy for discussions and comments. HAE was supported by NIH/NHLBI T32HL-07954 training grant. DED gratefully acknowledges discussions with D.A. Hammer and E.A. Evans. Funding was provided by NSF-MRSEC at the U. of Minnesota and U. of Pennsylvania, and NASA.


format: text

text:
# Comment on cond-mat/0107371: “Dynamical exponents of an even-parity-conserving contact process with diffusion”

## Acknowledgments

This work has initiated during a visit to UFMG. I thank R. Dickman introducing me to the field of nonequilibrium phase transitions. This work was partially supported by CNPq and FAPERJ.

marcio@if.uff
dickman@fisica.ufmg.br

FIGURE CAPTIONS

FIG. 1. Density of active sites $`\rho (t)`$ as a function of time (in Monte-Carlo steps). Results are shown for the diffusionless CP(2) (upper) and CP(2)d with d=0.05 (lower).

FIG. 2. Cummulative number of active sites, $`M(t)`$ as a function of time (in Monte Carlo steps) for the CP(2) (upper) and CP(2)d (lower) with d=0.05.

FIGURES


text:
# Self-organized branching process for a one-dimensional ricepile model

## 1 Introduction

Since the pioneering work of Bak, Tang and Wiesenfeld (BTW) ba\_ta\_wi\_87 ; ba\_ta\_wi\_88 , the sandpile model became one of prototype abstract models exhibiting self-organized criticality (SOC). The original BTW model as well as its variants (see e.g. ka\_na\_wu\_zho\_89 ; gra\_ma\_90 ; manna\_91 ; manna\_91a ; che\_sta\_ve\_za\_98 ) consists of a cellular automaton slowly driven by stochastic perturbation. The state of each site is described by number of grains on top of it. (Actually, this number is rather the slope than the height, if we like to interpret the model as a real sandpile. However, in 1D models, investigated here, the description through slope and height variables are strictly equivalent.) If the number of grains exceeds a threshold, the site becomes active, a toppling occurs and grains are transferred to neighbouring sites, which then may become active and the process continues. The driving consists in adding grains on randomly chosen sites. The critical state is reached asymptotically in the limit of infinitely slow driving so\_jo\_do\_95 . Fully deterministic versions were also studied, showing periodic wi\_the\_na\_90 ; ma\_ma\_92  or self-similar but non-random behaviour he\_cha\_de\_ro\_wa\_99 .

Even though experiments on real sandpiles did not confirm SOC behaviour, due to inertia effects ja\_li\_na\_89 ; rajchenbach\_90 ; he\_so\_ke\_ha\_ho\_gri\_90 ; nagel\_92 ; pra\_ola\_92 ; ba\_me\_96 ; ba\_wo\_96 , in the experiments using rice fre\_chri\_ma\_fe\_jo\_mea\_96 ; ma\_fe\_chri\_fre\_jo\_99  instead of sand it was found that large aspect ratio of the rice grains can lead to SOC behaviour fre\_chri\_ma\_fe\_jo\_mea\_96 , contrary to the case of sand, which has grains much closer to spherical.

Another difference between a typical sandpile and ricepile experiments is that the ricepiles used in the experiments are quasi one-dimensional fre\_chri\_ma\_fe\_jo\_mea\_96 ; ma\_fe\_chri\_fre\_jo\_99 . While the original BTW model in one dimension is trivial, there are several variants of 1D BTW model which exhibit non-trivial behaviour ka\_na\_wu\_zho\_89 ; he\_cha\_de\_ro\_wa\_99 ; me\_ba\_94 ; sorensen\_96 ; lu\_us\_98 ; di\_al\_mu\_pe\_ve\_za\_01 ; pri\_iva\_pov\_hu\_01 . Also the sandpiles on quasi one-dimensional stripes were investigated ma\_ta\_zha\_99 . Several one-dimensional models devised especially for modelling the ricepiles were studied ama\_la\_96 ; ama\_la\_96a ; ama\_la\_97 ; chri\_co\_fre\_fe\_jo\_96 ; pa\_bo\_96 ; sdzhang\_97 ; be\_be\_mhi\_zha\_99 ; be\_be\_ke\_lo\_mhi\_99 ; ma\_je\_la\_sne\_97 ; markosova\_00 ; sdzhang\_00 . The models taking into account a possible long-range rolling of grains are able to describe the transition from SOC behaviour typical for ricepiles to the inertia-dominated behaviour of sand heaps gle\_can\_tam\_zhe\_01 ; gleiser\_01 .

Besides numerous exact results and renormalisation-group calculations (to cite only a few items of a vast bibliography, see dhar\_89 ; dha\_ma\_90 ; markosova\_95 ; klo\_mas\_tan\_01 ; pie\_ve\_za\_94 ; ivashkevich\_96 ; zhang\_89a ), the mean-field approximation ta\_ba\_88a ; ka\_ko\_96 ; ve\_za\_97a  was very useful in clarifying the nature of the SOC state, even though it cannot give correct values of the exponents below the upper critical dimension.

It was realised soon that the mean-field approximation for sandpiles is related to critical branching processes alstrom\_88 ; garcia\_94 . This idea lead to the introduction of self-organized branching processes za\_la\_sta\_95 ; la\_za\_sta\_96 ; ca\_te\_ste\_96 ; ve\_ma\_ba\_97 ; chu\_ada\_99 ; chu\_ada\_99a , which describe the approach to the critical state. Similar approach consist in mapping the sandpile to the percolation on Bethe lattice so\_va\_an\_99 .

The approximation is based on the observation that in high dimension activity returns to the same site with very small probability. So, we can suppose that in each step the toppling occurs at a site, which never toppled before during the same avalanche. Each toppling is mapped to one branching. Statistical properties of avalanches are determined by the probability $`p`$ of branching. This probability is itself determined self-consistently. If the avalanche is sub-critical, it does not fall off the system and average number of grains, and thus $`p`$, increases. If, on the other hand, the avalanche is super-critical, it surely falls off the system, which leads to decrease of the average number of grains and decrease of $`p`$. It was shown za\_la\_sta\_95 , that this process sets the $`p`$ exactly to the critical value, where the avalanche sizes $`s`$ have power-law distribution $`P(s)s^\tau `$ with mean-field exponent $`\tau =\frac{3}{2}`$.

The purpose of this work is to modify the self-organized branching processes in order to describe one-dimensional ricepile models. Our model will be designed to comprise the one-dimensional BTW model as a special case. Clearly we cannot obtain correct values of the exponents. Our main question will be, whether there is a sharp transition from trivial 1D BTW behaviour to SOC behaviour or what is the nature of the crossover from the former to the latter.

The paper is organised as follows. In the next section we define our version of the branching process, suitable for treating the one-dimensional ricepile. We find the condition for the criticality and investigate the crossover from the trivial one-dimensional BTW behaviour to the critical branching process. The self-organization toward the critical state is investigated in the section 3. We first define the self-organized branching process, then find the fixed point of the dynamics and show that it exactly corresponds to critical branching process. We finally investigate the influence of finite size effects and find the finite-size scaling form. The section 4 concludes and summarises the work.

## 2 Branching process for one-dimensional model

### 2.1 Ricepile model

The ricepile models were already thoroughly investigated by numerical simulations. In fact, there are two variants of the one-dimensional ricepile model. The so-called “Oslo model” chri\_co\_fre\_fe\_jo\_96 ; pa\_bo\_96 ; sdzhang\_97 ; be\_be\_mhi\_zha\_99  supposes that the critical slope depends on space and time, and assumes new random value after each toppling event. Another approach ama\_la\_96 ; ama\_la\_96a ; ama\_la\_97  assumes that the toppling occurs with certain probability, which depends on actual slope. It is the second approach, which we will follow in this article. It may be also noted that a two-dimensional model which also implements stochastic topplings was studied before ou\_lu\_di\_93 .

We recall shortly the definition of the model. We consider a chain of $`L`$ sites. The state of site $`i`$, $`i=1,2,\mathrm{},L`$ is described by a slope $`z_i=h_ih_{i+1}`$ where the height $`h_i`$ is a non-negative integer, with boundary condition $`h_{L+1}=0`$. If the pile is in a stable state, a grain is dropped on the site $`i=1`$. The update then proceeds for all sites in parallel. We look for all sites which satisfy at least one of the two conditions (i) it just toppled, (ii) its right-hand or left-hand neighbour toppled ama\_la\_96 . If $`i`$ is such a site, it topples with probability 1 if $`z_i>2`$, with probability $`\alpha [0,1]`$ if $`z_i=2`$ and with probability 0 if $`z_i<2`$. A toppling at the site $`i`$ means that $`z_i`$ is decreased by 2 and $`z_{i1}`$ and $`z_{i+1}`$ are increased by 1.

For $`\alpha =0`$ or $`\alpha =1`$ we recover the standard one-dimensional BTW sandpile model with critical slope $`z_c=1`$ or $`z_c=2`$, respectively. In the intermediate region, $`0<\alpha <1`$, self-organized criticality was found in numerical simulations, with avalanche exponent $`\tau =1.55\pm 0.02`$ ama\_la\_97 . However, it is not clear, what is the behaviour of the model for $`\alpha `$ close to either 1 or 0. It seems, that for a finite system the behaviour is SOC (modified by finite size effects) only if $`\alpha `$ is not too close to 1 or 0 be\_be\_ke\_lo\_mhi\_99 ; markosova-unpublished . The behaviour of the system when the system size diverges and $`\alpha `$ stays close 0 or 1 was not clarified. We would like to study this question within the approximation provided by a self-organized branching process.

### 2.2 Characteristic functions

From the technical point of view we will use the method of characteristic function (discrete Laplace transform), defined for a function $`f(s)`$ on integer numbers $`s`$ as $`\widehat{f}(\zeta )=_{s=0}^{\mathrm{}}\zeta ^sf(s)`$.

We will see that the distribution of avalanches have generic form

$$P(s)s^\tau \mathrm{e}^{s/s_0}$$
(1)
for large $`s`$. In the mean-field approximation or in the branching process we have $`\tau =3/2`$, while in one-dimensional BTW sandpile the exponent is $`\tau =0`$. The process is critical, if the cutoff avalanche size $`s_0`$ diverges, $`s_0\mathrm{}`$.

In the language of characteristic functions the behaviour (1) translates in the properties of the singularity in $`\widehat{P}(\zeta )`$. Generally we have $`\widehat{P}(\zeta )(\zeta \zeta _0)^\eta +nonsingularpart`$. For the one-dimensional BTW process we have $`\eta =1`$, while true branching process has $`\eta =1/2`$. The cutoff is given by the distance of the singularity from the point $`\zeta =1`$, namely $`s_01/|\zeta _01|`$. The process is critical, if $`\zeta _0=1`$.

We will also see that the characteristic function for the branching process is typically the solution of a quadratic equation. The singular part of the characteristic function comes from the square root of the discriminant $`D(\zeta )`$ of the equation, i.e. $`\widehat{P}(\zeta )\sqrt{D(\zeta )}+nonsingularpart`$. Therefore, $`\eta =1/2`$ and the cutoff is given by the solution of the equation $`D(\zeta _0)=0`$. If $`D(1)=0`$, we have $`s_0=\mathrm{}`$ and the process is critical.

### 2.3 Branching process

Let us first recall how the branching process is used to describe the simplest case of the sandpile model, for which in each toppling event two grains are transferred to two randomly chosen nearest neighbours (Manna model manna\_91a ). There are $`N_0`$ sites in state $`z=0`$ and $`N_1`$ sites in state $`z=1`$. The branching process starts by dropping a grain to a randomly chosen site. The probability of becoming active (to topple) is $`p=\frac{N_1}{N_0+N_1}`$. Two new branches arise from an active site. Each of them is active with probability $`p`$ and a tree is created iteratively. The branching process stops, when no active sites are present at the end-points of the tree. The number of active sites, or number of branchings, corresponds to the size of the avalanche. The probability distribution of avalanche sizes can be easily obtained with the use of characteristic functions za\_la\_sta\_95 ; la\_za\_sta\_96 ; ca\_te\_ste\_96 ; ve\_ma\_ba\_97 ; chu\_ada\_99  and gives the mean-field value of the exponent $`\tau =3/2`$

Approximating the sand- or ricepile models by branching process is well justified in high dimensions, where the activity returns to the same point with very small probability. It seems, therefore, that the use of branching processes in the opposite limit, in one dimension, lacks sense, because the return of activity is very frequent. However, we can use a very simple property of the return of activity to make the approximation sensible. Indeed, the most frequent case when the activity returns to the same site is described by the following process.

If the site $`i`$ is active (it topples), a grain is transferred to site $`i+1`$ which can become active. If that happens, another grain is transferred back to site $`i`$ (and also to site $`i+2`$, but it is not important now) and thus the site $`i`$ may become active again. This observation leads to the modification of the branching process suitable for the one-dimensional case. We should take into account explicitly the return of the activity just in the next step. We will do it by setting different branching probabilities for site which was active just one step before (it is the site to the left) and for the site which did not have to be (the site to the right).

Because the grains are added only on the site $`i=1`$, we have $`z_i0i`$. The condition that the site topples with probability 1 if $`z>2`$ ensures that $`z_i2i`$. We denote $`N_a`$ number of sites with $`z=a`$. So, picking randomly a site, we have probability $`p_a=N_0/(N_0+N_1+N_2)`$ of having $`z=a`$, where $`a=0,1,2`$.

Let us now describe the construction of the branching process corresponding to the one-dimensional ricepile. There are three types of the points on the tree created by the branching process, according to the value of $`z\{0,1,2\}`$. We denote $`q_a`$ the probability that a point with $`z=a`$ do branch. The points with $`z=0`$ do not branch, i. e. $`q_0=0`$, while the points with $`z=2`$ always branch, so $`q_2=1`$. The points with $`z=1`$ branch with probability $`\alpha `$, i. e. $`q_1=\alpha `$. The approximation consists in supposing that if a site did not topple in the previous step, it has probability $`p_a`$ of having $`z=a`$, while if the site did topple in the last step, the probability of having $`z=a`$ is modified due to the previous toppling to the value

$$p_a^{}=\frac{q_{a+1}p_{a+1}}{_{b=0}^2q_{b+1}p_{b+1}}$$
(2)
where we used $`p_3=q_3=0`$ for convenience.

If a branching occurs at a site, two new branches (“left” and “right”) emanate from it. The probability that the right branch ends with a point with $`z=a`$ is $`p_a`$, while for the left branch the probability is $`p_a^{}`$. This way the tree corresponding to the branching process is created. The above described rules are illustrated in the Fig. 1.

The root of the tree should be treated separately. The reason is that in the ricepile model the avalanche starts by dropping a grain always on the left edge of the pile, i. e. on the site $`i=1`$. If it topples, it transfers a grain only to the right, while the grain going left falls off the system. If we translated this feature to the description of our branching process, the root would consist either of a single non-branching point, or a point with a single branch (the right one) emerging from it. However, we are interested in the regime of long trees, where the different behaviour of the root from the rest of the tree is irrelevant. So, we assume that in the branching process also the root obeys the same rules as all other points. Thus, all points, including the root, have either zero or two branches emanating from it.

The key quantity will be $`P_n^a(s)`$, the probability that a tree consisting of $`n`$ levels starting with a point of type $`a`$ contains $`s`$ branchings. The probability of having $`s`$ branchings (i. e. avalanche of size $`s`$) is then $`P_n(s)=_ap_aP_n^a(s)`$. We can easily derive the recurrence relation for $`P_n^a(s)`$ which becomes particularly simple if we use the characteristic function. We obtain

$$\widehat{P}_n^a(\zeta )=(1q_a)+q_a\zeta \underset{b,c=0}{\overset{2}{}}p_bp_c^{}\widehat{P}_{n1}^b(\zeta )\widehat{P}_{n1}^c(\zeta ).$$
(3)
A straightforward calculation leads to the following equations for the characteristic functions

$$\begin{array}{cc}\hfill \widehat{P}_n^0(\zeta )& =1\hfill \\ \hfill \widehat{P}_n^1(\zeta )& =1\alpha +\alpha \widehat{P}_n^2(\zeta )\hfill \end{array}$$
(4)
and

$$P_n(s)=(\alpha p_1+p_2)P_n^2(s)\mathrm{for}s>0.$$
(5)
Therefore the basic quantity of interest will be the characteristic function $`\widehat{P}_n^2(\zeta )`$. All properties of the branching process can be computed from it. The set of equations (3) thus represent a single recurrence equation for $`\widehat{P}_n^2(s)`$, which in the limit $`n\mathrm{}`$ leads to quadratic equation for the stationary distribution $`\widehat{P}^2(\zeta )=lim_n\mathrm{}\widehat{P}_n^2(\zeta )`$. We obtain explicitly

$$\begin{array}{cc}\hfill \frac{1}{\zeta }& \widehat{P}^2(\zeta )=\hfill \\ & \frac{\left(\alpha p_1+(1\alpha )p_2\right)\left(1\alpha p_1p_2\right)}{p_2+\alpha p_1}+\hfill \\ & +\frac{\alpha p_2+2\alpha (1\alpha )p_1p_2+(12\alpha )p_{2}^{}{}_{}{}^{2}+p_{1}^{}{}_{}{}^{2}\alpha ^2}{p_2+\alpha p_1}\widehat{P}^2(\zeta )\hfill \\ & +\alpha p_2(\widehat{P}^2(\zeta ))^2.\hfill \end{array}$$
(6)
### 2.4 Criticality

The discriminant $`D(\zeta )`$ of the equation (6) depends on the parameters $`p_1`$, $`p_2`$, and $`\alpha `$. The branching process is critical if $`D(1)=0`$. This implies the following relation

$$\alpha p_1(1\alpha )p_2+2\alpha p_1p_2+p_{1}^{}{}_{}{}^{2}\alpha ^2+p_{2}^{}{}_{}{}^{2}=0$$
(7)
which determines a surface in the parametric space. On this surface the process is critical and the distribution of avalanche sizes has a power-law tail with exponent $`\tau =3/2`$.

However, the latter statement is not strictly true in the sense that if the coefficient at the quadratic term in the equation (6) is zero, the process is not a true branching process, because each parent can have at most one offspring. It corresponds to a process with an exponential distribution of avalanche sizes, which we will call, in this work, a “one-dimensional BTW process”. The important feature which makes it different from a generic branching process is that there are no true branching points. Indeed, there may be a non-zero probability that the process stops at a given point, but there is zero probability to be split into more than one branch. Therefore, the process does not generate a tree-like structures, but linear chains of random length. Both one-dimensional BTW and branching processes have the same general form (1) of the distribution of avalanches for large $`s`$, but the one-dimensional BTW process is characterised by the exponents $`\tau =0,\eta =1`$. Therefore, together with checking the criticality condition (7) we must also look at the behaviour close to the singularity.

We will prove in the section 3.2 that in the thermodynamic limit our ricepile model self-organizes so that the parameters stabilise at values

$$\begin{array}{cc}\hfill p_1=& \mathrm{max}(0,\frac{2\alpha 1}{\alpha })\hfill \\ \hfill p_2=& 1\alpha .\hfill \end{array}$$
(8)
If we insert these values into the criticality condition (7), we find that it is satisfied for any value of $`\alpha `$, including the limit values of 0 and 1. At the same time we find that the singularity is always located at $`\zeta _0=1`$. (Indeed, as we discussed in section 2.2, the criticality of the process is equivalent to the condition $`\zeta _0=1`$.) However, we find that the type of the singularity corresponds to the exponents $`\eta =1/2`$, $`\tau =3/2`$ (critical branching process) only for $`\alpha `$’s within the open interval $`(0,1)`$, while at the points 0 and 1 the model corresponds to one-dimensional BTW process. This can be easily interpreted in the language of sand- and ricepiles. Indeed, for $`\alpha =0`$ and $`1`$ the system recovers the behaviour of one-dimensional BTW sandpile, which does not exhibit critical behaviour in the usual sense. (In fact, the avalanche distribution does exhibit a power-law distribution: all avalanche sizes have the same probability, which corresponds to the power with exponent 0. But this is not the situation we usually describe as critical behaviour.)

### 2.5 Crossover behaviour

The question arises, how the behaviour with exponent $`\tau =3/2`$ inside the interval $`[0,1]`$ crosses over to the exponent $`\tau =0`$ at the edges. As the critical behaviour is related to the singularities of the characteristic function, we will turn to the investigation of the function $`\widehat{P}^2(\zeta )`$ in more detail.

Indeed, we find that if we expand the solution of Eq. (6) for small values of the parameter $`\rho `$ defined as

$$\rho (\zeta )=\frac{2\alpha (1\alpha )}{\zeta ^11+2\alpha (1\alpha )}$$
(9)
we can express the solution in terms of $`\rho `$ and expand in the lowest order (for $`\rho ^21`$)

$$\widehat{P}^2(\zeta )=\frac{1}{\rho }\sqrt{\frac{1}{\rho ^2}1}\frac{\rho (\zeta )}{2}.$$
(10)
While, as noted earlier, the exact solution for $`\widehat{P}^2(\zeta )`$ has always the singularity of the type $`\eta =1/2`$ for $`\zeta \zeta _0=1`$, the approximate behaviour (10) has a singularity with $`\eta =1`$ located at the point $`\zeta _0^{}=(12\alpha (1\alpha ))^1>1`$. When $`\alpha `$ goes to either 0 or 1, the value of $`\zeta _0^{}`$ approaches 1. This suggests the following scenario. For large avalanches, i. e. $`1\zeta \zeta _0^{}1`$ the singularity at $`\zeta _0=1`$ is relevant and the avalanche size distribution has a power-law tail with exponent $`\tau =3/2`$.

For shorter avalanches, i. e. $`1\zeta `$ larger or comparable to $`\zeta _0^{}1`$ the singularity at $`\zeta _0^{}`$ becomes dominant. Therefore, for short avalanches we have one-dimensional BTW behaviour $`P(s)\mathrm{exp}(s/s_0)`$ with a cutoff

$$s_0=|1\zeta _0^{}|^1=\frac{12\alpha (1\alpha )}{2\alpha (1\alpha )}.$$
(11)
The next step is investigation of the behaviour of $`s_0`$ when $`\alpha `$ approaches either 0 or 1. We find it by expanding the expression for $`\zeta _0^{}`$ as a function of $`\alpha `$ around the points 0 and 1, respectively. To make the notation more compact, let us introduce the variable $`\mu \{0,1\}`$, which distinguishes the two limit points $`\alpha =0`$ and 1. We can see from (11) that the cutoff diverges as

$$s_0\frac{1}{2|\alpha \mu |}$$
(12)
for $`\alpha \mu `$.

On the other hand, sufficiently close to the singularity at $`\zeta \zeta _0=1`$ the exponent $`\eta =1/2`$ is relevant. The question is, how close to $`\zeta =1`$ one behaviour crosses over to the other. We have one-dimensional BTW behaviour for $`\rho ^21`$, while critical branching process behaviour for $`1\rho ^21`$. A typical crossover value $`\zeta _{\mathrm{cr}}`$ can be found by solving the equation

$$\rho (\zeta _{\mathrm{cr}})=\frac{1}{2}.$$
(13)
The avalanche size distribution will exhibit the crossover around $`s_{\mathrm{cr}}=1/|1\zeta _{\mathrm{cr}}|`$. For $`ss_{\mathrm{cr}}`$ the one-dimensional BTW behaviour with exponential cutoff, diverging to infinity for $`\alpha =0`$ and 1, will apply, while for $`ss_{\mathrm{cr}}`$ the distribution will have power-law tail with usual mean-field exponent $`3/2`$, and therefore exhibits self-organized criticality.

The point of the transition between SOC and one-dimensional BTW when $`\alpha `$ approaches to 1 or 0 lies in the diverging crossover value for the avalanche size. Similarly as in the case of $`s_0`$, by solving the equation (13) with defintition (9) we find the following limit behaviour

$$s_{\mathrm{cr}}\frac{1}{2|\alpha \mu |}s_0$$
(14)
for $`\alpha \mu `$.

We can see, comparing equations (12) and (14), that the cutoff for the one-dimensional BTW behaviour is asymptotically equal to the crossover at which the critical branching process behaviour sets on. This suggests the scaling form

$$P(s)\frac{1}{s_0(\alpha )}F(\frac{s}{s_0(\alpha )})$$
(15)
valid for $`s1`$ and $`\alpha `$ close to 0 and 1. The scaling function has the form $`F(x)\mathrm{e}^x`$ for $`x1`$ and $`F(x)x^{3/2}`$ for $`x1`$. Indeed, we can find the Laplace transform of the scaling function as

$$_0^{\mathrm{}}\mathrm{e}^{x(y1)}F(x)dx=y\sqrt{y^21}.$$
(16)
From here we obtain immediately the expression for the scaling function through the Bessel function of imaginary argument

$$F(x)=\frac{\mathrm{e}^x}{x}I_1(x).$$
(17)
The expected behaviour for $`x1`$ and $`x1`$ can be verified directly by inspecting the asymptotic behaviour of the Bessel function.

## 3 Self-organization

### 3.1 Self-organized branching process

In the basic setup of our branching process, all three parameters $`\alpha `$, $`p_1`$, $`p_2`$ are freely chosen. However, in the ricepile model the only free parameter is $`\alpha `$. The number of sites with given $`z`$ can change during an avalanche, so that also the probabilities $`p_1`$ and $`p_2`$ are modified. This defines a flow in the space of parameters $`p_1`$, $`p_2`$. Our task now is to establish stable fixed points of this dynamics and check whether they satisfy the condition (7). If that happens, we can conclude that the system is self-organized critical.

There are four types of events, which can happen during an avalanche. Let us denote them $`T2`$, $`T1`$, $`E1`$, $`E0`$. In the event $`T2`$, the point with $`z=2`$ receives a grain and topples. As a result, the number of sites with $`z=2`$ is decreased by 1, $`N_2N_21`$, and number of sites with $`z=1`$ is increased by 1, $`N_1N_1+1`$. Similarly, in the event $`T1`$ point with $`z=1`$ topples, $`N_1N_11`$ and $`N_0N_0+1`$, in event $`E1`$ site with $`z=1`$ receives a grain but does not topple, $`N_1N_11`$ and $`N_2N_2+1`$, and finally in event $`E0`$ site with $`z=0`$ receives a grain and does not topple, $`N_0N_01`$ and $`N_1N_1+1`$.

Using the variables $`y\{T,E\}`$ and $`a,b\{0,1,2\}`$, let us denote $`s_{ayb,n}`$ the number of events of the type $`yb`$ occurring at the level $`n`$ within the branching process, on condition that the very first site had $`z=a`$. There are $`s_{ayb}=_{n=0}^{\mathrm{}}s_{ayb,n}`$ such events in the entire realisation of the branching process. On average, there are $`s_{yb}=_ap_as_{ayb}`$ events of the type $`yb`$. The averages $`s_{yb}`$ are of central importance for the dynamics of the self-organization and can be easily obtained as follows.

For the characteristic function of the probability distributionof the number of events $`s_{ayb,n}`$ we obtain an equation analogical to (3). To study the self-organization, we will need only the average number of events, which is $`s_{ayb,n}`$, calculated as derivative of the characteristic function. Hence

$$s_{ayb,n}=q_a\underset{c}{}(p_c+p_c^{})s_{cyb,n1}.$$
(18)
This is a set of three recurrence relations, which may be reduced to one equation only, by considering the relations $`s_{0yb,n}=0`$ and $`s_{1yb,n}=\alpha s_{2yb,n}`$, valid for $`n>1`$. If we take as the basic quantity the average $`s_{2yb,n}`$, we get a recurrence relation determining a geometric sequence

$$s_{2yb,n+1}=\kappa s_{2yb,n}$$
(19)
with quotient

$$\kappa =\frac{\alpha p_2+(p_2+\alpha p_1)^2}{p_2+\alpha p_1}.$$
(20)
We recognise in the stationarity condition $`\kappa =1`$ the equation (7), implying the criticality of the branching process.

Summation of the infinite geometric series gives immediately

$$s_{yb}=\left(p_b+(p_b+p_b^{})\frac{\alpha p_1+p_2}{1\kappa }\right)s_{byb,1}$$
(21)
where the initial conditions are given by $`s_{bTb,1}=q_b`$ and $`s_{bEb,1}=1q_b`$.

The self-organization of the branching process is due to the changes in the numbers $`N_a`$, caused by the toppling (and non-toppling) events. These numbers determine the probabilities $`p_a`$. Therefore, for fixed $`\alpha `$ the self-organized branching process (SOBP) $`𝒮(\alpha )`$ consists of an (infinite) sequence of branching processes

$`𝒮(\alpha )=`$ (22)
$`=[(\alpha ,p_1^{(0)},p_2^{(0)}),(\alpha ,p_1^{(1)},p_2^{(1)}),(\alpha ,p_1^{(2)},p_2^{(2)}),\mathrm{}]`$
where $`(\alpha ,p_1,p_2)`$ is the branching process determined by fixed parameters $`\alpha ,p_1,p_2`$, defined above. The branching processes within the sequence differ only by the values of the parameters $`p_1`$, and $`p_2`$. Let us consider the $`t`$-th branching process in the sequence. When realised, it changes the original values of the numbers $`N_a`$, or, equivalently, the values of the parameters $`p_a`$. The average change is uniquely determined by the average number of events $`s_{yb}`$ . So, the SOBP is entirely determined by the transition relations connecting the values of the parameters in the $`t`$-th and $`t+1`$-st step

$$p_i^{(t+1)}p_i^{(t)}=T_i(p_1^{(t)},p_2^{(t)})$$
(23)
for $`i\{1,2\}`$. We find explicitly

$$\begin{array}{cc}& T_1(p_1,p_2)=\hfill \\ & \frac{\alpha p_1+(1\alpha )p_2\alpha (2\alpha )p_1^2p_2^22(1\alpha )p_1p_2}{\alpha p_1+(1\alpha )p_2+2\alpha p_1p_2+p_2^2+\alpha ^2p_1^2}\hfill \\ & T_2(p_1,p_2)=\hfill \\ & \frac{\alpha (1\alpha )p_1^2+(12\alpha )p_1p_2}{\alpha p_1+(1\alpha )p_2+2\alpha p_1p_2+p_2^2+\alpha ^2p_1^2}.\hfill \end{array}$$
(24)
### 3.2 Fixed point

The fixed point of the self-organization dynamics can be found immediately by equating the right-hand sides of equations (24) to zero. Direct solution of the two coupled equations gives three fixed points

$`p_1=`$ $`0,`$ $`p_2=0`$ (25)
$`p_1=`$ $`0,`$ $`p_2=1\alpha `$ (26)
$`p_1=`$ $`{\displaystyle \frac{2\alpha 1}{\alpha }},`$ $`p_2=1\alpha .`$ (27)
The correct solution is determined by stability considerations. The relations (24) are linearised around the fixed points and the eigenvalues of the resulting matrices of rank 2 are found. The result is that the fixed point (25) is always unstable, while (26) is stable for $`\alpha [0,1/2)`$ and (27) is stable for $`\alpha (1/2,1]`$. For $`\alpha =1/2`$ the fixed points (26) and (27) coincide and both of them are marginally stable (i. e. the eigenvalues have zero real part).

Therefore, we find that the fixed point corresponds to the values of the probabilities

$$\begin{array}{cc}\hfill p_1=& \mathrm{max}(0,\frac{2\alpha 1}{\alpha })\hfill \\ \hfill p_2=& 1\alpha \hfill \end{array}$$
(28)
which proves the already announced result of Eq.(8).

### 3.3 Finite-size effects

In the numerical simulations of the ricepile model be\_be\_mhi\_zha\_99 ; be\_be\_ke\_lo\_mhi\_99 ; markosova\_00 ; markosova-unpublished  attention is paid to the fact that the critical behaviour is observed only for large enough systems with $`\alpha `$ not too close to neither 0 nor 1. We have already shown how the crossover length blows up when $`\alpha `$ approaches the edge values 0 or 1. It is obvious then, that for small systems the crossover value of the avalanche size may not be accessible and the critical regime in the tail of the distribution is not observed at all. In this subsection we will investigate the consequences of finite length of the branching process. There are two phenomena where the finite size enters the problem. First, if the maximum number of generations in the branching process is $`L`$, instead of infinity, the distribution of the avalanche sizes will not extend to infinity either, but will be bounded by $`s<s_{\mathrm{max}}=2^L1`$. Moreover, if we take for example $`p_1=1`$, $`p_2=0`$, $`\alpha =1`$, then all avalanches will have size $`L`$, therefore a peak at $`s=L`$ will occur, $`P(s)=\delta (sL)`$. If we move slightly from this position by increasing $`p_2`$ and decreasing $`p_1`$ and $`\alpha `$, a structure of multiple peaks located at $`s=L,2L1,3L3,\mathrm{}`$ will appear. This makes the analysis very complicated.

Second consequence is the shift in the self-organized value of the parameters $`p_1`$ and $`p_2`$, which for finite $`L`$ will deviate from the critical values. Therefore, the avalanche-size distribution will develop an exponential cutoff in the form $`P(s)s^{3/2}\mathrm{exp}(s/s_1)`$.

As the first problem brings new particular difficulties, we will concentrate only on the second one. This makes the analysis less consistent, but feasible. Thus, we should stress that in the following we will suppose that the branching process in question has unbounded length, but the self-organization is made in such a way, that only the first $`L`$ generations of the branching process are taken into account.

Instead of working with finite-$`L`$ version of the equations (23) and (24), describing the approach to the fixed point, we can use the set of equations

$$\begin{array}{cc}\hfill s_{E1}=& s_{T2}\hfill \\ \hfill s_{E0}=& s_{T1}\hfill \end{array}$$
(29)
which determine the position of the fixed point. The only information lost in Eq. (29) is the stability of the fixed points. However, we suppose the stability will not be affected by the finite-size effects. Therefore, we will rely on the stability analysis performed for $`L=\mathrm{}`$ also in the case of finite $`L`$ and calculate the finite-size corrections starting with Eq. (29).

The point is that the equations (29) should hold also for finite $`L`$. In fact, the expression (21) for the averages $`s_{yb}`$ assume the same form, only the factor $`(\kappa 1)`$ arising in the $`L=\mathrm{}`$ version should be replaced by the factor $`K=(\kappa 1)/(\kappa ^{L1}1)`$. Assuming $`K`$ small for large $`L`$, we can find $`p_1`$ and $`p_2`$ in lowest order in $`K`$. Then, we return to the definition of $`K`$ and find that $`KL^1`$, confirming that our approach is consistent.

Hence, for finite $`L`$ we find, by solving the equations (29) to lowest order in $`1/L`$, for $`\alpha (0,1/2)`$

$$\begin{array}{cc}\hfill p_1& =\frac{1\alpha }{(2\alpha 1)^2}\frac{\mathrm{ln}(1\gamma )}{L}+O(\frac{1}{L^2})\hfill \\ \hfill p_2& =1\alpha \frac{1\alpha }{2\alpha 1}\frac{\mathrm{ln}(1\gamma )}{L}+O(\frac{1}{L^2})\hfill \end{array}$$
(30)
and for $`\alpha (1/2,1)`$

$$\begin{array}{cc}\hfill p_1& =\frac{2\alpha 1}{\alpha }+\frac{5\alpha ^25\alpha +1}{(2\alpha 1)^2\alpha }\frac{\mathrm{ln}(1\gamma )}{L}+O(\frac{1}{L^2})\hfill \\ \hfill p_2& =1\alpha +\frac{1\alpha }{2\alpha 1}\frac{\mathrm{ln}(1\gamma )}{L}+O(\frac{1}{L^2})\hfill \end{array}$$
(31)
where we denoted

$$\begin{array}{cc}\hfill \gamma & =\frac{1}{2}\frac{12\alpha }{1\alpha }\mathrm{for}\alpha (0,1/2)\hfill \\ \hfill \gamma & =\frac{1}{2}\frac{2\alpha 1}{\alpha }\mathrm{for}\alpha (1/2,1).\hfill \end{array}$$
(32)
The above formulae confirm that the explicit limit $`L\mathrm{}`$ gives the same result as obtained previously when working directly with $`L=\mathrm{}`$.

Using these results we can find the position of the square-root singularity in the characteristic function for the avalanche size distribution, solving the equation$`D(\zeta _0)=0`$. The distance from 1 then determines the exponential cutoff of the distribution. We find

$$1/s_1=|\zeta _01|=\frac{\sigma (\alpha )}{L^2}+O(\frac{1}{L^3})$$
(33)
where

$$\sigma (\alpha )=\frac{\mathrm{ln}^2(1\gamma )}{4\alpha (1\alpha )}\mathrm{for}\alpha (0,1)$$
(34)
and asymptotically for $`L\mathrm{}`$ and $`\alpha `$ fixed the avalanche distribution becomes the function of $`sL^2`$ only,

$$P(s;\alpha ,L)L^3G(sL^2\sigma (\alpha ))$$
(35)
and he scaling function has the form

$$G(x)=x^{3/2}\mathrm{e}^x.$$
(36)
This scaling holds well for all $`\alpha `$ with exception of the point $`\alpha =1/2`$, where we have $`\gamma =0`$ and hence $`\sigma (\alpha )=0`$. Then, the next order in $`1/L`$ takes over and the scaling changes.

Let us use again the variable $`\mu \{0,1\}`$, which distinguishes the two limit points $`\alpha =0`$ and 1. The factor $`\sigma (\alpha )`$ diverges as $`\sigma (\alpha )\sigma _0|\alpha \mu |^1`$ for $`\alpha \mu `$, where $`\sigma _0=(\mathrm{ln}2)^2/4`$. Therefore, we can write the following scaling form for the avalanche size distribution

$$P(s;\alpha ,L)L^3|\alpha \mu |^{\frac{3}{2}}G(sL^2|\alpha \mu |^1\sigma _0)$$
(37)
for $`\alpha \mu `$.

We can see that the power-law distribution holds only for avalanches shorter than $`L^2|\alpha \mu |`$. In other words, if the parameter $`\alpha `$ is close to the end-points of the interval $`[0,1]`$, we need to have systems of the size $`L1/\sqrt{|\alpha \mu |}`$ for being able to observe any sign of self-organized criticality.

In the above calculations we tacitly assumed that we beyond the regime we called “one-dimensional BTW” in the section 2.5. This means $`ss_{\mathrm{cr}}`$. In fact, we can always reach this regime by choosing $`L`$ large enough. Therefore the presence of the one-dimensional regime does not influence the scaling behaviour for large $`L`$. More precisely, we should have $`L^2|\alpha \mu |s_{\mathrm{cr}}`$. But because $`s_{\mathrm{cr}}`$ itself diverge for $`\alpha \mu `$ as $`|\alpha \mu |^1`$, we obtain a stronger condition for the scaling (37) to be valid, namely

$$L|\alpha \mu |^1$$
(38)
if $`\alpha \mu `$.

## 4 Conclusions

We investigated analytically the self-organized critical ricepile model. We defined a self-organized branching process, suitable for one-dimensional problems. The model is characterised by the parameter $`\alpha [0,1]`$, the probability of toppling at a sub-threshold site. For both limit values $`\alpha =0`$ and $`\alpha =1`$ the model is equivalent to the one-dimensional BTW model with trivial (uniform) distribution of avalanches.

We found that in the thermodynamic limit the system is self-organized critical for all values of $`\alpha `$ within the open interval $`(0,1)`$, with power-law tail in the distribution of avalanche sizes with mean-field value of the exponent, $`\tau =\frac{2}{3}`$. However, the power-law behaviour holds only for avalanches longer than a certain crossover value of the avalanche size. The crossover diverges when $`\alpha `$ approaches to either of the limit points of the interval $`[0,1]`$. We found the scaling as well as the exact form of the scaling function for avalanche distribution close to these limit points. This describes how the one-dimensional BTW behaviour develops when approaching the limit points.

The finite-size effects play important role in determining whether the model is self-organized critical or not. In our model the SOC behaviour starts to occur at the larger sizes the closer we are to the limit points $`\alpha =0`$ or $`1`$. We found the form of the finite size scaling in our self-organized branching process and determined the necessary condition for the the power-law regime in the avalanche distribution to be observable, when we approach to the limit points.

###### Acknowledgements.

I wish to thank Mária Markošová for numerous useful discussions which motivated me to this work. I am indebted to Petr Chvosta for clarifying comments regarding stochastic processes. This work was partially supported by the Grant Agency of the Czech Republic, project No. 202/00/1187.


text:
# A quasi-pure Bose-Einstein condensate immersed in a Fermi sea

## Abstract

We report the observation of co-existing Bose-Einstein condensate and Fermi gas in a magnetic trap. With a very small fraction of thermal atoms, the 7Li condensate is quasi-pure and in thermal contact with a 6Li Fermi gas. The lowest common temperature is $`0.28\mu `$K $`0.2(1)T_\mathrm{C}=0.2(1)T_\mathrm{F}`$ where $`T_\mathrm{C}`$ is the BEC critical temperature and $`T_\mathrm{F}`$ the Fermi temperature. Behaving as an ideal gas in the radial trap dimension, the condensate is one-dimensional.

Bose-Einstein condensation of atomic gases has been very actively studied in recent years . Because of the dilute character of the samples and the ability to control the atom-atom interactions, for the first time a detailed comparison with the theories of quantum gases has been made. Fermi gases, on the other hand, have only been investigated experimentally for two years . They are predicted to possess intriguing properties and may offer an interesting link with the behavior of electrons in metals and semi-conductors, and the possibility of Cooper pairing  such as in high temperature superconductors and neutron stars. Mixtures of bosonic and fermionic quantum systems, with the prominent example of 4He- 3He fluids, have also stimulated intense theoretical and experimental activity . This has led to new physical effects including phase separation, influence of the superfluidity of the Bose system on the Fermi degeneracy and to new applications such as the dilution refrigerator .

In this paper, we present a new mixture of bosonic and fermionic systems, a stable Bose-condensed gas immersed in a Fermi sea. Confined in the same magnetic trap, both atomic species are in thermal equilibrium with a temperature of $`0.2(1)T_F<<T_\mathrm{C}`$. The condensate is made of 7Li atoms in internal state $`|F=1,m_F=1`$ while the degenerate Fermi gas is made of 6Li atoms in $`|F=1/2,m_F=1/2`$ (fig. 1). All previous experiments performed with 7Li in $`|F=2,m_F=2`$ had condensate numbers limited to $`N1400`$ because of the negative scattering length, $`a=1.4`$nm, in this state . Our condensate is produced in a state which has a positive, but small, scattering length, $`a=+0.27`$nm . The number of condensate atoms is typically $`10^4`$, and BEC appears unambiguously both in the position distribution in the trap and in the standard time of flight images. An interesting feature is the one-dimensional character of this condensate, behaving as an ideal gas in the transverse direction of the trap.

Because of the symmetrization postulate, colliding fermions have no s-wave scattering at low energy and the p-wave contribution vanishes in the low temperature domain of interest. Our method for producing simultaneous quantum degeneracy for both isotopes

of lithium is sympathetic cooling ; s-wave collisions between two different atomic isotopes are allowed and RF evaporation selectively removes from the trap high energy atoms of one species. Elastic collisions subsequently restore thermal equilibrium of the two-component gas at a lower temperature. Our experimental setup has been described in detail in . A mixture of 6Li and 7Li atoms is loaded from a magneto-optical trap into a strongly confining Ioffe-Pritchard trap at a temperature of about $`2`$mK. As depicted in fig.1, this relatively high temperature precludes direct magnetic trapping of the atoms in their lower hyperfine state because of the shallow magnetic trap depth, 2.4 mK for 7Li in $`|F=1,m_F=1`$ and 0.2 mK for 6Li in $`|F=1/2,m_F=1/2`$. Therefore we proceed in two steps. Both isotopes are first trapped and cooled in their upper hyperfine states. There is no limitation for the trap depth and the confinement is better because the magnetic moment is higher. Evaporation is performed selectively on 7Li using a microwave field near $`803`$ MHz that couples $`|F=2,m_F=2`$ to $`|F=1,m_F=1`$. When both gases are cooled to a common temperature of about $`9\mu `$K, atoms are transferred using a combination of microwave and RF pulses into states $`|F=1,m_F=1`$ and $`|F=1/2,m_F=1/2`$ with an energy far below their respective trap depths. Evaporative cooling is then resumed until the BEC threshold is reached for 7Li.

In the first series of experiments, both Li isotopes are trapped in their higher HF states. 6Li is sympathetically cooled to Fermi degeneracy by performing 30 seconds of evaporative cooling on 7Li . Trap frequencies for 7Li are $`\omega _{\mathrm{rad}}=2\pi 4000(10)s^1`$ and $`\omega _{\mathrm{ax}}=2\pi 75.0(1)s^1`$ with a bias field of 2 G. Absorption images of both isotopes are recorded on a single CCD camera with a resolution of $`10\mu `$m. Images are taken quasi-simultaneously (only 1 ms apart) in the trap or after a time of flight expansion. Probe beams have an intensity below saturation and a common duration of $`30\mu `$s. Typical in-situ absorption images in the quantum regime can be seen in fig.2. Here the temperature $`T`$ is $`1.4(1)\mu `$K and $`T/T_\mathrm{F}=0.33(5)`$, where the Fermi temperature $`T_\mathrm{F}`$ is $`(\mathrm{}\overline{\omega }/k_\mathrm{B})(6N_\mathrm{F})^{1/3}`$, with $`\overline{\omega }`$ the geometric mean of the three oscillation frequencies in the trap and $`N_\mathrm{F}`$ the number of fermions. For images recorded in the magnetic trap, the common temperature is measured from the spatial extent of the bosonic cloud in the axial direction since the shape of the Fermi cloud is much less sensitive to temperature changes when $`T/T_\mathrm{F}<1`$. The spatial distributions of bosons and fermions are recorded after a 1 sec thermalization stage at the end of the evaporation. As the measured thermalization time constant between the two gases, $`0.15`$s, is much shorter than $`1`$ s, the two clouds are in thermal equilibrium. Both isotopes experience the same trapping potential. Thus the striking difference between the sizes of the Fermi and Bose gases  is a direct consequence of Fermi pressure. The measured axial profiles in fig.2 are in excellent agreement with the calculated ones (solid lines) for a Bose distribution at the critical temperature $`T_\mathrm{C}`$. In our steepest traps, Fermi temperatures as high as $`11\mu `$K with a degeneracy of $`T/T_F=0.36`$ are obtained. This $`T_\mathrm{F}`$ is a factor 3 larger than the single photon recoil temperature at 671 nm, opening interesting possibilities for light scattering experiments .

Our highest Fermi degeneracy in the 6Li $`F=3/2`$ state is $`T/T_\mathrm{F}=0.25(5)`$ with $`T_\mathrm{F}=4\mu `$K, very similar to ref.. This limit is set by the fact that the boson temperature cannot be lowered below $`T_C`$. Because of the negative scattering length in 7Li, $`|F=2,m=2`$, the number of condensed atoms cannot exceed $`300`$ in our trap . This, together with the condition that sympathetic cooling stops when the heat capacity of the bosons approximately equals that of the fermions, limits the Fermi degeneracy to about 0.3 .

In order to explore the behavior of a Fermi sea in the presence of a BEC with a temperature well below $`T_\mathrm{C}`$, we perform another series of experiments with both isotopes trapped in their lower HF state. As the 7Li scattering length is then positive (fig.1), a stable BEC with high atom numbers is now possible. To avoid large dipolar relaxation, 6Li must also be in its lower HF state . First, sympathetic cooling is performed on the initial stretched states to $`9\mu K`$.

Then, to facilitate state transfer, the trap is adiabatically opened to frequencies $`\omega _{\mathrm{rad}}=2\pi 100s^1`$ and $`\omega _{\mathrm{ax}}=2\pi 5s^1`$ (for 7Li, F=2). The transfer of each isotope uses two microwave $`\pi `$ pulses. The first pulse at 803 MHz for 7Li (228 MHz for 6Li) transfers the atoms from $`|2,2`$ to $`|1,1`$ ($`|3/2,3/2`$ to $`|1/2,1/2`$), a magnetically untrapped state (see fig.1). The second RF $`\pi `$ pulse at 1 MHz (1.3 MHz) transfers them to $`|1,1`$ ($`|1/2,1/2`$) a magnetically trapped state. Adiabatic opening of the trap cools the cloud, therefore decreasing the energy broadening of the resonance and giving more time for the passage through untrapped states. The duration of the $`\pi `$ pulses are $`17\mu s`$ and $`13\mu s`$ and more than $`70\%`$ of each isotope are transferred. Finally the trap is adiabatically recompressed to the steepest confinement giving $`\omega _{\mathrm{rad}}=2\pi 4970(10)s^1`$ and $`\omega _{\mathrm{ax}}=2\pi 83(1)s^1`$ for 7Li $`|F=1,m=1`$, compensating for the reduced magnetic moment.

Because of the very large reduction of the 7Li $`s`$-wave scattering cross section (by a factor of 28) from the F=2 to the F=1 state, we were unable to reach runaway evaporation with 7Li atoms alone in $`F=1`$. In contrast, the 6Li/7Li cross section is 28 times higher than the 7Li/7Li one. We therefore use 6Li atoms as a buffer gas to accelerate the thermalization rate of both gases. Two different methods were used to perform the evaporation. The first consists in using two RF ramps on the HF transitions of 6Li (from $`|1/2,1/2`$ to $`|3/2,3/2`$) and 7Li (from $`|1,1`$ to $`|2,2`$), which we balanced to maintain roughly equal numbers of both isotopes. After 10 s of evaporative cooling, Bose-Einstein condensation of 7Li occurs together with a 6Li degenerate Fermi gas (fig.3). Surprisingly, a single 25 s ramp performed only on 6Li achieved the same results. In this case the equal number condition was fulfilled because of the reduced lifetime of the 7Li cloud that we attribute to dipolar collisional loss .

The duration of the RF evaporation ramp was matched to this loss rate. In the following we concentrate on this second, and simpler, evaporation scheme, sympathetic cooling of 7Li by evaporative cooling of 6Li.

In fig.3 in-situ absorption images of bosons and fermions at the end of the evaporation are shown. The bosonic distribution shows the typical double structure: a strong and narrow peak forms the condensate at the center, surrounded by a much broader distribution, the thermal cloud. As the Fermi distribution is very insensitive to temperature, this thermal cloud is a very useful tool for the determination of the common temperature. Note that, as cooling was only performed on 6Li atoms, the temperature measured on 7Li cannot be lower than the temperature of the fermions. Measuring $`N_\mathrm{B}`$, $`N_\mathrm{F}`$, the condensate fraction $`N_0/N_\mathrm{B}`$, and $`\overline{\omega }`$, we determine the quantum degeneracy of the Bose and Fermi gases. In fig.3(top), the temperature is just below $`T_\mathrm{C}`$, $`T=1.6\mu `$K$`=0.87T_\mathrm{C}=0.57T_\mathrm{F}`$. In fig.3(bottom) on the contrary, the condensate is quasi-pure; $`N_0/N_\mathrm{B}=0.77`$; the thermal fraction is near our detectivity limit, indicating a temperature of $`0.28\mu `$K $`0.2T_\mathrm{C}=0.2T_\mathrm{F}`$ with $`N=\mathrm{8.2\hspace{0.17em}10}^3`$

bosons and $`\mathrm{4\hspace{0.17em}10}^3`$ fermions. Clearly a more sensitive thermal probe is required now to investigate the temperature domain $`T<0.2T_\mathrm{F}`$. An elegant method uses the measurement of thermalization rates with impurity atoms including Pauli blocking . The condensate fraction $`N_0/N_\mathrm{B}`$ as a function of $`T/T_\mathrm{C}`$ is shown in fig.4 (a), while the size of the fermi gas as a function of $`T/T_\mathrm{F}`$ is shown in fig.4 (b). With the strong anisotropy ($`\omega _{\mathrm{rad}}/\omega _{\mathrm{ax}}=59`$) of our trap, the theory including anisotropy and finite number effects differs significantly from the thermodynamic limit , in agreement with our measurements even though there is a $`20\%`$ systematic uncertainty on our determination of $`T_\mathrm{C}`$ and $`T_\mathrm{F}`$.

Because of the small scattering length, this 7Li condensate has interesting properties. The time of flight images, performed after expansion times of 0-10 ms with $`N_0=10^4`$ condensed atoms, reveal that the condensate is a one-dimensional (1D) condensate. In contrast to typical condensates in the Thomas-Fermi (TF) regime, where the release of interaction energy leads to a fast increase in radial size, our measurements agree to better than 5% with the time development of the radial ground state wave function in the harmonic magnetic trap (fig.5). This behavior is

expected when the chemical potential $`\mu `$ satisfies the condition $`\mu <\mathrm{}\omega _{\mathrm{rad}}`$. Searching for the ground state energy of the many-body system with a Gaussian ansatz radially and TF shape axially , we find that the mean-field interaction increases the size of the Gaussian by $`3\%`$. The calculated TF radius is $`28\mu `$m or 7 times the axial harmonic oscillator size and is in good agreement with the measured radius, $`30\mu `$m in fig. 3. Thus with $`\mu =0.45\mathrm{}\omega _{\mathrm{rad}}`$, the gas is described as an ideal gas radially but is in the TF regime axially. This 1D situation was also realized recently in sodium condensates . As $`\mu <\mathrm{}\omega _{\mathrm{rad}}`$ implies that the linear density of a 1D condensate is limited to $`1/a`$, the 1D regime is much easier to reach with 7Li (small $`a`$) than with Na or 87Rb which have much larger scattering lengths.

What are the limits of this BEC-Fermi gas cooling scheme? First, the $`1/e`$ condensate lifetime of about 3 s in this steep trap will limit the available BEC-Fermi gas interaction time. Second, the boson-fermion mean field interaction can induce a spatial phase separation  that prevents thermal contact between 7Li and 6Li. Using the method of  developed for $`T=0`$, we expect, for the parameters of fig.3 (top), that the density of fermions is only very slightly modified by the presence of the condensate in accordance with our observations. Third, because of the superfluidity of the condensate, impurity atoms (such as 6Li), which move through the BEC slower than the sound velocity $`v_\mathrm{c}`$, are no longer scattered . When the Fermi velocity $`v_\mathrm{F}`$ becomes smaller than $`v_\mathrm{c}`$, cooling occurs only through collisions with the bosonic thermal cloud, thus slowing down drastically. With $`10^4`$ condensed atoms, the velocity is $`0.9`$cm/s, corresponding to a limiting temperature of about 100 nK, a factor 3 lower than our currently measured temperature.

In summary, we have produced a new mixture of Bose and Fermi quantum gases. Future work will explore the degeneracy limits of this mixture with the sensitive temperature probe mentioned above . Phase fluctuations of the 1D 7Li condensate should also be detectable via density fluctuations in time of flight images, as recently reported in . The transfer of the BEC into $`|F=2,m=2`$ with negative $`a`$ should allow the production of bright solitons and of large unstable condensates where interesting and still unexplained dynamics has been recently observed . Finally, the large effective attractive interaction between 6Li $`|F=1/2,m_F=+1/2`$ and $`|F=1/2,m_F=1/2`$ makes this atom an attractive candidate for searching for BCS pairing if the temperature can be made sufficiently low .

We are grateful to Y. Castin, J. Dalibard, C. Cohen-Tannoudji, and G. Shlyapnikov for useful discussions. F. S., and K. C. were supported by a fellowship from the DAAD and by MENRT. This work was supported by CNRS and Collège de France. Laboratoire Kastler Brossel is Unité de recherche de l’Ecole Normale Supérieure et de l’Université Pierre et Marie Curie, associée au CNRS.


text:
# 1 Introduction

## 1 Introduction

Noncommutative field theories have recently received much attention after it was realized that in the presence of a background NS B-field, the gauge theory living on D-branes becomes noncommutative . Based on the existence of different regularization procedures in string theory, Seiberg and Witten argued that certain noncommutative gauge theories are equivalent to commutative ones and in particular that there exists a map from a commutative gauge field to a noncommutative one, which is compatible with the gauge structure of each. This map has become known as the Seiberg-Witten (SW) map. In this paper, we give a method for explicitly finding this map. We will consider gauge theories on the noncommutative space defined by

$$\left[x^i\stackrel{}{,}x^j\right]=i\theta ^{ij},$$
(1.1)
where $`\theta `$ is a constant Poisson tensor. Then the “$``$” operation is the associative Weyl-Moyal product

$$fg=fe^{\frac{i}{2}\theta ^{ij}\stackrel{}{_i}\stackrel{}{_j}}g.$$
(1.2)
We believe that our method is much more general, and can in fact be used even when $`\theta `$ is not constant.

In the next section, we review some previous work , which provides an essential starting point for our own. In Section $`3`$, we replace the gauge parameters appearing in the SW map with a ghost field, which makes explicit a cohomological structure underlying the SW map . We then discuss the ambiguities that appear, distinguishing the gauge and covariant types. In Section $`4`$, we define a homotopy operator, which can be used to explicitly write down the SW map order by order in $`\theta `$. In Section $`5`$, we discuss some complications that arise in this formalism and some ways to overcome them.

## 2 General Review

In this section, we review the formalism developed in , which provides an alternative method for obtaining an expression for the SW map.

The original equation which defines the SW map  arises from the requirement that gauge invariance be preserved in the following sense. Let $`a_i`$, $`\alpha `$ be the gauge field and gauge parameter of the commutative theory and similarly let $`A_i`$, $`\mathrm{\Lambda }`$ be the gauge field and gauge parameter of the noncommutative theory. Under an infinitesimal gauge transformation,

$$\delta _\alpha a_i=_i\alpha i[a_i,\alpha ],$$
(2.1)
$$\delta _\mathrm{\Lambda }A_i=_i\mathrm{\Lambda }i[A_i\stackrel{}{,}\mathrm{\Lambda }]_i\mathrm{\Lambda }i\left(A_i\mathrm{\Lambda }\mathrm{\Lambda }A_i\right).$$
(2.2)
Then, the SW map is found by requiring

$$A_i+\delta _\mathrm{\Lambda }A_i=A_i(a_j+\delta _\alpha a_j,\mathrm{}).$$
(2.3)
In order to satisfy (2.3) the noncommutative gauge field and gauge parameter must have the following functional dependence

$$\begin{array}{cc}A_i=A_i(a,a,^2a,\mathrm{})\hfill & \\ \mathrm{\Lambda }=\mathrm{\Lambda }(\alpha ,\alpha ,\mathrm{},a,a,\mathrm{}),\hfill & \end{array}$$
(2.4)
where the dots indicate higher derivatives. It must be emphasized that a SW map is not uniquely defined by condition (2.3). The ambiguities that arise  will be discussed shortly.

The condition (2.3) yields a simultaneous equation for $`A_i`$ and $`\mathrm{\Lambda }`$. For the constant $`\theta `$ case, explicit solutions of the Seiberg-Witten map have been found by various authors up to second order in $`\theta `$ . The solutions were found by writing the map as a linear combination of all possible terms allowed by index structure and dimensional constraints and then determining the coefficients by plugging this expression into the SW equation. The method we will describe in the rest of the paper provides a more systematic procedure for solving the SW map. For the special case of a U(1) gauge group, an exact solution in terms of the Kontsevich formality map is given in , while  present an inverse of the SW map to all orders in $`\theta `$.

An alternative characterization of the Seiberg-Witten map can be obtained following . In the commutative gauge theory, one may consider a field $`\psi `$ in the fundamental representation of the gauge group. If we assume that the SW map can be extended to include such fields, then there will be a field $`\mathrm{\Psi }`$ in the noncommutative theory with the following functional dependence

$$\mathrm{\Psi }=\mathrm{\Psi }(\psi ,\psi ,\mathrm{},a,a,\mathrm{}),$$
(2.5)
and with the corresponding infinitesimal gauge transformation

$$\delta _\alpha \psi =i\alpha \psi ,$$
(2.6)
$$\delta _\mathrm{\Lambda }\mathrm{\Psi }=i\mathrm{\Lambda }\mathrm{\Psi }.$$
(2.7)
An alternative to the SW condition (2.3) can now be given by

$$\mathrm{\Psi }+\delta _\mathrm{\Lambda }\mathrm{\Psi }=\mathrm{\Psi }(\psi +\delta _\alpha \psi ,\mathrm{},a_j+\delta _\alpha a_j,\mathrm{}).$$
(2.8)
More compactly, one writes

$$\delta _{\mathrm{\Lambda }_\alpha }\mathrm{\Psi }(\psi ,a_j,\mathrm{})=\delta _\alpha \mathrm{\Psi }(\psi ,a_j,\mathrm{}).$$
(2.9)
The dependence of $`\mathrm{\Lambda }`$ on $`\alpha `$ is shown explicitly on the left hand side, and on the right hand side $`\delta _\alpha `$ acts as a derivation on the function $`\mathrm{\Psi }`$, with an action on the variables $`\psi `$ and $`a_i`$ given by (2.6) and (2.1) respectively. Next, one considers the commutator of two infinitesimal gauge transformations

$$[\delta _{\mathrm{\Lambda }_\alpha },\delta _{\mathrm{\Lambda }_\beta }]\mathrm{\Psi }=[\delta _\alpha ,\delta _\beta ]\mathrm{\Psi }.$$
(2.10)
Since $`[\delta _\alpha ,\delta _\beta ]=\delta _{i[\alpha ,\beta ]}`$, the right hand side of (2.10) can be rewritten as

$$\delta _{i[\alpha ,\beta ]}\mathrm{\Psi }=\delta _{\mathrm{\Lambda }_{i[\alpha ,\beta ]}}\mathrm{\Psi }=i\mathrm{\Lambda }_{i[\alpha ,\beta ]}\mathrm{\Psi }=\mathrm{\Lambda }_{[\alpha ,\beta ]}\mathrm{\Psi }.$$

The last equality follows from the fact that $`\mathrm{\Lambda }`$ is linear in the ordinary gauge parameter, which is infinitesimal. As for the left hand side,

$$[\delta _{\mathrm{\Lambda }_\alpha },\delta _{\mathrm{\Lambda }_\beta }]\mathrm{\Psi }=\delta _{\mathrm{\Lambda }_\alpha }\left(i\mathrm{\Lambda }_\beta \mathrm{\Psi }\right)\delta _{\mathrm{\Lambda }_\beta }\left(i\mathrm{\Lambda }_\alpha \mathrm{\Psi }\right)$$

$$=i\left(\delta _\alpha \mathrm{\Lambda }_\beta \delta _\beta \mathrm{\Lambda }_\alpha \right)\mathrm{\Psi }+\left[\mathrm{\Lambda }_\alpha \stackrel{}{,}\mathrm{\Lambda }_\beta \right]\mathrm{\Psi }.$$

Equating the two expressions and dropping $`\mathrm{\Psi }`$ yields

$$\left(\delta _\alpha \mathrm{\Lambda }_\beta \delta _\beta \mathrm{\Lambda }_\alpha \right)i\left[\mathrm{\Lambda }_\alpha \stackrel{}{,}\mathrm{\Lambda }_\beta \right]+i\mathrm{\Lambda }_{[\alpha ,\beta ]}=0.$$
(2.11)
An advantage of this formulation is that (2.11) is an equation in $`\mathrm{\Lambda }`$ only, whereas (2.3) must be solved simultaneously in $`\mathrm{\Lambda }`$ and $`A_i`$. If (2.11) is solved, (2.2) with (2.3) then yields an equation for $`A_i`$ and (2.7) with (2.8) for $`\mathrm{\Psi }`$.

## 3 The Ghost Field and the Coboundary Operator

It is advantageous to rewrite equations (2.2), (2.7) and (2.11) in terms of a ghost field in order to make explicit an underlying cohomological structure. Specifically, we replace the gauge parameter $`\alpha `$ with a ghost $`v`$, which is an enveloping algebra valued, Grassmannian field111 In the $`U(1)`$ case, the introduction of a ghost has been considered by Okuyama .. We define a ghost number by assigning ghost number one to $`v`$ and zero to $`a_i`$ and $`\psi `$. The ghost number introduces a $`Z_2`$ grading, with even quantities commuting and odd quantities anticommuting. In our formalism, the gauge transformations (2.1) and (2.6) are replaced by the following BRST transformations:

$$\begin{array}{cc}\delta _vv=iv^2\hfill & \\ \delta _va_i=_ivi[a_i,v]\hfill & \\ \delta _v\psi =iv\psi .\hfill & \end{array}$$
(3.1)
We also take $`\delta _v`$ to commute with the partial derivatives,

$$[\delta _v,_i]=0.$$
(3.2)
The operator $`\delta _v`$ has ghost number one and obeys a graded Leibniz rule

$$\delta _v(f_1f_2)=(\delta _vf_1)f_2+(1)^{deg(f_1)}f_1(\delta _vf_2),$$
(3.3)
where $`deg(f)`$ gives the ghost number of the expression $`f`$. One can readily check that $`\delta _v`$ is nilpotent on the fields $`a_i`$, $`\psi `$ and $`v`$ and therefore, as a consequence of (3.3), we have

$$\delta _v^2=0.$$
(3.4)
Following the procedure outlined in the previous section, we characterize the SW map as follows. We introduce a matter field $`\mathrm{\Psi }(\psi ,\psi ,\mathrm{},a,a,\mathrm{})`$ and an odd gauge parameter $`\mathrm{\Lambda }(v,v,\mathrm{},a,a,\mathrm{})`$ corresponding to $`\psi `$ and $`v`$ in the commutative theory. $`\mathrm{\Lambda }`$ is linear in the infinitesimal parameter $`v`$ and hence has ghost number one. As before, we require that the SW map respect gauge invariance

$$\delta _\mathrm{\Lambda }\mathrm{\Psi }i\mathrm{\Lambda }\mathrm{\Psi }=\delta _v\mathrm{\Psi }.$$
(3.5)
The consistency condition (2.10) now takes the form

$$\delta _\mathrm{\Lambda }^2\mathrm{\Psi }=\delta _v^2\mathrm{\Psi }=0,$$
(3.6)
and again it yields an equation in $`\mathrm{\Lambda }`$ only. Since

$$0=\delta _\mathrm{\Lambda }^2\mathrm{\Psi }=\delta _\mathrm{\Lambda }(i\mathrm{\Lambda }\mathrm{\Psi })=i\delta _v\mathrm{\Lambda }\mathrm{\Psi }+\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Psi },$$

we can drop $`\mathrm{\Psi }`$ and obtain

$$\delta _v\mathrm{\Lambda }=i\mathrm{\Lambda }\mathrm{\Lambda }.$$
(3.7)
Once the solution of (3.7) is known, one can solve the following equations for $`\mathrm{\Psi }`$ and the gauge field

$$\delta _v\mathrm{\Psi }=i\mathrm{\Lambda }\mathrm{\Psi },\delta _vA_i=_i\mathrm{\Lambda }i\left[A_i\stackrel{}{,}\mathrm{\Lambda }\right].$$
(3.8)
It is natural to expand $`\mathrm{\Lambda }`$ and $`A_i`$ as power series in the deformation parameter $`\theta `$. We indicate the order in $`\theta `$ by an upper index in parentheses

$$\begin{array}{cc}\mathrm{\Lambda }=_{n=0}^{\mathrm{}}\mathrm{\Lambda }^{(n)}=v+_{n=1}^{\mathrm{}}\mathrm{\Lambda }^{(n)}\hfill & \\ A_i=_{n=0}^{\mathrm{}}A_i^{(n)}=a_i+_{n=1}^{\mathrm{}}A_i^{(n)}.\hfill & \end{array}$$
(3.9)
Note that the zeroth order terms are determined by requiring that the SW map reduce to the identity as $`\theta `$ goes to zero. Using this expansion we can rewrite equations (3.7) and (3.8) as

$$\begin{array}{cc}\delta _v\mathrm{\Lambda }^{(n)}i\{v,\mathrm{\Lambda }^{(n)}\}=M^{(n)}\hfill & \\ \delta _vA_i^{(n)}i[v,A_i^{(n)}]=U_i^{(n)},\hfill & \end{array}$$
(3.10)
where, in the first equation, $`M^{(n)}`$ collects all terms of order $`n`$ which do not contain $`\mathrm{\Lambda }^{(n)}`$, and similarly $`U_i^{(n)}`$ collects terms not involving $`A_i^{(n)}`$. We refer to the left hand side of each equation as its homogeneous part, and to $`M^{(n)}`$ and $`U_i^{(n)}`$ as the inhomogeneous terms of (3.10). Note that $`M^{(n)}`$ contains explicit factors of $`\theta `$, originating from the expansion of the Weyl-Moyal product (1.2). If the SW map for $`\mathrm{\Lambda }`$ is known up to order $`(n1)`$, then $`M^{(n)}`$ can be calculated explicitly as a function of $`v`$ and $`a_i`$. On the other hand, $`U_i^{(n)}`$ depends on both $`\mathrm{\Lambda }`$ and $`A_i`$, the former up to order $`n`$ and the latter up to order $`(n1)`$. Still, one can calculate it iteratively as a function of $`v`$ and $`a_i`$.

The structure of the homogeneous portions suggests the introduction of a new operator $`\mathrm{\Delta }`$

$$\mathrm{\Delta }=\{\begin{array}{cc}\delta _vi\{v,\}\hfill & \text{on odd quantities}\hfill \\ \delta _vi[v,]\hfill & \text{on even quantities}.\hfill \end{array}$$
(3.11)
In particular, $`\mathrm{\Delta }`$ acts on $`v`$ and $`a_i`$ as follows

$$\mathrm{\Delta }v=iv^2,\mathrm{\Delta }a_i=_iv.$$
(3.12)
As a consequence of its definition, $`\mathrm{\Delta }`$ is an anti-derivation with ghost-number one. It follows a graded Leibniz rule identical to the one for $`\delta _v`$ (3.3). Another consequence of the definition (3.11) is that $`\mathrm{\Delta }`$ is nilpotent

$$\mathrm{\Delta }^2=0.$$
(3.13)
The action of $`\mathrm{\Delta }`$ on expressions involving $`a_i`$ and its derivatives can also be characterized in geometric terms. Specifically, $`\mathrm{\Delta }`$ differs from $`\delta _v`$ in that it removes the covariant part of the gauge transformation. Therefore, $`\mathrm{\Delta }`$ acting on any covariant expression will give zero. For instance, if one constructs the field-strength, $`F_{ij}_ia_j_ja_ii[a_i,a_j]`$, one finds by explicit calculation

$$\mathrm{\Delta }F_{ij}=0.$$
(3.14)
It can also be checked that the covariant derivative, $`D_i=_ii[a_i,]`$, commutes with $`\mathrm{\Delta }`$

$$[\mathrm{\Delta },D_i]=0.$$
(3.15)
In terms of $`\mathrm{\Delta }`$ the equations (3.10) take the form

$$\begin{array}{cc}\mathrm{\Delta }\mathrm{\Lambda }^{(n)}=M^{(n)}\hfill & \\ \mathrm{\Delta }A_i^{(n)}=U_i^{(n)}.\hfill & \end{array}$$
(3.16)
In the next section, we will provide a method for solving these equations. Also note that since $`\mathrm{\Delta }^2=0`$, it must be true that

$$\begin{array}{cc}\mathrm{\Delta }M^{(n)}=0\hfill & \\ \mathrm{\Delta }U_i^{(n)}=0.\hfill & \end{array}$$
(3.17)
Indeed one should verify that (3.17) holds order by order. If (3.17) did not hold, this would signal an inconsistency in the SW map.

Many authors have commented on the ambiguities of the SW map . At any particular order, the ambiguities can be seen as an invariance of (3.16) when $`\mathrm{\Lambda }^{(n)}`$ is changed by an amount $`\mathrm{\Delta }S^{(n)}`$

$$\mathrm{\Lambda }^{(n)}\mathrm{\Lambda }^{(n)}+\mathrm{\Delta }S^{(n)},$$
(3.18)
which follows from the fact that $`\mathrm{\Delta }`$ is nilpotent. Then the corresponding change in the potential is

$$A_i^{(n)}A_i^{(n)}+D_iS^{(n)}.$$
(3.19)
This follows from the fact that the equation of order $`n`$ for the gauge field is always of the form

$$\mathrm{\Delta }A_i^{(n)}=D_i\mathrm{\Lambda }^{(n)}+\mathrm{},$$
(3.20)
where the ellipsis denotes terms which are explicitly $`\theta `$-dependent. Notice that (3.19) is a consequence of the fact that the coboundary operator $`\mathrm{\Delta }`$ commutes with the covariant derivative $`D_i`$. The ambiguities at order $`n`$ also affect the solutions at higher order.

These ambiguities can also be understood as an invariance of (3.7) and (3.8) under the following transformations

$`\mathrm{\Lambda }`$ $``$ $`G^1\mathrm{\Lambda }G+iG^1\delta _vG`$
$`A_i`$ $``$ $`G^1A_iG+iG^1_iG`$ (3.21)
$`\mathrm{\Psi }`$ $``$ $`G^1\mathrm{\Psi },`$
where all products are star products and $`G`$ is an arbitrary element of the enveloping algebra with ghost number zero. Notice that $`G`$ should also be unitary if we require that $`\mathrm{\Lambda }`$ and $`A_i`$ remain real.

To compare (3.21) with (3.18) and (3.19) it is useful to introduce the operators

$`\widehat{D}_i`$ $``$ $`_ii[A_i\stackrel{}{,}]`$ (3.22)
$`\widehat{\mathrm{\Delta }}`$ $``$ $`\{\begin{array}{cc}\delta _vi[\mathrm{\Lambda }\stackrel{}{,}]\hfill & \text{for even quantities}\hfill \\ \delta _vi\{\mathrm{\Lambda }\stackrel{}{,}\}\hfill & \text{for odd quantities},\hfill \end{array}`$ (3.25)
which satisfy

$$\widehat{\mathrm{\Delta }}^2=0,[\widehat{D}_i,\widehat{\mathrm{\Delta }}]=0$$
(3.26)
and which reduce to $`D_i`$ and $`\mathrm{\Delta }`$ in the limit of vanishing $`\theta `$. Then (3.21) can be rewritten as

$`\mathrm{\Lambda }`$ $``$ $`\mathrm{\Lambda }+iG^1\widehat{\mathrm{\Delta }}G`$ (3.27)
$`A_i`$ $``$ $`A_i+iG^1\widehat{D}_iG.`$ (3.28)
To recover (3.18) and (3.19) we set

$$G=1iS^{(n)},$$
(3.29)
and take (3.27) and (3.28) at order $`n`$. These ambiguities are of the form of a gauge transformation. Notice that in the particular case, $`\mathrm{\Delta }S^{(n)}=0`$, $`A^{(n)}`$ is modified while $`\mathrm{\Lambda }^{(n)}`$ is unaffected.

In  it has been observed that there are also other kinds of ambiguities, which don’t have the form of a gauge transformation, but are of a covariant type. To see this, we rewrite the SW equation for $`A`$ (3.8) using $`\widehat{\mathrm{\Delta }}`$

$$\widehat{\mathrm{\Delta }}A_i=_i\mathrm{\Lambda }.$$
(3.30)
It is then possible to add to the gauge potential a quantity $`S_i`$,

$$A_iA_i+S_i,\widehat{\mathrm{\Delta }}S_i=0,$$
(3.31)
while keeping $`\mathrm{\Lambda }`$ unchanged.

## 4 The Homotopy Operator

For simplicity, we begin by considering in detail the SW map for the case of the gauge parameter $`\mathrm{\Lambda }`$. Much of what we say actually applies to the other cases as well with minor modifications.

In the previous section, we have seen that order by order in an expansion in $`\theta `$, the SW map has the form:

$$\mathrm{\Delta }\mathrm{\Lambda }^{(n)}=M^{(n)},$$
(4.1)
where $`M^{(n)}`$ depends only on $`\mathrm{\Lambda }^{(i)}`$ with $`i<n`$. Clearly, if one could invert $`\mathrm{\Delta }`$ somehow, we could solve for $`\mathrm{\Lambda }^{(n)}`$. But $`\mathrm{\Delta }`$ is obviously not invertible, as $`\mathrm{\Delta }^2=0`$. In particular, the solutions of (4.1) are not unique, since if $`\mathrm{\Lambda }^{(n)}`$ is a solution so is $`\mathrm{\Lambda }^{(n)}+\mathrm{\Delta }S^{(n)}`$ for any $`S^{(n)}`$ of ghost number zero222These are precisely the ambiguities in the SW map that were first discussed in , where our operator $`\mathrm{\Delta }`$ was called $`\widehat{\delta ^{}}`$.. That is, $`\mathrm{\Delta }`$ acts like a coboundary operator in a cohomology theory, and the solutions that we are looking for are actually cohomology classes of solutions, unique only up to the addition of $`\mathrm{\Delta }`$-exact terms. The formal existence of the SW map is then equivalent to the statement that the cycle $`M^{(n)}`$ is actually $`\mathrm{\Delta }`$-exact for all $`n`$. Since we know that $`\mathrm{\Delta }^2=0`$, this fact would follow as a corollary of the stronger statement that there is no non-trivial $`\mathrm{\Delta }`$-cohomology in ghost number two. In other words, there are no $`\mathrm{\Delta }`$-closed, order $`n`$ polynomials with ghost number two which are not also $`\mathrm{\Delta }`$-exact. To prove this stronger claim, we could proceed as follows. Suppose that we could construct an operator $`K`$ such that

$$K\mathrm{\Delta }+\mathrm{\Delta }K=1.$$
(4.2)
Clearly, $`K`$ must reduce ghost number by one, and therefore must be odd. Consider its action on a cycle $`M`$, (so $`\mathrm{\Delta }M=0`$)

$$(K\mathrm{\Delta }+\mathrm{\Delta }K)M=\mathrm{\Delta }KM=M.$$
(4.3)
Therefore, $`M=\mathrm{\Delta }\mathrm{\Lambda }`$, with $`\mathrm{\Lambda }=KM`$, which not only shows that $`M`$ is exact, but also computes explicitly a solution to the SW map. We note that this method of solution is nearly identical to the method used by Stora and Zumino  to solve the Wess-Zumino consistency conditions for nonabelian anomalies. In fact, it was the parallels between these problems that motivated the current approach.

We now proceed to construct $`K`$. First we notice that $`M^{(n)}`$ depends on $`v`$ only through its derivative $`_iv`$, as one can see by looking at the explicit expressions. The same is true for $`U_i^{(n)}`$ since it depends on $`v`$ only through $`\mathrm{\Lambda }`$. It is convenient to define

$$b_i=_iv,$$
(4.4)
so that $`M`$ and $`U_i`$ can all be expressed as functions of $`a_i`$, $`b_i`$ and their derivatives only. Furthermore, we rewrite $`M^{(n)}`$ solely in terms of covariant derivatives, rather than ordinary ones. After these replacements, we may consider $`M^{(n)}`$ an element of the algebra generated by $`a_i`$, $`b_i`$, and $`D_i`$. As explained in the next section this algebra is not free, but for the moment we ignore this issue. The action of the operator $`\mathrm{\Delta }`$ takes on a particularly simple form in terms of these variables:

$$\mathrm{\Delta }a_i=b_i,\mathrm{\Delta }b_i=0,[\mathrm{\Delta },D_i]=0.$$
(4.5)
Let us first define an odd operator $`L`$, which obeys the super Leibniz rule, and satisfies

$$La_i=0,Lb_i=a_i,[L,D_i]=0.$$
(4.6)
Acting on either $`a`$ or $`b`$, we have $`L\mathrm{\Delta }+\mathrm{\Delta }L=1`$, but this is no longer true acting on monomials of higher order. The solution is to define

$$K=D^1L,$$
(4.7)
where $`D^1`$ is a linear operator which when acting on a monomial of total order $`d`$ in $`a`$ and $`b`$ multiplies that monomial by $`1/d`$. In can be proven that $`K`$ defined in this way satisfies (4.2) when acting on monomials of degree greater than or equal to one. Since $`L`$ satisfies the Leibniz rule, we see that $`L^2=0`$, by considering its action on the generators of the algebra (4.6). It then follows that

$$K^2=0.$$
(4.8)
Notice that this prescription requires that we rewrite any expression involving ordinary derivatives in terms of covariant derivatives and gauge fields only.

## 5 Constraints

We have so far only considered the free algebra, generated by $`a_i`$ , $`b_i`$ and $`D_i`$, where the construction of $`K`$ was relatively simple. To show that our algebra is not free consider the following

$$\begin{array}{cc}\mathrm{\Delta }F_{ij}\hfill & =\mathrm{\Delta }\left(D_ia_jD_ja_i+i[a_i,a_j]\right).\hfill \\ & =D_ib_jD_jb_i+i[b_i,a_j]+i[a_i,b_j].\hfill \end{array}$$
(5.1)
As an element of the free algebra, the right hand side is not zero, but according to (3.14), the left hand side should be. The problem becomes more serious when one rewrites $`M^{(n)}`$ in terms of the elements of the free algebra. Beyond first order, one finds that $`\mathrm{\Delta }M^{(n)}`$ is no longer zero in general, but vanishes only by using the following constraints

$$[F_{ij},]i[D_i,D_j]()=0,\mathrm{\Delta }F_{ij}=0.$$
(5.2)
If $`\mathrm{\Delta }M^{(n)}`$ is not zero identically, $`K`$ no longer inverts $`\mathrm{\Delta }`$ when acting on $`M^{(n)}`$, and we no longer have a method for solving (3.16) for $`\mathrm{\Lambda }^{(n)}`$. The origin of the constraints can be traced to the fact that partial derivatives commute

$$_i_j_j_i=0,_ib_j_jb_i=0,$$
(5.3)
since $`b_i=_iv`$. This is no longer manifest in our algebra. In fact, written in terms of covariant derivatives, (5.3) becomes (5.2). There seems to be no way to eliminate these constraints since $`K`$ is not defined on $`v`$, but only on $`b_i=_iv`$. One might expect that at higher orders one would have to use additional constraints to verify that $`\mathrm{\Delta }M^{(n)}`$ vanishes, but this is not the case. For example, when one rewrites

$$_i_kb_j_j_kb_i=0$$
(5.4)
in terms of covariant derivatives, the resulting expression is not an independent constraint, but can be written in terms of the two fundamental ones (5.2).

The reason why $`\mathrm{\Delta }M^{(n)}`$ is not zero in general is because the existence of the constraints allows us to write $`M^{(n)}`$ in terms of the algebra elements in many different ways. Our goal will then be to define a procedure for writing $`M^{(n)}`$ in terms of algebra elements so that $`\mathrm{\Delta }M^{(n)}=0`$, identically. We will describe two procedures.

The first is the method used in  to calculate some low order terms of the SW map. One begins by obtaining an expression for $`M^{(n)}`$ in terms of the algebra elements. Generically, $`\mathrm{\Delta }M^{(n)}`$ will be proportional to the constraints. At low orders, once $`\mathrm{\Delta }M^{(n)}`$ is calculated, it is easy to guess an expression $`m^{(n)}`$, which is proportional to the constraints, such that the combination $`M^{(n)}+m^{(n)}`$ is annihilated by $`\mathrm{\Delta }`$. Acting $`K`$ on this new combination then gives the solution $`\mathrm{\Lambda }^{(n)}`$. We believe this guessing method can be formalized, but at higher orders the second procedure which we will now describe seems to be more systematic.

First we introduce a new element of the algebra, $`f_{ij}`$, which is annihilated by all the operators defined in previous sections

$$\mathrm{\Delta }f_{ij}=Lf_{ij}=0.$$
(5.5)
We also introduce a new constraint

$$f_{ij}F_{ij}=0,$$
(5.6)
where $`F_{ij}`$ is considered a function of $`D_i`$ and $`a_i`$. We want to show that using this enlarged algebra and the constraints we can rewrite $`M^{(n)}`$ so that it has the following dependence

$$M^{(n)}=M^{(n)}(a,b,(D^ka)_s,(D^lb)_s,D^hf),$$
(5.7)
where the subscript $`s`$ indicates that all the indices within the parentheses should be totally symmetrized. It would then follow that $`\mathrm{\Delta }M^{(n)}`$ depends on the same variables. Since it is impossible that $`\mathrm{\Delta }M^{(n)}`$ contains any term antisymmetric in the indices of $`Da`$ or $`Db`$, the constraints (5.6) and (5.2) cannot be generated. However, we may find that $`\mathrm{\Delta }M^{(n)}`$ is proportional to the following constraints

$$[f_{ij},]i[D_i,D_j]()=0,D_if_{jk}+D_jf_{ki}+D_kf_{ij}=0.$$
(5.8)
Since these constraints commute with the action of both $`K`$ and $`\mathrm{\Delta }`$, if we add to $`M^{(n)}`$ a term proportional to (5.8), our result for $`\mathrm{\Lambda }^{(n)}=KM^{(n)}`$ is unchanged. To show that we can actually write $`M`$ in the form suggested above, we begin with an expression for $`M^{(n)}`$ as found by expanding the star product

$$M^{(n)}=M^{(n)}(a,(^k)_sa,(^l)_sv),$$
(5.9)
where we choose to explicitly write the derivatives in symmetric form. By replacing $`()D()+i[a,]`$, and $`vb`$ the expression takes the form

$$M^{(n)}=M^{(n)}(a,b,(D^k)_sa,(D^lb)_s).$$
(5.10)
The difference $`(D^ka)_sD^ka`$ contains terms that are proportional to the antisymmetric parts of $`DD`$ or $`Da`$. But using the constraints we can make the following substitutions

$$[D_i,D_j]()i[f_{ij},],D_ia_jD_ja_if_{ij}i[a_i,a_j].$$
(5.11)
This must be done recursively since the commutator term involving $`a`$’s above may again be acted on by $`D`$’s. But at each step, the number of possible $`D`$’s acting on $`a`$ is reduced by one. After carrying out this procedure $`M^{(n)}`$ will have the form (5.7).

## Acknowledgments

We are very indebted to Julius Wess for explaining to us some of his recent work . We are also grateful to Raymond Stora for the valuable comments we received from him. This paper is a modified version of the lecture that B. Zumino gave at the conference “2001: A Spacetime Odyssey”, Michigan Center for Theoretical Phyisics, Ann Arbor, May 21-25. This work was supported in part by the Director, Office of Science, Office of High Energy and Nuclear Physics, Division of High Energy Physics of the U.S. Department of Energy under Contract DE-AC03-76SF00098 and in part by the National Science Foundation under grant PHY-95-14797. B.L.C. is supported by the DFG (Deutsche Forschungsgemeinschaft) under grant CE 50/1-2.