In this paper, we study phase transitions for weakly interacting multiagent systems. By investigating the linear response of a system composed of a finite number of agents, we are able to probe the emergence in the thermodynamic limit of a singular behavior of the susceptibility. We find clear evidence of the loss of analyticity due to a pole crossing the real axis of frequencies. Such behavior has a degree of universality, as it does not depend on either the applied forcing or on the considered observable. We present results relevant for both equilibrium and nonequilibrium phase transitions by studying the Desai–Zwanzig and Bonilla–Casado–Morillo models.

Multiagent models feature in a very vast range of applications in natural sciences, social sciences, and engineering. We study here the Desai–Zwanzig (DZ) and Bonilla–Casado–Morillo (BCM) models, which are paradigmatic for equilibrium and nonequilibrium conditions, respectively. Phase transitions result from the coordination between the individual agents and are associated with the divergence of the linear response of the system. The occurrence of phase transitions is universal: it does not depend on the acting forcing and can be detected by looking at virtually any observable of the system. We showcase here how response theory is capable of providing a useful angle for understanding the universal properties of phase transitions in complex systems.

## I. INTRODUCTION

Agent based models are regularly employed to model various phenomena in the natural sciences, social sciences, and engineering.^{1,2} Multiagent systems are used to model diverse phenomena such as cooperation,^{3} synchronization,^{4} systemic risk,^{5} and consensus formation.^{6,7} They are fundamental in developing algorithms for sampling and optimization^{8} and they have also been used for the management of natural hazard^{9} and climate change impact.^{10}

Multiagent systems can often exhibit abrupt changes in their behavior, often corresponding to critical transitions that occur when a parameter, e.g., interaction strength or temperature, passes a certain threshold. Such transitions are often associated to cataclysmic events such as climate change, market crashes, etc.^{11,12} The importance of developing tools for predicting critical transitions has long been recognized. One of the main tools used in order to develop early warning signals for critical transitions is that of linear response theory.

Following the seminal contribution by Kubo,^{13} linear response theory represents a very powerful framework for studying the properties of statistical mechanical systems by investigating how they respond to external perturbations.^{14–16} Linear response theory has been successfully applied to classic problems of solid-state physics and optics^{17} as well as plasma physics and stellar dynamics (Chap. 5 of Ref. 18); see some examples of the application of the theory in both equilibrium and nonequilibrium systems.^{19–24} Rigorous mathematical foundations for linear response theory have been provided for the case of Axiom A systems^{25,26} (see, e.g., Ref. 27 for further developments in the context of deterministic systems) and for diffusion processes, both in finite and in infinite dimensions^{28,29} (see also the interesting contribution^{30} that bridges the deterministic and the stochastic viewpoints).

Critical transitions arise when the spectral gap of the transfer operator of the unperturbed system shrinks to zero^{31–34} as the Ruelle–Pollicott poles^{35,36} touch the real axis. Near criticality, the negative feedbacks of the system become increasingly ineffective, resulting in arbitrarily large, usually non-Gaussian, fluctuations and a divergence of correlation properties of the system.^{3,37,38}

In the thermodynamic limit, the multiagent system can also undergo a qualitative change of its properties through a different mathematical mechanism, namely, phase transitions,^{3,37} defined as the exchange of stability of nonunique stationary distributions as the parameters of the systems vary (see a detailed analysis in Ref. 39).

In a previous paper,^{40} we derived linear response formulas for a system of weakly interacting diffusions described by an $N$-particle Fokker–Planck equation and have explicitly identified two qualitatively different scenarios for the breakdown of the linear response, associated with the previously mentioned critical transitions and phase transitions. We focus here on the latter case.

Phase transitions are a genuine thermodynamic phenomenon, where the divergence of the response stems from the coordination taking place, in suitable conditions, because of the coupling between the infinite number of agents composing the total system. The coupling among the subsystems results in a memory effect that leads to obtaining the macroscopic response function of the system as a renormalized version of its microscopic counterpart,^{40} with formal similarities with the well-known Clausius–Mossotti relation.^{17,41,42} The role of memory in determining criticality due to endogenous processes has been emphasized in Refs. 12 and 43. The link between phase transitions and slow decay of correlations for interacting particle systems is well established (see, e.g., Ref. 44).

### A. This paper

In this paper, we focus in much greater detail on the relationship between the occurrence of phase transitions and the non-analyticity of the susceptibility of the system describing the frequency-dependent response of an observable to a given perturbation in the upper complex frequency plane. The singularity manifests itself as a pole that crosses the real axis of the frequency variable. We use a formalism that mirrors spectroscopic techniques that are used for investigating the frequency dependence of the optical properties of materials.^{17} By studying how the real and imaginary parts of the susceptibility of the systems depend on the number of agents, we are able to predict the position of the pole and the associated residue, which describe the emergence of the singularity in the thermodynamic limit. We verify that the position of the pole depends on the considered model, but, instead, that for a given model the loss of analyticity depends neither on the choice of the observable nor on the applied perturbation and is, in this sense, a universal feature of the system. Our numerical investigations are performed on the Desai–Zwanzig (DZ)^{45} and the Bonilla–Casado–Morillo (BCM)^{46} models. The DZ model exhibits a paradigmatic example of an equilibrium order–disorder phase transition, analogous to the Ising ferromagnetic transition,^{3,37} while the BCM model describes an out-of-equilibrium synchronization transition of an infinite collection of coupled nonlinear oscillators. As the transition point is crossed, the order parameter (magnetization) acquires a nonvanishing constant value for the DZ model and is periodically oscillating for the BCM model.

## II. THE GENERAL FRAMEWORK

We investigate a system composed of $N$ exchangeable interacting $M$ dimensional sub-systems whose dynamics is determined by the following Itô stochastic differential equation (SDE):

where $xk\u2208RM$ and $k=1,\u2026,N$. $F\alpha :RM\u2192RM$ is a smooth vector field, possibly depending on a parameter $\alpha $, and $W$ denotes a standard $P$-dimensional Brownian motion; $S:RM\u2192RM\xd7P$ is the volatility matrix and the parameter $\sigma >0$ controls the strength of the stochastic forcing, i.e., plays the role of the temperature. We consider a fully coupled system given by the quadratic (Curie–Weiss) interaction potential $U(y)=|y|22$. In this case, the order parameter is known and it is given by the first moment/magnetization.

The coefficient $\theta $ modulates the intensity of the coupling, which attempts at synchronizing all systems by attracting them to the center of mass. In the thermodynamic limit $N\u2192\u221e$, the one-particle distribution function converges to the distribution $\rho (x,t)$ that satisfies a nonlinear and nonlocal Fokker–Planck equation^{3,47–49}

where $D=SST$. This mean field partial differential equation might support multiple coexisting stationary measures at low temperatures/large interaction strengths. In particular, in a conservative system described by a confining potential $F\alpha (y)=\u2212\u2207V\alpha (y)$ with additive noise such that $S$ is the identity matrix, stationary solutions of Eq. (2) correspond to local minima of $V\alpha (y)$. In this case, the thermodynamic limit (2) can be written in the standard form as $limN\u2192+\u221e\u22121Nlog\u2061ZN=infF(\rho )$, where $ZN$ denotes the partition function of the $N$-particle system and $F(\rho )$ denotes the free energy of mean field system.^{50} A stationary state is characterized by the order parameter $\u27e8x\u27e90$ and the associated stationary distribution $\rho 0(x)$.

We now perturb the stationary state by setting $F\alpha (x)\u2192F\alpha (x)+\epsilon X(x)T(t)$ and we study the response of the system by expanding the distribution function as $\rho (x,t)=\rho 0(x)+\epsilon \rho 1(x,t)+O(\epsilon 2)$. Following a tedious calculation reported in Ref. 40, the response of the order parameter in the frequency domain is written in terms of a macroscopic (or renormalized) susceptibility $\Gamma ~i(\omega )$ as (repeated indices are summed),

where $Pij(\omega )=\delta ij\u2212\theta Yij(\omega )$ and the susceptibilities $\Gamma i(\omega )$ and $Yij(\omega )$ are, respectively, the Fourier transform of the microscopic response functions that can be written as correlation functions in the unperturbed state as^{40}

where $L\u27e8x\u27e90+$ represents the adjoint of $L\u27e8x\u27e90$ and $\u27e8\u22c5\u27e90$ is the expectation value on the unperturbed state $\rho 0(x)$. The Fokker–Planck operator $L\u27e8x\u27e90$ appears on the right hand side of (2) [evaluated at the stationary state $\rho 0(x)$] and its adjoint $L\u27e8x\u27e90+$ can be interpreted as the generator of the Koopman operator of the stationary dynamics.^{40,51} As such, correlation properties of the system are related to the spectrum of $L\u27e8x\u27e90+$.^{32,52,53} The renormalization of the susceptibility derives from the coupling among the subsystems; note that $\Gamma ~i(\omega )$ inherits the poles of both $\Gamma i(\omega )$ and of the matrix $Pij\u22121(\omega )$. Away from criticality, both the microscopic and macroscopic susceptibilities are analytic in the upper half of the $\omega $ complex plane.

As discussed in Ref. 40, the critical behavior of this class of multiagent systems, signified by the singular behavior of the susceptibility, originates from two distinct physical phenomena that are associated to either the poles of $\Gamma i(\omega )$ or $Pij\u22121(\omega )$.

The case where $\Gamma i(\omega )$ diverges pertains to the occurrence of critical transitions. It is of interest here the case where poles appear in the real $\omega $ axis for $\Gamma ~i(\omega )$ because the matrix $Pij\u22121(\omega )$ becomes singular. This corresponds to phase transitions originated by the coupling and do not show a divergence of correlation properties, because $\Gamma i(\omega )$ is, instead, analytic in the upper complex $\omega $ plane. Equivalently, the spectral gap of $L\u27e8x\u27e90+$ remains finite at a phase transition. However, at the transition point, the usual dispersion relations need to be modified.^{17,40} The conditions underpinning the breakdown of linear response theory do not depend on the perturbation field $X(x)$ or on the choice of the observable ($x$, in our case) and can be related to the spectral properties of a modified transfer operator.^{40}

## III. NUMERICAL RESULTS

Below, we present results for the Desai–Zwanzig (DZ)^{45} and the Bonilla–Casado–Morillo (BCM)^{46} models. These are composed by $N$ interacting agents evolving according to Eq. (1). Each individual agent of the DZ (BCM) model evolves in $R$ ($R2$). A detailed description of the two models is presented in the Appendix. We repeat our experiments for various choices of $N$ in order to detect the emergence of singularities for the combination of the parameters corresponding to phase transitions. Here, we keep fixed the values of the internal parameter $\alpha $. Both models undergo a phase transition at the transition line $\sigma ~=\sigma (\theta ;\alpha )$ in the parameter space $(\sigma ,\theta )$ (see the Appendix for the analytical evaluation of the transition line). Since one of the two parameters is redundant, we fix the coupling intensity $\theta $ and we vary, instead, the noise strength $\sigma $.

Following Ref. 54, we perform $n$ simulations where the initial conditions are chosen according to the unperturbed invariant measure $\rho 0(x)$ and where at time $\tau =0$, we apply a perturbation proportional to a Dirac $\delta $ function. The average of the response for the observable $x$ over the $n$ simulations gives an estimate $G~i(\tau ;N)$ of the renormalized response function. Details on the numerical simulations are also reported in the Appendix. Figure 1 shows the response functions $G~i(\tau ;N)$ for an additive perturbation $X(x)=1$ for the DZ model (left panel) and $X(x)=(1,0)$ for the BCM model (right panel). The two response functions are qualitatively different because, by and large, the one for the DZ model describes a monotonic decay, whereby the system relaxes toward the unperturbed state, while the one for the BCM combines the decay with an oscillatory behavior taking place at the natural frequency $\omega ~=1$.

In the DZ model, the response functions initially undergo a fast and substantial decay, both far from and at the phase transition, associated with a time scale of order $1$. However, at the phase transition, a new, much longer, timescale appears. This timescale increases monotonically with $N$. The same is observed in the case of the BCM model if one considers the envelope of the response function rather than the response function itself: at the transition, the decay of the oscillations becomes slower and slower as $N$ increases.

The origin of the new timescales resides in the appearance of simple pole at $\omega =\omega 0$ in the susceptibility $\Gamma ~i(\omega ;N)$, the Fourier transform of the response function. The pole is located at $\omega 0=0$ for the DZ model and at $\omega 0=\omega ~=1$ for the BCM model. When considering finite values of $N$, the susceptibilities describing the response of (virtually) any observable to (virtually) any external perturbation have a contribution of the form $\kappa \omega \u2212\omega 0+i\gamma (N)$, where $\gamma (N)\u21920+$ as $N\u2192+\u221e$ and $\kappa $ represents the residue of the pole, because $limN\u2192\u221e\kappa \omega \u2212\omega 0+i\gamma (N)=\u2212i\pi \kappa \delta (\omega \u2212\omega 0)+\kappa P(1/(\omega \u2212\omega 0))$. The quantity $\kappa \u2208C$ depends on the choice of observable and of the perturbation. We remark that the asymptotic property does not depend on how fast the function $\gamma (N)$ vanishes for increasing values of $N$. Following Ref. 38, one might conjecture that for the Desai–Zwanzig model and related models, the function $\gamma $ would scale as $1/N$. However, we have observed here that the behavior of $\gamma $ is different. This is an issue of fundamental importance that we will explore in future work, also in the case of nonequilibrium systems. Note also that, in the case of equilibrium systems, the mean field limit $N\u2192\u221e$ and the limit $T\u2192Tc$, where $Tc$ is the critical temperature, do not commute.^{55}

We next investigate the phase transitions by looking at the properties of the susceptibilities (see Fig. 2). When $\sigma \u2260\sigma ~$, the susceptibilities do not show any singularity or any remarkable dependence on $N$, thus indicating that the thermodynamic limit has been reached to a good approximation. As $N$ increases, for both the DZ model (left panel) and the BCM model (right panel), the resonance at $\omega =\omega 0$ of the real part of the susceptibility approaches the limiting Dirac function $\pi k\delta (\omega \u2212\omega 0)$ with coefficient $k>0$. This singular behavior is clear from the plot of the primitive function of the real part of the susceptibility (bottom inset) that tends to step function. For both models, $\kappa =ik$ is an imaginary number. Indeed, the imaginary part of the susceptibility behaves exactly as the Cauchy principal value distribution and can be used to get easily a quantitative estimate of $k$. The top insets of Fig. 2 show the function $(\omega \u2212\omega 0)Im{\Gamma ~i(\omega )}$. As $N\u2192\u221e$, this function converges to $k$ everywhere except for $\omega =\omega 0$.

An explicit expression for $k$ is known^{40} in the case of the DZ model:^{56} $k=\u27e8x2\u27e90\theta \u222b0+\u221edt\u27e8x(t)x(0)\u27e90$. Using the statistics of the unperturbed runs, we obtain $k\u22480.89$, which agrees within 2% with the one obtained from the limiting behavior of the susceptibility, thus validating our results. In the case of the BCM model, our procedure allows one to derive a direct estimate $k\u22480.44$; in this case, no expression for the residue is available in the literature, and following Refs. 46 and 40, its evaluation seems cumbersome. We here observe that, by evaluating the susceptibility for finite values of $N$, we are able to predict the residue of the pole at $\omega =\omega 0$, which appears, instead, only in the thermodynamic limit. The residue plays the role of a latent heat of phase change in classical thermodynamics. Our results, though, allow one to deal with the case of a dynamical latent heat, which is observed for perturbations occurring in a non-vanishing frequency. As discussed earlier, the singular behavior of the susceptibility has some degree of universality. By this, we mean that while for a given model, the value of the residue is forcing- and observable-dependent, its position is a fundamental property of the model itself (see the Appendix for additional examples).

## IV. CONCLUSIONS

The study of how a large network of identical agents responds to exogenous perturbations is of the uttermost importance in different areas of science. One might be interested not only in the smooth response of the system, where its properties change ever so slightly but also in the critical, nonsmooth regime, where small perturbations can lead to large and possibly undesired changes. Multiagent systems modeled as weakly interacting Itô diffusions represent a rich class of models exhibiting such critical behavior, and for which rigorous analysis and careful numerical investigations can be carried out. Usually, these phenomena are accompanied by a large spatial (among sub-systems) and temporal restructuring of the system where correlations get highly magnified. The critical behavior due to the emergence of a phase transition is a genuine thermodynamic phenomenon arising from the complex interactions among the infinite number of agents. Nevertheless, we have shown in this paper that linear response theory provides a powerful framework for detecting and anticipating phase transitions by investigating the response of a finite particle system to external perturbations. We have been able to predict the appearance of poles in the susceptibility, which describes the frequency-dependent response of the system, as well as to obtain a correct estimate of critical thermodynamic properties, such as the residue of the poles, based on the knowledge of the response for the finite particle system in two paradigmatic models describing equilibrium and nonequilibrium phase transitions. This is an encouraging starting point for improving our ability to understand and predict transitions in more complex multiagent systems.

## ACKNOWLEDGMENTS

V.L. acknowledges support received from the European Union’s Horizon 2020 program through the project TiPES (Grant Agreement No. 820970) and from the EPSRC (Project No. EP/T018178/1). The work of G.A.P. was partially funded by the EPSRC (Grant No. EP/P031587/1) and by J.P. Morgan Chase & Co. Any views or opinions expressed herein are solely those of the authors listed, and may differ from the views and opinions expressed by J.P. Morgan Chase & Co. or its affiliates. This material is not a product of the Research Department of J.P. Morgan Securities LLC. This material does not constitute a solicitation or offer in any jurisdiction. N.Z. was supported by an EPSRC studentship as part of the Centre for Doctoral Training in Mathematics of Planet Earth (Grant No. EP/L016613/1).

## DATA AVAILABILITY

The data that support the findings of this study are openly available in “Spectroscopy of Phase Transitions” at https://figshare.com/projects/Spectroscopy_of_phase_transitions/101846, Ref. 57.

### APPENDIX: THE MODELS

Weakly interacting diffusions represent a rich class of agent based models, describing a network of interacting subsystems. The local dynamics of each subsystem is determined by a smooth vector field $F\alpha (x)$. The local force is, in general, nonconservative, leading to irreversible and dissipative processes that can exhibit complex behaviors, such as deterministic chaos (see Fig. 3).

An all-to-all coupling between the subsystems is given by a matrix $Lij=\u2207U(xi\u2212xj)$, where $U(x)=U(\u2212x)$ represents the interaction potential and $xi$ is the state vector of the $i$th subsystem. Weakly interacting diffusions are characterized by a coupling strength that is inversely proportional to the number of subsystems $N$. As $N$ increases, the interaction structure gets more and more intricate, while the intensity becomes weaker and weaker (see Fig. 3). As mentioned in the main text, the DZ model has been introduced, and thereafter commonly used, as a paradigmatic example of an equilibrium continuous phase transition reminiscent of the Ising-like ferromagnetic transitions in spin systems.^{3,37} The DZ model describes a network of one-dimensional subsystems $xk$ whose dynamics is prescribed by the following equations (see main text for notation):

where $k=1,\u2026,N$ and the confining potential $V\alpha (x)=\u2212\alpha 2x2+x44$ has a double well shape for $\alpha >0$. Without loss of generality, we here consider $\alpha =1$. In the absence of coupling, $\theta =0$, the above equations describe the simple motion of a particle in a double well potential, subject to additive noise. The presence of the coupling allows for a long range coordination of the system that in the thermodynamic limit $N\u2192+\u221e$ results in a proper phase transition. In this regime, by varying the parameters $(\alpha ,\theta )$, the order parameter $\u27e8x\u27e9$ undergoes a continuous order–disorder transition, similar to the pitchfork bifurcation diagram for the Ising model. It is possible to show^{3} that the critical line is given by $D\u22123/2(\theta \u22121\sigma )D\u22121/2(\theta \u22121\sigma )=\sigma \theta $, where $D\nu (z)$ is a parabolic cylinder function. Here, the coupling $\theta $ is kept fixed ($\theta =0.55$) and we vary $\sigma $ to approach the transition point.

The BCM model describes an ensemble of bi-dimensional nonlinear oscillators undergoing an out of equilibrium self-synchronization transition. The time evolution of the network of oscillators is given by the following equations:

where the local force is not conservative, giving rise to the nonequilibrium features of the system, and reads $F\alpha (x)=(\alpha \u2212|x|2)x+x+$ where $x+=(\u2212x2,x1)$. The latter term corresponds to a rotation and makes the stationary state a nonequilibrium one. The parameter $\alpha >0$ controls the amplitude of the oscillations of the individual nonlinear oscillators. In fact, when $\theta =\sigma =0$, each subsystem oscillates as $xj(t)=\alpha (cos\u2061(t+\beta j),sin\u2061(t+\beta j))$, where $\beta j=tan\u2061(x2j(0)/x1j(0))$. The coupling tries to synchronize the subsystems by attracting them toward the center of mass $1N\u2211j=1Nxj$. In the thermodynamic limit and for sufficiently low values of the noise, the order parameter $\u27e8x\u27e9=\u27e8x\u27e9(t)$ exhibits a periodic time evolution, resulting from the subsystem oscillating in a coherent way. On the other hand, high values of the noise correspond to a nonsynchronized state where the order parameter vanishes. In particular, the transition happens at the surface of the $(\alpha ,\theta ,\sigma )$ parametric space defined by^{46}

where $A=\alpha \theta \u22121$ and $\delta =2\sigma 2\theta $. In the following, we have $\theta =\alpha =2$. The color code for the figures is given by

noncritical black: DZ $\sigma \u22481$ , BCM $\sigma \u22482$,

noncritical blue: DZ $\sigma \u22480.87$ , BCM $\sigma \u22481.8$, and

critical red: DZ $\sigma ~\u22480.75$ , BCM $\sigma ~\u22481.59$.

#### 1. Numerical linear response experiments

As mentioned in the main text, we perform $n$ simulations of the response given by Eqs. (A1) and (A2) where the initial conditions are chosen according to the respective unperturbed invariant measure $\rho 0(x)$ and where at time $\tau =0$, we apply a perturbation proportional to Dirac’s $\delta (\tau =0)$. The average of the response for the observable $x$ over the $n$ simulations gives an estimate of $G~i(\tau ;N)$. The response functions away from the transitions are estimated on an ensemble of $n=105$ simulations, while the critical response functions with $n=106$ for DZ and $n=7\xd7106$ for BCM. Furthermore, we investigate the response up to time $\tau =5\xd7103$. The corresponding susceptibility $\Gamma ~i(\omega ;N)$ is simply defined as the Fourier transform of $G~i(\tau ;N)$. In the main text, we show the results for an additive perturbation $X(x)\u2261X$. However, the critical behavior of the response does not depend on the type of perturbation, modulo, a potential degenerate class of perturbations that have zero projection on the invariant measure $\rho 0(x)$. We have thus decided to investigate the response of the DZ model for a spatially dependent perturbation $X(x)=x2$ (see Fig. 4). The response function $G~i(\tau ;N)$, both away and at the phase transition, has a rapid initial decay with a timescale that is different from the response function shown in the main text. As a matter of fact, the timescale associated to the dominant mode of the response function for $\tau \u21920$ does, in general, depend on the applied perturbation.^{40} As expected, the response function at the phase transition develops a much longer timescale that increases as the number of particles increases. A more accurate comparison with the result shown in the main text can only be performed in the frequency domain. Figure 4 (right panel) shows that, away from the transition, the susceptibilities have a smooth behavior and no evident dependence on $N$. At the phase transition, the susceptibility develops the expected singular behavior $\kappa \omega \u2212\omega 0+i\gamma (N)$, where $\gamma (N)\u21920+$ as $N\u2192+\u221e$ due to the appearance of a simple pole $\omega 0=0$. The residue $\kappa $ is purely imaginary and its magnitude $k$ can be inferred by visual inspection of the top inset representing the function$(\omega \u2212\omega 0)Im{\Gamma ~(\omega )}$ to be just less than 0.29. As mentioned in the main text, the residue depends both on the observable and on the perturbation $X(x)$.

## REFERENCES

*Interacting Multiagent Systems: Kinetic Equations and Monte Carlo Methods*(Oxford University Press, 2013).

*Critical Transitions in Nature and Society*, Princeton Studies in Complexity (Princeton University Press, Princeton, 2009).

*Extreme Events in Nature and Society*, edited by S. Albeverio, V. Jentsch, and H. Kantz (Springer, Berlin, Heidelberg, 2006) pp. 95–119.

*ICM Seoul 2014, Proceedings*(KYUNG MOON SA Co. Ltd., 2014), Vol. III, pp. 525–545.

*Ecole d’Eté de Probabilités de Saint-Flour XIX—1989*, Lecture Notes in Mathematics Vol. 1464, edited by P. L. Hennequin (Springer, Berlin, 1989).

*Semiclassical Analysis, Witten Laplacians, and Statistical Mechanics*, Series in Partial Differential Equations and Applications Vol. 1 (World Scientific Publishing Co., Inc., River Edge, NJ, 2002), p. x+179.

*Chaos, Fractals, and Noise*, 2nd ed., Applied Mathematical Sciences Vol. 97 (Springer-Verlag, New York, 1994), p. xiv+472.