Cellular automaton model for the simulation of laser dynamics

Cellular automaton model for the simulation of laser dynamics J. L. Guisado,1,*F. Jime ´nez-Morales,1,† and J. M. Guerra2,‡ 1Departamento de Fı ´sica de la Materia Condensada, Universidad de Sevilla, P.O. Box 1065, 41080 Sevilla, Spain 2Departamento de Optica, Facultad de CC. Fı ´sicas, Universidad Complutense de Madrid, 28040 Madrid, Spain 共Received 28 November 2002; published 20 June 2003兲 The classical modeling approach for laser study relies on the differential equations. In this paper, a cellular automaton model is proposed as an alternative for the simulation of population dynamics. Even though the model is simplified it captures the essence of laser phenomenology: 共i兲there is a threshold pumping rate that depends inversely on the decaying lifetime of the atoms and the photons; and 共ii兲depending on these lifetimes and on the pumping rate, a constant or an oscillatory behavior can be observed. More complex behaviors such as spiking and pattern formation can also be studied with the cellular automaton model. DOI: 10.1103/PhysRevE.67.066708 PACS number共s兲: 07.05.Tp, 42.55.Ah, 42.60.Jf, 42.65.Sf I. INTRODUCTION In a laser system, the interactions among simple atoms and the radiation they produce can give rise to cooperative phenomena 关1兴. However, the usual approach for its study is based on very detailed microscopical equations, which somehow mask the action of such cooperative properties. In this study, a simple cellular automaton model is presented, which reproduces much of the laser phenomenology, special attention being focused on these self-organizing cooperative effects. Our model is interesting in that it illustrates the emergence of laser properties as cooperative phenomena based on simple underlying rules. It can also be useful in calculating laser output in situations which are difficult to treat with the traditional approach based on the resolution of detailed differential equations 关2,3兴. One example is when dealing with complex boundary conditions or numerical difficulties in the integration of the equations. Cellular automata 关4,5兴are a class of spatially and temporally discrete mathematical systems, characterized by local interaction and synchronous dynamical evolution. They have the ability to generate a very complex behavior from sets of components that interact locally with relatively simple rules. They provide good models for a wide variety of physical systems 关6,7兴exhibiting cooperative phenomena, such as magnetization in solids 关8兴, reaction-diffusion processes 关9兴, fluid dynamics for complex situations 关10兴, growth phenomena 关11兴, etc. As an example of an application particularly close to laser dynamics, a cellular automaton model has been recently applied to successfully study the excited-state dynamics of atomic oxygen, which play a prominent role in creating the aurora borealis 关12兴. This paper is organized as follows. Section II presents a brief review of laser dynamics including the laser phenomenology that our model attempts to reproduce. Section III describes the cellular automaton model. The simulations are carried out and their results are presented in Sec. IV. Finally, the conclusions of this study are explained in Sec. V. II. LASER DYNAMICS A laser is a device that generates or amplifies electromagnetic radiation based on the stimulated emission phenomenon. The basic components of a laser system are 共1兲alaser medium—an appropriate collection of atoms, molecules, ions or a semiconductor crystal 共these elements will generally be referred to as ‘‘atoms’’兲;共2兲apumping process that excites electrons from these atoms to upper energy levels, due to some external electrical, optical, or chemical energy source; 共3兲optical feedback elements that reflect repeatedly the radiation beam into the laser medium 共in a laser oscillator兲,or allow it to pass only once through it 共in a laser amplifier兲. The working principle of laser is stimulated emission, i.e., an excited atom can decay to a lower state stimulated by the presence of a photon with energy equal to the difference between the two energy levels, emitting a second photon with the same frequency and propagating in the same direction. The process of absorption has the same probability, so stimulated emission dominates only when a population inversion is induced in the material by some pumping mechanism. A simplified but yet realistic model of many real lasers is the four-level laser system shown in Fig. 1. The population dynamics of a laser 共the variation with time in the number of laser photons and in the population inversion, or number of electrons in the upper laser level minus the number of electrons in the lower laser level兲is usually described as a system of coupled differential equations called laser rate equations. The rate equations can be put into their simplest forms when the lifetimes of E1and E3level electrons are negligible as compared to the lifetime of E2level electrons. Therefore the E1level electron population is N1⯝0, and thus the population inversion is approximately equal to the upper laser level population N(t)⫽N2(t)⫺N1(t)⯝N2(t); and the absorption of laser photons by electrons at level E1is negligible. Additionally, the pumping into level E2can be described by a simple pumping rate R. Under these assumptions, which are quite realistic, only *Present address: Centro Universitario de Me ´rida, Universidad de Extremadura, 06800 Me ´rida 共Badajoz兲, Spain. Electronic address: [email protected] †Electronic address: [email protected] ‡Electronic address: [email protected] PHYSICAL REVIEW E 67, 066708 共2003兲 1063-651X/2003/67共6兲/066708共8兲/$20.00 ©2003 The American Physical Society67 066708-1 levels E0and E2can be considered, and the rate equations are a system of two coupled differential equations 关2,3兴: dn共t兲 dt ⫽KN共t兲n共t兲⫺n共t兲 ␶ c,共1兲 dN共t兲 dt ⫽R⫺N共t兲 ␶ a ⫺KN共t兲n共t兲.共2兲 Here Ris the pumping rate, ␶ ais the decay time of the upper laser level (E2), ␶ cis the decay time of the laser photons in the cavity, and Kis a constant called ‘‘coupling constant.’’ The first equation reflects the variation with time of the number of laser photons, n(t), which is related to the laser beam intensity. The second equation represents the temporal variation of the population inversion N(t). This is a set of two coupled nonlinear equations due to the presence of the product term KN(t)n(t). Laser operation does not generally produce a smooth and continuous response, but exhibits different kinds of characteristic transient or modulation behaviors such as spiking, relaxation oscillation, gain switching. Laser spiking, in particular, refers to the pulses that typically occur during the initial turn-on phase of many lasers. In these cases, the laser signal presents a sequence of sharp, large-amplitude narrow pulses or ‘‘spikes.’’ Each spike is typically a fraction of a microsecond wide, and they are separated by a few microseconds. In some lasers, the spiking behavior jumps in an erratic way. In other cases, or under more stable conditions, it is possible to obtain a more regular spiking behavior, in which the amplitude of the spikes gradually damps out with time into a relaxation process. These are sometimes called relaxation oscillations, and include lasers with largeamplitude oscillations when enough time has elapsed for the amplitude to relax. This kind of behavior varies for the different types of lasers. The spiking oscillatory response is characteristic of most solid-state and semiconductor lasers, which have a substantially longer decay time of the upper laser level than the decay time of the photons in the laser cavity: ␶ aⰇ ␶ c. Laser spiking and relaxation oscillations are not observed, on the other hand, in most gas lasers, which usually have decay times of the same order of magnitude, ␶ a⬇ ␶ c. The steady-state solution to the rate equations above the threshold 共i.e., with laser emission兲is Ns⫽1 K ␶ c,共3兲 ns⫽ ␶ c 冉 R⫺1 K ␶ a ␶ c 冊 .共4兲 In order for this solution to have a physical meaning, the number of photons, n, has to be greater than zero. The minimum value of Rproducing a laser emission in the cavity, known as threshold pumping rate Rtis found as Rt⫽1 K ␶ a ␶ c.共5兲 For the case when the spiking behavior in the laser has relaxed to small-amplitude fluctuations, a linearized smallsignal analysis of the rate equations about the steady-state solution can be carried out. Following two different situations arise. 共a兲Nonspiking lasers: For lasers in which ␶ a⬇ ␶ c, the solutions for n(t) and N(t) are real exponentials. The system is overdamped and the response of the laser to any perturbation dies out exponentially towards the steady state. The solution is a stable node. 共b兲Relaxation oscillations: When ␶ aⰇ ␶ cthe solutions have an exponentially damped sinusoidal form towards the steady-state values. The system responds to any perturbation exhibiting relaxation oscillations. The solution is a stable focus. The relaxation-oscillation frequency, often referred to as spiking frequency, is found as ␻ sp⫽ 冑 冉 1 2 ␶ a R Rt 冊 2 ⫺1 ␶ a ␶ c 冉 R Rt ⫺1 冊 .共6兲 The necessary condition for the appearance of the relaxation oscillations is that the quantity inside the square root be positive, yielding ␶ a ␶ c ⬎ 冉 R Rt 冊 2 4 冉 R Rt ⫺1 冊 .共7兲 III. CELLULAR AUTOMATA MODEL Cellular automata 共CA兲are fully discrete dynamical systems, where the states are chosen in a finite set and distributed on a discrete regular lattice, and the time evolution is run synchronously in all the sites of this lattice and each site FIG. 1. Schematic view of a four-level laser system. GUISADO, JIME ´NEZ-MORALES, AND GUERRA PHYSICAL REVIEW E 67, 066708 共2003兲 066708-2 changes its state according to a local rule that only depends upon its neighboring values. Despite the simplicity of their construction, CA are found to be capable of a diverse and complex behavior and are often used as a prototype for the analysis of the spontaneous emergence of ordered behavior in spatially extended systems that are locally coupled. A. States of the cells We consider a system composed of a square lattice of Nc⫽L⫻Lcells with periodic boundary conditions. Two variables ai(t) and ci(t) are associated to each node of this lattice. The first one, ai(t), represents the state of the electron in node iat a given time t. An electron in the laser ground state takes the value ai(t)⫽0, while an electron in the upper laser state takes the value ai(t)⫽1. A temporal variable a ˜ i(t)苸 兵 0,1,2,..., ␶ a 其 is also introduced in order to take into account the finite time ␶ a, for which an electron can remain in the upper state. If the electron is in the base state, then a ˜ i(t)⫽0, otherwise a ˜ i(t⫹1)⫽a ˜ i(t)⫹1 until the maximum value ␶ ais reached and then a ˜ i(t⫹1)⫽0. The second variable ci(t)苸 兵 0,1,2,...,M 其 represents the number of photons in node iat time t. A large enough upper value of Mis taken to avoid saturation of the system. There is also another temporal variable c ˜ i j(t)苸 兵 0,1,2,..., ␶ c 其 , which measures the amount of time since a photon j 苸 兵 1,2,...,M 其 was created at node i. ␶ cis the lifetime of each photon. For a given photon j,c ˜ i j(t⫹1)⫽c ˜ i j(t)⫹1 until the lifetime ␶ cis reached and then c ˜ i j(t⫹1)⫽0. B. Neighborhood The neighborhood considered is the Moore neighborhood, each cell having nine neighbors: the cell itself, its four nearest neighbors 共situated in the positions north, south, east, and west兲, and the four next neighbors 共in the positions northeast, southeast, northwest, and southwest兲. C. Transition rules The time evolution of the CA is given by a set of rules which determines the state of any particular cell of the system at time t⫹1, depending on the state of the cells included in its neighborhood at time t. The evolution of the temporal variables a ˜ i(t) and c ˜ i j(t) was described beforehand. Here we describe only the evolution of ai(t) and ci(t). 共i兲If ai(t)⫽0, then ai共t⫹1兲⫽1 with probability ␭. 共ii兲If ai(t)⫽1 and ⌫i(t)⫽兺neighborsci(t)⬎ ␦ , then ci共t⫹1兲⫽ci共t兲⫹1, ai共t⫹1兲⫽0, where ⌫iis the number of photons included in the neighborhood of the cell iand ␦ is a threshold value, which has been taken to be 1 in our simulations. 共iii兲If ci(t)⬎0 and there is one photon jfor which c ˜ i j(t)⫽ ␶ c, then ci共t⫹1兲⫽ci共t兲⫺1. 共iv兲If ai(t)⫽1 and a ˜ i(t)⫽ ␶ a, then ai共t⫹1兲⫽0. These transition rules represent the different physical processes that function at the microscopical level in a laser system. Rule 共i兲represents the pumping process. Rule 共ii兲models the stimulated emission—if the electronic state of a cell has a value of ai(t)⫽1 and the number of laser photons in the neighborhood is greater than a certain threshold, then at the time t⫹1 a new photon will be created in that cell and the electron will decay to the ground level. Rule 共iii兲represents the photon decay. Rule 共iv兲represents the electron decay in a way similar to rule 共iii兲—after time ␶ aof the excited electrons, those electrons will decay to the ground level. To simplify the model as much as possible, we consider this decay to be entirely nonradiative, i.e., spontaneous emission is not taken into account. Also, as in an ideal four-level laser the population of level E1is negligible, stimulated absorption has not been considered. Additionally, in order to represent the noise level observed in practice 共typically of the order of one noise photon per cavity mode兲, we introduce a small continuous noise level of random photons in the laser mode at every time step. This is done by setting ci(t⫹1)⫽ci(t)⫹1 for a number of cells smaller than 0.01% of the total number of system cells, whose positions are randomly chosen. This noise 共along with the population inversion induced by pumping兲is responsible for the initial start up of the laser process, and can prevent the eventual extinction of this process, in case the number of laser photons drops down to negligible values. IV. SIMULATIONS The simulations were carried out using lattices of 200 ⫻200 and 300⫻300 cells. Three parameters determine the response of the system: the pumping probability ␭, the lifetime of photons ( ␶ c), and the lifetime of excited electrons ( ␶ a). The pumping probability ␭introduced in the CA model is linearly related to the constant pumping rate Rincluded in the rate equations, as explained in the Appendix. Thus the laser rate equations’ prediction 关Eq. 共5兲兴 for the dependence of threshold pumping rate Rtwith the lifetimes ␶ aand ␶ ccan be expressed for the threshold pumping probability ␭tas ␭t⫽1 K⬘ ␶ a ␶ c ,共8兲 where K⬘is a constant. Figure 2 shows the threshold pumping probability ␭tversus the lifetime of the upper laser level ( ␶ a) and the cavity lifetime ␶ c. Each point has been obtained in the following manner. After a transient time, the system evolves for 200 time steps and the average number of laser photons (n ¯ )is CELLULAR AUTOMATON MODEL FOR THE SIMULATION . . . PHYSICAL REVIEW E 67, 066708 共2003兲 066708-3 recorded. For each pair of values ␶ cand ␶ a, the procedure is repeated for values of the pumping probabilities ␭ranging from 0.0001 to 0.1. If the laser action does not initiate in the system, n ¯ is approximately equal to the number of noise photons (nnp) introduced. When the laser action starts, photons are produced by stimulated emission 关via rule 共ii兲兴, and the average number of photons is greater than nnp . Thus, ␭t can be computed as the minimum value of the pumping probability for which the average number of photons is higher than a given reference value nmin . In our simulations, we have taken nmin⫽1.25nnp . The number of noise photons (nnp) is the number of noise photons introduced in each time step multiplied by the lifetime of each laser photon. Figure 3 shows the threshold pumping probability ␭tversus ␶ a ␶ con a logarithmic scale. It can be seen that all the different curves of Fig. 2 collapse in a unique straight line, having a slope close to ⫺1, which is in agreement with the behavior predicted by the laser rate equations, i.e., Eq. 共8兲. To analyze the different dynamic behaviors of the system, we allow it to evolve for 500 time steps. Initially all the electrons are in the ground level 关ai(0)⫽0兴and there are no laser photons present 关ci(0)⫽0兴, except a small fraction (0.01%) of noise photons. In each experiment, the total number of laser photons, n(t)⫽兺i⫽1 Ncci(t), and the total number of electrons in the upper laser state 共population inversion兲,N(t)⫽兺i⫽1 Ncai(t) are measured. In order to classify the range of values of the parameters for which oscillations or a constant regime in n(t), and N(t) appear, the Shannon entropy Sof the distribution of values taken by n(t) and N(t) is calculated. This is carried out by dividing the range of values taken in 103intervals, and computing the frequency fiat which the magnitude value lies inside every particular nonvoid interval i. The Shannon entropy is expressed as S共␭, ␶ c, ␶ a兲⫽⫺兺 i⫽1 m filog2fi,共9兲 where mis the number of nonvoid intervals. Figure 4 shows the contour plot of the Shannon entropy of the values distribution of the number of laser photons, Sc, for a fixed value of ␶ c⫽10. Similar plots are obtained for other values of ␶ c. The black line is the theoretical stability curve: ␶ a ␶ c ⫽ 冉 R Rt 冊 2 4 冉 R Rt ⫺1 冊 , which follows from Eq. 共7兲. Here R/Rtis used instead of ␭/␭t, because Eq. 共A5兲from the Appendix, where (␭/␭t) ⫽(R/Rt), was taken into account. The Shannon entropy measures the dispersion in the distribution of the number of laser photons 共or in the population inversion兲. If this number is approximately constant, Swill tend toward zero; if oscillations appear, Stakes higher values. The maximum value would result from an equiprobable distribution. Sis thus a FIG. 2. Threshold pumping probability ␭tfrom the CA laser model: 共a兲␭ tdependence on the upper laser level lifetime ␶ afor different values of the cavity lifetime ␶ c;共b兲␭ tdependence on ␶ c for different values of ␶ a. Both figures are plotted on a logarithmic scale. ␶ aand ␶ care measured in time steps. FIG. 3. Dependence of the threshold pumping probability ␭t from the CA laser model versus ␶ a ␶ c, plotted on a logarithmic scale. The dots correspond to the results of the simulations and the continuous line to the theoretical prediction of the laser rate equations, i.e., Eq. 共8兲. ␶ aand ␶ care measured in time steps. GUISADO, JIME ´NEZ-MORALES, AND GUERRA PHYSICAL REVIEW E 67, 066708 共2003兲 066708-4 good indicator of the presence of oscillations in the system, as it must be low when the number of photons 共or the population inversion兲is essentially constant and high when this number oscillates. Figure 4 shows the presence of oscillations above and to the right of the theoretical stability curve 共dark zones兲, where ␶ a⬎ ␶ cand R⬎Rt. On the other hand, in the bright zones, where ␶ a⬇ ␶ cor R⭐Rt, the number of laser photons does not oscillate. This is in good agreement with the theoretical predictions of Eq. 共7兲, based on the linearization of the laser rate equations, and with the behavior observed in real systems. Laser spiking is characteristic of most solid-state and semiconductor lasers, which have a substantially longer decay time of the upper laser level than that of the photons in the laser cavity. This is not generally observed in gas lasers, whose decay times usually have the same order of magnitude. Figures 5 and 6 show the temporal evolution of the number of laser photons and the population inversion for two different values of the system parameters, characteristic of the two different regimes in which it can operate. In this case, a lattice of 300⫻300 cells has been used in order to define the behavior more clearly. In both figures, the temporal evolution is shown at the left, and the evolution in a phase space with the number of laser photons versus the population inversion is shown at the right. Figure 5 corresponds to the parameter values for which the Shannon entropy is low 共bright zones in Fig. 4兲and, after an initial gain switching peak, the system shows a constant behavior. In phase space, after an initial transient the system goes to a fixed point. Figure 6 corresponds to the parameters values for which the Shannon entropy is high 共dark zones in Fig. 4兲and the system exhibits oscillations. In this particular case, the system shows damped oscillations which are associated to a spiral trajectory towards a steady state in phase space. For the values of the parameters corresponding to Fig. 6, a threshold pumping probability of ␭t⫽0.0015 is found. According to Eq. 共A5兲, this means (R/Rt)⫽(␭/␭t)⫽(0.01/0.0015), which leads to a spiking period value predicted theoretically by Eq. 共6兲of Tsp th⫽2 ␲ / ␻ sp th⫽137.3 time steps. The average spiking period of the oscillations resulting from our simulations is Tsp sim⫽(127⫾10) time steps, which is in good agreement with the theoretical value predicted by Eq. 共6兲. The CA model is also very useful for observing the evolution of spatiotemporal patterns. Recently, quasiinstantaneous transverse patterns in a broad aperture laser have been measured for the first time 关13兴. In the case of a CO2laser, a disordered distribution of spots, nonreproducible from shot to shot, which yields a boundary-determined regular structure when integrated in progressively longer time windows, has been found. Figure 7 shows some snapshots of the instantaneous transverse patterns shown by the laser photon population in the lattice. Black color represents a cell with no photon present, gray represents cells with one photon, and white represents cells with two or more photons. Figures 7共a兲and 7共b兲correFIG. 4. Contour plot of the Shannon entropy of the distribution of the number of laser photons for a fixed value of ␶ c⫽10. Low values of Sc共bright zones兲indicate that the response of the system is nonoscillatory, while high values 共dark zones兲indicate an oscillatory response. The black line is the theoretical stability curve. FIG. 5. Evolution of the system for the parameters: ␭⫽0.1, ␶ c⫽8, ␶ a⫽30, where ␶ cand ␶ aare measured in time steps. Lattice: 300 ⫻300 cells. 共a兲Number of laser photons and population inversion versus time. 共b兲Evolution in a phase space formed by the number of laser photons versus the population inversion. CELLULAR AUTOMATON MODEL FOR THE SIMULATION . . . PHYSICAL REVIEW E 67, 066708 共2003兲 066708-5 spond to the oscillatory behavior after a transient time. In this case, the system shows spatial structures that evolve in time resembling wave fronts that interact destructively with each other. The typical scale of these patterns is dependent on the value of the parameters, as can be seen by comparing Figs. 7共a兲and 7共b兲. Figures 7共c兲and 7共d兲correspond to the nonoscillatory behavior and both figures have the same parameters values. Figure 7共c兲was taken at time t⫽50 time FIG. 6. Evolution of the system for the parameters: ␭⫽0.01, ␶ c⫽14, ␶ a⫽160, where ␶ cand ␶ aare measured in time steps. Lattice: 300⫻300 cells. 共a兲Number of laser photons and population inversion versus time. 共b兲Evolution in phase space. FIG. 7. Instantaneous transverse spatiotemporal patterns shown by the laser photon population in the system. 共a兲and 共b兲correspond to oscillatory behavior, whereas 共c兲and 共d兲correspond to nonoscillatory behavior. The values of the parameters are 共a兲␭ ⫽0.01, ␶ c⫽14, ␶ a⫽160; 共b兲␭ ⫽0.03, ␶ c⫽6, ␶ a⫽150; 共c兲␭ ⫽0.1, ␶ c⫽8, ␶ a⫽30, t⫽50; 共d兲 ␭⫽0.1, ␶ c⫽8, ␶ a⫽30, t⫽200. ␶ c, ␶ a, and tare measured in time steps. GUISADO, JIME ´NEZ-MORALES, AND GUERRA PHYSICAL REVIEW E 67, 066708 共2003兲 066708-6 steps, whereas Fig. 7共d兲at time t⫽200 time steps. In this case, it can be observed that the laser process is activated initially with expanding avalanches and then reaches a stationary state at which the laser photons distribution is homogeneous. V. DISCUSSION Depending on the values of the three parameters ␶ a, ␶ c,␭ of the system, our simulations show that two distinct characteristic behaviors appear. The Shannon entropy of the laser photons distribution and of the population inversion allows us to obtain the parameter ranges in which each behavior takes place. A high Shannon entropy 共dark zones in Fig. 4兲indicates a greater dispersion in the populations than in other areas of the parameters space. For these values of the parameters, the CA model reproduces the typical kind of behavior known as laser spiking—the temporal evolution of the populations of laser photons and inversion electrons oscillate in a correlated way, as shown in Fig. 6共a兲. A series of sharp, narrow ‘‘spikes’’ can be observed, whose amplitude decreases with time. In the phase space description 关Fig. 6共b兲兴, the system follows a spiral trajectory which converges towards a steadystate limit point. On the other hand, a low Shannon entropy 共bright zones in Fig. 4兲indicates a smaller dispersion in the populations. For these parameter values, the behavior shown by the CA model corresponds to a constant regime—after an initial transient, the total number of laser photons and inversion electrons remain approximately constant, as is shown in Fig. 5共a兲. In the phase space description, Fig. 5共b兲, after a short transient the system reaches a fixed point. Our model is successful in showing either one or the other of these two regimes, depending on the characteristic decay times involved. Laser spiking appears generally in situations in which the lifetime of the electrons in the excited state is substantially longer than the lifetime of the laser photons, whereas the constant regime mostly appears when both lifetimes are of the same order of magnitude. The CA model can be applied during the interesting period of strong spiking, for which there is no simple analytic solution for the laser rate equations and the traditional approach relies on the numerical integration of the system of differential equations. In Fig. 7, we can observe that the kind of spatiotemporal transverse patterns shown by the laser photons population depend on the population dynamics regime. For the laser spiking regime, spatiotemporal structures appear in the form of propagating wave fronts, Figs. 7共a兲and 7共b兲. The characteristic scale of these structures depends on the three parameters of the system. It is remarkable that our rather simple and generic laser model produces spatiotemporal structures similar to those observed experimentally and with the numerical integration of very detailed Maxwell-Bloch equations models 关13,14兴. On the other hand, when the system is in the constant regime, there is a homogeneous distribution of laser photons, Fig. 7共d兲. This is in agreement with the approximately constant values shown in the total number of laser photons. The laser process is initially activated with expanding avalanches, as shown in Fig. 7共c兲. When the avalanches have spread throughout the whole system, the constant steady state is reached. In conclusion, we have presented here a CA model for the simulation of the population dynamics in a laser cavity. While simplified, the model is able to capture the essential features and phenomenology encountered in lasers such as relaxation oscillations, spiking behavior, and pattern formation. The CA approach for lasers is a very promising tool that can supplement the traditional approach based on the solution of sets of coupled differential equations. Its performance is very good, so that the simulation can be carried out on large lattices. The full parallelism of the CA can be implemented on special purpose computers or in computers with full parallel architectures, and therefore an important saving in computer time can be obtained. A CA model can represent an advantage in cases in which the system of differential equations have convergence problems, for example, stiff differential equations. The Maxwell-Bloch equations that describe many laser systems are of this type. Our model is free from this kind of convergence problems, so it can be very useful for this kind of lasers. However, work must be done in order to test this potential application, including a quantitative comparison of the computing times of the CA model versus the numerical integration of the full set of MaxwellBloch equations. In the current state of our work, we see no obstacle to fine tune the proposed CA model in order to meet the characteristics of any particular laser system. This must be addressed in the near future. Additionally, we hope that similar models can be used to study other problems of current interest, such as chaotic laser dynamics 关15兴or the dynamics of a twophoton laser 关16兴. ACKNOWLEDGMENT J. L. Guisado would like to thank Ramon Risco for helpful suggestions and discussions. APPENDIX: RELATIONSHIP BETWEEN THE PUMPING RATE RAND THE PUMPING PROBABILITY ␭ Let us see how the pumping rate R, included in the laser rate equations, is related to the pumping probability ␭, introduced in the CA model. Ris defined as the number of electrons that are excited from the ground state E0to the state E3 per unit time, considered as a constant in the rate equations: R⫽dn3 dt ⇒n3⫽Rt.共A1兲 In the CA model, pumping is introduced as a probabilistic process, where ␭is the pumping probability per unit time of an electron from state E0to E3. This means that if there are n0electrons in E0, then the number of excitations in a time interval dt will be n␭dt. So the population of state E0decreases as CELLULAR AUTOMATON MODEL FOR THE SIMULATION . . . PHYSICAL REVIEW E 67, 066708 共2003兲 066708-7 dn0⫽⫺␭dtn0⇒n0⫽n0ie⫺␭t.共A2兲 If we assume that n3is initially unpopulated, then n0i ⫽n0⫹n3and the population of E3will increase with pumping as n3⫽n0i⫺n0⫽n0i共1⫺e⫺␭t兲.共A3兲 Then, the pumping rate Rwill be R⫽dn3 dt ⫽␭n0ie⫺␭t.共A4兲 In general, if the only process in the system is the pumping, Rwill decrease exponentially with time as the population of electrons in E0decreases. But the electrons pumped to level E3decay down again to level E0repopulating the ground level, and 共as long as the pumping probability is small兲after a transient period, Rstabilizes to a constant value: R⫽␭n0c⯝const, 共A5兲 where n0cstands for the approximately constant value toward which the population of level E0tends. So after a transient period, there is a linear relationship between the pumping probability considered in our model and the constant pumping rate considered in the laser rate equations. This means that the pumping probabilities from the CA model must follow Eq. 共5兲with the exception of the value of the constant K: ␭t⫽1 K⬘ ␶ a ␶ c .共A6兲 关1兴H. Haken, Synergetics: An Introduction. Nonequilibrium Phase Transitions and Self-organization in Physics, Chemistry and Biology 共Springer-Verlag, Berlin, 1983兲. 关2兴A. Siegman, Lasers 共University Science Books, Mill Valley, CA 1986兲. 关3兴O. Svelto, Principles of Lasers 共Plenum Press, New York, 1989兲. 关4兴S. Wolfram, Physica D 10,1共1984兲. 关5兴S. Wolfram, Cellular Automata and Complexity 共AddisonWesley, New York, 1994兲. 关6兴B. Chopard and M. Droz, Cellular Automata Modeling of Physical Systems 共Cambridge University Press, Cambridge, 1998兲. 关7兴T. Toffoli and N. Margolus, Cellular Automata Machines: A New Environment for Modelling 共MIT, Cambridge, MA, 1987兲. 关8兴G. Vichniac, Physica D 10,96共1984兲. 关9兴B. Chopard, P. Luthi, and M. Droz, Phys. Rev. Lett. 72, 1384 共1994兲. 关10兴A.K. Gunstensen, D.H. Rothman, S. Zaleski, and G. Zanetti, Phys. Rev. A 43, 4320 共1991兲. 关11兴J. Gravner and D. Griffeath, Adv. Appl. Math. 21, 241 共1998兲. 关12兴P. Seybold, L. Kier, and C.-K. Cheng, Int. J. Quantum Chem. 75, 751 共1999兲. 关13兴F. Encinas-Sanz, I. Leyva, and J.M. Guerra, Phys. Rev. A 62, 043821 共2000兲. 关14兴F. Encinas-Sanz, I. Leyva, and J.M. Guerra, Phys. Rev. Lett. 84, 883 共2000兲. 关15兴G. Huyet and J. Tredicce, Physica D 96, 209 共1996兲. 关16兴O. Pfister, W.J. Brown, M.D. Stenner, and D.J. Gauthier, Phys. Rev. Lett. 86, 4512 共2001兲. GUISADO, JIME ´NEZ-MORALES, AND GUERRA PHYSICAL REVIEW E 67, 066708 共2003兲 066708-8 View publication statsView publication stats