text:
# Cellular Models for River Networks

## I Introduction

Through experimental studies, it has become evident in the past few years that the geometrical and topological structures of river basins are characterized by the absence of a single well-defined length scale. This is reflected in the appearance of power laws in the distribution of several quantities, chiefly total contributing area at a point  and stream lengths , and by the clear experimental assessment of scaling properties (yielding either self-similarity or self-affinity ) for many geometrical descriptors of the river basin . The discovery of the general underlying mechanisms yielding scale-free features is the present theoretical challenge.

The network associated to a given natural terrain pertaining to a river basin can be experimentally analyzed by using the so-called Digital Elevation Map (DEM) technique  which allows to determine the average height of areas (pixels) of the order $`10^2\text{Km}^2`$. Thus a fluvial basin is represented in a objective manner often over 4 log scales of linear size. Lower bounds are imposed by channel initiation processes at O($`10100`$) m. Crossovers of geological nature provide altered aggregation processes and thus an upper cutoff, usually beyond scales of O($`10^5`$ \- $`10^6`$) m. Thus the observational evidence yields a much reliable framework over many scales for comparison with dynamical models aimed at the origin of scale-free features.

Much interest has been recently attracted by landscape evolution models. Chief among those are the detailed deterministic models which address the description of the detailed dynamics acting on the landscapes . The reductionist approach, where a precise description of the details of the dynamics is sought, is successful, and much interesting, in the pursuit of the description and the classification of landforms. Nevertheless, as standard in critical phenomena, the mechanism producing scale-free structures is expected to depend only on a few key features common to all the networks rather than on the details of the particular system under study. Hence in this work, centered on the dynamic origin of fractal river networks, we follow a nonreductionist approach based on the simplest possible, parameter-free models capable of allowing the emergence of complexity.

A flowrate unit is associated with each pixel and the flow contributing to any pixel follows the steepest descent path through drainage directions whose collection defines the planar structure under consideration. The resulting network is therefore the two-dimensional projection of the three-dimensional tree-like structure of the steepest descent paths draining a given basin. The planar patterns of network aggregation are obtained by employing the cellular model for the evolution of a fluvial landscape originally introduced by  and further studied in  It is aimed at describing in a crucially simple manner the sole fluvial component of landscape evolution. Although such component must be coupled to other - chiefly hillslope - transport processes to yield a comprehensive dynamical description , it rules the planar imprinting of the network. Hence the detailed study of the model is deemed significant.

Starting from a three dimensional landscape, evolution occurs according to a threshold dynamics similar to the one proposed in self organized critical (SOC) models . The main idea of the erosion dynamics is that whenever the local shear stress exceeds a given threshold, erosion starts an ‘avalanche’ and a related rearrangement of the network patterns takes place. The model of self-organizing fluvial structures may be seen as a modification of the sandpile model developed as a paradigm of the dynamics of open, dissipative systems with many degrees of freedom. It may be thought of as belonging to the set of models in which the threshold for activity, say $`\tau _c`$, rather than being a constant value, depends on non-local properties of the self-organizing structure. In the fluvial case, the non-local character of the threshold value follows from the fact that the threshold at the arbitrary site equals a shear stress, i.e. $`\tau h\sqrt{a}`$, where $`h`$ is the local landscape elevation and $`a`$ is total contributing area surrogating total flow collected from a distributed rainfall event. As such, the exceedence of $`\tau _c`$ depends not only on local conditions (i.e., a critical value of $`h`$), but also on non-local conditions defined by the contributing area $`a`$ computed through drainage directions, i.e. it depends on the entire state of the system which is self-organizing. Notice that the physical rationale for the nonlocal dependence lies in the fact that the system is open, i.e. injected from outside, allowing flow rates to be proportional to total contributing drainage area. The long-range nature of the threshold dynamics tends to hide from the observer the temporal fluctuations which take place in the evolutionary time scale. In this sense the above model was classified  as one of spatial self-organized criticality.

Whether or not river self-organization qualifies as a more general framework of self-organized criticality remains to be seen. If SOC must necessarily refer to the occurrence of a critical state in the sense of critical phenomena, where a small local perturbation can cause a significant change in the configuration of the whole system and thus the system shows both spatial and temporal scaling, then the time dynamics should be specifically considered. One way to do this is through the oscillation of the threshold in time, i.e. $`\tau _c(t)`$, simulating climatic fluctuations (see ), through which indeed temporal evolution appear, or through perturbations of random location and strength in the evolution of the landscape. This is also true, as we will discuss later, if the landscape-forming rainfall events are described as nonuniform in space, leading to patches of activity randomly scattered spatially (in such a case the outflow response of the system becomes a $`1/f`$ signal). However, regardless of any additional features, we believe that the central scope of SOC is the dynamic explanation of the growth of fractal structures of the type appearing in nature, i.e. the physics of fractals. As such we feel that our classification of the model as a particular case of SOC is a suitable one regardless of the description of the embedded temporal activity because the system always reaches a fractal state. Moreover, questioning on this basis  the self-organized critical nature of the model by Ref.  is irrelevant because it has been shown on thermodynamics grounds that scaling properties of energy and entropy yield limit states which, depending on the constraints, are temporally frozen or active . Furthermore, ref.  shows that optimal states like the ones dynamically accessed by the above model may exist in temporally active states precisely at the edge of a chaotic behaviour.

A question, indeed more interesting than the semantics of SOC, is whether the constraints in the model may be relaxed to produce a ’hot’ fluvial landscape more closely resembling an ordinary sandpile. This question is addressed in  and in more detatil in .

Our main goal is twofold. On one hand we will extend previous investigations both in accuracy and in statistics by performing simulations at much larger scales. On the other hand, we will consider important issues such as the effect of the boundaries and of the initial conditions bearing much significance on geologic influences. In particular we will show that both the aforementioned effects play an important role in the results previously obtained. We also study the effects of disorder, say through the presence of small, uncorrelated inhomogeneities in the initial conditions, in particular with regard to the robustness to single/multiple outlet arrangements.

The paper is organized as follows. In the first section 1 the model is recalled. Section 2 presents the results with emphasis on scaling analyses, while the following section focuses on the important effects of quenched pinning on the system. A set of conclusions closes then the paper.

## II The model

We consider a lattice model of a real landscape. Let $`h_x`$ be the height of the landscape associated to every site $`x`$ of a square lattice of size $`L\times L`$. The lattice is tilted at an angle $`\theta `$ with respect to a given axis to mimic the effects of gravity. Two possibilities will be analyzed:

* All the sites on the lowest side (kept at height $`h=0`$) are possible outlets (i.e. the multiple outlet arrangement) of an ensemble of rivers which are competing to drain the whole $`L\times L`$ basin;
* Only one site is kept at $`h=0`$ and it is the outlet of a single river in the $`L\times L`$ basin.

In addition, for both the above cases, two types of initial conditions will be considered: (a) a regular initial landscape, e.g. flat, and (b) an irregular surface obtained by superposing to a smooth sloping surface a suitable noise.

Each site collects an unit amount of water from a distributed injection (here a constant rainfall rate as in the original approach) in addition to the flow which drains into it from the upstream sites. A unit of water mass is assigned to each pixel of drainage area so that the total area drained into a site is also a measure of the total water mass collected at that site. From each site water flows to one of the eight sites, four nearest neighbors and four next nearest neighbors, having the lowest height (i.e. the steepest descent path). We shall indicate all these eight neighbouring sites as nearest-neighbours ($`nn`$). This construction allows the assignment of drainage directions to an arbitrary landscape. The drained area $`a_x`$ is associated to each site $`x`$ according to the equation

$`a_x`$ $`=`$ $`{\displaystyle \underset{y(x)}{}}a_y+1`$ (1)
where the sum runs over the subset $`y(x)nn(x)`$ of neighbor sites whose area is actually drained by $`x`$. The second term in eq.(1) represents the uniform injection.

Also, the (up)stream-length $`l_x`$ from site $`x`$ to the source is computed according to the following procedure. At a given site $`x`$ the areas of all $`nn(x)`$ of that sites are checked. Following the ordinary meaning of downstream and upstream sites (i.e. downstream is the site one finds following the river to the outlet, upstream is the site following which one reach the source from which the largest incoming river enters the site.) Following Ref., the $`nn`$ with largest value leads to the outlet and is defined to be a downstream site. The $`nn`$ with the second largest value indicates the longest path toward the source and is defined to be the upstream site. The sum of all the upstream sites from site $`x`$ to the source is $`l_x`$. The downstream length could be defined through an analogou procedure. Experimental measures are available for both $`a_x`$ and $`l_x`$ .

The time evolution of the model follows the following steps:

* The shear stress $`\tau _x`$ acting at every site is computed according to

$$\tau _x=\mathrm{\Delta }h_x\sqrt{a_x}$$
(2)
where $`\mathrm{\Delta }h_x`$ is the local gradient along the drainage direction.
* If the shear stress at a site exceeds a threshold value, $`\tau _c`$, then the corresponding height $`h_x`$ is reduced (i.e. by erosion) in order to decrease the local gradient. The shear stress is set just at the threshold value. This produces a rearrangement of the network followed by a reupdating of the whole pattern as in step 1.
* When all sites have shear stress below threshold the system is in a dynamically steady state. Since this situation is not necessarily the most stable, a perturbation is applied to the network with the aim of increasing the stability of a new steady state. A site is thus chosen at random and its height is increased in such a way that no lakes, i.e. sites whose height is lower than that of their eight neighbours, are formed. Steps 1 and 2 then follow as before.

After a suitable number of the perturbations (step 3), the system reaches a steady state which is unsensitive to further perturbations and where all statistics of the networks are stable. This resulting state is scale-free, i.e. it is characterized by power-law distributions of the physical quantities of interest.

## III Results

### A Landscape evolutions

Our numerical calculations were carried out on a bidimensional square lattice (where each site has eight nearest-neighbours) for sizes up to $`L=200`$ with reflecting boundary conditions in the direction transversal to the flow and open boundary condition in the parallel one. We considered the following initial conditions:

Model A A comb-like structure with a single outlet This was the situation originally studied in  and our results are in agreement with theirs.

Model B An inclined plane with all sites at the bottom of the plane allowed to be possible outlets. This choice was selected with the aim of investigating the differences arising when arranging the boundary conditions with multiple outlets versus single outlet. The former allows for competition for drainage area among rivers.

Model C,D The two previously considered situations with the addition of a random, uncorrelated noise (whose strenght - i.e. variance - is less than $`10\%`$ of the average height). That is, on top of the height computed according to the rules of model A and B respectively (a comb-like lattice and an inclined plane), we added a random $`Dh`$ that is extracted in the interval $`[<h>/10,<h>/10]`$ where $`<h>`$ represents the mean altitude of the landscape.

An average over a few (up to five) configurations was taken. This choice, especially when coupled to large sizes of the system, proves sufficient for a statistical descriptions sought in view of the self-averaging nature of the random perturbation.

In Figure 1 typical landscapes sculpted by the above dynamical process and the corresponding networks drawn through the steepest descent construction are shown for models A,B. The same picture for models C,D is shown in Figure 2.

Two features can be grasped from these pictures. First, in the results of both model A and B there is a strong memory of the initial configuration despite the fact that the dynamics of the erosion process was somewhat expected to be sufficiently strong to soon lose the imprinting of its initial condition. Secondly, the single outlet restriction imposed in model A appears to be a severe constraint because it increasingly affects the wandering of the main river towards the lowest part of the basin. Our results suggest that this is indeed the case for flat initial conditions (A and B) while for noisy initial conditions boundary effects are of lesser importance.

### B Area and length exponents

Let us define $`P(a,L)`$ and $`\mathrm{\Pi }(l,L)`$ as the exceedence (cumulative) probability distributions of the drainage area $`a`$ and stream length $`l`$ respectively arising in a domain of linear size $`L`$. The following scaling forms are expected to hold :

$`P(a,L)`$ $`=`$ $`a^{1\tau }F({\displaystyle \frac{a}{L^{1+H}}})`$ (3)
$`\mathrm{\Pi }(l,L)`$ $`=`$ $`l^{1\gamma }G({\displaystyle \frac{l}{L^{d_l}}})`$ (4)
Here $`H`$ is the Hurst exponent  and $`d_l`$ is the stream-length (or chemical distance) fractal exponent.

As it was already noted , for self-affine river networks ($`H<1`$, $`d_l=1`$) the scaling relations relate all exponents in terms of $`H`$. For self-similar river networks ($`H=1`$, $`d_l>1`$), the same happens in terms of $`d_l`$.

Another important indicator of basin morphology is the relation between the mean total contributing area $`a`$ and the length of the main stream $`l_{\text{max}}L^{d_l}`$ , which is commonly known as Hack’s law :

$`a`$ $``$ $`l_{max}^{1/h}L^{d_l/h}`$ (5)
The related exponent has been studied in all simulations. A summary of the scaling relations between the various exponents involved is reported in Table II.

Experimental values of $`\tau `$ and $`\gamma `$ are available from earlier analyses of DTMs from basins of different size, geology, exposed lithology, climate and vegetation . It was observed that, while a majority of basins tend to seemingly universal values $`\tau =1.43\pm 0.02`$ and $`\gamma =1.8\pm 0.1`$, exceptions are observed where altered values are observed although always in a concerted manner. Since it was suggested  that scaling laws for river networks are related, e.g. $`\gamma =1+(\tau 1)/h`$, it was concluded there that no universal exponents are expected in nature. Rather, the roles of geology and tectonics concert a coordinated scaling structure which strives for fractality yet adapted to its geological environment. The results of the model described here, revisited in the above light, conform to this view.

The results for the four models A,B,C and D for the area distributions are shown in Figure 3. It is apparent that, due to the pathological initial conditions, the scaling behaviour for models A and B is somewhat more noisy than for models C and D. Figure 4 contains the collapse plot for all the cases. Figure 5 shows the stream-length distribution for the four models. For this picture the same remarks of Figure 3 apply. In Figure 6 we show the collapse plot corresponding for the stream-length distributions.

A summary of the scaling exponents obtained is included in Table II where we observe a consistent picture of related scaling exponents as theoretically expected - see Table II.

### C Energy Dissipation and Optimal Channel Networks

During the evolution of the landscape we also monitored the change in total energy dissipation of the system, defined as $`E=_xa_x^{0.5}`$ (where $`x`$ spans all sites of the lattice) . The reason of this name come from the computation at any site of the gravitational energy lost by the falling of the water. In any point $`x`$ one can expect a gravitational energy loss of the order of $`a_x\mathrm{\Delta }h_x`$ where $`\mathrm{\Delta }h_x`$ represents the local gradient along the drainage direction. By using the observed scaling $`\mathrm{\Delta }h_xa_{x}^{}{}_{}{}^{0.5}`$ one obtains the above formula. The interest in this quantity comes from the fact that an extensively class of models known as Optimal Channel Network (OCN) models, assume this quantity is minimized by natural landscape evolution. By using this principle OCN describe evolution from random spanning graphs to network more similar to the real ones. It is interesting to note that in this model where no hypothesis is made on $`E`$, we still observe an almost monotonical decrease of $`E`$ associated with the dynamical evolutions, and a stabilization on different plateaus of values of $`E`$. The actual figures for a sample $`30`$ are as follows: $`E`$ starts from an initial value of $`7600`$ and decreases towards a plateau of $`68006700`$ where this monotonic decrease becomes slower ($`1\%`$ decrease in 50000 steps). This behavior, also observed in other models , bears important consequences in the light of the suggested connection of fractality and optimality

## IV Geological constraints and Quenched Random Pinning

This section presents a detailed study on the effects on landcape evolution of quenched randomness, simulated by a random choice of sites unable to evolve regardless of the threshold value developed therein. It is found that this form of disorder tends to favour aggregation patterns characterized by values of $`\tau =1.43\pm 0.02`$ for both models, say, A and B (i.e. with single outlets or open boundary conditions). This results suggests that the origin of the recurrent values observed in nature could be related to the ubiquity of heterogeneity in surface properties characterizing locally the critical shear threshold.

Within the river basin morphological and geological constraints play a definite role in the dynamical evolution of landforms. The effects of quenched constraints, simulating any heterogeneity is the distribution of surficial properties affecting erosion properties, is to favour some sites for the flowpaths, thus excluding other sites from the capture of the developing network which ultimately shapes the evolution process. In order to mimic such effects we analyzed the effect of a random pinning of a small region of the total surface (typically $`59\%`$) where the evolution is frozen, that is, the height is pinned to its initial value. We find that this constraint tends to favour aggregation even in the presence of random initial conditions.

As regards the effect of the pinning, figure 7 (on the left) shows a sample whose dimension is $`100`$ with multiple outlet and random initial noise (model D). Figure 7 (on the right) shows the same configuration (evolved from the same initial conditions) but with a $`5\%`$ dilution pinning. It is evident that some of the smaller streams on the left have increased their size thus leading to a bigger aggregation.

In this case we found $`\tau =1.43\pm 0.02`$. Moreover, all the other exponents verify the correct scalings predicted in Table I. Purely for comparison purposes with the prevous results, we also report the plot of the area distribution in Figure 8 on the left. We also found that model D reproduces the same results, i.e. $`\tau =1.43\pm 0.04`$ and $`\gamma =1.60\pm 0.04`$. We are confident that this result at least for model C is quite robust with respect to changes of the pinning dilution. In the case of a $`9\%`$ dilution for a smaller number of simulation we found quite similar results $`\tau =1.44\pm 0.05`$ and a cumulative plot of $`P(a^{}>a,L)`$ is shown on the right part of Fig.8.

This result suggests that the origin of recurrent values observed in nature could indeed be related to the ubiquity of geological and morphological constraints in the surface properties locally characterizing the critical shear stress.

## V Conclusions

In this paper we revisited the model originally introduced in  which we extended both in accuracy and goals. Specifically, we analyzed the stability of the universality class of the original model with respect to the initial conditions and to the change from single to multiple outlets. We found that if one starts with structured initial conditions, critical exponents are sensible to a change from single to multiple outlets. On the other hand, upon starting from disordered initial conditions, we found critical exponents belonging to a new class which appears to be robust to the change from single to multiple outlet. Thus this simple model, under controlled conditions, yields somewhat different yet internally consistent scale-free fluvial landforms depending on the dominant conditions affecting evolution.

The above results conform to the experimental observation  suggesting that the relevant scaling exponents for river networks are not universal. Rather, the fractal nature of river networks adjusts to the constraints imposed by the geological environment in a coordinated manner. It is interesting to observe that the final state of all simulation yields indeed fractal structures, as observed in nature, though characterized by different aggregation properties. The exponents characaterizing these different aggregates, nevertheless follow in a rather good agreement the scaling relations shown in Table I.

We suggest that the lack of robustness in the value of the scaling exponents with respect to boundary and initial conditions is related to the non-local character of the shear-based threshold, differently from what is observed in classical sandpile models of self-organized criticality . The intrinsic interest of the different aggregation properties of the steady states of the dynamics is related to their optimality with respect of total energy dissipation . In fact, depending on external conditions, the same dynamical process may indeed get trapped in steady state configurations yielding local minima of the total energy dissipation functional, in what we may define as a feasible optimality process .

Finally, we have found that quenched disorder, modeled by random pinning, has a profound effect on the robustness of the resulting planar patterns by favouring aggregation and by locking the planar landforms into modes quite similar to the ones observed in nature.

The remarkable success obtained by such a simple model in enlightnening some crucial features of the real basins is promising for a future success in a general characterization of the dynamics of fractal growth.

###### Acknowledgements.

Enlightening discussions with J. Banavar, F. Colaiori A. Flammini, A. Maritan, and A. Rinaldo are acknowledged. Suggestion by A. Giacometti have been particularly helpful. Author acknowledge support from EU contract FMRXCT980183.


format: text

text:
# Large N limit of Calogero–Moser models and Conformal Field Theories

## Abstract

We discuss the large $`N`$ limit of Calogero–Moser models for the classical infinite families of simple Lie algebras $`A_N`$, $`B_N`$, $`C_N`$ and $`D_N`$. We show that the limit defines two different Conformal Field Theories with central charge $`c>1`$. The value of $`c`$ and the dimension of the primary field are dictated by the underlying algebraic symmetries of the model.

preprint:    INFNCA-TH00018
IFUM-670-FT

The classical and quantum integrable model of $`N`$ interacting bodies in one spatial dimension introduced by Calogero  has been studied intensively and still attracts much attention. Surprisingly, the model and its various generalizations  have been found relevant to describe properties of theories ranging from spin chains, quantum Hall effect, two dimensional gravity and random matrices. Recently, in the high energy and gravity context the interest was renewed since the proposal  that the model (in a supersymmetric extension) could be relevant to provide a microscopic description of the extremal Reissner–Nordström black hole.

The model is described by the Hamiltonian $`H=_{i=1}^Np_i^2/2+U(q)`$, where $`p=(p_1,\mathrm{},p_N)`$ and $`q=(q_1,\mathrm{},q_N)`$ are momentum and coordinate operators in the $`N`$ dimensional space $``$ with $`p_l=i/q_l`$ (particle masses and the Planck constant are set equal to unity). In the most general formulation of the model the two–body potential $`U(q)`$ is defined by means of the subsystem $`R_+`$ of positive roots $`\{\alpha \}`$ of a root system defining a Coxeter group.

We will be concerned with the Calogero–Moser (CM) case,

$$U(q)=\underset{\alpha R_+}{}g_\alpha ^2q_\alpha ^2+\frac{\omega ^2}{2}\underset{\alpha R_+}{}q_\alpha ^2,$$
(1)
where $`q_\alpha =(q,\alpha )`$ and the couplings obey the stability condition

$$g_\alpha ^2\frac{|\alpha |^2}{8},$$
(2)
$`|\alpha |`$ being the length of the root $`\alpha `$. Since we are interested in the large $`N`$ limit we consider only the infinite (classical) root systems $`A_N`$, $`B_N`$, $`C_N`$ and $`D_N`$, corresponding, respectively, to the simple Lie algebras $`su(N+1)`$, $`so(N+1)`$, $`sp(N)`$ and $`so(2N)`$. The only, nontrivial, nonreduced case, $`BC_N`$, is equivalent to $`A_{2N}`$ with proper symmetry conditions preserved by the time evolution . Every other reducible case corresponds to the union of noninteracting systems.

We aim to study the $`N\mathrm{}`$ limit of the model (1) pointing out a precise, yet very simple, relation with Conformal Field Theory (CFT) universality classes, or, more exactly, with, say, the holomorphic sector of them. We consider a CFT in $`D=1+1`$ defined on a strip or on a cylinder. We will show that in the large $`N`$ limit the CM model is equivalent to two different CFTs, depending on the root system considered.

We stress that the large $`N`$ limit we consider here is not the thermodynamic limit, where the particle density is kept fixed. The latter limit was considered in previous papers dealing with similar models. In Ref.  for instance, using finite–size scaling relations, it was found that the low–energy physics of a Calogero–Sutherland (CS) model (see Eq. (18) below) is described by a $`c=1`$ CFT. In Ref. , a second quantized ($`D=1+1`$) Hamiltonian has been studied, which corresponds to the CS mechanical model. It was shown that that if one takes into account the leading ($`1/N`$) and subleading ($`1/N^2`$) terms in the thermodynamic limit, the dynamics of small fluctuations around the Fermi surface is described by an effective field theory endowed with an extended $`W_{1+\mathrm{}}`$ conformal symmetry. Our approach is rather different: we do not consider a thermodynamic limit and will not eventually write down an effective field theory; on the contrary we use the *entire* spectrum, not only the lowest–lying exited states, to reconstruct a CFT2 via a large $`N`$ limit. Needless to say that we rely on the exact integrability of the CM model, while, for instance, in Ref.  no use is made of the exact form of the spectrum.

The starting point is the very simple form of the spectrum of the CM model. It is well known that the spectrum is generated by $`N`$ decoupled harmonic oscillators . Moreover, it is possible to define similarity or unitary transformations that map the rational CM model (1) into a set of decoupled harmonic oscillators . In this sense the equivalence of the CM model with a CFT in the large $`N`$ limit is self-evident. We want to identify the CFT precisely.

Let us briefly recall the spectral properties of the CM model (see Ref.  for a review). Whenever the stability condition (2) is satisfied the Hamiltonian $`H`$ is self–adjoint and the energy level $`E_k`$ ($`H\psi _k=E_k\psi _k`$) is characterized by $`N`$ integers. If the couplings $`g_\alpha `$ are written as usual as $`g_\alpha ^2=\stackrel{~}{\mu }_\alpha (\stackrel{~}{\mu }_\alpha 1)|\alpha |^2/2`$, with $`\stackrel{~}{\mu }_\alpha =\mu _\alpha +\mu _{2\alpha }`$, the ground state energy $`E_0`$ is given by

$$E_0=\omega \left(\frac{N}{2}+\mu \right),\mu =\underset{\alpha R_+}{}\mu _\alpha .$$
(3)
The energies of the excited states read

$$E_{n_1,\mathrm{},n_N}=E_0+\omega (\nu _1n_1+\mathrm{}+\nu _Nn_N),$$
(4)
where the positive integers $`\nu _1,\mathrm{},\nu _N`$ are the degrees of the $`N`$ polynomial invariants of the Coxeter group, i.e. the orders of the Casimir operators of the Lie algebra under consideration. Given the order of the invariants, it is immediate to compute the quantum partition function $`Z(\beta )`$:

$`Z(\beta )=`$ (6)
$`e^{\beta E_0}[(1e^{\beta \omega \nu _1})(1e^{\beta \omega \nu _2})\mathrm{}(1e^{\beta \omega \nu _N})]^1.`$
The orders of the invariants for the classical root systems are known from the classic work of Racah , and are quoted in table I. Neglecting as usual the zero point energy, which is divergent in the large $`N`$ limit, see Eq. (3), the partition functions have a very simple form in all cases.

Let us now consider the CFT side. The Hilbert space of a generic two-dimensional conformal field theory is given by (see for instance Ref. )

$$\underset{h,\overline{h}}{}V(c,h)\overline{V}(c,\overline{h}),$$
(7)
where the sum is taken over the conformal dimensions of the theory and $`c`$ is the central charge. $`V(c,h)`$ is the Verma module generated from a vacuum $`|h`$ with conformal weight $`h`$, $`L_0|h=h|h`$, by the action of the holomorphic Virasoro generators $`L_n`$. $`\overline{V}(\overline{h},c)`$ describes the antiholomorphic sector (the holomorphic and antiholomorphic weights are in general independent quantities). The Verma modules need not to be irreducible representations of the conformal group. Indeed in general they are reducible. That means that the states

$$L_{k_1}L_{k_2}\mathrm{}L_{k_n}|h,1k_1\mathrm{}k_n,$$
(8)
are not necessarily all independent. It may happen that a combination of the states (8) defines a null vector $`|\chi `$, that is a vector annihilated by all of the generators $`L_n`$ with $`n>0`$ (we are considering only the holomorphic sector. For the antiholomorphic one the same statements hold). A null vector is orthogonal to the whole Verma module and in particular to himself, $`\chi |\chi =0`$. This is of course also true for all of its descendants. The null submodule arising from a null state does not contribute to the partition function, which is computed tracing over an irreducible representation of the Virasoro algebra. Given the module $`V(c,h)`$, an irreducible representation $`M(c,h)`$ is obtained by quotienting out of it all the null modules, i.e. identifying two states differing by a null vector.

The null vectors are in correspondence with the zeros of the Kac determinant of the Gram matrix $`G`$, defined as $`G_{ij}=i|j`$, where $`|i`$ and $`|j`$ are the states (8) of the Verma module $`V(c,h)`$. The Gram matrix is block diagonal, with blocks $`G^{(l)}`$ corresponding to the level $`l`$, i.e. the $`L_0h`$ eigenvalue of the states (8). The expression for the Kac determinant is (see e.g. Ref. )

$$detG^{(l)}=\alpha _l\underset{\begin{array}{c}r,sl\\ rsl\end{array}}{}[hh_{r,s}]^{p(lrs)},$$
(9)
where $`p(lrs)`$ is the number of partitions of the integer $`lrs`$ (not to be confused with the product $`p\times (lrs)`$), and $`\alpha _l`$ denotes a positive constant depending on $`l`$. $`h_{r,s}`$ are functions of the central charge $`c`$ and can be written as

$`c=136\left(t+{\displaystyle \frac{1}{t}}\right),`$ (10)
$`h_{r,s}={\displaystyle \frac{1}{4}}(r^21)t+{\displaystyle \frac{1}{4}}(s^21){\displaystyle \frac{1}{t}}{\displaystyle \frac{1}{2}}(rs1),`$ (11)
where $`t`$ is a function of $`c`$. Once we have identified the null vectors and constructed the irreducible representation $`M(c,h)`$, the partition function for the holomorphic sector with conformal weight $`h`$ can be defined as

$$Z_h(\beta )=\text{Tr}_{M(c,h)}q^{(L_0h)},$$
(12)
where $`qe^\beta `$. Notice that the trace is taken in $`M(c,h)`$. We bear in mind that $`L_0`$ is, in proper units, the energy operator $`H`$ for $`M(c,h)`$. In our definition we have again neglected the zero point energy (basically the vacuum conformal weight). We note that $`Z_h(\beta )`$ depends on $`h`$ since $`M(c,h)`$ does. Now the question is if one can define a CFT such that $`Z_h(\beta )`$ defined by Eq. (12) equals $`Z(\beta )`$ of Eq. (6) in the limit $`N\mathrm{}`$. The answer is simple and affirmative for the four cases $`A_N`$, $`B_N`$, $`C_N`$ and $`D_N`$.

$`A_N`$ is the case of major interest, both from the physical and mathematical point of view. We set $`\omega =1`$. From the table I it follows that all the modes but the first are present. In order to have the same for the CFT partition function, we simply need a null vector at the first level. This implies $`h=h_{1,1}=0`$, because at level one the only state is $`L_1|h`$. The corresponding field $`\varphi _{(1,1)}(z)`$ is a constant since it satisfy the differential equation

$$\frac{}{z}\varphi _{(1,1)}(z)X=0,$$

for every string of fields $`X`$. $`\varphi _{(1,1)}`$ is called the identity operator and denoted by $`𝕀`$. If we choose $`c`$ such that no more null vectors are present, then $`Z_h(\beta )`$ of Eq. (12) equals $`Z(\beta )`$ given by Eq. (6) since the vectors

$$L_1\left(L_{k_1}L_{k_2}\mathrm{}L_{k_n}\right)|0,1k_1\mathrm{}k_n,$$

do not contribute to the partition function and have the proper degeneracy (see also Eq. (9)). It is also evident that the Hilbert spaces of the two theories have the same structure. The condition that $`h_{r,s}0`$ for $`r1`$ and $`s1`$ restricts the possible values of $`c`$. A general elementary unitary bound on $`c`$ is $`c0`$ (we consider only unitary theories). From Eq. (10) we immediately see that $`c`$ cannot be equal to 1, since in this case $`h_{r,r}=0`$ $`r>0`$. On the other hand if $`c<1`$ the unitary theories are given by the well known minimal models whose Hilbert space have a more complicated structure. Indeed the unitary minimal models with $`c<1`$ can be labeled by $`(p=m+1,p^{}=m+1)`$ with $`m=2,3,\mathrm{}`$ and satisfy the periodicity property $`h_{r,s}=h_{r+p,s+p^{}}`$ which is sufficient to exclude them. We conclude that $`c>1`$ is the only possibility. From Eq. (10) it follows that $`t`$ as a function of $`c`$ has two branches and it is given by

$$t=1+\frac{1}{12}\left[1c\pm \sqrt{(1c)(25c)}\right].$$
(13)
From the previous equation we see that $`h_{r,s}`$ never vanishes if $`1<c<25`$. If $`c25`$ we can resort to another well-known expression of the roots of the Kac determinant, namely

$`c=1{\displaystyle \frac{6}{m(m+1)}},`$ (14)
$`h_{r,s}={\displaystyle \frac{[(m+1)rms]^21}{4m(m+1)}}.`$ (15)
If $`c25`$ then $`1<m<0`$, but this implies that $`h_{r,s}`$ vanishes only for $`r=s=1`$. We conclude that the limit $`N\mathrm{}`$ of the $`A_N`$ systems defines a two-dimensional CFT with $`c>1`$, $`h=0`$ and the conformal family given by the descendants of $`𝕀`$.

The cases $`B_N`$, $`C_N`$ and $`D_N`$ are less interesting. Setting $`\omega =1/2`$ we see from table I that the all the modes $`1,2,\mathrm{}`$ contribute to $`Z(\beta )`$ giving the same large $`N`$ limit for the three cases. The partition function basically is the character of the $`\widehat{u}(1)`$ algebra, the affine extension of $`u(1)`$, known as the Heisenberg algebra. This means that in the limiting CFT the Kac determinant never vanishes. This condition is ensured taking $`h>0`$ and $`c>1`$. It follows that we have one conformal family descendent from the primary field $`|h=\varphi (0)|0`$.

It is interesting to note that while the CM model can be defined even for the noncrystallographic dihedral group $`I_2(N)`$, with $`N=5`$, $`N7`$, the large $`N`$ limit is trivial in this case, since there are only two Casimirs of order $`2`$ and $`N`$.

For high temperatures $`\beta \omega 1`$, we obtain from the partition function (6) the energy–temperature relation in the large $`N`$ limit,

$$E=E_0+\frac{\pi ^2T^2}{6\omega }.$$
(16)
This can be compared with the generic CFT behavior $`EE_0=\frac{c}{12}\pi LT^2`$, which is valid if the characteristic size $`L`$ of the system is much larger than the thermal wavelength of the particles. In this way, we get the equation

$$\frac{2\pi }{\omega }=cL,$$
(17)
which relates the frequency $`\omega `$ of the CM model and the product $`cL`$ of the central charge and the characteristic size of our two-dimensional CFT on the strip or cylinder.

As a final consideration we note that the procedure used in this letter cannot be employed for the CS model,

$$H=\frac{1}{2}\underset{j=1}{\overset{N}{}}\frac{^2}{x_j^2}+\gamma (\gamma 1)\frac{\pi ^2}{L^2}\underset{i<j}{}\frac{1}{\mathrm{sin}^2(\pi (x_ix_j)/L)},$$
(18)
where $`\gamma `$ is a positive constant, often taken to be an integer. Indeed in this case the difference between two successive energy levels is proportional to $`N`$, so that a large $`N`$ limit is meaningless.

Summarizing, in this letter we have considered a large $`N`$ limit for the CM model in order to reconstruct a CFT2 from the exactly known spectrum of exited states of the model. We have found two different CFTs both with central charge $`c>1`$. The value of the central charge as well as the dimension of the primary field are dictated only by the algebraic underlying Lie symmetries of the CM model and are independent of the couplings strength (provided the stability condition (2) is satisfied).


text:
# Automatic Debugging Support for UML Designs 11footnote 1In M. Ducassé (ed), proceedings of the Fourth International Workshop on Automated Debugging (AADEBUG 2000), August 2000, Munich. COmputer Research Repository (http://www.acm.org/corr/), cs.SE/0011017; whole proceedings: cs.SE/0010035.

## 1 Introduction

Size and complexity of software systems has increased tremendously. Therefore, the development of high-quality software requires rigorous application of sophisticated software engineering methods. One such method which has become very popular is the Unified Modeling Language. UML  has been developed by the “three amigos” Booch, Jacobson, and Rumbaugh as a common framework for designing and implementing object-oriented software. UML contains many different notations to describe the static and dynamic behavior of a system on all different levels and phases of the software design process.

Although UML provides a common notational framework for requirements and design, UML, as any other language, does not eliminate bugs and errors. These bugs must be found and fixed in order to end up with a correctly working and reliable system. It is well known, that debugging a large software system is a critical issue and can be a major cost-driving factor. Changes which have to be applied to the system (e.g., to fix a bug) are becoming substantially more expensive, the later they are detected (Figure 1). When an error is detected early during the definition phase, its cost is relatively low, because it only influences the requirements definition. Bugfixes in a product already shipped can be up to 60–100 times more expensive .

Therefore, it is mandatory to start with debugging as early in the project as possible. In this paper, we will discuss an approach which supports debugging of scenarios (more precisely UML sequence diagrams) with respect to given domain knowledge. This is done as a part of an algorithm  which can synthesize UML statecharts from a number of sequence diagrams. This synthesis step can be seen as a transformation from requirements to system design. It does not only facilitate fast and justifiable design from requirements (sequence diagrams), but also substantially helps to debug the generated designs. Because sequence diagrams usually cover only parts of the system’s intended behavior, the generated statecharts need to be refined and modified manually. By applying the synthesis algorithm in a “backward” way, the refined statechart can be checked against the requirements. Each conflict is reported to the user and indicates a bug.

For practical applicability of any debugging aid, the presentation of the bug, its cause and effect is of major importance. In our approach, we rely on logic-based explanation technology: all conflicts correspond to failure in logical reasoning about sequence diagrams, statecharts, and domain knowledge. Ongoing work, as discussed in the conclusions, uses methods from automated deduction to point the user to the exact place where the conflict occurred and which parts of the models and specification are affected.

This paper is organized as follows: Section 2 gives an overview of major UML notations and a typcial iterative software design process. Then we will describe how sequence diagrams are annotated for a justified synthesis of statecharts (Section 4). Based on this algorithm we discuss methods for debugging a sequence diagram and a synthesized statechart. In Section 7 we discuss future work and conclude.

Throughout this paper, we will use one example to illustrate our approach. The example concerns the interaction between an espresso vending machine and a user who is trying to obtain a cup of coffee. This example (based on the ATM example discussed in ) is rather small, yet complex enough to illustrate the main issues. The requirements presented here are typical scenarios for user interaction with the machine (e.g., inserting a coin, selecting the type of coffee the user wants, reaction on invalid choices, and pressing the cancel button). More details of the requirements will be discussed when the corresponding UML notations have been introduced.

## 2 UML

The Unified Modeling Language is the result of an effort to bring together several different object-oriented software design methods. UML has been developed by Booch, Jacobson and Rumbaugh  and has gained wide-spread acceptance. A variety of tools support the development in UML; among them are Rhapsody , Rational’s Rose , or Argo/UML .

On the top-level, requirements are usually given in the form of use cases, describing goals for the user and system interactions. For more detail and refinement, UML contains three major groups of notations: class diagrams for describing the static structure, interaction diagrams for requirements, and state diagrams and activity diagrams for defining dynamic system behavior. Below, we will illustrate the notations which are important for our approach to debugging of UML designs.

### 2.1 Software Development with UML

Although no explicit development process is prescribed for UML, UML design usually follows the steps of Inception, Elaboration, Construction, and Transition, used in an iterative manner. In this paper, we will not elaborate on the process model. For details, cf., e.g., . The importance of support for debugging of UML designs on the level of sequence diagrams (requirements), and statecharts becomes evident, when we look at a graphical representation of an iterative development process (Figure 2). The design starts by analyzing the (physical) process at the lower left part of the figure. The result of the analysis comprises the requirements (e.g., as a set of sequence diagrams), and knowledge about the domain (henceforth called domain theory). Based on these, a model of the system is developed, consisting of class diagrams, statecharts and activity diagrams. This model must now be implemented. Modern software engineering tools provide automatic code-generation (or at least support) for this step. Finally, the produced system must be verified against the physical process, and its performance tuned.

Traditionally, the way to get a working system is simulation (process–requirements–model), and testing (requirements–model–system). Here, errors and bugs have to be found and removed. Within an iterative design process, these steps are performed over and over again, depicted by the circular arcs. To keep these iterations fast (and thus cost-effective), powerful techniques for debugging requirements against domain knowledge, and models against requirements are vital. Our approach supports this kind of debugging and it will be discussed in the next section, following a short description of the basic concepts of class diagrams, sequence diagrams, and statecharts.

### 2.2 Class Diagram

A *class diagram* is a notation for modeling the static structure of a system. It describes the classes in a system and the relationships between them. Figure 3 shows an example of a class diagram for our coffee-vending machine example. In an object-oriented fashion, the main class (here “coffee machine”) is broken down into sub-classes. The aggregation relation ($`\text{}\mathrm{}`$) shows when one class is *part of* another one. The generalization relation ($`\text{}`$) shows when one class is *an instance of* another. For further details, see e.g., .

### 2.3 Statecharts

*Statecharts* , are finite state machines extended with hierarchy and orthogonality. They allow a complex system to be expressed in a compact and elegant way. Figure 4 shows a simple example of a statechart. Nodes can either be simple nodes (A1, A2, A3, B, and C), or composite nodes (node A in the figure) which themselves contain other statecharts. The initial node in a statechart is marked by $``$. Transitions between states have labels of the form $`e[c]/a`$. If event $`e`$ occurs and guard $`c`$ holds, then the transition may be selected to fire which results in action $`a`$ being taken and a state change occurring. This behavior is extended in a natural way to handle composite nodes.

### 2.4 Sequence Diagrams

Scenarios describe concrete examples of the system’s intended behavior. In UML scenarios can be expressed as *sequence diagrams*. A sequence diagram (SD) shows the interaction between objects of a system over time. The SD in Figure 6 is an example for interactions between the objects “User”, the user interface of the coffee machine (“Coffee-UI”), and the machine (“Control”) itself. The vertical lines represent the time-line for the given object, defining the object’s life during the interaction. Messages (like “Insert coin”) are exchanged between the objects. Figure 6 is a different scenario for our coffee-machine. It describes an invalid selection by the user (e.g., choosing sugar and sweetener at the same time).

## 3 Extending Sequence Diagrams

The simplicity of sequence diagrams makes them suitable for expressing requirements as they can be easily understood by customers, requirements engineers and software developers alike. Unfortunately, the lack of semantic content in sequence diagrams makes them ambiguous and therefore difficult to interpret. Let us assume that in our example, there exists an additional sequence diagram, SD0, identical to SD1 in Figure 6 except that there are two “Insert coins” messages adjacent to each other. There are three possible ways to interpret the conjunction of the two SDs — either a cup of coffee costs one or two coins (ridiculous!), or it costs just one coin, in which case SD0 is incorrect. The other case (two coins needed) invalidates SD1. In practice, such ambiguities are often resolved by examining the informal requirements documentation but, in some cases, ambiguities may go undetected leading to costly software errors.

For the automatic generation of (conflict-free) designs, such documents are usually too informal. On the other hand, the need to provide a full formal domain theory containing all semantic information is clearly too much a burden for the designer and thus not acceptable in practice.

Our approach allows for a compromise: the user can annotate messages in a sequence diagram with a pre/post-condition style specification expressed in OCL, UML’s logic-based specification and constraint language. For successful conflict detection (and statechart synthesis), only a small percentage of messages need to be annotated at all. This specifications should include the declaration of global *state variables*, where a state variable represents some important aspect of the system, e.g., whether or not a coin is in the coffee-vending machine. Pre- and post-conditions should then include references to those variables. Our experience with the case studies carried out so far (see Conclusions) is that the state variables and their data types usually directly “fall out” from the class diagram. Note that not every message needs to be given a specification, although, clearly, the more semantic information that is supplied, the better the quality of the conflict detection. Currently, our algorithm only exploits constraints of the form $`var=value`$, but there may be something to be gained from reasoning about other constraints using an automated theorem prover, e.g.,  or constraint solving techniques.

Fig. 7 gives specifications for selected messages in our coffee-machine example. Here, the state variables are the boolean variables CoinInMachine, CoinInReturnSlot, CoffeeTypeSelected, the variable Coin reflecting the number of coins in the machine (0, or 1), and SelectedCoffeeType. In order to talk about all values of the state variables at a given point, we use the notion of a *state vector*. This is a vector of values of the state variables. In our example, the state vector has the following form:

```
           <CoinInMachine^, CoinInReturnSlot^, CoffeeTypeSelected^,
            Coin^, SelectedCoffeeType^>

```

The notation $`var\widehat{}`$ extends the possible value for a state variable by an undetermined value, denoted by a “?”, i.e., $`var\widehat{}Dom(var)\{\mathrm{?}\}`$. For use with our algorithm, we will annotate each message of a sequence diagram with a statevector where the values of the state variables are determined by the algorithm described below.

## 4 Automatic Synthesis of Statecharts from Sequence Diagrams

The framework for debugging UML designs is based upon an algorithm for automatic synthesis of statecharts from sequence diagrams and a domain theory . The process to convert a number of SDs into a structured statechart consists of several steps: in the first step, each SD is annotated and conflicts between the SD and the domain theory (and hence, other SDs) are detected and reported to the user. Then, a statechart for each of the objects in the SD is generated; and all statecharts for an object are merged into a single statechart. The final step of the synthesis introduces hierarchy by grouping nodes into composite nodes, thus enhancing readability. In this paper, we are only concerned with the first, conflict detection part (as a basis for debugging), and the final result, the statechart. For details on the algorithm see .

There are two kinds of constraints imposed on a sequence diagram: constraints on the state vector given by the OCL specification, and constraints on the ordering of messages given by the SD itself. These constraints must be solved and arising conflicts be reported to the user. More formally, the process of conflict detection can be written as follows. An annotated sequence diagram is a sequence of messages $`m_1,\mathrm{},m_n`$, with

$$s_0^{\text{pre}}\stackrel{m_1}{}s_0^{\text{post}},s_1^{\text{pre}}\stackrel{m_2}{}\mathrm{}\stackrel{m_{r1}}{}s_{r1}^{\text{post}},s_r^{\text{pre}}\stackrel{m_r}{}s_r^{\text{post}}$$
(1)
where the $`s_i^{\text{pre}}`$, $`s_i^{\text{post}}`$ are the state vectors immediately before and after message $`m_i`$ is being sent. $`S_i`$ will be used to denote either $`s_i^{\text{pre}}`$ or $`s_i^{\text{post}}`$; $`s_i^{\text{pre}}[j]`$ denotes the element at position $`j`$ in $`s_i^{\text{pre}}`$ (similarly for $`s_i^{\text{post}}`$).

In the first step of the synthesis process, we assign values to the variables in the state vectors as shown in Figure 8. The variable instantiations of the initial state vectors are obtained directly from the message specifications (lines 1,2): if message $`m_i`$ assigns a value $`y`$ to a variable of the state vector in its pre- or post-condition, then this variable assignment is used. Otherwise, the variable in the state vector is set to an undetermined value $`\mathrm{?}`$. Since each message is specified independently, the initial state vectors will contain a lot of unknown values. Most (but not all) of these can be given a value in one of two ways: two state vectors, $`S_k`$ and $`S_l`$ ($`kl`$), are considered the same if they are unifiable (line 6). This means that there exists a variable assignment $`\varphi `$ such that $`\varphi (S_k)=\varphi (S_l)`$. This situation indicates a potential loop within a SD. The second means for assigning values to variables is the application of the frame axiom (lines 8,9), i.e., we can assign unknown variables of a pre-condition with the value from the preceeding post-condition, and vice versa. This means that values of state variables are propagated as long as they are not changed by a specific pre- or post-condition. This also assumes that there are no hidden side-effects between messages.

A conflict (line 11) is detected and reported if the state vector immediately following a message and the state vector immediately preceding the next message differ.

Example. Let us consider how this algorithm operates on the first few messages of SD1 from Figure 6. When annotating the first message (“Display Ready Light”), we obtain the following state vector on the side of the user-interface: $`S_1=\text{<F,F,?,?,?>}`$. The values of the first two state variables are determined by the message’s pre-condition in the domain theory. The state-vector $`S_2`$ on the receiving side of our message only consists of “?”. As a pre-condition for the message “Insert coin” we have CoinInMachine = F. Thus we have $`S_3=\text{<F,?,?,?,?>}`$ as the state vector. All other messages in SD1 are annotated in a similar way. Now, our algorithm (lines 4–12) tries to unify state vectors and propagate the variable assignments. In our case, the attempt to unify $`S_2`$ with $`S_3`$ would assign the value F to the first variable in $`S_2`$, yielding $`S_2=\text{<F,?,?,?,?>}`$. Now, both state vectors are equal. Then, variable values are propagated using the frame axiom. In our case, we can propagate the value of CoinInReturnSlot = F (from $`S_1`$) into $`S_2`$ and $`S_3`$, because the domain theory does not prescribe specific values of this state variable at these messages. Hence, its current value F can be used in the other state vectors, finally yielding $`S_2=S_3=\text{<F,F,?,?,?>}`$. After performing all unification and propagation steps, we obtain an annotated sequence diagram as shown in Figure 9. The conflict indicated there will be discussed in the next section.

## 5 Debugging a Sequence Diagram

The algorithm from the previous section detects conflicts of a SD with the domain theory (and thus with other sequence diagrams). Any such conflict which is detected corresponds to a bug which needs to be fixed. The bug can be in the sequence diagrams, which means that one or more sequences of actions are not compatible with the domain theory, and henceforth with other SDs. Such a situation often occurs when sequence diagrams and domain theory for a large system are developed by different requirements engineers. Our algorithm is capable of directly pointing to the location where the conflict with the domain theory occurs. The respective message, together with the instantiated pre- and post-conditions, as well as the required state vector values are displayed. This feature allows to easily debug the sequence diagram. Of course, the error could be in the domain theory instead. For example, one designer could have set up pre- or post-conditions which are too restrictive to be applicable for scenarios, specified by other designers. In that case, the domain theory must be debugged and modified. Our algorithm can also provide substantial support here, because it is able to display the exact location where the conflicting state variables have been instantiated. Especially in long sequence diagrams the place where a state variable is instantiated and the place where the conflict occurs can be far apart. The current version of our algorithm provides only rudimentary feed-back as demonstrated in the example below. Future work (which also allows richer OCL constructs to be used) requires more elaborate, human-readable descriptions of the error trace. Automated theorem provers and work on proof presentation, like the ILF system  will be used for that purpose. Such a system will not only explain  the possible reasons for a conflict, but can also give (heuristics-driven) hints to the user on how to fix the problem.

Example. The following example shows, how conflict detection can be used for debugging: Figure 9 shows SD1 from Figure 6 after the state vectors have been extended by our algorithm of Figure 8. Our procedure has detected a conflict with the domain theory. As an output it provides the messages and state vectors which are involved in the conflict:

```
Conflict in SD1: Object Coffee-UI
 statevector after  "Insert coin"        = <T,F,T,1,none> [Msg 2]
 statevector before "Request Selection"  = <T,F,F,1,none> [Msg 3]
  conflict in variable "CoffeeTypeSelected"
  conflict occurred as consequence of unification of
   statevector after "Display Ready Light" = <F,F,T,0,none> [Msg 1]
   statevector after "Display Ready Light" = <F,F,T,0,none> [Msg 11]
   statevector after  "Take coin"          = <F,F,T,0,none> [Msg 10]

```

This arises because state vectors SV1 (state vector before “Display Ready Light”) and SV2 (after “Take coin”) are unified (Figure 9 shows the instantiations of the vectors after unification). This corresponds to the fact that the coffee machine returns to its initial state after “Take coin” is executed. The state vectors tell us that there is a potential loop at this point. A second execution of this loop causes the state variable “CoffeeTypeSelected” to true, when the system asks for a selection. However, the domain theory tells us that this variable must be false as a pre-condition of the “Request Selection” message. Hence, there is a conflict, which represents the fact that the developer probably did not account for the loop when designing the domain theory.

The user must now decide on a resolution of this conflict — i.e., to debug this situation. The user either

* can tell the system that the loop is not possible, in which case the unifier that detected the loop is discarded. This amounts to modifying the annotated sequence diagram (by restricting possible interpretations). The user can
* modify the sequence diagram at some other point, e.g., by adding messages; or
* modify the domain theory. In our example, the action taken might be that the domain theory is updated by giving “Release coin” the additional postcondition CoffeeTypeSelected = false. This extra post-condition resets the value of the variable (i.e., the selection) when the user is asked to remove the coin. The position of the change has been obtained by systematically going backwards from SV2. Although possible locations are automatically given by the system, the decision where to fix the bug ( at “Release coin” or at “Take coin”) must be made by the user. Here, the second possibility was chosen, because the specification for that message modified a state variable which is related to the variable which caused the conflict.

## 6 Debugging a Synthesized Statechart

When the statechart synthesis algorithm successfully terminates, it has generated a human-readable, hierarchically structured statechart, reflecting the information contained in the SDs and the domain theory. In general, however, sequence diagrams usually describe only parts of the intended dynamic behavior of a system. Therefore, the generated statechart can only be a skeleton rather than a full-fledged system design. Thus, the designer usually will extend, refine, and modify the resulting statechart manually. Our approach takes this into account by generating a well structured, human-readable statechart which facilitates manual refinement and modification.

However, these manual actions can be sources of errors which will have to be found and removed from the design. In the following, we describe two approaches, addressing this problem.

### 6.1 Classical Debugging

The traditional way to find bugs in a statechart is to run simulations and large numbers of test cases. Most commercial tools for statecharts, like Betterstate, Statemate, or Rhapsody support these techniques. Some tools also provide more advanced means for analysis, like detection of deadlocks, dead branches, non-deterministic choices, or even model checking for proving more elaborate properties. In this paper, we will not discuss these techniques.

### 6.2 Debugging w.r.t. Requirements

Whenever a design (in our case the statechart) is modified, care must be taken that all requirements specifications are still met, or that an appropriate update is made. Traditionally, this is done manually by updating the requirements document (if it is done at all). Bugs are usually not detected (and not even searched for) until the finished implementation is tested. Thereby, late detection of bugs leads to increased costs. By considering the “reverse” direction of our synthesis algorithm, we are able to

* check that all sequence diagrams are still valid, i.e., that they represent a possible sequence of events and actions of the system
* detect conflicts between the current design (statechart) and one or more SDs, and
* detect inconsistencies with respect to the domain theory.

The basic principle of that technique is that we take one sequence diagram after the other, together with the domain theory, and check if that sequence of messages is a possible execution sequence in the given statechart. Here again we use logic-based techniques, similar to those described above (unification of state vectors, value propagation with the frame axiom). An inconsistency between the (modified) statechart and the SD indicates a bug (in the SD or SC). By successively applying patches to the SD (by removing or adding messages to the SD) the algorithm searches for possible ways to obtain an updated and consistent SD. Since in general more than one possible fix for an inconsistency exists, we perform an iterative deepening search resulting in a solution with the fewest modifications to the sequence diagram. We are aiming to extend this search by applying heuristics to select “good” fixes.

Here again, the form of feed-back to the user is of major importance. We are envisioning that the system can update the requirements and provide explanations for conflicts in a similar way as described above.

Example. The statechart in Figure 10 has been refined. The transition between $`N_2`$ and $`N_3`$ has been extended in such a way that first event $`e_2`$, then $`e_3`$ with action $`a_3`$ has to occur before the state $`N_3`$ is reached. The original statechart has been generated from a sequence diagram as shown on the right-hand side of Fig. 10. The modification of the statechart is propagated back to the sequence diagrams where the change is clearly marked. In this example, the extension could be made without causing a conflict. However, it is advisable for the designer and/or the requirements engineer to carefully observe these changes in order to make sure that these modified requirements still meet the original intended system behavior.

## 7 Future Work and Conclusions

We have presented a method for debugging UML sequence diagrams and statecharts during early stages in the software development process. Based on an algorithm, designed for justified synthesis of statecharts, we have identified two points where conflicts (as a basis for debugging) can be detected: during extending the annotations of a SD (conflicts w.r.t. the domain theory), and updating of sequence diagrams based upon a refined or modified statechart.

The algorithm which is described in  has been implemented in Java and has been used for several smaller case studies in the area of object-oriented systems, user interfaces, and agent-based systems . Current work on this part include integration this algorithm into a commercial UML tool (MagicDraw). Currently we are extending our synthesis algorithm to provide the debugging facilities described in this paper. Future work will mainly focus on integrating and extending explanation technology into our system.

Debugging large designs with lengthy and complex domain theories vitally depends upon an elaborate way of providing feed-back to the user. Starting from the basic information about a conflict (i.e., a failed unification), we will use theorem proving techniques of abduction and counter-example generation to provide as much feed-back as possible on where the bug might be, and how to fix the problem<sup>2</sup><sup>2</sup>2 This problem essentially is equivalent to finding which hypotheses are missing or wrong when a conjecture cannot be proven valid. . These techniques will be combined with tools capable of presenting a logic statement in human-readable, problem-specific way (e.g., ILF ). Only, if debugging feedback can be given in the notation of the engineering domain rather than in some logic framework, such debugging aids will be accepted in practice.

It is believed that UML (and tools based upon this notation) will have a substantial impact on how software development is made. By providing techniques which do not only facilitate design by synthesis, but also provide powerful means to debug requirements and designs in early stages we are able to contribute to tools which are useful in design of large software systems.


text:
# A young stellar group associated with HD 199143 (𝑑 = 48 pc) Based on observations collected at the European Southern Observatory, La Silla, Chile, observations made by IUE at NASA-GSFC, and with the Isaac Newton Telescope operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos

## 1 Introduction

Zuckerman & Webb (2000) sketch a picture of the recent star formation history of the solar neighbourhood in which 10–40 million years ago an ensemble of molecular clouds were forming stars at a modest rate near the present position of the Sun. About 10 Myrs ago, the most massive of these newly formed stars exploded as a supernova, terminating the star formation episode and generating the very low density region seen in most directions from the present Sun (Welsh et al. 1998). This scenario can not only explain the presence of young stellar groups close to the earth, but also explains how the $`\beta `$ Pic moving group can be so young (20 Myr; Barrado y Navascués et al. 1999), and yet so close. However, currently this scenario is largely speculative.

HD 199143 is a poorly studied bright ($`V`$ = 7$`\stackrel{m}{.}`$27) star in the constellation of Capricornius. It has been classified as F8V in the Michigan Spectral Survey (Houk & Smith-Moore 1988), after an initial classification of G0 by Cannon & Mayall (1949). The star would be completely inconspicuous, if it hadn’t been detected as a bright extreme-ultraviolet source by the ROSAT (2RE J205547$``$170622) and Extreme Ultraviolet Explorer (2EUVE J2055$``$17.1) missions (Pounds et al. 1993; Malina et al. 1994; Pye et al. 1995; Bowyer et al. 1996).

In this Letter we present new optical and ultraviolet spectroscopy of HD 199143 and show that it is a variable and rapidly rotating F8V star. We infer that all characteristics of the HD 199143 system can be explained by assuming that it is a binary system, in which the primary has been spun up by accretion of mass from a low-mass companion. Its association with a previously studied T Tauri-like system (BD$``$176128) suggests that these two stars could be the first two members of a close (48 pc) new region of recent star formation and may provide compelling support for the star formation history of the solar neighbourhood outlined in the first paragraph.

## 2 Optical Observations

Low-resolution (0.53 Å pix-1) spectra of HD 199143 in the wavelength region of 5700–6800 Å were taken with the 2.5 m Isaac Newton Telescope (INT) at La Palma during the nights of August 29 (JD 2451023.072), 30 (JD 2451024.122) and 31 (JD 2451025.076), 1998. The spectra were reduced with the usual steps of bias subtraction, flatfielding, background subtraction and spectral extraction, and wavelength and flux calibration. Apart from a multiplicative factor, due to the fact that the nights in which the spectra were taken were not of photometric quality, the spectra taken in the different nights were identical. In Fig. 1 we show the spectrum of Aug. 30, 1998. For comparison we also show the spectrum of HD 6111 (F8V), obtained from the spectral database by Jacoby et al. (1984). The resolution of this spectrum is about three times lower than that of the INT spectrum. Apart from the differences expected because of the differences in spectral resolution, the two spectra are identical, confirming the F8V spectral classification of HD 199143.

High-resolution (0.05 Å pix-1) spectra of HD 199143 in the H$`\alpha `$ (6536–6591 Å) and Na i D (5858–5910 Å) wavelength ranges were obtained with the Coudé Auxiliary Telescope (CAT) at La Silla, Chile, on Dec. 16 (at JD 2450432.016) and Dec. 15 (JD 2450431.016), 1996. The spectra were reduced in a standard fashion, after which the continuum was normalized to unity. They are shown in Figs. 2 and 3, together with the spectra of HR 963 (F8V), obtained during the same night as the HD 199143 spectra. Apart from a number of very sharp absorption features due to water vapour in the earth’s atmosphere, a number of highly broadened (FWHM $``$ 250 km s-1) photospheric absorption lines (most prominently H$`\alpha `$ and the Na i doublet) are visible in the HD 199143 spectra. The same lines are present in HR 963, but much narrower, again confirming the spectral classification of F8. From the photospheric lines we measure a radial velocity of $``$9 $`\pm `$ 16 km s-1 for HD 199143. The wings of the H$`\alpha `$ profile appear identical in HD 199143 and HR 963, demonstrating that the broadening of the lines in HD 199143 is not due to a luminosity classification smaller than V, but must be caused by a high (a few hundred km s-1) value of $`v\mathrm{sin}i`$.

## 3 Ultraviolet Observations

HD 199143 was observed with IUE on several occasions in 1995, under a discretionary program (OD89Z), and the observatory program (USSBS). Archival data taken by the NC119 program, in November 1992, are also included in our analysis. A total of 19 archived images are available, both in low- (1.68 Å pix-1 for SWP, 2.66 Å pix-1 for LWP) and high-dispersion (25 km s-1, $`\lambda 0.2`$ Å resolution), which were secured through the large aperture (oval-shaped: 10″ $`\times `$ 20″). There are 11 SWP camera (1150–2000 Å) exposures all in low-dispersion, and eight LWP (1800–3200 Å) exposures, of which six are in high-dispersion.

A preliminary inspection of the short-wavelength data reveals an emission spectrum typical of T Tauri, Herbig Ae/Be (HAeBe) stars and planetary nebulae. The SWP low-dispersion images present clear emissions in N v (1240 Å), C ii (1335 Å), Si iv (1394 Å), C iv (1550 Å), C iii (1577 Å), He ii (1640 Å) and C i (1657 Å). With the exception of a few lines such as N v and He ii, the emission spectrum of HD 199143 resembles the spectra of T Tauri stars such as RW Aur and GW Ori (Imhoff & Appenzeller 1987). In Fig. 4 we present the emission spectrum corresponding to the well-exposed (80 min.) image SWP 55183. This spectrum is compared with the spectrum of HR 963 (SWP 52921), which is a typical F8V, high-proper motion star.

In the long-wavelength spectrum the most conspicuous emission is the Mg ii doublet at 2798.0 and 2802.7 Å, which is only detectable in high-dispersion exposures. This normally P-Cygni emission has been shown to be typical of chromospherically active stars and other pre-main sequence (PMS) stars such as T Tauri and HAeBe objects. The Mg ii doublet from the well-exposed (110 min.) image LWP 30968 is presented in Fig. 5 along with the comparison spectrum of HR 963. Note the self-absorption feature on top of the h & k emissions. The Mg ii spectrum of HD 199143 resembles that of the T Tauri star GW Ori (Imhoff & Appenzeller 1987), or that of the “double emission peak” HAeBe stars classification described by Imhoff (1994).

## 4 Analysis

The ultraviolet emission spectrum of HD 199143 is somewhat peculiar because emission lines such as N v, which corresponds to a temperature regime of about $`2\times 10^5`$ K, is rarely present in PMS stars. Furthermore, lines such as He ii, which has a complex origin, are commonly present in planetary nebulae and in only a handful young objects. The ratio of the emission fluxes for C iv to Si iv in T Tauri stars ranges from 2–3 (typical of chromospheric activity) to less than unity (Imhoff & Appenzeller 1987). In HD 199143 this range is from 1.2 to 1.4, compatible with an origin in a T Tauri star.

Repeated observations with IUE allow us to address the issue of ultraviolet variability, beyond the instrument flux repeatability (3%). Comparisons of the well-exposed section of the IUE images indicate a variability of 10–20% in the continuum and 20-50% in the emission lines of C ii, C iv, N v and He ii. From the repeated SWP images taken closely in time, we detected that flux variability (continuum and lines) was found to be random and not associated with the period of 1.6 days recently suggested by Handler (1999), whom classified HD 199143 as a $`\gamma `$ Doradus candidate.

We have searched for infrared emission by checking the raw IRAS scans at the position of HD 199143 using routines from the Groningen Image Processing System (GIPSY). In the 12 $`\mu `$m band, a point-like source is clearly present at the position of HD 199143. No source was detected at longer wavelengths. From these data we derive a flux of 0.24 $`\pm `$ 0.04 Jy in the IRAS 12 $`\mu `$m band, and upper limits of 0.12, 0.12 and 0.30 Jy for the fluxes at 25, 60 and 100 $`\mu `$m, respectively.

Using the newly determined IRAS fluxes and the optical photometry of HD 199143 by Olsen (1983) and Cutispoto et al. (1999), we constructed a Spectral Energy Distribution (SED) of HD 199143, shown in Fig. 6. In this plot we also show UV fluxes of HD 199143 from archive IUE data. Also plotted is a Kurucz (1991) model for the photosphere of a F8V star, fitted to extinction-corrected optical photometry of HD 199143. In the SED we can see that both the IRAS 12 $`\mu `$m flux and the IUE fluxes below 2250 Å are significantly higher than that expected from the stellar photosphere. One explanation for the infrared excess could be the presence of circumstellar dust in the system, similar to that found in Vega-type systems. However, the tight 25 and 60 $`\mu `$m upper limits show that if this is the case, only a very warm ($`>`$ 1000 K) dust component must be present, which is unlikely. A more likely source of the observed infrared excess might be a late-type companion to HD 199143, or to infer a modification of the photospheric structure of HD 199143 due to its rapid rotation (Sect. 2).

The Hipparcos catalogue (ESA 1997) lists a parallax of 21.0 $`\pm `$ 1.0 milliarcseconds for HD 199143 (HIP 103311). Using a distance $`d`$ of 47.7 $`\pm `$ 2.4 pc obtained from this parallax we computed the luminosity of HD 199143 by integrating the flux of the Kurucz (1991) stellar photosphere model fitted to the optical photometry, and multiplying by $`4\pi d^2`$ to correct for spatial dilution. Note that if the photospheric structure of HD 199143 has been altered by its rapid rotation, this procedure might not be completely correct since the emitted flux will be non-spherically symmetric distributed. However, the resulting stellar luminosity of 2.4 $`\pm `$ 0.2 L is in agreement with that expected of a F8V star (2.1 L; Schmidt-Kaler 1982), showing that this is not a big effect. It also confirms our earlier conclusion that HD 199143 belongs to luminosity class V.

## 5 Discussion and Conclusions

The presence of a normal late-type companion could not explain the ultraviolet excess, or the EUVE and ROSAT detections of HD 199143. However, the presence of an accretion disk around our hypothetical companion, such as that found in LMXB or T Tauri systems, might easily explain those properties, as well as the infrared excess, the presence of emission lines and the variability. In such a scenario, the high rotational velocity of HD 199143 could be due to a spin-up in its past by accretion from the companion.

At first glance, a scenario in which a nearby main-sequence star like HD 199143 would have a T Tauri-like companion would seem far-fetched. However, Mathioudakis et al. (1995) report the presence a strongly flaring K7e–M0e dwarf with a high Li abundance only 5 arcminutes from HD 199143. The optical spectrum of this star, BD$``$176128, is identical to that of many T Tauri stars. From Digital Sky Survey images we identify BD$``$176128 with HD 358623. An inspection of the Tycho-2 Catalogue (Høg et al. 2000) shows that this star has a proper motion of 59 $`\pm `$ 3 and $``$63 $`\pm `$ 3 mas yr-1 in $`\mu _\alpha `$ and $`\mu _\delta `$, identical to that of HD 199143. From the fact that HD 358623 is the only star within a 5 degree radius for which this is the case, we exclude the possibility that this could be a coincidence and conclude that the two stars form a genuine proper motion pair. Using the data by Mathioudakis et al. (1995), and the newly determined distance, we compute the absolute luminosity of BD$``$176128 to be 0.34 $`\pm `$ 0.06 L, employing a similar procedure to that followed for HD 199143. Comparison with the pre-main sequence evolutionary tracks by D’Antona & Mazzitelli (1997) yields an age of 107 years for BD$``$176128, consistent with a T Tauri nature of this star.

Using the radial velocity of HD 199143 determined in Section 2, and the parallax and proper motions listed in the Hipparcos catalogue, we compute the galactic space velocity components $`(U,V,W)`$ of HD 199143 to be $`(10\pm 13,13\pm 6,13\pm 6)`$ km s-1. This space motion is similar to that of many stars in the vicinity of the Tucanae and TW Hydra associations (Zuckerman & Webb 2000), suggesting that these stars might have formed from the same cloud complex. We conclude that HD 199143 and BD$``$176128 could very well be the first two members of a region of recent star formation similar to the TW Hydrae Association and the newly identified Tucanae Association (Kastner et al. 1997; Zuckerman & Webb 2000). If confirmed, a further study of these two enigmatic stars could lead to a better understanding of the star formation history in the solar neighbourhood.

###### Acknowledgements.

The authors would like to thank the referee, Ben Zuckerman, for valuable suggestions for improving the manuscript. Bruce McCollum and Mario Pérez thank Yoji Kondo for Discretionary Time to observe HD 199143 with IUE. This research has made use of the Simbad data base, operated at CDS, Strasbourg, France.


text:
# An entropy–driven noise–induced phase transition

## Abstract

We introduce a class of exactly solvable models which exhibit an ordering noise-induced phase transition driven by an entropic mechanism. In contrast with previous studies, order does not appear in this case as a result of an instability of the disordered phase produced by noise, but of a balance between the relaxing deterministic dynamics and the randomizing (entropic) character of the fluctuations. A finite-size scaling analysis of the phase transition reveals that it belongs to the universality class of the equilibrium Ising model. All these results are analyzed in the light of the nonequilibrium probability distribution of the system, which can be derived analytically. The relevance of our result as a possible scenario of the so-called Lower Critical Solution Temperature (LCST) transitions is analysed.

PACS: 05.40.-a, 64.60.-i, 05.10.Gg

A surprising discovery in the field of stochastic processes is that the interaction between noise and spatial degrees of freedom in extended nonlinear systems can lead, under appropriate conditions, to noise-induced order . This seemingly counterintuitive phenomenon arises when a distributed medium develops some kind of spatiotemporal self-organization which requires noise to exist. The resulting ordered state can take many forms, such as a static spatial pattern , traveling excitable waves , propagating pulses , spatiotemporal chaotic states , or a coherent response to an external periodic forcing , among others. From a fundamental point of view, however, the most interesting example is that of noise–induced ordering phase transitions (NIOPTs) . In this case, transitions between true extended phases (in a thermodynamic sense) are produced when the noise intensity is used as control parameter. These non-equilibrium phase transitions can therefore be characterized using standard techniques of critical phenomena, such as dynamic renormalization group, mean-field approximations, as well as finite-size scaling analyses using, for instance, results obtained from numerical simulations of stochastic partial differential equations. Up to now, all the arguments presented to account for the occurrence of NIOPTs have been dynamical ones. This is mainly because it has been impossible so far to find a model whose nonequilibrium steady-state probability distribution and the associated effective potential are known. In particular, noise-induced phase transitions have been systematically explained in terms of a short-time instability of the local dynamics, which becomes “frozen” at larger times by the spatial coupling .

In this Letter we introduce a class of systems exhibiting NIOPTs for which the steady-state probability distribution can be obtained exactly, so that one can define the corresponding nonequilibrium free energy or potential. As a consequence, the NIOPT can be studied in the steady state, with no dynamical reference. It turns out that the noise–induced phase transition that appears in these model systems is not a consequence of an instability of the disordered homogeneous state. Rather, the situation is similar to what happens for noise–induced transitions in 0–dimensional systems, where the disordered state is linearly stable, but the effective non-equilibrium potential in the steady state is bimodal . A phase transition is obtained by coupling adequately many of those 0-dimensional systems. For a strong enough spatial coupling, the neighboring variables tend to a common value and a macroscopic ordered phase appears as a result of ergodicity breaking. This situation contrasts vividly with NIOPTs reported so far, where the need of a short-time dynamical instability prevents systems that undergo noise-induced transitions in 0-d from exhibiting NIOPTs in the presence of spatial coupling .

The phase transition that we report in what follows shares some features with the so-called Lower Critical Solution Temperature (LCST) transition which appears in some polymer blends . In these very special systems, one can go from a disordered phase (homogeneous mixture) to an ordered phase (separated phases) by increasing the temperature (inverted phase diagram). This is an equilibrium phenomenology which goes in the opposite direction of what is expected intuitively from statistical physics. Our findings present a possible scenario, in terms of a Langevin description, in which these transitions can take place. To reach our goal we consider two important theoretical ingredients: the presence of a field-dependent kinetic coefficient and the existence of a fluctuation-dissipation relation accordingly.

Let us consider the following generic deterministic model for a real field $`\varphi (\stackrel{}{r},t)`$:

$$\frac{\varphi (\stackrel{}{r},t)}{t}=\mathrm{\Gamma }(\varphi (\stackrel{}{r},t))\frac{\delta V}{\delta \varphi (\stackrel{}{r},t)},$$
(1)
which corresponds to a relaxational flow  in a potential $`V(\{\varphi \})`$ with a field-dependent kinetic coefficient $`\mathrm{\Gamma }(\varphi )`$. Field-dependent kinetic coefficients can appear in coarse-grained, master-equation derivations of macroscopic (deterministic or stochastic) field equations . Their explicit form depends on the assumptions in their derivation, but the general trend is that the coefficient is large in the disordered phase and small in the ordered phase . These coefficients have been introduced in the context of phase separation in binary mixtures, either to explore ad hoc their influence, or by necessity to take into account the effect of a gravitational field . In those cases, the kinetic coefficients are also a function of the spatial derivatives of the field. in order to reflect the conserved nature of the concentration field. Here, on the other hand, we consider for simplicity a non-conserved situation.

When stochastic terms are introduced in this description a fluctuation-dissipation relation is used . Following this approach we define that the Langevin equation (in the Stratonovich interpretation) corresponding to the deterministic equation (1) to be:

$$\frac{\varphi (\stackrel{}{r},t)}{t}=\mathrm{\Gamma }(\varphi (\stackrel{}{r},t))\frac{\delta V}{\delta \varphi (\stackrel{}{r},t)}+\mathrm{\Gamma }(\varphi (\stackrel{}{r},t))^{1/2}\xi (\stackrel{}{r},t).$$
(2)
The noise is taken to be Gaussian, with zero mean and correlation

$$\xi (\stackrel{}{r},t)\xi (\stackrel{}{r}^{},t^{})=2\epsilon \delta (\stackrel{}{r}\stackrel{}{r}^{})\delta (tt^{}),$$
(3)
where $`\epsilon `$ is the noise intensity. Under these conditions, the stationary solution of the probability distribution for the field, $`P_{\mathrm{st}}(\{\varphi \})`$ is of the Boltzmann’s type:

$$P_{st}(\{\varphi \})\mathrm{e}^{V_{\mathrm{eff}}/\epsilon },$$
(4)
in terms of an effective potential

$$V_{\mathrm{eff}}(\{\varphi \})V(\{\varphi \})+\frac{\epsilon _0}{2}𝑑\stackrel{}{r}\mathrm{ln}\mathrm{\Gamma }(\varphi (\stackrel{}{r})),$$
(5)
where $`\epsilon _0`$ is a renormalized parameter, proportional to $`\epsilon `$, which includes an ultraviolet cutoff . The previous expression defines the exact nonequilibrium free energy of model (2). We will now show that this class of systems can exhibit a noise-induced phase transition.

So far, the discussion has been made in very general terms. In order to proceed further, we assume that the deterministic free-energy potential has the usual Landau form

$$V(\{\varphi \})=𝑑\stackrel{}{r}\left[V_0(\varphi (\stackrel{}{r}))+\frac{K}{2}|\stackrel{}{}\varphi (\stackrel{}{r})|^2\right],$$
(6)
and we choose a monostable local potential $`V_0(\varphi )=\frac{a}{2}\varphi ^2`$, with $`a>0`$. We also consider a non-conserved order parameter system in which the kinetic coefficient $`\mathrm{\Gamma }`$ is only a function of the field and not of its spatial derivatives. According to the discussion above, we adopt the functional form

$$\mathrm{\Gamma }(\varphi )=1/(1+c\varphi ^2),$$
(7)
which means that fluctuations are larger in dilute regions and smaller in dense ones. In this particular case, and according to Eqs. (4) and (5), the local effective potential can be seen to become bistable for $`\epsilon _0>a/c`$. One can then expect that for a strong enough coupling $`K`$, a true phase transition towards an ordered state controlled by the noise intensity $`\epsilon `$ exists. In order to check this expectation, we perform a Weiss mean-field analysis of the corresponding local effective potential . The analysis is more easily done by considering a spatial discretization in a regular d–dimensional lattice with spacing $`\mathrm{\Delta }x=1`$, i.e. $`\varphi (\stackrel{}{r}_i)\varphi _i`$, with $`i`$ the cell index. The gradient term is then approximated by the sum over nearest neighbors on the lattice in the standard way, $`|\stackrel{}{}\varphi |^2_{j=1}^d(\varphi _j\varphi _i)^2`$. The mean-field approximation is tantamount to replacing the exact value of the neighbors by a mean-field common value $`M`$. In this way, the local effective potential becomes

$$V_{\mathrm{eff}}(\varphi ,M)=V_0(\varphi )+\frac{\epsilon _0}{2}\mathrm{ln}\mathrm{\Gamma }(\varphi )+dK(\varphi M)^2,$$
(8)
and depends on the yet unknown mean–field value $`M`$. This can be computed by using the consistency relation:

$$M=\frac{_{\mathrm{}}^{\mathrm{}}\varphi \mathrm{e}^{V_{\mathrm{eff}}(\varphi ,M)/\epsilon }𝑑\varphi }{_{\mathrm{}}^{\mathrm{}}\mathrm{e}^{V_{\mathrm{eff}}(\varphi ,M)/\epsilon }𝑑\varphi },$$
(9)
which, in general, has to be solved numerically. The solution of this equation is shown as a solid line in Fig. 1 for the two-dimensional case. The result predicts a phase transition from disorder ($`M=0`$) to order ($`M0`$) as the noise intensity increases. This theoretical result has been confirmed by numerical simulations of a discretized version of model (2) in a two-dimensional square lattice. Following standard procedures of equilibrium statistical mechanics, we have characterized the phase transition by means of two quantities: the intensive order parameter $`M=m`$, where $`m=|_i\varphi _i|/L^2`$ is the absolute value of the spatial average of the field, and the generalized susceptibility

$$\chi =L^2\frac{m^2m^2}{\epsilon },$$
(10)
which measures the relative fluctuations of the field in its steady state ($`\epsilon `$ plays the role of $`K_BT`$ in equilibrium phase transitions and $`L`$ is the linear system size). The numerical results are plotted in Figs. 1 and 2, which show the behavior of the order parameter $`M`$ and the susceptibility $`\chi `$, respectively, as the noise intensity $`\epsilon `$ increases, for four different system sizes. The figures clearly denote the existence of a noise–induced ordering phase transition with the usual features of a second–order equilibrium phase transition, in agreement with the theoretical analysis. Monte Carlo simulations of the system making use of the nonequilibrium effective free energy give completely compatible results. The position of the critical point has been computed both by extrapolating to infinite size the location of the maxima of the susceptibility and by means of a finite-size scaling analysis (see below). The two approaches lead to similar results.

As mentioned above, noise–induced phase transitions have been related so far to an instability of the homogeneous phase at short times. This had been clearly manifested by a linear stability analysis of the equation for the first moment, and was understood as a renormalization of the linear forces due to the noise term. Let us now perform in our model a similar linear stability analysis of the homogeneous disordered state $`\varphi =0`$  and we will be able to show that this is not the mechanism that induces the phase transition in our model. The dynamical equation for the first statistical moment, in the linear approximation, reads

$$\frac{\varphi }{t}=(a+\epsilon _0c)\varphi +K^2\varphi .$$
(11)
where Novikov’s theorem has been applied to the term $`\mathrm{\Gamma }(\varphi )^{1/2}\xi `$.

The contribution of the noise is in the same direction as the one from the local potential, and thus the homogeneous state $`\varphi =0`$ is linearly stable. Hence the physical mechanism responsible for the appearance of a new phase must be different from the one of previously known cases of noise induced phase transitions. A simple interpretation can be found in terms of a balance between the role of the deterministic monostable local potential, which tends to take the system towards the disordered ($`\varphi =0`$) phase, and the stochastic motion, which is due to the fact that fluctuations are more intense in the disordered phase, and thus push the system away from it. The presence of a spatial coupling $`K`$ helps to break the symmetry of the homogeneous state.

A natural question arises concerning the universality class of the noise-induced phase transition described above. The availability of an exact nonequilibrium potential for this model, given by Eq. (5), allows us in a very simple way to analyze its universal behavior at the critical point. If an expansion of $`\mathrm{ln}\mathrm{\Gamma }(\varphi )`$, as chosen in Eq. (7), is performed on (5), the nonequilibrium effective potential can be seen to have a Ginzburg-Landau form with infinite higher-order (irrelevant) terms. Since the universality class of the Ginzburg-Landau model is that of the Ising model, we predict that this transition will belong to the d-dimensional Ising universality class. This is indeed confirmed by a finite-size scaling analysis of the numerical results. We make use of the following scaling laws, generalized from equilibrium statistical mechanics :

$`M`$ $`=`$ $`L^{\beta /\nu }\stackrel{~}{M}(\stackrel{~}{\epsilon }L^{1/\nu }),`$ (12)
$`\chi `$ $`=`$ $`L^{\gamma /\nu }\stackrel{~}{\chi }(\stackrel{~}{\epsilon }L^{1/\nu })`$ (13)
where $`\stackrel{~}{M}`$ and $`\stackrel{~}{\chi }`$ are the corresponding scaling functions, and $`\stackrel{~}{\epsilon }`$ is the reduced control parameter, $`\stackrel{~}{\epsilon }=(1\epsilon /\epsilon _c)`$. The values of the critical exponents $`\beta `$, $`\gamma `$ and $`\nu `$ define the universality class. We have chosen the corresponding values of the 2-d Ising universality class, and plotted expressions (12) and (13) in Figs. 3 and 4. All data are seen to scale fairly well into the same curve in the region close to the critical point. As expected, those data corresponding to the smallest system size ($`L=16`$) show the major discrepancies. Hence, we can conclude that this new kind of noise-induced phase transitions belongs to the 2d-Ising universality class.

In conclusion, we have introduced a dynamical system that exhibits a nonequilibrium phase transition from disorder to order as the intensity of the noise is increased. As a consequence, the system is disordered below a certain amount of noise and can be considered a possible scenario, within the context of stochastic differential equations, for the appearance of a lower critical temperature in physical systems. At variance with other noise-induced nonequilibrium phase transitions reported before, the physical mechanism inducing the transition is not a noise–induced short time instability of the disordered state. Rather, we are able to show that the effective local potential, which can be obtained analytically, exhibits bistability as the intensity of the noise increases, whereas the deterministic one is monostable. Hence, the mechanism is seen to be given by a balance between deterministic forces and stochastic motion. Our model shows that it is possible to obtain a noise–induced phase transition by coupling adequately 0–dimensional systems undergoing a noise–induced transition. We have argued that this type of coupling is consistent with a mesoscopic description in terms of fields, kinetic coefficients and the fluctuation–dissipation relation. Our findings have been sustained by numerical simulations in a two-dimensional lattice model. Finally, by studying the exact stationary probability distribution, we are able to show that the model belongs to the same universality class as the Ising model.

We acknowledge the Dirección General de Enseñanza Superior (Spain) for financial supports under projects PB94-1167, PB96-0241, PB97-0141-C02-01, and PB98-0935.