ARTICLE IN PRESS
Physica A 376 (2007) 199–207 www.elsevier.com/locate/physa
Optimal range of a nonlocal kernel for stochastic resonance in extended systems B. von Haeftena, H.S. Wiob, a
Departamento de Fı´sica, FCEyN, Universidad Nacional de Mar del Plata, Dea´n Funes 3350, 7600 Mar del Plata, Argentina b Instituto de Fı´sica de Cantabria, Universidad de Cantabria and CSIC, E-39005 Santander, Spain Received 22 September 2006; received in revised form 23 October 2006 Available online 21 November 2006
Abstract Here we study stochastic resonance (SR) in a spatially extended system described by a reaction–diffusion equation for a scalar (activator-like) field including a nonlocal contribution. We assume that such a contribution arises from an effective adiabatic elimination of an auxiliary (inhibitor-like) field. Our aim is to study the role played by the range of the nonlocal kernel on the SR phenomena. We have found that increasing the nonlocal coupling reduces the system’s response and that, similar to the so-called system size SR, there is an ‘‘optimal’’ value of the kernel’s range, yielding a maximum in the system’s response. r 2006 Elsevier B.V. All rights reserved. Keywords: Stochastic resonance; Reaction–diffusion systems; Lyapunov functional; Nonlocal interaction
1. Introduction The phenomenon of stochastic resonance (SR) has attracted the attention of an increasing number of researchers for more than 20 years. The range of phenomena for which this mechanism can offer an explanation is very broad, and can be appreciated in recent reviews [1,2] as well as in the references therein. Many of the phenomena that could be described within an SR framework occur in extended systems [3,4], for example, experiments carried out to explore the role of SR in sensory and other biological functions [5] or in chemical systems [6]. These, together with the possible technological applications, were the main motivation to many recent studies showing the possibility of achieving an enhancement of the system’s response by means of the coupling of several units in what conforms an extended medium [3,4,7–11], as well as analyzing the possibility of making the system’s response less dependent on a fine tuning of the noise intensity, as well as studying different ways to control the phenomenon [12,13]. Regarding the problematic of SR in extended systems, in a series of papers [7–11,19], we have studied such a phenomenon for the transitions between two different patterns, exploiting the concept of nonequilibrium potential (NEP) [14–17]. The NEP is a special Lyapunov functional of the associated deterministic system Corresponding author.
E-mail address:
[email protected] (H.S. Wio). 0378-4371/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2006.10.091
ARTICLE IN PRESS B. von Haeften, H.S. Wio / Physica A 376 (2007) 199–207
200
which for nonequilibrium systems plays a role similar to that played by a thermodynamic potential in equilibrium thermodynamics [14]. Such an NEP, closely related to the stationary solution of the system’s Fokker–Planck equation, characterizes the global properties of the dynamics, that is, attractors, linear and relative stability of these attractors, height of the barriers separating attraction basins and, in addition, it allows us to evaluate the transition rates among the different attractors [14–18]. The SR between the attractors was theoretically investigated in terms of the two-state approximation [1]. It was shown that the knowledge of the NEP allows us to obtain a rather complete picture of the behavior of the output signal-to-noise ratio (SNR). In another recent paper we have explored the characteristics of this SR phenomenon in an extended system composed by an ensemble of noise-induced nonlinear oscillators coupled by a nonhomogeneous, densitydependent diffusion, externally forced and perturbed by a multiplicative noise, showing an effective noiseinduced bistable dynamics [18]. Also, in Ref. [19] we have reviewed the results of these studies of SR in extended systems done exploiting the knowledge of the NEP. Other recent studies refer to the so-called system size SR. For instance, studies in biological models of the Hodgkin–Huxley type [20,21] have shown that ion concentrations along biological cell membranes present intrinsic SR-like phenomena as the number of ion channels is varied; that the regularity of the collective firing of a set of coupled excitable FitzHugh–Nagumo units—even in the absence of external forcing—results optimal for a given value of the number of elements [22]; and similar phenomena found in a set of globally coupled units described by a f4 theory [23] and even in an opinion formation model [24]. As system size SR arises in extended systems, it was natural to try to include it within an NEP approach, a fact studied in very recent papers where we have shown that this is the most natural framework to understand such a phenomenon [25,26]. In this work we analyze an activator–inhibitor-like system, where the inhibitor-like field has been adiabatically eliminated yielding an effective, nonlocal, one-field reaction–diffusion equation. Our study is based on a previous one [16], where the nonlocal kernel Gðx; x0 Þ comes from the adiabatic elimination of a fast inhibitor with a diffusive transport behavior. Here we change the form of the kernel, arguing that the transport mechanism of the adiabatically eliminated inhibitor-like field has a nondiffusive character, but is such that it yields a kernel more localized in space and with a controllable interaction range. Exploiting the knowledge of the NEP for this system, we have studied the dependence of the SR phenomenon on the range (or nonlocality degree) of the kernel entering into the nonlocal term. We have analyzed the comparative weight of the local (diffusive) and nonlocal terms in contributing to the SR response, and have also obtained that the range of such nonlocal kernel has an optimum value yielding a maximum for the SNR. The paper is organized as follows. In Section 2 we present the model and the formalism to be used. After that we discuss in Section 3 the SR between the stationary solutions, and in Section 4 we show some results. Finally, in Section 5 we present some conclusions. 2. The model 2.1. The new kernel In Ref. [16] the studied system was quðx; tÞ q2 uðx; tÞ ¼ Du þ f ðuðx; tÞÞ vðx; tÞ, qt qx2
(1)
qvðx; tÞ q2 vðx; tÞ ¼ Dv þ buðx; tÞ gvðx; tÞ, (2) qt qx2 where Du and Dv are the diffusion coefficients of the activator uðx; tÞ and the inhibitor vðx; tÞ, respectively. Both u and v are real fields, and f ðuÞ is a nonlinear (autocatalytic) term. The system is confined to the domain x 2 ½L; L, with Dirichlet boundary conditions at both boundaries ðuðL; tÞ ¼ vðL; tÞ ¼ 0Þ. The timescale associated with each field enters through the relation Z ¼ tv =tu . This allows us to perform an adiabatic approximation and obtain a particular form of the NEP for this system. Assuming Z51, we can Z
ARTICLE IN PRESS B. von Haeften, H.S. Wio / Physica A 376 (2007) 199–207
201
eliminate the inhibitor (which is now slaved to the activator) by solving the second equation using the Green function method. This slaving procedure reduces our system to a nonlocal equation for the activator only, which has the form Z L quðx; tÞ q2 uðx; tÞ ¼ Du þ f ðuÞ b Gðx; x0 Þuðx0 Þ dx0 . (3) qt qx2 L From this equation, and taking into account the symmetry of the Green function Gðx; x0 Þ, we can obtain the Lyapunov functional for this system, which has the form ) Z Z L ( 2 Z u Du qu b L 0 0 0 F½u ¼ dx f ðwÞ dw þ dx Gðx; x Þuðx ÞuðxÞ . (4) 2 L 2 qx L This spatial nonlocal term in the NEP takes into account the repulsion between activated zones. When two activated zones come near each other, the exponential tails of the inhibitor concentration overlap, increasing its concentration between both activated zones and creating an effective repulsion between them, the Green function playing the role of an exponential screening between the activated zones. In Ref. [8] the knowledge of such an NEP was exploited to study SR on the system indicated by Eq. (3). According to the discussion in the Introduction, the starting point of our analysis will be the effective, nonlocal and stochastic, reaction–diffusion equation for the real (activator-like) field fðx; tÞ, analogous to Eq. (3), defined in the one-dimensional domain x 2 ½L; L by Z L qf q2 f ¼ D 2 þ f ðfÞ b Kðx; x0 Þfðx0 Þ dx0 þ xðx; tÞ, (5) qt qx L where the diffusivity D is constant, and we assume a nonlinear cubic term f ðfÞ ¼ fðf bÞð2 fÞ. Here xðx; tÞ is an additive Gaussian white noise, of zero mean value and correlation hxðx; tÞxðx0 ; t0 Þi ¼ 2dðt t0 Þdðx x0 Þ, where is the noise intensity. We also assume that the system is subject to Dirichlet boundary conditions, i.e., fðL; tÞ ¼ 0. Finally, Kðx; x0 Þ is the new nonlocal kernel to be defined as follows. Similar to Eq. (3), this system could be written in a variational form, with the functional F½f given by Eq. (4). As indicated above, in Ref. [16] the kernel Gðx; x0 Þ comes from the adiabatic elimination of a fast inhibitor with a diffusive behavior. Here we argue that the transport mechanism of the adiabatically eliminated inhibitor-like field is different from diffusion, but still is one that through an adiabatical elimination yields a nonlocal contribution with a kernel. Moreover, we assume that such a kernel is more compact in space and has a controllable interaction range. In order to keep our analysis simple, we propose the following form: 8 < 1 if jx x0 jpl; Kðx; x0 Þ ¼ 2l (6) : 0 if jx x0 j4l; and we can do the analysis just varying the interaction range 2l. With the indicated form of the kernel, the new equation adopts the form Z qf q2 f b xþl ¼ D 2 þ f ðfÞ fðx0 Þ dx0 þ xðx; tÞ, qt qx 2l xl and the NEP can be written in a functional form similar to Eq. (4): ) Z L ( 2 Z f Z D qf b L 0 0 0 F½f ¼ dx f ðjÞ dj þ dx fðx ÞKðx; x ÞfðxÞ 2 qx 2 L L 0 ) Z L ( 2 Z f Z D qf b xþl 0 0 ¼ dx f ðjÞ dj þ dx fðx ÞfðxÞ . 2 qx 4l xl L 0
(7)
ð8Þ
The new effective reaction–diffusion equation contains local and nonlocal couplings corresponding to the diffusive and the nonlocal contribution, respectively. The last one contains the nonlocal kernel given by Eq. (6), with a variable range 2l, that corresponds to the interaction of the field at points x0 belonging to the
ARTICLE IN PRESS B. von Haeften, H.S. Wio / Physica A 376 (2007) 199–207
202
region ½x l; x þ l. However, such points will contribute if and only if they are inside the domain ½L; L. As was discussed in previous works, the NEP functional allows us to analyze the changes in behavior associated with the variation of the system’s model parameters (like the system length or the kernel range), as those parameters determine the form of the fields that enter into the definition of the functional (Eq. (8)). We are now in a position to study the role played by the nonlocal kernel, particularly its range 2l, on the SR phenomenon. 2.2. Stationary solutions The above indicated system, given by Eq. (7), can be numerically solved, founding the stable stationary patterns. However, it was not possible, for arbitrary values of b, to find the unstable pattern (the one corresponding to the saddle separating the attraction basins of both stable stationary solutions). To overcome this difficulty, we have applied a (first order) perturbation approach, with the nonlocal coupling parameter b playing the role of the small parameter. Such a perturbation method consists in writing the field fðx; tÞ as an expansion in powers of b, as f ¼ fð0Þ þ bfð1Þ þ , and, in order to obtain the stationary solutions, replacing the indicated expansion in Eq. (7), with qf=qt ¼ 0, obtaining two expressions that allow us to evaluate the nonperturbed contribution fð0Þ and the correction fð1Þ . These two contributions should satisfy the Dirichlet boundary conditions. This perturbation expansion allows us to obtain the unstable stationary solution fu ðxÞ. The stable stationary solutions are the homogeneous f ¼ f0 ¼ 0, which exists for all values of the kernel range, and the inhomogeneous fs . The form of these patterns is illustrated in Fig. 1. 2.3. Stochastic resonance As indicated in the Introduction, the SR between the stationary solutions was investigated within the twostate approach. All details about the procedure and the evaluation of the SNR can be found in Refs. [10,19]. Our system is subject to a weak (in amplitude and frequency) external signal b ¼ b0 þ AðtÞ ¼ b0 þ Db cosðo0 tÞ, which corresponds to a modulation in the nonlinear term f ðfÞ, rocking the NEP [19]. In order to have a subthreshold signal, the amplitude Db satisfies Db5b. We have chosen the value of b0 corresponding to the value of b, at which we have F½fs ¼ F½f0 ¼ 0, when b ¼ 0. Up to the first order in the small amplitude Db the transition rates W i take the form W i ðtÞ ¼ mi ai Db cosðo0 tÞ,
(9)
Fig. 1. Stationary patterns: dashed lines correspond to the unperturbed system with b ¼ 0, and solid lines, to the perturbed one. The two attractors are shown, the nonhomogeneous fs and homogeneous f ¼ f0 ¼ 0, that coincide with the horizontal axis, as well as the unstable pattern fu . We have set D ¼ 0:6, 2L ¼ 6:25, b ¼ 0:02, b0 ¼ 0:719123 and 2l ¼ 2.
ARTICLE IN PRESS B. von Haeften, H.S. Wio / Physica A 376 (2007) 199–207
203
where the constants m1;2 and a1;2 are obtained from the Kramers-like formula for the transition rates [27] sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi lþ det F½fi Fmod ½fi Fmod ½fu (10) Wi ¼ exp 2p j det F½fu j between the attractors fi . Here Fmod is the NEP affected by the modulation, lþ is the unstable eigenvalue of the deterministic flux at the saddle point fu and m1 ¼ W 1 jAðtÞ¼0 ; dW 1 a1 ¼ dAðtÞ
m2 ¼ W 2 jAðtÞ¼0 , dW 2 ; a2 ¼ dAðtÞ
AðtÞ¼0
.
AðtÞ¼0
It is worth remarking here that in Ref. [28]—for the case of soliton’s nucleation in finite strings and assuming different boundary conditions, including Dirichlet ones—an expression similar to Eq. (10) was obtained. In this way we are able to calculate the SNR (denoted by R) through the equation R¼
pða2 m1 þ a1 m2 Þ2 p mm ¼ 2 1 2 C, 4 m1 þ m2 4m1 m2 ðm1 þ m2 Þ
where C is a factor depending on the nonhomogeneous attractor fs as Z L 2 fs 2 C¼ fs 1 dx . 3 L
(11)
(12)
After fixing the length of the system L and the kernel range 2l, we can use the above-indicated expressions to find the SNR that, as function of the noise intensity, shows the usual bell-shaped form of SR. In the next section we present some results. It is worth remarking that here we indicate the adopted value of 2l by a number of sites in our discrete lattice. Hence, the true value will be 2l Dx, with Dx the lattice size. 3. Results In Fig. 2 we show the dependence of the SNR on the kernel range 2l for a fixed value of 2L. There is a nonmonotonic behavior in the system’s response against variation of 2l. For fixed 2L, the transition rates are decreasing functions of the nonlocal kernel range 2l. Therefore, the ratio m1 m2 =ðm1 þ m2 Þ in the expression (11) also reflects this behavior. The other factor in this expression has a maximum for a kernel range that corresponds only to first neighbors sites. The maximum in the system’s
Fig. 2. SNR as a function of 2l for 2L fixed and other parameters as in Fig. 1.
ARTICLE IN PRESS 204
B. von Haeften, H.S. Wio / Physica A 376 (2007) 199–207
response as a function of the kernel range is due to the interplay between these two factors. The same fact can be physically interpreted analyzing the behavior of the NEP (computed from Eq. (8) for both field configurations: f0 and fs ) for varying 2l. In Fig. 3 we depict the dependence of the NEP on 2L. The fast initial increase of the barrier, followed by a slower one, is apparent explaining the result shown in Fig. 2. In Fig. 4 we show the dependence of the SNR on the length 2L for a fixed range 2l (100 sites). There is a monotonic behavior in the system’s response against variation of 2L. For a fixed range kernel 2l, both the transition rates and C are increasing functions of the length 2L. Therefore, the expression in Eq. (11) also reflects this behavior. As in the previous case, we can understand this fact through the analysis of the behavior of the NEP as function of 2L. In Fig. 5 we depict such a dependence for a fixed value of 2l. From the figure it is apparent that increasing 2L, the barrier height slowly increases. However, such a barrier increase is compensated by the rate’s prefactors, implying an increase of the SNR with 2L for fixed 2l. For not too large values of L, this behavior is approximately linear, as shown in Fig. 4.
Fig. 3. NEP (computed from Eq. (8) for f0 and fs ) as a function of 2l, with 2L fixed. The solid lines correspond to b ¼ 0:02 and the dashed ones to b ¼ 0:06. In both cases the top line corresponds to the nonhomogeneous unstable pattern, while the lower corresponds to the stable one. The NEP for f0 ¼ 0 coincides with the horizontal axis. Other parameters are as in Fig. 1.
Fig. 4. SNR as a function of 2L, with 2l fixed and other parameters as in Fig. 1.
ARTICLE IN PRESS B. von Haeften, H.S. Wio / Physica A 376 (2007) 199–207
205
Fig. 5. NEP (computed from Eq. (8) for f0 and fs ) as a function of 2L, with 2l fixed and the other parameters as in Fig. 1.
Fig. 6. NEP (computed from Eq. (8) for f0 and fs ) as a function of 2L, with 2l fixed. The solid lines corresponds to 2l ¼ 4 and the dashed ones to 2l ¼ 100. In both cases the top line corresponds to the nonhomogeneous unstable pattern, while the lower corresponds to the stable one. Other parameters are as in Fig. 1.
Finally, in Fig. 6 we depict the dependence of the NEP on b, for fixed 2L, and a couple of values of 2l. From the figure it is apparent that an increase in b implies a departure from the symmetric situation, together with a slow approach between the stable and unstable nonhomogeneous patterns that, for a finite value of b, finally coalesce. This result indicates that an arbitrary increment of b will not be necessarily followed by an increase in the response, but just the opposite will happen. 4. Conclusions Following the results in Ref. [16], we have analyzed an activator–inhibitor-like system, where the inhibitorlike field has been adiabatically eliminated yielding an effective, nonlocal, scalar reaction–diffusion equation. Arguing that the transport mechanism of the adiabatically eliminated inhibitor-like field has a nondiffusive character, we have chosen the form of the nonlocal kernel in such a way that it is more localized in space and has a controllable interaction range. By exploiting the knowledge of the NEP for the indicated case, we have analyzed the dependence of the SNR on the nonlocal interaction kernel range. We have found that the system’s response is reduced when the nonlocal interaction parameter b is increased. Also, the results clearly
ARTICLE IN PRESS B. von Haeften, H.S. Wio / Physica A 376 (2007) 199–207
206
show that an optimum value for the size of such kernel range exists, and that it corresponds to a very localized interaction. All the results could be physically understood through the analysis of the NEP’s behavior. In order to compare with the present theoretical results, detailed numerical simulations of the system defined in Eq. (7) are under way. This is a new and clear example of the usefulness of the NEP, when it is completely or even only partially known [9], in order to physically understand the origin of some noise-induced phenomena in extended systems. Acknowledgments The authors thank G. Izu´s for his help with the perturbation approach. H.S.W. thanks the European Commission for the award of a Marie Curie Chair at the Universidad de Cantabria, Spain.
References [1] [2] [3] [4]
[5]
[6]
[7]
[8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
[18]
[19] [20]
L. Gammaitoni, P. Ha¨nggi, P. Jung, F. Marchesoni, Rev. Mod. Phys. 70 (1998) 223. T. Wellens, V. Shatokin, A. Buchleitner, Rep. Prog. Phys. 67 (2004) 45. J.F. Lindner, B.K. Meadows, W.L. Ditto, M.E. Inchiosa, A.R. Bulsara, Phys. Rev. E 53 (1996) 2081. A.R. Bulsara, G. Schmera, Phys. Rev. E 47 (1993) 3734; P. Jung, U. Behn, E. Pantazelou, F. Moss, Phys. Rev. A 46 (1992) R1709; P. Jung, G. Mayer-Kress, Phys. Rev. Lett. 74 (1995) 2130; J.F. Lindner, B.K. Meadows, W.L. Ditto, M.E. Inchiosa, A.R. Bulsara, Phys. Rev. Lett. 75 (1995) 3; F. Marchesoni, L. Gammaitoni, A.R. Bulsara, Phys. Rev. Lett. 76 (1996) 2609. J.K. Douglas, et al., Nature 365 (1993) 337; J.J. Collins, et al., Nature 376 (1995) 236; S.M. Bezrukov, I. Vodyanoy, Nature 378 (1995) 362. A. Guderian, G. Dechert, K. Zeyer, F. Schneider, J. Phys. Chem. 100 (1996) 4437; A. Fo¨rster, M. Merget, F. Schneider, J. Phys. Chem. 100 (1996) 4442; W. Hohmann, J. Mu¨ller, F.W. Schneider, J. Phys. Chem. 100 (1996) 5388. H.S. Wio, Phys. Rev. E 54 (1996) R3075; H.S. Wio, F. Castelpoggi, Unsolved problems of noise, in: C.R. Doering, L.B. Kiss, M. Schlesinger (Eds.), Proceedings of Conference UPoN’96, World Scientific, Singapore, 1997, p. 229; F. Castelpoggi, H.S. Wio, Europhys. Lett. 38 (1997) 91. F. Castelpoggi, H.S. Wio, Phys. Rev. E 57 (1998) 5112. H.S. Wio, M.N. Kuperman, F. Castelpoggi, G. Izu´s, R. Deza, Physica A 257 (1998) 275; M. Kuperman, H.S. Wio, G. Izu´s, R. Deza, Phys. Rev. E 57 (1998) 5122. S. Bouzat, H.S. Wio, Phys. Rev. E 59 (1999) 5142. B. von Haeften, R. Deza, H.S. Wio, Phys. Rev. Lett. 84 (2000) 404. C.J. Tessone, H.S. Wio, P. Ha¨nggi, Phys. Rev. E 62 (2000) 4623. M.A. Fuentes, R. Toral, H.S. Wio, Physica A 295 (2001) 114. R. Graham, T. Tel, in: E. Tirapegui, W. Zeller (Eds.), Instabilities and Non-equilibrium Structures III, Kluwert, 1991; H.S. Wio, in: P. Garrido, J. Marro (Eds.), Fourth Granada Seminar in Computational Physics, Springer, Berlin, 1997, p. 135. G. Izu´s, et al., Phys. Rev. E 52 (1995) 129; G. Izu´s, et al., Int. J. Mod. Physics B 10 (1996) 1273. G. Drazer, H.S. Wio, Physica A 240 (1997) 571. D.H. Zanette, H.S. Wio, R. Deza, Phys. Rev. E 53 (1996) 353; F. Castelpoggi, H.S. Wio, D.H. Zanette, Int. J. Mod. Phys. B 11 (1997) 1717; G. Izu´s, R. Deza, H.S. Wio, Phys. Rev. E 58 (1998) 93; S. Bouzat, H.S. Wio, Phys. Lett. A 247 (1998) 297; G.G. Izu´s, R.R. Deza, H.S. Wio, Computer Phys. Comm. 121–122 (1999) 406. B. von Haeften, et al., Phys. Rev. E 69 (2004) 021107; B. von Haeften, et al., in: Z. Gingl, J.M. Sancho, L. Schimansky-Geier, J. Kerstez (Eds.), Noise in Complex Systems and Stochastic Dynamics, Proceedings of SPIE, vol. 5471, 2004, pp. 258–265. H.S. Wio, S. Bouzat, B. von Haeften, in: A. Robledo, M. Barbosa (Eds.), Proceedings of the 21st IUPAP International Conference on Statistical Physics, STATPHYS21, Physica A 306C (2002) 140–156. G. Schmid, I. Goychuk, P. Ha¨nggi, Europhys. Lett. 56 (2001) 22; G. Schmid, I. Goychuk, P. Ha¨nggi, Phys. Biol. 1 (2004) 61.
ARTICLE IN PRESS B. von Haeften, H.S. Wio / Physica A 376 (2007) 199–207
207
[21] P. Jung, J.W. Shuai, Europhys. Lett. 56 (2001) 29; J.W. Shuai, P. Jung, Phys. Rev. Lett. 88 (2003) 068102. [22] R. Toral, C. Mirasso, J. Gunton, Europhys. Lett. 61 (2003) 162. [23] A. Pikovsky, A. Zaikin, M.A. de la Casa, Phys. Rev. Lett. 88 (2002) 050601. [24] C.J. Tessone, R. Toral, Physica A 351 (2005) 106–116. [25] B. von Haeften, G.G. Izu´s, H.S. Wio, Phys. Rev. E 72 (2005) 021101. [26] H.S. Wio, in: T. Gonzalez, J. Mateos, D. Pardo (Eds.), Noise and fluctuations, Proceedings of the 18th International Conference on Noise and Fluctuations-ICNF2005, vol. 780, AIP, 2005, p. 55. [27] P. Ha¨nggi, P. Talkner, M. Borkovec, Rev. Mod. Phys. 62 (1990) 251. [28] F. Marchesoni, C. Cattuto, G. Constantini, Phys. Rev. B 57 (1998) 7930.