Computational simulation of laser dynamics as a cooperative phenomenon

Computational Simulation of Laser Dynamics as a Cooperative Phenomenon J. L. Guisado*1, F. Jim´enez-Morales†2 and J. M. Guerra‡3 1Centro Universitario de M´erida, Universidad de Extremadura, 06800 M´erida (Badajoz), Spain 2Departamento de F´ısica de la Materia Condensada, Universidad de Sevilla, P.O. Box 1065, 41080 Sevilla, Spain 3Departamento de Optica, Facultad de CC. F´ısicas, Universidad Complutense de Madrid, 28040 Madrid, Spain Abstract The different kinds of behavior exhibited by the system in a laser dynamics simulation using a cellular automata model are analyzed. Three distinct types of behavior have been found: laser constant operation, laser spiking and a complex behavior showing irregular oscillations. In the last case, the power spectrum follows a power law of the type 1/f −with exponent close to =2. In the laser spiking regime, the dependence of the decay rate of the oscillations is found to be in good agreement with the predictions of the theoretical laser rate equations and the experimental phenomenology. In our model the system components evolve under local rules which reproduce the physics of the laser system at the microscopic level, and the laser properties appear as cooperative emergent phenomena associated to these rules. 1. Introduction The standard approach to model laser systems is based on coupled differential equations (Maxwell-Bloch equations) [1, 2]. Nevertheless, under certain circumstances, it is not possible to obtain satisfactory results with this kind of description. In the one hand, there are many lasers for which the corresponding Maxwell-Bloch equations are stiff differential equations, very difficult to integrate. In the other hand, the most recent devices usually have physical dimensions comparable to or even smaller than the wavelength of light, with an active medium of arbitrary geometry, and working with a small number of photons. Under these conditions, the standard approach (which usually involves approximations) may not be entirely applicable and therefore can give rise to not very accurate results. A different approach for laser dynamics, based on a cellular automata (CA) model, has been recently presented [3]. This CA approach can be an alternative to the standard treatment based on differential equations for the cases involving the difficulties mentioned above. Cellular automata are a class of spatially and temporally discrete mathematical systems, characterized by local interaction and synchronous dynamical evolution [4, 5]. Thanks to their capability to generate complex behavior from sets of components which interact locally with relatively simple rules, CA provide an excellent modeling approach for complex systems. In the last two decades, CA have been used to build models of a wide variety of physical systems, such as reactiondiffusion processes, fluid dynamics, magnetization in solids, growth phenomena, molecular excited-state dynamics, etc. [6, 7]. Additionally, CA are intrinsic parallel systems, very suitable to be naturally and efficiently implemented in parallel computers to carry out high performance simulations [8, 9]. *E-mail: [email protected] †E-mail: [email protected] ‡E-mail: jmguerra@fis.ucm.es In this paper we supplement the results presented in [3] by analysing in more detail the distinct kinds of behavior that are exhibited by the cellular automata model of laser dynamics that was introduced there. The outline of the paper is the following: in Section 2 the cellular automata model is presented; in Section 3 the results of the simulations are shown and the different kinds of behavior exhibited by the system are analyzed; finally, the conclusions are summarized in Section 4. 2. Cellular automata model The CA model that has been used in our simulations has been presented in detail in reference [3] and thus will only be summarized here. The system is modeled by a cellular automaton defined on a two-dimensional square lattice of Nc=400 ×400 cells with periodic boundary conditions. Two variables ai(t) and ci(t) are associated to each node. ai(t) represents the state of the electron in node iat time t:ifai(t)=0 the electron is in the laser ground state and if ai(t)=1 it is in the upper laser state. 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. The neighborhood considered is the Moore neighborhood, each cell having nine neighbors: the cell itself, its four nearest neighbors and the four next neighbors. The time evolution of the CA is given by a set of transition rules which determine the state of a cell at time t+1 depending on the state of the cells included in its neighborhood at time t. These rules represent the different physical processes that work at the microscopic level in a laser system: rR1. Pumping: If the electronic state of a cell has a value of ai(t)=0 in time t, then in time t+1 that state will have a value of ai(t+1) =1 with a probability . rR2. Stimulated emission: If, in time t, the electronic state of a cell has a value of ai(t)=1 and the sum of the values of the laser photons states in the nine neighbor cells is greater than a certain threshold (which in our simulations has been taken to be 1), then in time t+1 a new photon will be created in that cell: ci(t+1) =ci(t)+1 and the electron will decay to the ground level: ai(t+1) =0. rR3. Photon decay: A finite life time cis assigned to each photon when it is created. The photon will be destroyed c time steps after it is created. rR4. Electron decay: A finite life time ais assigned to each electron that is promoted from the ground level to the upper laser level. That electron will decay to the ground level again atime steps after it was promoted, if it has not yet decayed by stimulated emission. The alifetime include spontaneous radiative and eventually non-radiative processes.As in an ideal four level laser the population of the lower laser level is negligible, stimulated absorption has not been considered. Spontaneous emission as well as thermal contributions, are simulated by a continuous noise of random photons introduced at every time step in the laser mode, being responsible of the initial start-up as occurs in real lasers. This is done by making ci(t+1) =ci(t)+1 for a small number of cells (<0.01% of total) with randomly chosen positions. 3. Simulation results The CA model is dependent on three parameters: the pumping probability (), the life time of photons (c) and the life time of excited electrons (a). In each simulation, we provide an initial state (ai(0) =0, ci(0) =0, ∀i, except a small fraction 0.01% of noise photons present) and let the system evolve for a number of time steps. In each step, two macroscopic magnitudes are measured: the total number of laser photons n(t)=Nc i=1ci(t), and the total number of electrons in the upper laser state (population inversion) N(t)=Nc i=1ai(t). It has been shown that it is possible to classify the main kinds of behavior exhibited by the system in the parameters space by calculating the Shannon’s entropy of the distribution of values taken by nor N, after running the simulation for a time interval [10]. In order to calculate the Shannon’s entropy, the range of values taken by nor Nin the time interval is divided into 103equally spaced bins, and the frequency (fi) at which the magnitude value lies inside every particular non-void bin iis computed. The Shannon entropy is then calculated as: S(,c,a)=− m  i=1 filog2fi(1) where mis the number of non-void bins. This magnitude vanishes if the outcome of the system is constant, so that all the bins except one are void, and it increases as the probability distribution of values of the outcome is wider. The resulting Shannon’s entropy is shown in figure 1, where R is the laser pumping rate and Rtis the threshold laser pumping rate, which are linearly related to the pumping probability and the threshold pumping probability tthat appear in the CA model, so that R Rt= t.tis calculated as the smallest value of the pumping probability for which after a transient time the number of laser photons is clearly greater than the number of noise photons introduced. The system exhibits two main characteristic types of behavior: a constant response for values of the parameters for which the Shannon’s entropy Scis low, such as point ain figure 1, and an oscillatory behavior when Scis higher, such as point b. This result must be compared with the predictions of the standard approach based on differential equations. They can be put in their simplest form as a system of two coupled differential equations (laser rate equations) [1, 2]: one of them giving the temporal variation of the number of laser photons n(t) and the other one the temporal variation of the population inversion N(t): dn(t) dt=KN(t)n(t)−n(t) c , (2) dN(t) dt=R−N(t) a −KN(t)n(t). (3) Here Kis a constant called “coupling constant”. From a linearization of these equations for the case of small amplitude fluctuations, two different behaviors are expected, depending on the values of the laser parameters [3]: Damped oscillations Fig. 1. Contour plot of the Shannon’s entropy of the distribution of the number of laser photons for a fixed value of c=10 time steps. Low values of Sc(bright zones) indicate that the response of the system is non-oscillatory, while high values (dark zones) correspond to an oscillatory response. The black line is the theoretical stability curve. Points a,band ccorrespond to the values of the parameters for which the time evolution is shown in Figures 2, 3 and 4. a: R/Rt=16, a/c=3. b: R/Rt=5, a/c=18. c: R/Rt=12.4, a/c=18. A 200 ×200 lattice has been used for this figure. Fig. 2. Evolution of the system for the values of the parameters marked as point ain Figure 1. Up: number of laser photons and population inversion versus time. Down: evolution in a phase space with the number of laser photons versus the population inversion. Parameters: =0.192, c=10, a=30. appear for values of the parameters obeying a c >R Rt2 4R Rt−1(4) and a constant behavior appears if this condition is not satisfied. The black line in figure 1 is the theoretical stability curve from equation (4). This line separates the area of damped oscillations (above and to the right) and the area of a constant behavior (below and to the left) in the parameters space. As seen in figure 1, the dependence on the parameters of the type of behavior exhibited by the system (classified by the Shannon’s entropy), is in a good qualitative agreement with the theoretical stability curve. The time evolution of the system for values of the parameters characteristic of these two regimes is shown in figures 2 and 3. In figure 2, after an initial gain switching peak, the system reaches an stationary constant behavior. In phase space, after an initial transient, the systems goes to a fixed point. In figure 3, the system shows the typical behavior known as laser relaxation oscillations or laser spiking. Correlated large amplitude damped oscillations are shown in the number of laser photons and the population inversion. In phase space, the system follows a spiral trajectory which converges toward a steady-state limit point. As well as the fixed point and the relaxation regular oscillations, the CA model can display another kind of complex behavior. Fig. 4. Evolution of the system (number of laser photons and population inversion versus time) after a transient of 500 time steps, for the values of the parameters marked as cpoint in Figure 1. Fluctuations on a wide range of time scales are exhibited. Parameters: =0.031, c=10, a=180. Figure 4 shows the time evolution of the system for the values of the parameters corresponding to point cin figure 1. After a transient time, in this regime the system exhibits irregular oscillations involving fluctuations on a wide range of time scales. It is important to distinguish these oscillations from those that could be induced by the introduction of noise photons in the system. In order to do that, for this case noise photons have only been introduced in the first 100 time steps to activate the laser process, and the populations values have been recorded after a transient time of 500 time steps. Figure 5 shows the corresponding Fig. 5. Power spectrum of the population inversion (up) and the number of laser photons (down) for the time series depicted in Figure 4, showing a power law behavior. A linear fit of the spectrum (dotted line) varies as f−with =2.15 for the population inversion and =2.30 for the number of laser photons. Fig. 3. Evolution of the system for the values of the parameters marked as point b in Figure 1. Up: number of laser photons and population inversion versus time. Down: evolution in a phase space with the number of laser photons versus the population inversion. Parameters:  = 0.0125, c = 10, a = 180. power spectrum, which follows a power law of the kind 1/f − with exponent close to =2. Since the value of is actually greater than 2.0, the origin of the observed behavior doesn’t seem to be a process of self-organized criticality. Rather, in this regime the system can be in a chaotic state, which would be in agreement with the frequent finding of a chaotic dynamics in real lasers for sufficiently high values of the pumping. A full study of the dynamics of the system in this regime must be carried out in the future. In order to validate the feasibility of the CA model to reproduce the competitive processes between the electrons in the laser levels and the laser photons, the dependence of the laser spiking behavior on some laser parameters can be investigated. The decay rate of the spiking oscillations is found experimentally to depend on the life time of excited electrons (a). Its dependence can be found from a linearization of the laser rate equations [1, 2]: sp =1 2a R Rt . (5) In figure 6 the behavior of the system in the laser spiking regime is represented for different values of akeeping the other two parameters constant. The decay rate of the spiking oscillations increase as adecreases, in good agreement with the laser rate equations predictions (i.e. equation (5)) and the experimental results. Fig. 6. Dependence of the decay rate of the relaxation oscillations on the upper laser level lifetime a. The values of the other two parameters are: =0.0125, c=10. Fig. 7. Dependence of the decay rate of the relaxation oscillations on the number of noise photons () introduced at every time step. The values of the parameters are: =0.0125, c=10, a=180. The effect of the number of noise photons () introduced in the system at every time step on the decay rate of the oscillations is shown in figure 7. When the noise level is increased, the decay rate of the oscillations increases. In our CA model we have not considered explicitly the spontaneous emission process. However, in addition to representing the photons noise level observed experimentally, the noise photons introduced in random positions at every time step can represent too a spontaneous emission process, and the finite life time of excited electrons (a) can represent electron decaying by spontaneous emission. In this way, is linearly related to the spontaneous emission probability, and as shown in figure 7, the decay rate of the spiking oscillations increases as this probability is increased. This behavior is in agreement with the laser rate equations predictions (equation (5)) and with the experimental results [1, 2]. Additionally, figure 7 shows that the amplitude of the oscillations decreases as the number of noise photons increases. This is a logic result, as increasing the noise level intensifies the stochastic character of the system, and therefore the oscillations are smoothed. 4. Conclusions and future prospects The distinct kinds of behavior exhibited by a cellular automata model of laser dynamics (alternative to the standard treatment based on differential equations) have been analyzed. The results have been found to be in good agreement with the predictions of the laser rate equations and with the experimental results. In addition to the two kinds of behavior that had been previously reported (constant response and oscillatory behavior), the system has been found to exhibit a new kind of complex behavior showing 1/f −2noise in its power spectrum, which can be a fingerprint of chaotic dynamics. The CA model is robust under rule modifications which doesn’t change the physical mechanisms of the model. One example is the modification of rules R3 and R4, using a probabilistic (instead of deterministic) decay. Another example is the modification of rule R2, using a threshold value different from 1. In these cases, the qualitative behavior of the system doesn’t change. However, the rule modification can induce some consequence in the dynamics of the system. This happens in the last case, in which increasing the threshold of rule R2 is equivalent to decreasing the stimulated emission cross section and therefore decreasing constant Kin the laser rate equations. As the threshold pumping rate Rtis linearly related to 1/K (see [3]), this should give rise to an increase in Rt, and this is in fact the effect that can be obtained in the CA simulations. On the other hand, the CA model has enough flexibility to easily introduce modifications in the CA rules, involving different physical mechanisms, in order to reproduce specific phenomenology of real lasers. For example, the CA model could be used to model pulsed lasers if rule R1 is modified by using a pulsed pumping rule. After having verified that a very simple CA model can reproduce much of the laser phenomenology, we expect that more sophisticated CA models can be very useful in this area. In particular, we believe that a three dimensional model with more realistic boundary conditions could successfully simulate the characteristics of specific real optoelectronic devices. In addition, they can be used to study problems of current interest, such as cooperative phenomena in lasers, chaotic lasers or two-photon lasers. Acknowledgments We are grateful to the referee for helpful suggestions. References 1. Siegman, A. E., “Lasers”, (University Science Books, 1986). 2. Svelto, O., “Principles of lasers”, (Plenum Press, 1989). 3. Guisado, J. L., Jim´enez-Morales, F. and Guerra, J. M., Phys. Rev. E 67, 066708 (2003). 4. von Neumann, J., “Theory of Self-Reproducing Automata”, (University of Illinois Press, Urbana, 1966). 5. Wolfram, S., “Cellular automata and complexity”, (Addison-Wesley, 1994). 6. Chopard, B. and Droz, M., “Cellular automata modeling of physical systems”, (Cambridge University Press, 1998). 7. Toffoli, T. and Margolus, N., “Cellular automata machines: a new environment for modelling”, (The MIT Press, 1987). 8. Sloot, P. M. A., Kaandorp, J. A., Hoekstra, A. G. and Overeinder, B. J., Distributed simulation with cellular automata: architecture and applications. In (J. Pavelka, G. Tel, and M. Bartoˇsek, editors), SOFSEM’99: Theory and Practice of Informatics, volume 1725 of Lecture Notes on Computer Science, pages 203–248, (1999). 9. Talia, D., IEEE Computer 33, 44 September (2000). 10. Guisado, J. L., Jim´enez-Morales, F. and Guerra, J. M., Application of shannon’s entropy to classify emergent behaviors in a simulation of laser dynamics. In (T. E. Simos, editor), Computational Methods in Sciences and Engineering 2003, pages 213–216, Singapore, (2003). World Scientific.