text:
# Comment on “Deuterium–tritium fusion reactors without external fusion breeding” by Eliezer et al.

(July 30, 1998 \[Rev. August 18, 1998\])

## Abstract

Inclusion of inverse Compton effects in the calculation of deuterium-deuterium burn under the extreme conditions considered by Eliezer et al. \[Phys. Lett. A 243 (1998) 298\] are shown to decrease the maximum burn temperature from about 300 keV to only 100–150 keV. This decrease is such that tritium breeding by the DD $``$ T + p reaction is not sufficient to replace the small amount of tritium that is initially added to the deuterium plasma in order to trigger ignition at less than 10 keV.

Published in Phys. Lett. A, Vol. 253 (1999) 119-121

Simplified models of thermonuclear burn are well known to lead to first order results that are generally in reasonable agreement with more sophisticated models. However, the prerequisite for such agreement is that, in the parameter range under consideration, all physical phenomena that are relevant to this order of approximation be consistently taken into account.

In a recent Letter, the results of space-time independent calculations of very-high temperature ($`T_i300`$ keV, $`T_e150`$ keV) thermonuclear burn of a very-high density (compression of $``$ 30’000 times solid density) DD plasma containing a small initial addition of tritium are presented . From these calculations, the authors conclude that the initial amount of about 2% of tritium needed to trigger ignition at less that 10 keV is regenerated by the DD reaction if the burning temperature is higher that 200 kev. If this is correct, it would be possible to conceive DT fusion reactors which would not need external tritium breeding for their operation.

However, in these calculations, a number of effects and processes which become important at extremely high temperatures and densities seem not to have been fully taken into account. These phenomena comprise various relativistic and degeneracy effects, as well as several radiation processes which include inverse bremsstrahlung, electron-electron bremsstrahlung, Compton scattering, etc. Of these phenomena, the most important is electron energy loss by inverse Compton scattering. The purpose of this Comment is to show that its inclusion in the simulation totally changes the main conclusion of Ref. .

The energy lost by electrons of temperature $`kT_e`$ in a photon gas by the inverse Compton effect is

$$W_C\frac{dE_C}{dt}=\frac{32}{3}\pi r_{e}^{}{}_{}{}^{2}c\frac{E_r}{mc^2}N_e(kT_ekT_r),$$
(1)
where $`E_r=\frac{4}{c}\sigma T_{r}^{}{}_{}{}^{4}`$ is the radiation energy density and $`N_e`$ the electron number density . An estimate of $`E_r`$ can be made by assuming that the average traveling time of the bremsstrahlung photons through a thin plasma of radius $`R`$ is $`\mathrm{\Delta }tR/c`$. Therefore, if $`W_B`$ is the bremsstrahlung radiation loss rate, $`E_rW_BR/c`$. Using this estimate, it is readily seen that $`W_C>W_B`$ for $`T_e>60`$ keV in the parameter range considered in Ref. , so that Compton effects cannot be neglected.

To make this statement more quantitative, and in order to calculate how much of the initial tritium can effectively be regenerated during thermonuclear DD burn, we have used our simulation program ISRINEX  to reproduce the results of Ref. 1 and to show the consequences of taking radiation effects more fully into account.

In our program, inverse bremsstrahlung and Compton effects are simulated by a model attributed to H. Hurwitz. This model was first published in  and is used in both small (e.g., ) and large (e.g., ) simulation programs. However, since there is some discrepancy in the published literature on the details of this model, we take this Comment as an opportunity to review it briefly.

Hurwitz’s model is a three-temperature model in which the electron and photon energy equations are coupled as follows:

$`{\displaystyle \frac{dE_e}{dt}}`$ $`=`$ $`{\displaystyle \frac{dE_B}{dt}}{\displaystyle \frac{dE_C}{dt}},`$ (2)
$`{\displaystyle \frac{dE_r}{dt}}`$ $`=`$ $`+{\displaystyle \frac{dE_B}{dt}}+{\displaystyle \frac{dE_C}{dt}}R_{loss}.`$ (3)
The first term corresponds to bremsstrahlung effects. It is obtained by integrating over frequency the difference between emission by the electrons of bremsstrahlung radiation and reabsorption of the local radiation field (inverse bremsstrahlung). With the notation of Ref. , this is

$$\frac{dE_B}{dt}=_0^{\mathrm{}}c\kappa _\nu (U_{\nu p}U_\nu )𝑑\nu =W_BG(\gamma )\frac{T_eT_r}{T_e},$$
(4)
where $`\gamma =T_r/T_e`$ and

$$G(\gamma )=\frac{1}{1\gamma }_0^{\mathrm{}}\frac{1}{2}\mathrm{exp}(x/2)\mathrm{K}_0(x/2)\frac{1\mathrm{exp}(x/\gamma x)}{1\mathrm{exp}(x/\gamma )}𝑑x,$$
(5)
where $`\mathrm{K}_0`$ is the Bessel function characteristic of the quantum-mechanical bremsstrahlung emission spectrum. For $`\gamma 0`$, $`G(\gamma )1`$, making $`dE_B/dt`$ equal to the pure bremsstrahlung rate $`W_B`$; for $`\gamma \mathrm{}`$, $`G(\gamma )\pi ^2/4`$, giving pure inverse bremsstrahlung. Except for integer values of $`\gamma `$, $`G(\gamma )`$ cannot be calculated analytically . A simplified fit, $`G(\gamma )\pi ^2/4(\gamma +1.03)/(\gamma +2.54)`$ is therefore used in the calculations.

The second term in (2,3) corresponds to Compton effects as given by (1). The third term in (3) is the photon energy loss rate which for the purpose of this Comment is approximated by:

$$R_{loss}=\frac{3}{R}\sigma T_{r}^{}{}_{}{}^{4},$$
(6)
or

$$R_{loss}=\mathrm{min}[\frac{3}{R}\sigma T_{r}^{}{}_{}{}^{4},W_B].$$
(7)
Expression (6) is the maximum possible photon energy loss, so that its use tends to underestimate the value of $`T_r`$ and thus the effect of Compton scattering. On the other hand, expression (7), which gives the correct limit in the case of an infinitely thin plasma (and which gives results in close agreement with simulations presented in ), tends to minimize radiation losses and therefore to overestimate Compton effects.

The results of our calculations for the simulation shown in Figs. 2–4 of Ref.  are given in Table 1. The first column gives our results when the only radiation effect included is bremsstrahlung losses, i.e., $`dE_e/dt=W_B`$, as in Ref. . They are in reasonable agreement with Ref. , except possibly for the final tritium content which is only 0.83 instead of slightly more than 1. This may partly be explained by the fact that in our calculation the maximum ion temperature is only 270 keV instead of about 300 keV. Otherwise, the ratio of tritium content goes through a minimum at a burn time of about 2 ps, and through a maximum at about 6 ps, as in Fig. 2 of Ref. .

The second column corresponds to Hurwitz’s three-temperature model with inverse bremsstrahlung and radiation losses included according to (4) and (6). The effects on $`T_i`$ and $`T_e`$ are small, and the final fuel burn-up fraction is unchanged. However, the maximum radiation temperature, i.e., 17 keV, is large enough to imply that Compton effects cannot be neglected.

The third column gives our results when both bremsstrahlung and Compton effects are fully included. The maximum ion temperature is only half the $``$ 300 keV found in . Consequently, $`T_i`$ is not large enough to breed enough tritium to replace the initial tritium content. The final tritium content is only about 1/4 of its initial value.

Finally, the last column shows what happens when the photon losses are reduced, e.g., when the thermonuclear fuel is surrounded by a some material opaque to low energy photons but relatively transparent to bremsstrahlung. $`T_i`$ is further reduced to $``$ 100 keV, the final tritium content ratio is less than 1/5, and the final burn-up ratio is only 0.34.

In conclusion, we have shown that by taking radiation effects such as inverse Compton scattering fully into account the maximum burn temperature of a highly compressed DD pellet is reduced from about 300 keV to only 100 keV. This decrease in ion temperature leads to a substantial reduction in the fuel burn-up ratio. Moreover, this decrease is such that tritium breeding by the DD reaction is not sufficient to replace the amount of tritium that is added to the deuterium plasma in order to trigger ignition at less than 10 keV. These conclusions would be reinforced by taking relativistic and degeneracy effects into account. Finally, if the size of the pellet were larger, or if the fuel were surrounded by various materials, as is the case in a realistic pellet design, the importance of radiation effects would be further enhanced. These problems illustrate the difficulty of burning advanced thermonuclear fuels, such as DTx, D3He, or p11B , relative to burning those containing a large tritium fraction, such as DT or 6LiDT.


format: text

text:
# Modeling of Spiking-Bursting Neural Behavior Using Two-Dimensional Map

## I Introduction

Understanding dynamical principles and mechanisms behind the control of activity, signal and information processing that occur in neurobiological networks is hardly possible without numerical studies of collective dynamics of large networks of neurons. These simulations need to take into account the complex architecture of couplings among individual neurons that is suggested by data from biological experiments. One of the complicating factors in understanding the simulation results is the complexity of temporal behavior of individual biological neurons. This complexity is due to the large number of ionic currents involved in the nonlinear dynamics of neurons. As a result, realistic channel-based models proposed for a single neuron is usually a system of many nonlinear differential equations (see, for example  and the review of the models in ). The strong nonlinearity and high dimensionality of the phase space is a significant obstacle in understanding the collective behavior of such dynamical systems. The dynamical mechanisms responsible for restructuring collective behavior of a network of channel-based neuron models are difficult or impossible to analyze because these mechanisms are well hidden behind the complexity of the equations. If there is a chance to identify possible dynamical mechanisms behind a particular behavior of the network without the use of complex models, then this knowledge can be used to guide through numerical study of this behavior with the high-dimensional channel-based models. Based on these results one can specify the actual dynamical mechanism occurred in the network and understand the biological processes which contribute to it.

One way to reveal the dynamical mechanisms is to use a simplified or phenomenological model of the neuron. However, many experiments indicate the restructuring of collective behavior utilizes a variety of dynamical regimes generated by individual dynamics of a neuron. When neural dynamics include the regime chaotic spiking-bursting oscillations the system of differential equations describing each neuron should be at least a three-dimensional system (see for example ). Further simplification of phenomenological models for complex dynamics of the neuron can be obtain using dynamical systems in the form of a map. Despite the low-dimensional phase space of this nonlinear map, it is able to demonstrate a large variety of complex dynamical regimes.

The use of low-dimensional model maps can be useful for understanding the dynamical mechanisms if they mimic the dynamics of oscillations observed in real neurons, show correct restructuring of collective behavior, and are simple enough to study the reasons behind such restructuring. The development of a model map which is capable of describing various types of neural activity, including generation of tonic spikes, irregular spiking, and both regular and irregular bursts of spikes, is the goal of this paper.

It is known that constructing a low-dimensional system of differential equation which is capable of generating fast spikes bursts excited on top of the slow oscillations, one needs to consider a system which has both slow and fast dynamics (see for example ). Using the same approach one can construct a two-dimensional map, which can be written in the form

$`x_{n+1}`$ $`=`$ $`f(x_n,y_n+\beta _n),`$ (1a)
$`y_{n+1}`$ $`=`$ $`y_n\mu (x_n+1)+\mu \sigma _n,`$ (1b)
where $`x_n`$ is the fast and $`y_n`$ is the slow dynamical variable. Slow time evolution of $`y_n`$ is due to small values of the parameter $`\mu =0.001`$. Terms $`\beta _n`$ and $`\sigma _n`$ describe external influences applied to the map. These terms model the dynamics of the neuron under the action of the external DC bias current $`I_{DC}`$ and synaptic inputs $`I_n^{syn}`$. The term $`\sigma _n`$ can also be used as the control parameter to select the regime of individual behavior.

A model in the form of two-dimensional map, whose form is similar to (1) was used in the study of dynamical mechanisms behind the emergence and regularization of chaotic bursts in a group of synchronously bursting cells coupled through a mean field. The effect of anti-phase regularization was also modeled before with one-dimensional maps . In both these models the oscillations during the burst were described by a chaotic trajectory. Map (1) improves these models by adding a feature that enables one to mimic the dynamics of individual spikes within the burst. This is achieved using a modification of the shape of the nonlinear function $`f(x,y)`$ which is now a discontinuous function of the form

$`f(x,y)=\{\begin{array}{cc}\alpha /(1x)+y,\hfill & x0\hfill \\ \alpha +y,\hfill & 0<x<\alpha +y\hfill \\ 1,\hfill & x\alpha +y\hfill \end{array}`$ (2)
where $`\alpha `$ is a control parameter of the map. The dependence of $`f(x,y)`$ on $`x`$ computed for a fixed value of $`y`$ is shown in Fig.1. In this plot the values of $`\alpha `$ and $`y`$ are set to illustrate the possibility of coexistence of limit cycle, $`P_k`$, corresponding to spiking oscillation in (1a), and fixed points $`x_{ps}`$ and $`x_{pu}`$. Note that when $`y`$ increases or decreases the graph of $`f(x,y)`$ moves up or down, respectively, except for the third interval $`x\alpha +y`$, where the values of $`f(x,y)`$ always remain equal to -1.

The paper is organized as follows: Section II considers the features of the fast and slow dynamics which explain the formation of various types of behavior in the isolated map. It also discusses the bifurcations responsible for the qualitative change of activity. The dynamics of response generated by the map after it is excited or inhibited by an external pulse is discussed in Section III. The goal of this section is to illustrate dynamical mechanisms that control the response, and to show how the relation between the two inputs of the map influence the properties of the response. Section IV presents the results of a synchronization study in two coupled maps.

## II Individual Dynamics of the Map

Consider the regimes of oscillations produced by the individual dynamics of the map. In this case the inputs $`\beta _n=\beta `$ and $`\sigma _n=\sigma `$ are constants. Note that if $`\beta `$ is a constant, then it can be omitted from the equations using the change of variable $`y_n+\beta y_n^{new}`$. Therefore, the individual dynamics of the map depends only on the control parameters $`\alpha `$ and $`\sigma `$, and the map can be rewritten in the form

$`x_{n+1}`$ $`=`$ $`f(x_n,y_n),`$ (3a)
$`y_{n+1}`$ $`=`$ $`y_n\mu (x_n+1)+\mu \sigma .`$ (3b)
Typical regimes of temporal behavior of the map are shown in Figs 2 and 3. When the value of $`\alpha `$ is less then 4.0 then, depending on the value of parameter $`\sigma `$, the map generates spikes or stays in a steady state (see Fig 2). The frequency of the spikes increases as the value of parameter $`\sigma `$ is increased (see Fig 2).

For $`\alpha >4`$ the map dynamics are capable of producing bursts of spikes. The spiking-bursting regimes are found in the intermediate region of parameter $`\sigma `$ between the regimes of continuous tonic spiking and steady state (silence). The spiking-bursting regimes include both periodic and chaotic bursting. A few typical bursting-spiking regimes computed for different values of $`\sigma `$ are presented in Fig 3.

Due to the two different time scales involved in the dynamics of the model, the mechanisms behind the restructuring of the dynamical behavior can be understood through the analysis of the fast and slow dynamics separately. In this case, time evolution of the fast variable, $`x`$, is studied with the one-dimensional map (3a) where the slow variable, $`y`$, is treated as a control parameter whose value drifts slowly in accordance with equation (3b). It follows from (3b) that the value of $`y`$ remains unchanged only if $`x=x_s`$ given by

$`x_s=1+\sigma .`$ (4)
If $`x<x_s`$, then the value of $`y`$ slowly increases. If $`x>x_s`$, then $`y`$ decreases.

From (3a) one can find the equation for the coordinate of fixed points $`x_p`$ of the fast map. This equation is of form

$`y=x_p{\displaystyle \frac{\alpha }{1x_p}},`$ (5)
where $`x_p0`$. Equation (5) defines the branches of slow motion in the two-dimensional phase space $`(x_n,y_n)`$ (see Fig 4). The stable branch $`S_{ps}(y)`$ exists for $`x_p<1\sqrt{\alpha }`$ and the unstable branch $`S_{pu}(y)`$ exists within $`1\sqrt{\alpha }x_p0`$

Considering the fast and slow dynamics together, one can see that, if $`x_s`$ is in the stable branch $`S_{ps}(y)`$, then the map (3) has a stable fixed point. The stable fixed point corresponds to the regime of silence in the neural dynamics. The oscillations in the map dynamics will appear when $`x_s>1\sqrt{\alpha }`$. This is the threshold of excitation, which corresponds to the bifurcation values of $`\sigma `$ given by

$`\sigma _{th}=2\sqrt{\alpha }.`$ (6)
To understand the dynamics of excited neurons within this model one needs to consider branches that correspond to the spiking regime. To evaluate the location of the spiking branches, consider the mean value of $`x_n`$ computed for the periodic trajectory of the fast map (3a) as a function of $`y`$. Here $`y`$ is treated as a parameter. Use of such an approximation is quite typical for an analysis involving fast and slow dynamics. It works well for small values of parameter $`\mu `$.

It follows from the shape of $`f(x,y)`$ that for any value of $`y`$, map (3a) generates no more than one periodic trajectory $`P_k`$, where $`k`$ is the period (i.e. $`x_n=x_{n+k}`$). The cycles $`P_k`$ always contain the point $`x_n=1`$ and the index $`k`$ increases stepwise ($`kk+1`$) as $`y`$ decreases. Since the trajectory $`P_k`$ always has a point in the flat interval of $`f(x,y)`$, all these cycles are super-stable except for bifurcation values of $`y`$ for which the trajectory contains the point $`x_n=0`$ (see Fig 1). The location of the spiking branch, $`S_{spikes}`$, of ”slow” motion in the phase plane ($`x_n,y_n`$) can be estimated as the mean value of $`x`$ computed for the period of cycle $`P_k`$,

$`x_{mean}={\displaystyle \frac{1}{k}}{\displaystyle \underset{m=1}{\overset{k}{}}}f^{(m)}(1,y),`$ (7)
where $`k`$ is the period of $`P_k`$, and $`f^{(m)}(x,y)`$ is the $`m`$-th iterate of (3a), started at point $`x`$ and computed for fixed $`y`$. The spiking branch of ”slow” dynamics evaluated with (7) is shown in Fig4. One can see that this branch has many discontinuities caused by the bifurcations of the super-stable cycles $`P_k`$.

To complete the picture of fast and slow dynamics of the model for $`\alpha >4`$, one needs to consider the fast map bifurcation associated with the formation of homoclinic orbit $`h_{pu}`$ originating from the unstable fixed point, $`x_{pu}`$. This homoclinic orbit occurs when the coordinate of $`x_{pu}`$ become equal to -1. It can be easily shown that such a situation can take place only if $`\alpha >4`$. The homoclinic orbit forms at the value of $`y`$ where the unstable branch $`S_{pu}(y)`$ crosses the line $`x=1`$, see Fig.4b. When the map is firing spikes and the value of $`y`$ gets to the bifurcation point, the cycle $`P_k`$ merges into the homoclinic orbit, disappears, and then the trajectory of the map jumps to the stable fixed point $`x_{ps}`$.

Typical phase portraits of the model, obtained under the assumptions made above, are presented in Fig.4. Fig.4a shows the typical behavior for $`2<\alpha 4`$. Here, only two regimes are generated. The first regime is the state of silence, when the operating point (OP), given by the intersection of $`x_s=1+\sigma `$ and one of the branches, is on the branch $`S_{ps}(y)`$. The second regime is the regime of tonic spiking, when the operating point is selected on the spiking branch $`S_{spikes}`$.

When $`\alpha >4`$ the picture changes qualitatively (see Fig 4b). Now, stable branches $`S_{ps}(y)`$ and $`S_{spikes}`$ are separated by the unstable branch $`S_{pu}(y)`$. If the operating point is selected on $`S_{pu}(y)`$, then the phase of silence, corresponding to slow motion along $`S_{ps}(y)`$, and spiking, when system moves along $`S_{spikes}`$, alternate forming the regime of spiking-bursting oscillations. The beginning of a burst of spikes corresponds to the bifurcation state of the fast map where fixed points $`x_{ps}`$ and $`x_{pu}`$ merge together and disappear. Before this bifurcation the system is in $`x_{ps}`$ and, therefore, $`y`$ increases. The termination of the burst is due to the bifurcation of the fast map associated with the formation of the homoclinic orbit $`h_{pu}`$. Here, the limit cycle of the spiking mode merges into the homoclinic orbit and disappears. After that the fast subsystem flips to the stable fixed point $`x_{ps}`$. Then the process repeats (see arrows in Fig 4b).

It is clear from Fig.4a that when operating point is set on $`S_{ps}(y)`$ the model will be in the regime of silence. When the operating point is set on the branch $`S_{spikes}`$ the model produces tonic spiking, unless the point is set close to the formation of homoclinic orbit $`h_{pu}`$. One can see that, at the vicinity of this bifurcation, the branch $`S_{spikes}`$ becomes densely folded. As a result, the behavior of $`y`$, which is governed by the mean value of $`x_n`$, can become extremely sensitive to small perturbations and even lead to instability caused by high-gain feedback. This is one of the reasons for the irregular, chaotic spiking-bursting behavior which occurs in the map when the operating point is set close to the area of the branch $`S_{spikes}`$ where this branch is densely folded. The detailed and rigorous analysis of chaotic dynamics cannot be done within the approximations made above and require more precise computation of $`S_{spikes}`$ which is beyond the scope of this paper.

The results of the analysis presented above are summarized in the sketch of the bifurcation diagram plotted on the parameter plain ($`\sigma ,\alpha `$) (see Fig 5). The bifurcation curve $`\sigma _{th}`$ corresponds to the excitation threshold (6) where the fixed point of the 2-d map becomes unstable and the map starts generating spikes. Curve $`L_{ts}`$ shows the approximate location of the border between spiking and spiking-bursting regimes, obtained in the numerical simulations. Note that separation of spiking and bursting regimes is not always obvious, especially in the regime of chaotic spiking. The regime of spiking-bursting oscillations takes place within the upper triangle formed by curves $`\sigma _{th}`$ and $`L_{ts}`$. The regimes of chaotic spiking or spiking-bursting behavior are found in the relatively narrow region of the parameters located around $`L_{ts}`$. This region contains complex structure of bifurcations, associated with multi-stable regimes, and is not shown in Fig 5.

The bifurcation diagram shows the role of control parameters in the selection of dynamical features of the considered neuron model. Both parameters can be used to mimic a particular type of neural behavior. Parameter $`\sigma `$ can be used to model the external DC current injection that depolarizes or hyperpolarizes the neuron. In this case $`\sigma `$ can be written as

$`\sigma =\sigma _u+I_{DC},`$ (8)
where $`\sigma _u`$ is the parameter that selects the dynamics of isolated neuron, and $`I_{DC}`$ is the parameter that models the DC current injected into the cell. The changes in behavior caused by $`I_{DC}`$ are similar to the action of the parameter $`I`$ in the well known Hindmarsh-Rose model .

If the individual dynamics of the modelled neuron are capable of generating only regimes of silence and tonic spiking, and do not support the regime of bursts of spikes, then the value of $`\alpha `$ should be set below 4. In this case, no matter what the level of $`I_{DC}`$ injected is, the system will not show bursts of spikes (see Fig 5).

## III Modeling of Response to the Injection of Current

To study the dynamical regimes of neuronal behavior in experiments, biologists change the type of neural activity using the injection of electrical current into the cell though an electrode. It was shown above that the injection of DC current can be modeled in the map (1) using parameter $`\sigma _n=\sigma `$ (see equation (8)). In this case, since the external influence does not vary in time, the role of parameter $`\beta _n=\beta `$ is not important, because the behavior of the map after the transient is independent of the value of $`\beta `$. However, when the injected current changes in time, it may be useful to consider the dynamics of parameter $`\beta _n`$ in order to provide more realistic modeled behavior during the transient. Taking this into account, the input of the model can be considered in the form:

$`\beta _n=\beta ^eI_n,\sigma _n=\sigma ^eI_n,`$ (9)
where $`I_n`$ is injected current, and coefficients $`\beta ^e`$ and $`\sigma ^e`$ are selected to achieve the desired properties of response behavior.

This Section briefly illustrates how the relation between these coefficients effects the dynamics of response to the pulse of $`I_n`$. To be specific, consider the map in the regime of tonic spiking with $`\alpha =5.0`$, $`\sigma =0.33`$, $`\beta =0`$. To study the response behavior, positive and negative pulses of amplitude 0.8 and duration of 100 iterations were applied to the continuously spiking map. In the simulations presented below the coefficient $`\sigma ^e`$ was selected to be equal to one.

Figures 6 and 7 show the response to the positive and negative pulse, respectively, computed with $`\beta ^e=0`$. In this case, during the action of a positive pulse, the value of $`y_n`$ increases monotonically because of the increased value of $`\sigma _n`$ (see (1b)). The increase of $`y_n`$ pushes the fast map up, see Fig 1. This leads to an increase in the frequency of spiking. After the action of the pulse ends the value of $`y_n`$ monotonically decreases back to the original state (see Fig. 6).

When a negative pulse is applied, it pushes the operating point down and, if the amplitude of the pulse is sufficiently large, shuts off the regime of spiking (see Fig.7). This happens because the trajectory of the fast map gets to the perturbed operating point which is now on the stable branch $`S_{ps}(y)`$, (see Fig 4b). After the pulse is over the spiking does not re-appear immediately because the system spends some time drifting along the stable branch of slow motions $`S_{ps}(y)`$, during which the variable $`y`$ overshoots its original value for the spiking regime. As a result, after the system switches to the spiking branch $`S_{spikes}`$, $`y_n`$ monotonically drifts down to the unperturbed operating point. The dynamics of the slow evolution are clearly seen in the lower panel of Fig.7.

Figures 8 and 9 illustrate how response to the pulse of $`I_n`$ changes when coefficient $`\beta _n`$ is not equal to zero. To be brief and specific consider the case of $`\beta ^e=1.0`$ and $`\sigma ^e=1.0`$. Comparing Fig.8 with Fig.6 one can see that the response dynamics to the positive pulse changes qualitatively. Indeed, when the pulse current is applied, then, acting through $`\beta _n`$, it immediately forces the fast map to shift up. As a result, the rate of spiking and mean value of $`x_n`$ increases sharply. Variable $`y_n`$ reacts to this change and decreases its value to compensate for the sudden change of the mean value of $`x_n`$. When the action of the pulse is over, the value of $`\beta _n`$ returns to its original value and, due to the updated levels of $`y_n`$, the fast map overshoots its original state. As a result, the trajectory of the fast map reaches the stable fixed point. To return to the original regime of spiking the system has to go all the way along the branches of slow motion $`S_{ps}(y)`$ and $`S_{spikes}`$ back to the original operating point (see Fig. 4b). This type of response is observed in real neurobiological experiments, see for example .

Using the same analysis one can understand the new effects in the response to a negative pulse caused by the action of $`\beta _n`$. These new effects can be clearly seen comparing Fig.9 with Fig. 7.

The presented results illustrate how 2-d map (1) along with equation (9) can model a large variety of transient neural behavior induced by injected current. Due to the simplicity of the model one can clearly see the nonlinear mechanisms behind the response behavior and apply them to select the desired balance between $`\sigma ^e`$ and $`\beta ^e`$ to model a particular type of response.

## IV Regimes of synchronization in two coupled maps

This section presents the results of studies of synchronization regimes in coupled chaotically bursting 2-d maps. The goal of this study is to reproduce the main regimes of synchronous behavior found in a real neurobiological experiment . The experiment was carried out on two electrically coupled neurons (the pyloric dilators, PD) from the pyloric CPG of the lobster stomatogastric ganglion . The regimes found in the experiment were also reproduced in numerical simulations using Hindmarsh-Rose model , one-dimensional map model , and in experiments with electronic neurons .

The equations used in the numerical simulations of the coupled maps are of form

$`x_{i,n+1}`$ $`=`$ $`f(x_{i,n},y_{i,n}+\beta _{i,n}),`$ (10)
$`y_{i,n+1}`$ $`=`$ $`y_{i,n}\mu (x_{i,n}+1)+\mu \sigma _i+\mu \sigma _{i,n},`$ (11)
where index $`i`$ specifies the cell, and $`\sigma _i`$ is the parameter that defines the dynamics of the uncoupled cell. The coupling between the cells is provided by the current flowing from one cell to the other. This coupling is modeled by

$`\beta _{i,n}`$ $`=`$ $`g_{ji}\beta ^e(x_{j,n}x_{i,n})`$ (12)
$`\sigma _{i,n}`$ $`=`$ $`g_{ji}\sigma ^e(x_{j,n}x_{i,n})`$ (13)
where $`ij`$, and $`g_{ji}`$ is the parameter characterizing the strength of the coupling. The coefficients $`\beta ^e`$, $`\sigma ^e`$ set the balance between the couplings for the fast and slow processes in the cells, respectively. In the numerical simulations the values of the coefficients are set to be equal: $`\beta ^e=1.0`$ and $`\sigma ^e=1.0`$. The other parameters of the coupled maps (11) that remain unchanged in the simulations have the following values: $`\mu =0.001`$, $`\alpha _1=4.9`$, $`\alpha _2=5.0`$. The coupling between the maps is symmetrical, $`g_{ji}=g_{ij}=g`$.

First, consider the main regimes of synchronization between the maps generating irregular, chaotic bursts of spikes. To set the individual dynamics of the maps to this regime, the parameters $`\sigma _i`$, that take into account the DC bias current injected into the neurons (see  for details), are tuned to the following values: $`\sigma _1=0.240`$, $`\sigma _2=0.245`$. When these maps are uncoupled ($`g=0`$) they produce chaotic spiking-bursting oscillations shown in Fig 10a. When the coupling becomes sufficiently large the slow components of the bursts synchronize while spikes within the bursts remain asynchronous see the waveforms in Fig. 10b obtained with $`g=0.043`$. This regime of synchronization is typical for naturally coupled PD neurons (see Figure 2a in ).

Introduction of negative coupling, $`g<0`$, also leads to synchronization of bursts, but in this regime of synchronization the systems burst in antiphase. Typical waveforms produced in this regime are presented in Fig. 10c. This regime of synchronization of chaotic bursts is also observed in the experimental study of coupled PD neurons (see Figure 2c in ). It is important to emphasize that, both in this simulation and in the experiment, the regime of antiphase synchronization characterized by the onset of regular bursts.

The simplicity of this model enables one to understand the possible cause for the onset of regular bursting in the antiphase synchronization. It is shown in Section II that chaotic dynamics of bursts occurs when the operating point of the system appears close to the leftmost area of the spiking branch $`S_{spikes}`$. In this region, $`S_{spikes}`$ is densely folded and, if system slows down in this area, the timing for the end of the burst becomes very sensitive to infinitesimal perturbations.

This fact is illustrated in Fig. 11a. This figure shows the attractor computed from the waveforms $`x_{1,n}`$ and $`y_{1,n}`$ of the first cell when the cells are uncoupled (this regime is shown in Fig. 10a). To plot the attractor, the waveforms of $`x_{1,n}`$ and $`y_{1,n}`$ were filtered by a fourth order low-pass filter with cutoff frequency 0.6. As a result, the attractor does not contain sharp spikes and the approximate location $`S_{spikes}`$, which is close to the center of filtered spiking oscillations, is easy to see. The level of coordinate $`x_f`$, corresponding to the operating point, is equal to $`x_s=1+\sigma _1=0.76`$. One can see from Fig. 11a that when the center of oscillations gets close to that level the trajectory becomes rather complex, as shown by the dense set of trajectories in the left part of the attractor.

The numerical simulations show that this effect of regularization takes place even when $`\sigma ^e=0`$. Therefore, the vertical shift of operating point, $`x_s`$, caused by the slow coupling current is not very important for this effect. The coupling term $`\beta _{i,n}`$ directly influences the fast dynamics by shifting the graph of the 1-d map (see Fig. 1) up or down depending on the sign of $`\beta _{i,n}`$. In the regime of antiphase synchronization $`\beta _{i,n}`$ is positive, when the $`i`$-th cell is spiking, and negative, when it is silent. Therefore, from the viewpoint of fast dynamics, the coupling forces the cells to stay on the current branch of slow motion. One can understand these dynamics from the graph of $`f(x,y)`$ and (1a). This effect is also seen as the formation of extended shape of attractor plotted for the regime of antiphase synchronization (see Fig. 11b). Note that axes in Fig. 11a and Fig. 11b have different scales.

In the regime of chaotic bursts the level of $`x_s`$ is close to the stable branch of spiking $`S_{spikes}`$, and, therefore, the slow evolution along the branch of silence $`S_{pu}`$ is faster then along $`S_{spikes}`$. This means that duration of the phase of silence in the cell is shorter than duration of the burst of spikes. When the cell switches from silence to spiking it sharply changes the value and the sign of $`\beta _{i,n}`$. This change pushes $`S_{spikes}`$ of the spiking cell down to its original location, and, as a result, quickly drags the trajectory through the area of complex behavior toward the branch $`S_{pu}`$. Due to this fast change the duration of bursting is determined only by the dynamics of the silent phase, which is regular.

It was observed in the experiment that, in the regime of tonic spiking, spikes in PD neurons synchronize at low levels of coupling (i.e. without added artificial coupling, see  for details). This property of synchronization in the regime of tonic spiking is also typical for the dynamics of the maps considered here (see Fig.12). In these numerical simulations, the uncoupled maps were tuned to generate spikes with slightly different rates. One can clearly see the beats between the waveforms plotted in Fig.12a. When small coupling, $`g=0.008`$, is introduced the beats disappear as soon as the spikes get synchronized (see Fig. 12b).

Although the synchronization between continuously spiking maps is easily achieved, it is important to understand that the dynamics of such synchronization in the maps is a much more complex process than it is in the models with continuous time. Due to the discrete time, the periodic spikes can lock with different frequency ratios. This results in the existence of a complex multistable structure of synchronization regimes. This structure becomes more noticeable in the dynamics of synchronization as the number of iterations in the period of spikes decreases.

## V Discussion and outlook

The simple phenomenological model of complex dynamics of spiking-bursting neural activity is proposed. The model is given by two-dimensional map (1). The first equation (1a) of the map describes the fast dynamics. Its isolated dynamics is capable of generating stable limit cycles, which mimics the spiking activity of a neuron, and a stable fixed point, which corresponds to phase of silence. This 1-d map has a region of parameters were both these stable regimes coexist. Existence of such multistable regimes is the reason for generation of bursts, when the operating point, defined by the dynamics of the second equation (1b), is properly set (see Fig 4b).

The shape of the nonlinear function $`f(x,y)`$ used in (1a) is selected in the form (2) based on the following considerations: (i) The fast map should generate limit cycles whose waveforms mimics those of the spikes. (ii) Each spike generated by the map always has a single iteration residing on the most right interval $`x\alpha +y`$. Therefore, the moment of time corresponding to the appearance of a trajectory on this interval can be used to define the time of a spike. This feature is important for modeling the dynamics of chemical synapses in a group of coupled neurons. (iii) The analytic expression of the nonlinear function in interval $`x<0`$ should be simple enough to allow rigorous analysis of bifurcations of fixed points of the fast map. (iv) Use of the fixed level, $`1`$, in the rightmost interval of the function simplifies the analysis of dynamics at the end of the bursts. The end of a burst is associated with the formation of a homoclinic orbit, which corresponds to the case when the trajectory from this interval maps into the unstable fixed point (see Section II).

It is clear that the shape of function (2) can be modified to take into account other dynamical features that need to be modeled. Due to the low-dimensionality of the model, the dynamical mechanisms behind its behavior are easy to understand using phase plane analysis. This allows one to see how introduced modifications influence the dynamics.

Some modification of the fast map is required to enhance the region of parameters where the model can still be used to mimic neural dynamics properties. For example, from (2) and Fig. 1 one can see that, due to the dynamics of coupling terms, the trajectory of the fast system can stay in the middle interval of $`f(x,y)`$ for several iterations. This can happen when the external influence monotonically pushes the function up while the trajectory is located in the middle interval. In this case the trajectory will map to the middle interval again and again, increasing the duration of a spike. This artifact can be removed using one of the following modifications. Introduction of a sufficient gap between the right end of this interval and the diagonal (see Fig. 1) will help the map to terminate the spike despite the monotonic elevation of the function. Alternatively one can introduce an additional condition to the fast map (1a) that forces the map always to iterate its trajectory from the middle interval to the rightmost one, despite the dynamics of $`y`$ (see (2)). This can be achieved using the function $`f(x_n,y)`$ of the following form

$`f(x_n,y)=\{\begin{array}{cc}\alpha /(1x_n)+y,\hfill & x_n0\hfill \\ \alpha +y,\hfill & 0<x_n<\alpha +y\hfill \\ 1,\hfill & x_n\alpha +y\text{ or }x_{n1}>0\hfill \end{array}`$ (14)
An important feature of the model discussed here is that one can use two inputs, $`\beta _n`$ and $`\sigma _n`$, to achieve the desired response dynamics. Although these inputs are not directly related to dynamics of specific ionic currents, they can be used to capture the collective dynamics of these currents. Selecting a proper balance between these inputs, one can model a large variety of the responses that are seen in different neurons. Again, the simplicity of the model helps one to understand the dynamical properties of each input and set the proper balance between them.

## VI Acknowledgment

The author is grateful to M.I. Rabinovich, R. Elson, A. Selverston, H.D.I. Abarbanel, P. Abbott, V.S. Afraimovich and A.R. Volkovskii for helpful discussions. This work was supported in part by U.S. Department of Energy (grant DE-FG03-95ER14516), the U.S. Army Research Office (MURI grant DAAG55-98-1-0269), and by a grant from the University of California Institute for Mexico and the United States (UC MEXUS) and the Consejo Nacional de Ciencia y Tecnologia de México (CONACYT).


text:
# Prospects of the Zee model

## 1 INTRODUCTION

The Zee model  is one of the attractive models for neutrino mass matrix. The mass matrix form is given by

$$M_\nu =\left(\begin{array}{ccc}0& a& c\\ a& 0& b\\ c& b& 0\end{array}\right),$$
$`(1.1)`$
where

$$a=f_{e\mu }(m_\mu ^2m_e^2)K,$$

$$b=f_{\mu \tau }(m_\tau ^2m_\mu ^2)K,$$
$`(1.2)`$
$$c=f_{\tau e}(m_e^2m_\tau ^2)K,$$

and $`K`$ is a common factor. The model has only 3 free parameters and it can naturally lead to a large neutrino mixing . Especially, for the case $`a=cb`$, it leads to a bi-maximal mixing

$$U\left(\begin{array}{ccc}\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& 0\\ \frac{1}{2}& \frac{1}{2}& \frac{1}{\sqrt{2}}\\ \frac{1}{2}& \frac{1}{2}& \frac{1}{\sqrt{2}}\end{array}\right),$$
$`(1.3)`$
with $`\mathrm{\Delta }m_{12}^2/\mathrm{\Delta }m_{23}^2\sqrt{2}b/a`$.

However, from the standpoint of the unification of quarks and leptons, there are some problems to be overcome.

(1) The masses of the quarks and charged leptons are generated by the Higgs mechanism, so that they depend on the Yukawa coupling constants $`y_{ij}`$, while only neutrino masses are generated radiatively, so that they depend on the Zee coupling constants $`f_{ij}`$. In general, $`f_{ij}`$ are independent of $`y_{ij}`$. We must seek for a principle which relates $`f_{ij}`$ to $`y_{ij}`$. (For such an attempt, for example, see Ref.,)

(2) The Zee model is not embedded in a Grand Unification Theory (GUT). We must seek for a room for the Zee scalar $`h^+`$ in a GUT scenario.

(3) The Zee model leads to a severe constraint  $`\mathrm{sin}^22\theta _{solar}>0.99`$ under the observed ratio $`\mathrm{\Delta }m_{solar}^2/\mathrm{\Delta }m_{atm}^21`$. The present solar neutrino data  show $`\mathrm{sin}^22\theta _{solar}0.8`$. We must investigate a modified Zee model which can be free from the severe constraint $`\mathrm{sin}^22\theta _{solar}>0.99`$.

The problem (3) has recently pointed out by the author. A parameter independent investigation leads to a severe constraint on the value of $`\mathrm{sin}^22\theta _{solar}`$

$$\mathrm{sin}^22\theta _{solar}1\frac{1}{16}\left(\frac{\mathrm{\Delta }m_{solar}^2}{\mathrm{\Delta }m_{atm}^2}\right)^2.$$
$`(1.4)`$
The conclusion cannot be loosened even if we take RGE effects into consideration: Also, we can show that the two-loop effect is negligibly small.

The simple ways to escape from the constraint (1.8) will by as follows: One is to consider  that the Yukawa vertices of the charged leptons can couple to both scalars $`\varphi _1`$ and $`\varphi _2`$. Another one  is to introduce a single right-handed neutrino $`\nu _R`$ and a second singlet Zee scalar $`S^+`$. Also, a model with a new doubly charged scalar $`k^{++}`$ is interesting because the two loop effects in such a model can give non-negligible contributions to the neutrino masses . As another attractive model, there is an idea  that in an $`R`$-parity violating SUSY model we identify Zee scalar $`h^+`$ as slepton $`\stackrel{~}{e}_R`$. Then, we can obtain additional contributions from $`d`$-quark loops to the neutrino masses. However, if we want to extend the model to a GUT scenario, we will meet a new trouble “proton decay” as stated in the next section.

Anyhow, these models are not connected to GUT scenarios. We want to embed the Zee model into a GUT scenario. For an extended Zee model based on a GUT scenario, there is, for example, the Haba-Matsuda-Tanimoto model . They have regarded the Zee scalar $`h^+`$ as a member of the messenger field $`M_{10}+\overline{M_{10}}`$ of SUSY-breaking on the basis of an SU(5) SUSY GUT. However, their model cannot escape from the constraint (1.8) as we discuss in the next section.

## 2 EXTENDED ZEE MODEL WITH GUT

We identify the Zee scalar $`h^+`$ as a member of SU(5) 10-plet scalar (including a case of $`\overline{5}+10`$ sfermion). Then, we must pay attention to the following items:

(a) Proton decay is safely forbidden.

(b) The model is free from the severe constraint $`\mathrm{sin}^22\theta _{solar}>0.99`$ in the Zee model.

### 2.1 How to avoid proton decay

Generally, the $`\overline{5}`$-plet and 10-plet scalars $`\varphi _{\overline{5}}`$ and $`\varphi _{10}`$ can couple to the $`\overline{5}+10`$-plet fermions as follows.

$$\begin{array}{c}A_5(\overline{\psi }_{\overline{5}}^c)^A(\psi _{10})_{AB}(\varphi _{\overline{5}})^B,\hfill \\ A_{10}(\overline{\psi }_{\overline{5}}^c)^A(\psi _{\overline{5}})^B(\varphi _{10})_{AB},\hfill \\ B_5\epsilon ^{ABCDE}(\overline{\psi }_{10}^c)_{AB}(\psi _{10})_{CD}(\overline{\varphi }_{\overline{5}})_E.\hfill \end{array}$$
$`(2.1)`$
The terms $`A_{10}`$ can contribute to the radiative neutrino mass, while the terms $`A_5`$ and $`B_5`$ induce the proton decay. We want the terms $`A_{10}`$, but do not want terms $`A_5`$ and $`B_5`$.

The conventional $`R`$-parity violating SUSY model  contains terms $`A_{10}`$, but also contains terms $`A_5`$ and $`B_5`$, so that the model is ruled out from the candidates of the extended Zee model based on a GUT scenario.

For example, in the Haba-Matsuda-Tanimoto mode , there is no SU(5) 5-plet scalar (except for the conventional Higgs scalars), so that there is no proton decay due to $`A_5`$ and $`B_5`$ terms. In their model, in principle, the down-quark loop diagrams can contribute to the neutrino masses in addition to the charged lepton loop diagrams. Such the down-quark loop diagrams can contribute to the non-diagonal elements of the neutrino mass matrix. However, the down-quark loop diagrams include a colored Higgs scalar (the triplet component of the SU(5) 5-plet Higgs scalar). We suppose that the colored Higgs scalar has a mass of the order of the grand unification scale in order to suppress the proton decay due to the colored Higgs scalar. Therefore, the contributions are negligibly small, so that the model cannot still be free from the severe constraint $`\mathrm{sin}^22\theta _{solar}>0.99`$.

### 2.2 $`R`$-parity violating model

Note that if we assume that $`R`$-parity violating interactions are allowed only for special generations (families), then we can choose such the $`A_5`$ and $`B_5`$ interactions as they do not contribute to the proton decay without breaking the SU(5) GUT.

For example, we assume

$$\begin{array}{c}A_5\lambda _{ij}^3(\overline{\psi }_{\overline{5}}^c)_i^A(\psi _{10})_{3AB}(\stackrel{~}{\psi }_{\overline{5}})_j^B,\hfill \\ A_{10}\lambda _{ij}^3(\overline{\psi }_{\overline{5}}^c)_i^A(\psi _{\overline{5}})_j^B(\stackrel{~}{\psi }_{10})_{3AB},\hfill \end{array}$$
$`(2.4)`$
i.e., only the interactions with the SU(5) 10-plet superfield of the third generation violate the $`R`$-parity. Besides, we must assume this expression (5.4) is true on the basis on which the up-quark mass matrix is diagonal, because if it is not so, the proton decay still occurs through the CKM mixing $`V_{13}^{}`$. In order to forbid the $`u`$-$`t`$ mixing at energy scale, we must assume additional flavor symmetries (a discrete symmetry, and so on).

Although the assumption (2.4) is somewhat unnatural, anyhow, the model can give nonvanishing diagonal elements in the neutrino mass matrix. The contributions give the form

$$M_{ij}^{e,dloop}=m_0f_if_j\mathrm{with}f_i=\lambda _{i3}^3.$$
$`(2.5)`$
This matrix form (2.5) corresponds to the Drees-Pakvasa-Tata-Veldhuis model .

### 2.3 Model with a SU(2) triplet scalar

Except for the “restricted” R-parity violating SUSY GUT model, we have failed to build a GUT scenario in which the original Zee scalar $`h^+`$ is embedded into an SU(5) 10 scalar (or SO(10) 120 scalar), as far as we adhere to the model-building without the constraint $`\mathrm{sin}^22\theta _{solar}>0.99`$.

Alternatively, we consider an SU(2) triplet scalar $`\varphi =(\varphi ^{++},\varphi ^+,\varphi ^0)`$. Usually, such a Higgs scalar can directly give Majorana masses for the left-handed neutrinos $`\nu _L`$. We assume that the potential for $`\varphi `$ does not have a minimum. Then, we obtain the following radiative masses

$$M_{ij}=f_{ij}\left((m_i^e)^2+(m_j^e)^2\right)K,$$
$`(2.6)`$
instead of the Zee mass matrix (1.1) with (1.2).

We regard the triplet scalar $`\varphi `$ as a member of the SO(10) 126-plet scalar

$$\begin{array}{c}126=1_{10}+\overline{5}_2+10_6+\overline{15}_{+6}+45_{+2}+\overline{50}_2,\hfill \end{array}$$
$`(2.7)`$
of $`SU(5)\times U(1)`$,where the triplet scalar $`\varphi `$ belongs to the component $`\overline{15}_{+6}`$ and the $`45_{+2}`$ component plays a role in the unified description of quark and charged lepton mass matrices . Of course, there are no dangerous terms for the proton decay in the coupling $`\overline{\psi }_{\overline{5}}^c\psi _{\overline{5}}\varphi _{15}`$.

However, the present model with a triplet scalar $`\varphi `$ causes another problem, because the neutral component $`\varphi ^0`$ generally takes non-vanishing value or the vacuum expectation value $`\varphi ^0`$ .

## 3 SUMMARY

The embedding of the $`R`$-parity violating SUSY models into a GUT scenario, except for a special case, is ruled out because of the proton decay. The SU(5) 10$`+\overline{10}`$ model , which is safe for the proton decay, is ruled out if the solar neutrino data establish $`\mathrm{sin}^22\theta _{solar}0.8`$. On the other hand, an SO(10) 126 model is promising if neutrino data show $`\mathrm{sin}^22\theta _{solar}0.8`$. However, such a model with a triplet scalar will cause another problems. More study is required. An SUSY GUT model with restricted $`R`$-parity violating which was discussed in 2.2 is also interesting. The details will be given elsewhere.


text:
# On the quality of random number generators with taps

## 1 Introduction

The large-scale simulations in statistical physics use some deterministic procedures to generate a sequence of uniformly distributed (pseudo)random numbers. It is possible to generate $`10^8`$ numbers per second and $`10^{15}`$ numbers in 100 days on the best processors. The widely known linear congruential methods cannot be used in a such simulations on $`32`$-bit computers because of the short period $`T<2^{32}10^9`$ of the generated sequence Knuth2 . The next most used class of generators are based on the linear recursion

$$X_n=(X_{nq}\mathrm{𝖮𝖯}X_{np})\mathrm{mod}m,$$
(1)
where $`m`$ is the word length and $`p`$ and $`q`$ are feedback taps $`p>q`$. Here OP stands for the exlusive-or operation for shift registers, $`+`$ operation for the lagged Fibonacci (LF) and $``$ operation for the substract-with-carry (SWC) recipes. In the latter case one has to add unity if the result of the previous substraction was negative. The correlations in 2-tap generators were found to be crucial in large-scale simulations of Ising model, self-avoiding random walks Grassberger  , percolation ziff-4tap  and cluster algorithms FLW  (see, also, many ). It was shown recently by Heringa, Blöte and author that these deviations can be rigorously estimated numerically and even analyticallyljh . Using this idea, Blöte and author discovered the scaling of systematic deviations for Ising modellh .

The method gives us some rigorous and quantitative way to estimate the possible level of deviations in the results of simulations. So, we propose some practical definition of goodness of a concrete recipe for generating random numbers. This definition is based on the concept of observability: the recipe should be considered as a good one if it passes all known tests using the sequence of the pseudo-random numbers generated in a couple of month on the nowadays computers. So, for today we need those recipes which do not produce observable deviations when $`10^{15}`$ random numbers are used.

## 2 Random walk test

The test is based on the one-dimensional directed random walk ljh : a walker starts at some site of an one-dimensional lattice and, at discrete times $`i`$, either he takes a step in a fixed direction with a probability $`\mu _i`$ or he stops with a probability $`\nu _i1\mu _i`$. In the case of unbiased probabilities $`\mu _i=\mu `$ and the probability of a walk with length $`n`$ is then $`P(n)=\mu ^{n1}(1\mu )`$. We could apply some particular recipe of random number generation to calculation of the probabilities $`\mu _i`$. Clearly, they would differ from $`\mu `$ because of the correlations in the recipe (1). Therefore, it was proposed in ljh ; lh  to calculate the relative deviations

$$\delta P(n)\frac{P^{}(n)}{P(n)}1$$
(2)
of the probability of the walk $`P^{}(n)=\nu _i^{}_{i=1}^n(\mu _i^{})^{i1}`$ from the unbiased probability $`P(n)`$.

Some results for SR anf LF generators from ljh  and for SWC one from lp  are presented in the Table 1. The first interesting fact is that the lagged Fibonacci and substract-with-carry leads to the same deviations up to the order of $`2^m`$, which is usually of the order of $`10^9`$. Next, the deviations for shift registers are of the same nature and and of the same order as for LF and SWC. Both these statements contradict to the widespread belief that LF and SWC produce ”better” randomness.

## 3 Scaling of the deviations

We performed Wolff simulations of the 2D Ising model at criticality, using SR with feed-back positions $`(p,q)`$=(36,11), (89,38), (127,64) and (250,103) for a variety of the lattice sizes $`L`$ lh . We found that the data could be collapsed on a single curve with the appropriate rescaling of deviations, say, of the energy $`\delta \stackrel{~}{E}p^{0.88}\delta E`$ and the system size $`\stackrel{~}{L}p^{0.43(5)}L`$. This enables us to bound deviations of energy $`E`$

$$\delta E<\mathrm{\hspace{0.17em}0.3}L^{0.84}p^{0.52},$$
(3)
of the specific heat $`C`$

$$\delta C<\mathrm{\hspace{0.17em}0.85}L^{0.21}p^{0.42},$$
(4)
and of the universal ratio $`Q=m^2^2/m^4`$, where $`m`$ is the magnetisation,

$$\delta Q<\mathrm{\hspace{0.17em}0.244}L^{0.45}p^{0.41}.$$
(5)
The analysis has been performed also for 3D Ising model, and it yields

$$\delta CL^{0.24}p^{0.19}$$
(6)
and

$$\delta QL^{1.15}p^{0.09}.$$
(7)
These estimations are very useful in the planning of computer experiments spp3 . Their knowledge is equivalent to the knowledge of the accuracy of the experimental equipment.

## 4 Consclusion

Using the method developed in ljh  we found some support for the previously proposed modifications of the 2-tapped methods. Correlations would be sharply diminished by one of the following tricks:

1. concept of luxuries lusher ; lp , i.e., generate $`r>p`$ random numbers but use only $`p`$ of them;
2. by decimation of the random number sequence, or, by using a combination of the two or more generators At ; lh ; spp3 ;
3. by using four-tap generators At ; ziff-4tap .

Our results are shown to be important for the simulations of Ising model, the percolation problem, for random walks and cluster algorithms.

So, instead of the qualifications ”bad” or ”good” we propose to apply to random number generators a more detailed description of their features lp , for example, that no deviations would be observed in the simulation of the above mentioned problems when up to $`10^{15}`$ random numbers are used.

Discussions and collaboration with H. Blöte, J. Herringa and P. Butera are kindly acknowledged.


text:
# On the origin of the broad, relativistic iron line of MCG–6-30-15 observed by XMM-Newton

## 1 Introduction

Wilms et al. (2001; cited as W01 hereafter) recently presented and discussed an extremely broad and red-shifted iron K$`\alpha `$ feature detected in the 06/11-12/2000 100 ksec XMM-Newton observation of MCG–6-30-15. This Seyfert 1 galaxy is well known for possessing one of the best studied broad iron lines, whose profile is explained by relativistic effects (Tanaka et al., 1995, Guainazzi et al., 1998; see Fabian et al., 2000, for a review). The Fe K$`\alpha `$ profile observed by XMM-Newton’s EPIC-pn camera is similar to the one observed by Iwasawa et al. (1996) using ASCA data during a short ($`15.2`$ ksec) period of low X-ray flux. W01 may have caught the source in a similar “deep minimum state”, i.e. a state in which the primary flux is lower ($`F_{210\mathrm{keV}}=2.3\times 10^{11}`$ erg s-1 cm-2) and the line Equivalent Width (EW) higher (up to $`300÷400`$ eV) than the time-averaged values.

The line profile indicates that a large fraction of the emission comes from $`r<6r_\mathrm{g}`$ ($`r_\mathrm{g}=m=\frac{GM}{c^2}`$). This implies either that the central Black Hole (BH) is rotating – thus the radius of the disc innermost stable orbit, $`r_{\mathrm{ms}}`$, lies between $`6r_\mathrm{g}`$ (=$`r_{\mathrm{ms}}`$ for a static BH) and $`1.23r_\mathrm{g}`$ ($`r_{\mathrm{ms}}`$ of a canonically spinning BH: Thorne, 1974) – or that the fluorescent line emission originates from matter falling freely below $`r_{\mathrm{ms}}`$ (Reynolds & Begelman, 1997). Sako et al. (2001) proposed that also some spectral features observed at lower energies in this source, as well as in Mkn 766, are Ly$`\alpha `$ lines of carbon, nitrogen, and oxygen, affected by relativistic broadening in the spacetime of a rotating BH.

In most works on the subject, a simple power law parameterization of the disc emissivity $`ϵ(r)r^\beta `$ is usually adopted. This is done also in the analysis of W01: letting $`\beta `$ be a free fitting parameter, they find $`\beta 4`$, a value much larger than usually found in Seyfert galaxies (Nandra et al., 1997).

In order to provide a physical picture of a so steep emissivity, W01 invoke strong magnetic stresses acting in the innermost part of the system, which dissipate a considerable amount of energy in the disc at very small radii. If the magnetic field lines thread the BH horizon, this would imply magnetic extraction of the BH rotational energy – the so called Blandford-Znajek effect (Blandford & Znajek, 1977; cited as BZ hereafter). However, the efficiency of the BZ effect has been questioned in recent years by e.g. Ghosh & Abramowicz (1997) and Livio, Ogilvie & Pringle (1999). These works argue that the electromagnetic output from the inner disc regions should in general dominate over that due to the BH. Thus the BH spin would probably be irrelevant to the expected electromagnetic power output from the system.

Krolik (1999), Agol & Krolik (2000) and Li (2000) proposed that MHD processes play a dominant role if magnetic field lines connect downfalling plasma near the hole with more distant regions: high efficiency of energy extraction can be achieved in this way even if the magnetic field does not thread the horizon itself. This magnetized accretion offers an alternative to the original BZ process and to its follow-up generalizations (e.g., Phinney 1983); however, the mechanism is violently non-stationary and such situations have not been quantitatively modelled yet (cf. Koide et al., 2000, and Tomimatsu & Takahashi, 2001, for the first attempts of such modelling).

In this paper we show that the required steep emissivity law, as well as the line EW and the amount of Compton reflection, may be reproduced with a phenomenological model in which a X-ray illuminating source is located on the BH symmetry axis (lamp-post model: Martocchia & Matt, 1996, Petrucci & Henri, 1997, Bao, Wiita & Hadrava, 1998, Reynolds et al., 1999, Dabrowski & Lasenby, 2001). This can be considered as a simplified scheme, appropriate for various physical scenarios, including the mentioned MHD energy extraction. Indeed, Agol & Krolik (2000) state that magnetized accretion may also lead to enhanced coronal activity immediately above the plunging region. “If so, this would provide a physical realization for models (…) which call for a source of hard X-rays on the system axis a few gravitational radii above the disc plane”. Alternatively, shock waves in an aborted jet close to the BH axis have been proposed as a source of the central irradiation by Henri & Petrucci (1997). This model assumes that a point source of relativistic leptons (e+,e-) illuminates the accretion disk by Inverse Compton process; the resulting angular and spectral distribution of soft and hard radiation has been derived.

## 2 A centrally illuminated disc in Kerr metric

The primary source is schematically supposed to be pointlike and located at a height $`h`$ on the system symmetry axis. This allows to estimate the main effects resulting from various degrees of anisotropy of the illumination just by varying $`h`$.

If the primary source, assumed to be isotropic in its own reference frame, is very close to the BH, a fraction of photons emitted towards infinity are deflected by the BH gravitational field and illuminate the accretion disc, at the same time increasing the number of X-rays able to produce line emission and reducing the primary radiation observed at infinity (Figure 1; see Martocchia & Matt, 1996, and Martocchia, 2000, for details). For static BHs this effect is counterbalanced by the loss of solid angle subtended by the matter to the source when the latter is very close to the BH (because $`r_{\mathrm{ms}}=6r_\mathrm{g}`$). Martocchia & Matt (1996) showed that in the case of a spinning BH the increase in the line intensity can be up to a few times the value for a static BH, and that the increase in the equivalent width can be even more dramatic. Due to the high (relativistic) orbital speed of the disc medium at very low radii, the impinging photons as seen in the matter’s frame arrive with high incident angles, which further increases the local emissivity of fluorescence and reflection. The local emissivity of fluorescence is further increased by the blue-shift of these photons.

Using the lamp-post framework, iron line profiles with the underlying Compton-reflected continuum have been presented in Martocchia, Karas & Matt (2000). Monte Carlo simulations have been used to calculate the photon transfer within disc matter on the base of Compton scattering, self-consistently taking into account GR effects such as those on the illuminating photons’ impinging angle.

### 2.1 Emissivity law

In the lamp-post as well as in more realistic models, the actual profile of the disc emissivity – which depends on the geometry of the illuminating matter and is affected by many factors, including general-relativistic (GR) effects – is more complex than a simple power-law. Appropriate emissivity laws $`ϵ(r)`$ have been presented in Martocchia (2000) and Martocchia, Karas & Matt (2000).

In particular, $`ϵ(r)`$ steepens when $`h`$ decreases, because of the enhanced anisotropy of the primary emission (see Figure 2). It can be approximated by functions of the form $`ϵ(r)=Ar^B+Cr^D`$. The best-fit coefficients for this formula, obtained by least square approximation, have been presented in Martocchia, Karas & Matt (2000).

With decreasing $`h`$, the effect of light bending is enhanced (Martocchia & Matt, 1996) and the fraction of (primary) photons impinging onto the innermost regions of the disc increases. It is easy to recognize (cf. Figure 2) that an emissivity law with $`\beta 4`$, i.e. the one derived by W01, may be produced in our model with a small height of the primary source ($`h3÷4r_\mathrm{g}`$). For an emissivity law of this kind, the line profile comes out to be very broad and red-shifted (cf. Figure 3; see Martocchia, Karas & Matt, 2000, for details). The differences between this profile and that obtained by a powerlaw emissivity with $`\beta 4`$ are clearly too small to be measurable by XMM.

### 2.2 Line EW and reflection continuum

In the previous section we have shown that an emissivity law similar to the one derived by W01 in the XMM-Newton spectrum of MCG–6-30-15 may be obtained in the framework of the simple lamp-post model. Here we would like to stress that the large observed EW and $`R`$ may be also self-consistently reproduced.

Martocchia & Matt (1996) showed that the model predicts an anti-correlation between the intensity of the reflected features (line and continuum) and the intensity of the primary flux, assuming that the latter is due to a change in the average height of the emitting matter. Indeed, for a low primary source ($`h3r_\mathrm{g}`$) values of EW$`250`$ eV or more may be easily obtained (e.g. Martocchia, Karas & Matt, 2000.) When allowing the source to be located off the axis of rotation, an even stronger enhancement can be obtained (Dabrowski & Lasenby, 2001).

The lamp-post picture predicts that the intensity of the Compton-reflected continuum correlates with the intensity of the fluorescent emission line. The best-fit value for $`R`$ found by W01 is not well constrained, but may be as high as $`7.6`$. In our model (see Figure 1), a primary source located at $`h3r_\mathrm{g}`$ produces an anisotropy ratio $`1.6`$ (thus $`\mathrm{\Omega }_{\mathrm{disc}}3.2\pi `$), i.e. we find

$$R=\frac{\mathrm{\Omega }_{\mathrm{disc}}}{4\pi \mathrm{\Omega }_{\mathrm{disc}}}4.$$

Let us finally determine the order of magnitude of the intrinsic luminosity of the point source in these assumptions. From the observed flux we can estimate the intrinsic (rest frame) value of the illumination in the 2–10 keV band by use of the anisotropy and gravitational redshift factors (the latter given by $`g_{h=3m}0.63`$, cp. Martocchia & Matt 1996), and assuming an energy index $`\alpha =1.3`$, $`H_070`$ km s-1 Mpc-1 and $`z0.008`$, we get the following primary source luminosity:

$$L_\mathrm{s}4\pi D^2F_{210\mathrm{keV}}\times g_h^{1\alpha }\times 1.65.3\times 10^{42}\mathrm{erg}\mathrm{s}^1.$$

In conclusion, we have shown that a simple lamp-post model is able to describe the results presented by W01 equally well and in a consistent way. It reproduces the observed emissivity and explains the large amount of line flux and reflection at the same time, providing that the primary X-ray source is located at $`h3r_\mathrm{g}`$.

###### Acknowledgements.

AM and GM acknowledge financial support from MURST under grant cofin-00-02-36, GM also from ASI, VK from GACR 205/00/1685. We thank the referee, J. Wilms, for his useful suggestions.


text:
# Theoretical models for Bump Cepheids

## 1 Introduction

Classical Cepheids are key objects for estimating stellar distances and for calibrating secondary distance indicators necessary to evaluate the Hubble constant (Ferrarese et al. 2000; Saha et al. 1999). From a theoretical point of view Cepheids can also be regarded as fundamental physical laboratories to assess both the accuracy of the input physics (Simon 1989) and the reliability of the physical assumptions adopted in modeling radial pulsations (Bono, Marconi, & Stellingwerf 1999, hereinafter BMS). Thus, any new test concerning the accuracy of theoretical predictions appears of paramount relevance. Moreover, this effort is mandatory to improve our confidence on Cepheids as a robust rung of the cosmic distance ladder and to disentangle relevant open questions such as the dependence of the Cepheid distance scale either on metal content (Bono et al. 1999; Bono, Castellani & Marconi 2000a; Caputo et al. 2000; Laney 2000; Storm et al. 2000) or on blending (Mochejska et al. 2000).

In this context, one has to remind that the current theoretical framework based on nonlinear, convective pulsation models, is facing with a fundamental problem, i.e. the calibration of the turbulent convection (TC) model adopted to account for the coupling between pulsation and convection (Castor 1968; Stellingwerf 1982). The approach originally suggested by Bono & Stellingwerf (1993, 1994, hereinafter BS) relies on the combination of leading physical arguments and empirical constraints on the topology of the RR Lyrae instability strip in Galactic globulars. Such a theoretical scenario accounts for some relevant properties of radial variables in the Cepheid instability strip (Bono et al. 1997a,b; BMS; Bono, Caputo, & Marconi 2001). However, the predicted RR Lyrae light curves close to the low-temperature instability edge are somehow at variance with empirical data (Kovacs & Kanbur 1998). This appears as a disturbing evidence, since empirical light curves supply tighter constraints on the accuracy of theoretical predictions than mean magnitudes and/or pulsational periods.

In a recent paper (Bono, Castellani,& Marconi 2000b, hereinafter BCM) we have approached such a problem and successfully reproduced the luminosity variation over a full pulsation cycle of the field, first overtone RR Lyrae U Com. At the same time, we also constrained the adopted TC model, and indeed in that paper we found that the Bump located just before the luminosity maximum can be reproduced by nonlinear models only by assuming a vanishing overshooting efficiency at the boundaries of the convective unstable regions. To explore the general reliability of this finding, we undertook a further investigation aimed at checking whether this TC model can also account for Cepheid light curves. To this purpose we selected two Bump Cepheids from the OGLE database (Udalski et al. 1999) of the LMC, with the Bump either along the decreasing (OGLE194103, short-period) or along the rising (OGLE56087, long-period) branch of the light curve. Table 1 lists the empirical data of the selected objects.

The reason for this choice follows from the evidence that both the amplitude and the phase of the Bump are crucial features to nail down the predicting power of theoretical models. Moreover, both stars present quite similar apparent magnitudes, and therefore the period difference suggests that the short-period Cepheid is located close to the blue (hot) edge, while the long-period one close to the red (cool) edge of the instability strip. When moving from hotter to cooler effective temperatures the thickness of the convective unstable region, and in turn the efficiency of the convective transport increases. This means that the envelopes of the two Bump Cepheids present substantially different physical structures. Finally, the intimate nature of the Hertzsprung progression is still controversial (Bono, Marconi, & Stellingwerf 2000c), and a reliable fit to empirical light curves can supply new insights into the physical mechanisms that govern the occurrence of this phenomenon.

## 2 Theoretical best fit models

To fit the empirical light curve of U Com, BCM computed several iso-period sequences of nonlinear models at fixed chemical composition. For each given assumption on the stellar mass, the iso-period sequence was constructed by changing a single input parameter, namely the effective temperature, since the luminosity was constrained by the observed period. However, the fit to Bump Cepheid light curves will be performed by adopting a different and stronger constraint. At first we undertake the fit of one of the two variables according to suitable Mass-Luminosity (ML) relations based either on canonical or on mild overshooting evolutionary models (BMS; Girardi et al. 2000), both neglecting mass-loss. This approach obviously decreases the degrees of freedom of the problem. In fact, for each selected luminosity we adopt the stellar mass as given by the ML relation, thus only one effective temperature accounts for the constraint on the observed period. If and when the best fit of the light curve of the former Cepheid is accomplished, we derive the stellar parameters (L, M, Te), and in turn an estimate of the LMC distance modulus. According to this result, the fit to the light curve of the latter object will be accomplished under the even stronger constraint to adopt the same ML relation and the same LMC distance modulus already derived for the former object.

For both stars we adopted a chemical composition typical of LMC Cepheids (Y=0.25, Z=0.008). We considered at first the short-period Cepheid OGLE194103 and constructed various iso-period ($`P=8.7d`$) sequences at different luminosity levels by adopting the canonical ML relation used in previous papers (BMS; Bono et al. 1999; Bono et al. 2000d). A glance at these iso-period sequences disclosed that the shape of the light curve is mainly governed by the stellar parameters L, M, and Te, whereas the luminosity amplitude sensitively depends on the adopted mixing length parameter. However, we found that under the quoted assumptions theoretical models appear unable to reasonably fit observational constraints. Several numerical experiments were performed to test the dependence of both pulsation amplitude and shape of the light curve on the free parameters adopted in the TC model (BS), but the shape of the predicted light curve was still at variance with the observed one.

Therefore we decided to test models with a noncanonical ML relation. To mimic the luminosity predicted by evolutionary models that account for mild convective core overshooting (Girardi et al. 2000), we increased the luminosity of canonical models, for each stellar mass, by 0.25 dex (Chiosi et al. 1993). We repeated once again the quoted experiments and now we found that our basic TC model gives a reasonable fit to the OGLE194103 light curve for the couple of parameters log $`L/L_{}`$=3.55 ($`M/M_{}`$=5.57), and $`T_e=5660`$ K, and our canonical mixing length value i.e. l=1.5 Hp (BS94). Note that the value of the effective temperature refers to the static linear model. The temperature variation along the pulsation cycle of the best fit nonlinear model is equal to 850 K. The agreement between theory and observations appears quite satisfactory over the entire pulsation cycle (see Fig. 1) and accounts for the phase and the amplitude of the Bump.

The same approach has been applied to fit the light curve of the long-period Bump Cepheid (OGLE56087, P=13.6 d). In particular, we adopted the same ML relation and the same distance modulus (DM=18.48) we found for the short-period Cepheid. The bottom panel of Fig. 1 does show that a remarkable agreement between theory and observations is obtained by adopting, log $`L/L_{}`$=3.61 ($`M/M_{}`$=5.70), and $`T_e=5160`$ K, but now we were forced to assume l=1.8 Hp. In fact, the model constructed by adopting l=1.5 Hp (dotted line) presents the same secondary features but the luminosity amplitude is systematically larger, i.e. $`\mathrm{\Delta }I`$=0.88 against the observed 0.51 mag. A glance at the data plotted in the bottom panel clearly shows the accuracy with which the best fit model does reproduce the well-defined dip of the observed light curve at $`\varphi 0.4`$, soon after the Bump secondary maximum. Present computations thus suggest that the mixing length parameter l becomes larger when a star becomes cooler, and the convective unstable region becomes thicker, i.e. the efficiency of the convection increases more than predicted by the assumption of a constant mixing length.

The use of a slightly higher metal abundance Z=0.01, according to recent empirical evidence (Luck et al. 2000) that LMC Cepheids seem to cover a wide metallicity range ($`0.006Z0.013`$), does not improve the fit. The same outcome applies to plausible changes in the He content, and indeed a change from 0.25 to 0.26 does affect neither the shape of the light curve nor the luminosity amplitude.

To investigate the robustness of current best fit solutions we also explored whether plausible changes in the input parameters still supply a reasonable fit to the observed light curves. The top panel of Fig. 2 displays that a nonlinear model constructed by adopting log $`L/L_{}`$=3.586 ($`M/M_{}`$=5.7), $`T_e`$=5780 K, and a mixing length l=1.5 Hp properly fits the luminosity variation of the short-period Cepheid (OGLE194103). This finding implies a LMC distance modulus of 18.58 mag and the fit to the long-period Cepheid (OGLE56087) now requires: log $`L/L_{}`$=3.65 ($`M/M_{}`$=5.95), $`T_e=5225`$ K, and a mixing length l=1.85 Hp (see Fig. 2). No agreement was found outside this range of stellar luminosities.

As a whole, previous numerical experiments support the evidence that the pulsational LMC distance modulus can be safely bracketed between 18.48 and 18.58 mag. Moreover, a decrease in the evolutionary mass between canonical and best fit models of the order of 15% is, in any case, required. Finally, the variation of a single free parameter, the mixing length, allows us to account for the shape of empirical Bump Cepheid light curves across the instability strip.

## 3 Discussion

According to the results given in the previous section, our pulsational computations require a larger mixing length value when the star becomes cooler. We note that this is not a very surprising result, since a mixing length roughly equal to 1.9 Hp is also required by the evolutionary models that best fit the red giant branch of Galactic globular clusters (Castellani 1999). This notwithstanding, as already discussed in BCM, a substantial variation in the mixing length parameter causes a sizable change in the efficiency of the convective transport across the driving regions, and in turn in the temperature width of the instability strip. To estimate the dependence of the instability edges on l we constructed new models with the stellar mass of the long-period best fit model. We find that when moving from l=1.5 to 1.8 Hp the temperature width decreases from 1100 K to 800 K. This change is mainly due to a shift toward hotter effective temperatures of the red boundary and could be marginally at odds with empirical estimates. In fact, estimates by Pel & Lub (1978, see their Fig. 3) do suggest that in this period range the strip is slightly larger than 800 K. Unfortunately, we are not aware of any recent empirical estimates that can allow us to supply firm constraints on this observable.

The pulsational scenario we are dealing with accounts for the empirical light curves of the selected Bump Cepheids provided that their stellar masses are 15% smaller than predicted by canonical evolutionary models. This theoretical evidence, if taken at face value, could be due either to the occurrence of mass-loss among intermediate-mass stars, as originally suggested by Cox (1980), or to convective core overshooting, or both of them. Unfortunately, the pulsational approach adopted in this investigation does not allow us to discriminate between the two different hypothesis. Current results supply some support to the evidence recently brought out by Beaulieu et al. (2001) concerning the discrepancy between evolutionary and pulsational masses among Magellanic Cepheids. On the basis of a new approach they found that evolutionary masses are up to $`\mathrm{\Delta }\mathrm{log}M/M_{}=0.1`$ larger than pulsational ones. On the other hand, we found that for the two selected Bump Cepheids the discrepancy is of the order of 0.07 and 0.08 dex for the short and the long-period respectively.

Note that Bono et al. (2001) performed a detailed analysis of both pulsation and evolutionary masses of Galactic Cepheids for which are available both accurate empirical estimates of mean radii, distances, and individual reddenings. Interestingly enough, they found that pulsational masses are approximately 10% smaller than evolutionary ones. A similar result was found by Wood, Arnold, & Sebo (1997) by fitting the light curve of a LMC Bump Cepheid (HV905). However, they found that the luminosity of this object was substantially higher than predicted by ML relations based either on canonical or on mild overshooting evolutionary models.

To assess whether our best fit models support either the occurrence of a resonance between fundamental and second overtone or the Christy’s echo mechanism (Bono et al. 2000c, and references therein) we finally constructed two new models by adopting the same input parameters of the models plotted in Fig. 1, but by perturbing the linear second overtone radial eigenfunctions. We find that both of them, after a transient phase lasting $``$200 (long-period) and $``$700 (short-period) cycles, undergo a mode switch from the second to the fundamental mode. The period ratios between second overtone and fundamental period are $`P_2/P_0(OGLE194103)`$=0.523-0.528 and $`P_2/P_0(OGLE56087)`$=0.495-0.503 respectively. For each object the former value was derived by adopting linear $`P_2`$ and $`P_0`$ periods, while the latter one with linear $`P_2`$ and nonlinear $`P_0`$ periods. The values of these ratios still fall in the so-called resonance region (Simon & Schmidt 1976). However, we note that the physical structure of the best fit models does not show any nodal line (Bono et al. 1997c) along the pulsation cycle. This finding together with the good agreement between theoretical predictions and empirical data further strengthens the plausibility that the Bumps are triggered by the Christy’s echo mechanism.

This work was supported by MURST-Cofin 2000, under the scientific project ”Stellar Observables of Cosmological Relevance”.

| TABLE 1. Empirical dataa | | |
| --- | --- | --- |
|  | OGLE194103 | OGLE56087 |
| R.A. (J2000) | $`05^h06^m33.96^s`$ | $`05^h00^m48.3^s`$ |
| DEC. (J2000) | $`68^{}25^{\mathrm{}}12.8^{\mathrm{}}`$ | $`69^{}31^{\mathrm{}}54.8^{\mathrm{}}`$ |
| Period (d) | 8.71306 | 13.62635 |
| $`m_V`$ (mag) | $`14.68\pm 0.01`$ | $`14.66\pm 0.01`$ |
| $`m_I`$ (mag) | $`13.90\pm 0.01`$ | $`13.68\pm 0.01`$ |
| $`\mathrm{\Delta }_V`$ (mag) | $`0.80`$ | 0.86 |
| $`\mathrm{\Delta }_I`$ (mag) | 0.46 | 0.51 |

a Empirical data for the two selected Bump Cepheids according to Udalski et al. (1999).