text:
# Velocity distributions in dilute granular systems

## I Introduction

Granular materials consisting of macroscopic particles or grains can exhibit behavior reminiscent of conventional phases of matter. Sand, for instance, can flow like a liquid under some conditions. Dilute granular systems, or gases, have been extensively studied both experimentally and theoretically, in large part as simple model systems exhibiting nonequilibrium and dissipative behavior. These systems are intrinsically dissipative and out of equilibrium, even though it is tempting to apply such equilibrium notions as temperature. Since the collisions in such a gas are inelastic, it is necessary to supply energy, or to drive them in order to maintain a gas-like steady state. Otherwise, inelastic collapse can occur, in which all motion ceases after only a finite timemcnamara ; kadanoff1 . In principle, it is possible to heat or drive the system uniformly throughout the container (uniform heating), as has been done in simulationsmackintosh ; moon , and as is assumed in analytic theoriesernst . In experiments, however, one usually drives a granular gas by shaking or vibrating the walls of the container. Such boundary heating means that the energy is inserted in a spatially inhomogeneous wayluding ; losert ; kadanoff2 ; rouyer ; kudrolli . As a consequence, the gas will develop a gradient in density and mean kinetic energy.

Even for uniform heating, however, significant deviations from equilibrium gases, e.g., in density correlations, are observedmackintosh . It is often assumed that in the bulk, where the boundaries are far away and spatial gradients are small on the scale of the mean free path of the particles, both heating methods should give the same behavior. Here, we show that there are striking differences in the behavior of granular gases heated uniformly or through the boundary.

There has been much interest recently in another aspect of these non-equilibrium gases: namely, the velocity distributions. This is in part because they deviate from the Gaussian distributions that one would expect if the collisions were elastic. It has been suggested that the velocity distributions can be described by a sort-of stretched Gaussian, of the form $`P(v)=C\mathrm{exp}[\beta (v/\sigma )^\alpha ]`$, where $`\sigma =v^2^{\frac{1}{2}}`$ is sometimes called the granular temperature. Rouyer and Menon rouyer  have reported such a distribution with exponent $`\alpha =1.5`$ over the whole observed range of velocities, which was unaffected by changes in amplitude and frequency of driving. This observation was particularly intriguing, given the theoretical predictions obtained before by Van Noije and Ernst ernst . They developed a kinetic theory that predicts a high-velocity tail described by a distribution with an exponent $`\alpha =\frac{3}{2}`$.

Unfortunately, results for velocity distributions are never unambiguous. Both in simulation and experiment, different setups and driving mechanisms usually give different behavior of the velocity distribution. For a setup where particles on a horizontal plate were driven in the vertical direction, Olafsen and Urbach olafsen  found a crossover from exponential to Gaussian distributions as the amplitude of the driving was increased. The result of Rouyer and Menon rouyer  was obtained for a different configuration where particles were confined between two vertical plates and driven in the vertical direction. Although the exponent $`\alpha =1.5`$ they find is reminiscent of theoretical results obtained by Van Noije and Ernst ernst , both results are actually inconsistent in that the exponent $`\alpha =\frac{3}{2}`$ can only describe the high-velocity tail of the distribution. Such a high-velocity tail was actually obtained by Moon et al. moon . Using a simulation with a variation on uniform heating they find a velocity distribution that has an exponent of $`\alpha =2.0`$ for low velocities, but crosses over to an exponent of $`\alpha =1.5`$ for the high-velocity tail. The results of Rouyer and Menon are also limited in the sense that they were obtained for almost elastic particles with $`\eta =0.93`$. Blair and Kudrolli kudrolli  use a different setup where particles move along an inclined plane. Friction with the plane during collisions reduces the coefficient of restitution to $`\eta 0.5`$. They find the distribution with exponent $`\alpha =1.5`$ only in the very dilute case. Otherwise, the distributions deviate strongly from both Gaussian and the distribution obtained by Rouyer and Menon. It is interesting to notice that for the denser case Blair and Kudrolli find a distribution with a crossover much like the type observed by Moon et al.

At present there is no agreement on what the velocity distribution of a granular gas looks like exactly nor is it clear what causes the velocity distributions to deviate from a Gaussian. Puglisi et al. puglisi  have suggested that the deviations are caused by the spatial correlations in the gas. They propose that for a region of uniform density the velocity distribution of the gas actually is Gaussian, with a density dependent width. The spatial correlations cause density fluctuations and they claim that the non-Gaussian distributions arise as an average over the velocity distributions over these regions of different density. It is true that an average over Gaussian distributions with different widths in general yields a stretched Gaussian, but experiments performed by Olafsen and Urbach olafsen  showed that this does not happen in granular gases. They find that even within small windows of uniform density the velocity distributions remained non-Gaussian. Here, we show that non-Gaussian behavior arises even in a simple model with no spatial degrees of freedom.

We study behavior of the velocity distributions of the granular gas as a function of $`\varphi `$, the area fraction, and $`\eta `$, the coefficient of restitution. Specifically, we consider the effect on the velocity distribution of driving the gas by heating uniformly, as is assumed in theory and many prior simulations, and by heating through a boundary, as is done in most experiments. We will show that there exists clear qualitative difference between the velocity distributions for uniform and boundary heating. Furthermore, we will show that there is no evidence for a universal velocity distribution with a constant exponent $`\alpha =1.5`$. We demonstrate instead that a family of distributions showing apparent exponents covering a wide range of values $`\alpha <2`$ is expected, depending on both material and experimental conditions. Furthermore, we show that the velocity distribution is governed primarily by the relative importance of collisions to heating, *i.e.*, the way in which energy flows through the system of particles. Specifically, we introduce a new parameter $`q=N_H/N_C`$, which measures the ratio between mean numbers of heating events $`N_H`$ and mean number of collisions $`N_C`$ experienced by a typical particle. These theoretical observations can explain both the observed non-Gaussian behavior, as well as the ambiguities in the experimental and theoretical literature on dissipative gases to date. We also demonstrate that the behavior of the velocity distributions can be captured quantitatively by a simple model that takes only $`\eta `$ and $`q`$ into account, with no spatial degrees of freedom.

## II Numerical Simulation

We use an event-driven algorithm to simulate $`N`$ particles of radius $`r`$ moving in a two-dimensional box. Particles gain energy by heating and lose energy through inelastic collisions. When two particles $`i`$ and $`j`$ collide their final velocities depend on their initial velocities in the following way:

$$𝐯_i^{}=𝐯_i\frac{1+\eta }{2}(𝐯_i\widehat{𝐫}_{ij}𝐯_j\widehat{𝐫}_{ij})\widehat{𝐫}_{ij},$$
(1)
where $`0\eta <1`$ is the coefficient of restitution and $`\widehat{𝐫}_{ij}`$ is the unit vector connecting the centers of particles $`i`$ and $`j`$.

For uniform heating we adapted an one-dimensional algorithm described in mackintosh . When heating uniformly, each individual particle is heated by adding a random amount to the velocity of each particle during a time step $`\mathrm{\Delta }t`$:

$$𝐯_i(t+\mathrm{\Delta }t)=𝐯_i+\sqrt{h\mathrm{\Delta }t}𝐟(t),$$
(2)
where $`𝐟(t)`$ is a vector whose components are uniformly distributed between $`\frac{1}{2}`$ and $`\frac{1}{2}`$ and $`h`$ is proportional to the heating rate. After heating the system is transferred to the center-of-mass frame. Particles move in a box with sides $`L=2`$. We use periodic boundary conditions to simulate bulk behavior. The time step $`\mathrm{\Delta }t`$ is chosen in such a way that on average the number of collisions per time step is less than one. It should be noted that this heating mechanism is significantly different from the spatially homogeneous heating used in some experiments olafsen . In the experiment all particles feel the same forcing, so the motion of the neighboring particles is strongly correlated in space and time. In uniform heating however, all particles are independently driven by a stochastic source and as a consequence correlations are very weak.

When heating through the boundaries, particles gain velocity upon collision with the boundary. For simplicity, we assume that the collision with the boundary is elastic. In that case, a collision occurs by reflecting $`𝐯_{}`$, the component of the velocity perpendicular to the boundary. Heating occurs by adding a random amount of velocity to $`𝐯_{}`$. Then after collision with the boundary one has:

$$𝐯_i^{}=𝐯2𝐯_{}+\sqrt{h}𝐟(t).$$
(3)
Particles move in a circular box of radius $`R=1`$. A symmetrical container has the advantage that it allows us to examine density and granular temperature gradients along a single coordinate $`r`$, the distance from the center of the box, as in the one-dimensional case kadanoff2 .

We start the simulation by distributing the particles uniformly over the box. When using boundary heating, we give each particle a small, uniformly distributed velocity to enable particles to reach the boundary. Then particles are heated and we allow the system to reach a steady state before taking data. For both uniform heating and boundary heating, data is taken periodically every time step $`\mathrm{\Delta }t`$. For uniform heating, data is taken when the particles are heated, so $`\mathrm{\Delta }t`$ equals the time between heating events.

## III Simulation results: clustering

One of the first striking differences between uniform heating and boundary heating is clustering. When heating through the boundary, a stable liquid-like cluster surrounded by a hot gaseous state will form for low coefficients of restitution $`\eta `$ or high area fraction $`\varphi `$. This occurs as particles are compressed in the center of the box by the pressure of particles moving in from the boundary. As the cluster grows in size, it can no longer be destroyed by the impact of high velocity particles and the cluster remains stable. A typical example is shown in Fig. 1. We do not find these clusters when heating uniformly. This is because all particles are heated all the time, which prevents the collapse to a cluster. In experiments that find clusters in homogeneously driven systems olafsen  this is different, because here the driven motion of neighboring particles is highly correlated. The formation of a cluster spoils measurement of the velocity distribution because the system is no longer in a pure gas-like state. To avoid values of $`\varphi `$ and $`\eta `$ corresponding to the formation of clusters in our simulation, we constructed a phase diagram. We did this by counting for every particle the number $`N_{6r}`$ of neighbors with their center within a distance $`6r`$ from that particle. When the gas is in a hexagonal close packed state $`N_{6r}=32`$. We obtained the distribution $`P(N_{6r})`$ for different values of $`\varphi `$ and $`\eta `$. An example for $`N=350`$ and $`\varphi =0.1`$ is shown in figure 2.

For $`\eta =0.9`$ the distribution corresponds to a state with the particles uniformly distributed over the box and the peak of the distribution at the mean value $`\overline{N}_{6r}=3.6`$. For $`\eta =0.7`$ the distribution becomes bimodal, with a broad peak at high $`N_{6r}`$ corresponding to the densely-packed cluster and a peak at $`N_{6r}=1`$ corresponding to the surrounding dilute gas. The distribution shows a continuous transition for $`\eta `$ in between, which makes it hard to pinpoint an exact value of $`\eta `$ for which the gas enters the clustered state. Still, by looking at the shape of the distributions, it can be argued that the transition occurs somewhere between $`\eta =0.75`$ and $`\eta =0.85`$. This was repeated for different values of $`\varphi `$, which allowed us to determine a sort of phase diagram. Specifically, we determined the limit of a pure gas-like phase, and all results presented below were obtained in this state.

## IV Simulation results: velocity distributions

The velocity distributions $`P(v_x)`$ for uniform heating are shown in Figs. 3 and 4. The velocity component $`v_x`$ is scaled by $`\sigma _x=v_x^2^{\frac{1}{2}}`$ and the maximum of the distribution $`P(v_x/\sigma _x)`$ is scaled to be unity. For a broad range of the parameters $`\varphi `$ and $`\eta `$ the velocity distributions are very close to Gaussian. For $`\eta =0.8`$ the velocity distributions can be fitted by a distribution with $`\alpha =2.0`$. This decreases only slightly for $`\eta =0.1`$, which can be fitted by a distribution with $`\alpha =1.9`$. Values of $`\alpha `$ are found to be independent of $`\varphi `$. These exponents are constant over the entire observed range of velocities and we see no high-velocity tails with $`\alpha =1.5`$ for the range of $`\varphi `$ and $`\eta `$ we examined.

For boundary heating the gas develops a gradient in both density and mean kinetic energy as shown in figures 5 and 6. Ideally, we want to measure velocity distributions in a region where the gradient is small. To this end we divided the box in five rings of width $`0.2`$. These rings are indicated in figures 5 and 6. Only for values of $`\varphi `$ and $`\eta `$ close to the clustering state, does the density within a ring vary by more than $`10\%`$. The velocity distributions $`P(v_x)`$ for particles within the different rings are shown in figure 7. Figure 7(a) shows $`P(v_x)`$ with the velocity component $`v_x`$ scaled by $`\sigma _x=v_x^2^{\frac{1}{2}}`$ and the maximum of the distribution scaled to be unity. All distributions seem to have the same shape. This is remarkable, since naively one would expect the dynamics of the particles to be different for the different rings. In the center the density is high and particles have many collisions, causing strong spatial correlations. Along the boundary the gas is dilute. Here, many particles have experienced recent collisions with the boundary, and spatial correlations should be much weaker. It is has been suggested that these correlations cause the velocity distributions to deviate from Gaussian. But contrary to the expectations, within the container the shape of the distributions seems to be determined by the granular temperature $`\sigma _x`$ only. This is a feature that is observed for all values of $`\varphi `$ and $`\eta `$, even close to the cluster state.

Figure 7(b) shows the behavior of the exponent $`\alpha `$. This behavior is very different from the case of uniform heating. For uniform heating $`\alpha `$ has the same value over the entire observed range of velocities. For boundary heating, on the other hand, $`\alpha `$ has a constant value $`\alpha _1`$ over the low-velocity range but crosses over to different value $`\alpha _2`$ when above a critical velocity $`v_c`$. For the two innermost rings $`\alpha _1=1.7`$ and $`\alpha _2=1.0`$. For the third ring $`\alpha _1=1.8`$ and the outer ring has $`\alpha =1.5`$. For the outer two rings the velocity distributions actually undergo two crossovers. First it crosses over to $`\alpha _21.0`$ and then the exponent increases again to $`\alpha _31.5`$. In figure 8 we show the effect of a change in $`\varphi `$ and $`\eta `$ on the shape of the velocity distributions. Here we focus on the velocity distribution as measured in the ring with $`0.4<r0.6`$. This has the advantage of good statistics, but for values of $`\varphi `$ and $`\eta `$ close to a cluster, we might see effects due to the density gradient in the gas. As shown in figure 8(a) the exponent $`\alpha _1=1.8`$ except for $`\eta =0.4`$, where $`\alpha _1=1.6`$. For $`\eta =0.9`$ there is no crossover in the observed range of velocities. In the other distributions one does observe a crossover and the point where it occurs shifts down to lower velocities as $`\eta `$ is decreased. It is clear that the distribution for velocities above the crossover cannot be described by a single exponent. For low enough $`\eta `$, the distribution seems to approach a constant exponent for high velocities. This exponent decreases from $`\alpha _2=1.3`$ for $`\eta =0.7`$ to $`\alpha _2=1.0`$ for $`\eta =0.4`$. In figure 8(b) $`\alpha _1=1.8`$ for $`\varphi =0.01`$ and $`\varphi =0.02`$, $`\alpha _1=1.9`$ for $`\varphi =0.03`$ and $`\alpha _1=1.7`$ for $`\varphi =0.05`$. For every $`\varphi `$ does one observe a crossover and the velocity at which the crossover occurs shifts only a bit as $`\varphi `$ is varied. The distributions approach a constant exponent for high velocities. This exponent goes down from $`\alpha _2=1.5`$ for $`\varphi =0.01`$ to $`\alpha _2=1.0`$ for $`\varphi =0.05`$. In general, the deviations from Gaussian become more pronounced as dissipation increases, i.e. as $`\varphi `$ increases or as $`\eta `$ decreases.

To test whether the velocity distributions we find here are only observed for this specific driving mechanism of heating through a circular boundary, we constructed different systems that drive through boundaries in a different way. For instance, we constructed a box with periodic boundary conditions that includes a small circular region around the center. Within this circular region particles are uniformly driven but outside of the region they are not heated at all. For particles within the circular region we observe velocity distributions that are Gaussian. On the other hand, for particles outside of the circular region we observe the same non-Gaussian velocity distributions as seen in the case of a circular boundary.

For uniform heating the distribution of velocities of particles that are heated is Gaussian, because the velocities undergo a random walk. This is not the case for the particles heated at the boundary. To see if velocity distributions are sensitive to the precise way of heating at the boundary, we changed our heating algorithm so that, when a particle hits a boundary, it’s new velocity is drawn from a Gaussian distribution. This has a modest effect on the far end of the high-velocity tail, but leaves all major differences between uniform and boundary heating intact.

Finally, we studied the behavior of the velocity distribution for different particle numbers $`N`$. Figure 9(a) shows that, as $`N`$ increases, the velocity distributions fall off more rapidly, but the high-velocity tails grow longer. This is shown more clearly in figure 9(b). For the high-velocity tail we find $`\alpha =1.7`$ for $`N=100`$ decreasing to $`\alpha =0.7`$ for $`N=1000`$. The crossover shifts to lower velocity and becomes sharper as $`N`$ is increased. As we will demonstrate in more detail below, for boundary heating there is no thermodynamic limit. Instead of approaching a limiting velocity distribution as $`N`$ is increased, we find that the shape of the velocity distribution depends not only on $`\eta `$ and $`\varphi `$, but also on $`N`$.

For their experiments Rouyer and Menon used $`N`$ particles with $`\eta 0.9`$, where $`100<N<500`$ and $`0.05<\varphi <0.25`$ rouyer . In figure 10 we plotted the velocity distribution for $`\eta =0.9`$, $`\varphi =0.05`$ and several values of $`N`$. The solid black lines indicated the fit with $`\alpha =1.52`$ as made in rouyer . This line clearly coincides with the high-velocity tail of the velocity distribution found by simulation. This suggests that it is possible that instead of a universal distribution with $`\alpha =1.5`$, they only observed a part of a more complex velocity distribution, with two exponents and a crossover. For higher $`\varphi `$ the high-velocity tails have exponents that are smaller than the $`\alpha =1.5`$ observed for $`\varphi =0.05`$.

The main difference between uniform and boundary heating is that in the first case heating takes place homogeneously throughout the box, whereas in the latter case energy is injected inhomogeneously at the boundaries. This is not the direct cause for the difference in velocity distributions. This is demonstrated by the occurrence of non-Gaussian velocity distributions in a system studied by Moon et. al. moon . They study a uniformly driven gas, but nevertheless find velocity distributions with a clear crossover. The main difference is the way they implement uniform heating. Every time $`\mathrm{\Delta }t`$, they select at random two particles and add to these particles a random but opposite velocity to conserve the total momentum. On average heating is spatially homogeneous and in the limit of very small $`\mathrm{\Delta }t`$ the behavior of both their and our uniform heating algorithm is the same. For small enough $`\mathrm{\Delta }t`$ the number of collisions between heating events is smaller than one and the heating dominates over the dissipative behavior of the gas. When $`\mathrm{\Delta }t`$ becomes larger, many collisions might happen between two heating events. Apperently this changes the velocity distributions. In figure 11 we show the velocity distribution for $`N=350`$, $`\varphi =0.02`$ and $`\eta =0.4`$. The gas is heated using the two-point heating algorithm described above, varying the time between heatings, $`\mathrm{\Delta }t`$. For $`\mathrm{\Delta }t=0.01`$ the distribution has a exponent $`\alpha =1.7`$ that is constant approximately over the observed range. When $`\mathrm{\Delta }t`$ is reduced a clear crossover develops. The behavior of the velocity distribution for velocities higher than the crossover velocity is more complicated than in boundary heating. There is also a sharp kink at the high-velocity end, that we have been unable to explain so far.

This can explain the differences between velocity distributions for uniform and boundary heating. When heating uniformly all the particles are heated every timestep $`\mathrm{\Delta }t`$, while only a few collisions occur in that time. The system is dominated by heating and we see a velocity distribution with a constant exponent over the entire observed velocity range. If we are heating through a boundary, there are many collisions in the dense region in the center of the box, but only a few particles escape to the boundary to get heated. In this case, the system is dominated by dissipation and we see in general a velocity distribution that strongly deviates from a Gaussian. This difference between uniform and boundary heating can be quantified by the ratio $`q=N_H/N_C`$, where $`N_H`$ and $`N_C`$ are the mean number of heatings and the number of collisions a particle has experienced, respectively. For a uniformly heated gas of $`N=350`$, $`\varphi =0.02`$ and $`\eta =0.4`$ we find that $`q=120`$. For the same parameters but heating through the boundary we find $`q=0.08`$. For again the same parameters but using two-point heating with $`\mathrm{\Delta }t=1.0`$ we find $`q=0.012`$.

We tested the effect of changing $`q`$ on velocity distributions obtained with boundary heating. When increasing the number of particles $`N`$ or the area fraction $`\varphi `$, the average number of collisions increases. One can show in a mean field approximation that $`q(N\varphi )^{1/2}`$. The average distance a particle travels between collisions is given by $`l_{coll}1/\varphi `$. The average distance a particles travels between heatings $`l_{heat}`$ is given by the dimensions of the container. For a box of area $`A`$ the average distance between boundaries is given by $`l_{heat}A^{1/2}=(N/\varphi )^{1/2}`$. Finally, we know $`N_H/N_Cl_{coll}/l_{heat}`$. Our simulation obeys this approximation very well. In Fig. 12 we show velocity distributions for $`\eta =0.8`$ and different combinations of $`N`$ and $`\varphi `$. We measure the heating-dissipation ratio $`q`$ in the simulation and show velocity distributions with the same $`q`$ on top of each other. For $`q=1.3`$ and $`q=0.13`$ we find excellent collapse for different $`N`$ and $`\varphi `$, even when we scale the system by a factor $`8`$. For $`q=0.013`$, where spatial correlations become very strong, we still find reasonable collapse. As we increase $`q`$ we observe the usual pattern, where a crossover appears in a distribution that was initially close to a Gaussian. This means that the velocity distributions are almost exclusively controlled by two parameters; the coefficient of restitution $`\eta `$, which is a material parameters, and the heating-dissipation ratio $`q`$, which depends on experimental conditions.

## V a simple model without spatial degrees of freedom

To study the dependence of the velocity distributions on the ratio $`q`$ we constructed a model of the inelastic gas inspired by work done by Ulam ulam  on elastic gases. Ulam rederived the Maxwell-Bolzmann distribution for the velocities of particles in a perfect gas using the following very simple model: he considered $`N`$ particles with random initial velocities. At every time step a pair of particles was selected and the velocities of the particles were changed as if they had collided elastically. Ulam found that, regardless of the initial velocities, the system quickly evolved to a state where on average all energy was distributed equally over all particles and the deviation of the velocities of the particles around the average velocity of the system was given by the Maxwell-Bolzmann distribution. This approximation of random collisions is only justified when the gas is sufficiently dilute. The same approximation can also be made for the granular gas by replacing the elastic collisions between particles with inelastic collisions. Since the gas cools down without any further energy input, heating has to be included as well. This allows us to study the dissipative behavior of the granular gas while neglecting all spatial correlations that develop in the more detailed simulations described before. In particularly we can study the effect on the velocity distributions of changing the fraction $`q=N_H/N_C`$ described in the previous section by varying the ratio of heatings and collisions.

We adapted Ulam’s procedure in order to model a two-dimensional granular gas consisting of $`N`$ particles. Equation 1 can be cast into the following form:

$$𝐯_i^{}=𝐯_i\frac{1+\eta }{2}\left(\begin{array}{cc}\mathrm{cos}^2\theta & \mathrm{sin}\theta \mathrm{cos}\theta \\ \mathrm{sin}\theta \mathrm{cos}\theta & \mathrm{sin}^2\theta \end{array}\right)(𝐯_i𝐯_j),$$
(4)
where $`𝐯_i`$ and $`𝐯_i`$ are the velocities of particles $`i`$ and $`j`$, $`\eta `$ is the coefficient of restitution and $`\theta `$ is the angle between the separation vector $`𝐫_{ij}`$ and a reference axis. Collisions in our model occur by selecting at random particles $`i`$ and $`j`$ and an uniformly distributed impact parameter $`2R<b<2R`$. We then use the above collision rule with $`\theta =\mathrm{arcsin}(b/2R)+\mathrm{arccos}(𝐯\widehat{𝐬}/v)`$, where $`𝐯=(𝐯_j𝐯_i)/2`$ is the velocity in the center-of-mass frame and $`\widehat{𝐬}`$ is a unit vector along the reference axis. We discard values of $`\theta `$ corresponding to $`(𝐯_j𝐯_i)𝐫_{ij}<0`$ as these represent unphysical collisions. We heat the gas by selecting a random particle $`i`$ and adding a random amount of velocity according to equation 2. To prevent the velocities from running away, we subtract the center-of-mass velocity after heating. In a single time step, we have let $`2C`$ particles collide and we heat $`H`$ particles. This gives us $`q=H/(2C)`$.

When we use inelastic instead of elastic collisions and drive the system, we find that the gas heats up until it reaches a steady state, where on average the energy lost in collisions is compensated by the heat inserted into the system through our driving mechanism. As $`H`$ is chosen increasingly bigger in comparison to $`C`$, the granular temperature $`\sigma `$ of the gas increases accordingly.

In figure 13 we plotted the result for an inelastic gas with $`\eta =0.4`$. We varied the number of heatings and the number of collisions in a single time step from $`H=100`$ and $`C=1`$ to $`H=1`$ and $`C=1000`$. As $`q`$ is lowered, the velocity distributions develop a crossover and for $`q1`$ the distributions are strongly non-Gaussian, similar to the velocity distributions obtained for $`\mathrm{\Delta }t1`$ in two-point heating. In figure 14 we keep $`q=0.025`$ fixed and vary $`\eta `$. We see that the crossover point shifts down as $`\eta `$ is lowered and that the kink in the velocity distribution shifts down.

We also compared velocity distributions found in the model with those acquired by simulation. To this end we measured $`q`$ in simulations using different ways of heating. For a uniformly heated gas of $`N=350`$, $`\varphi =0.02`$ and $`\eta =0.4`$ we found $`q=120`$. For the same parameters but heating through the boundary we found $`q=0.08`$. For two-point heating with $`\mathrm{\Delta }t=1.0`$ we found $`q=0.012`$. We ran the model with the corresponding set of parameters and measured velocity distributions. The results are shown in figure 15.

The velocity distributions agree qualitatively, but especially for the case of q=0.08 there is considerable quantitative difference. Still, the model is able to generate non-Gaussian velocity distributions, while neglecting all spatial correlations. This suggests that the non-Gaussian behavior of the velocity distributions is not caused by spatial fluctuations and correlations in the gas. Instead it is the flow of energy through the system, mediated by the inelastic collisions, that determines the shape of the velocity distribution. This also means that there is no a priori difference between inhomogeneous and homogeneous heating: the difference in the shape of the velocity distribution, particularly the occurrence of the crossover, is only a consequence of the fact that in inhomogeneous heating in general the number of collisions between particles exceeds the number of particles being heated.

## VI relation to experiments

Due to the idealized nature of our system, it is not possible to do a direct comparison between our simulation and experiments. In the experiment of Rouyer and Menon the driving is through the boundary, but there are some significant differences between their heating mechanism and the one we use in simulations with boundary heating. Due to gravity and the geometry of the setup the injection of energy in the experiment of Rouyer and Menon is mainly in the vertical direction. This energy is transferred into the horizontal direction by collisions between particles. Another difference is that in the experiment the frequency of driving is relatively low. Because of this the dynamics of the gas close to the driving boundary is strongly dependend on the phase of the driving cycle. In fact, it has been shown in simulation that for a system similar to the experiment by Rouyer and Menon, a shockwave propagates up through the gas bougie . At a certain distance from the boundary the time dependence has decayed away and the gas enters a steady state. It is in this steady state that the velocity distributions are measured.

It is not yet established how this time dependence and the occurence of a shock wave influences the velocity distributions in the steady state. A priori it is not clear if it is possible to compare velocity distributions in systems with a strong time-dependence, like the experiments, with those that have no time dependence, as is the case in our simulations. There are, however, reasons to assume this is possible. The velocity distributions in the experiment of Rouyer and Menon are measured only in the direction orthogonal to the driving direction. Simulations bougie  show that the effect of the shock in the direction orthogonal to the shock are usually relatively weak and decay rapidly in height. In the steady state, influence of the shock is absent in the orthogonal direction, even while it still may be apparent in the direction perpendicular to the driving direction. So, if we only look at velocity distributions in the orthogonal direction and in the steady state, a comparison between the experiment and our simulations is be valid.

It is also not clear how in the experiments the dynamics of the gas shape the velocity distribution and whether it is controlled by the parameter $`q`$. We speculate that in the steady state the system behaves in fact like a one dimensional inelastic gas. Fast upward moving particles inject energy in the orthogonal direction when colliding with particles in the steady state, effectively functioning as a heat source. In this picture the mean number of collisions between fast upwards moving particles and particles in the steady state would be $`N_H`$, the average number of heatings, and collisions between the particles in the steady state mutually would be $`N_C`$, the average number of collisions. One way of changing the shape of the velocity distribution would be changing the fraction of particles in the steady state. More particles in the steady state would lower $`N_H`$ and increase $`N_C`$ leading to more non-Gaussian velocity distributions.

The above considerations not only apply to the experiment of Rouyer and Menon but to most of the other experiments as well. In the experiments of Blair and Kudrolli and those of Olafsen and Urbach velocity distributions are measured orthogonal to the driving direction. In both cases it is not so much the collisions with the bottom plate that drive the gas in the ortogonal directions, but mainly collisions between fast upward moving particles with particles that have low velocities in the orthogonal directions. It is in these experiments rather than those of Rouyer and Menon that we find a similar dependence on $`\eta `$ and $`q`$ as we describe in this article.

In the setup of Olafsen and Urbach olafsen  velocity distributions go from non-Gaussian to Gaussian when a rough plate is used instead of a flat plate. On a flat plate, energy is only injected in the in-plane directions by off-angle collisions between neighboring particles. With a rough plate, energy is injected directly into the directions parallel to the plate every time a particle collides with the plate, effectively increasing the number of heatings over collisions. Baxter and Olafsen baxter  observe the same behavior in a system where a layer of heavy particles is inserted between the other layer of particles and a flat bottom plate. Particles from the upper layer have off-angle collisions with the layer of heavy particles, injecting energy in the in-plane directions every cycle. Particles in the upper layer show Gaussian velocity distributions, whereas particles in the lower layer have non-Gaussian velocity distributions.

Most convincing is the experiment by Blair and Kudrollli kudrolli . Here the number of collisions is increased by adding more particles. As a result, their velocity distributions develop the same crossover we see both in our simulations and model. The reason why these transitions are not visible in the experiment by Rouyer and Menon as they increase number of particles, is that in the first case the effective coefficient of restitution is much lower, $`\eta 0.5`$, due to friction with the inclined plane.

Again, one of the main problems considering the velocity distributions in granular gases is that different setups and experiments usually find different velocity distributions. As we have shown in this section, our finding of the controlling parameter $`q`$ could ultimately explain these seemingly inconsistent results.

## VII conclusion

We compared the velocity distributions of a granular gas that was driven by uniform heating and by heating through the boundary. Although it is usually assumed that uniform heating yields the same behavior as boundary heating when the boundary is far away and the spatial gradients are small, we find that there are clear qualitative differences. When driven through the boundary, for instance, the gas can form coexisting cool liquid-like clusters surrounded by a hot gaseous state for certain values of $`\varphi `$ and $`\eta `$. Such clusters do not occur in our simulations with uniform heating.

The difference between uniform heating and boundary heating also extends to the velocity distributions. For uniform heating, we find velocity distributions that are close to Gaussian over the observed range of velocities and for a wide range of $`\varphi `$ and $`\eta `$. When heating through the boundaries, on the other hand, we find that velocity distributions are stretched Gaussians that usually cross over from one exponent $`\alpha `$ to another for the high-velocity tail. For the high-velocity tail we find that the exponent $`\alpha `$ varies over a wide range. For strong dissipation this tail cannot be described by a single exponent.

This means that there is no universal velocity distribution with $`\alpha =1.5`$, as was proposed by Rouyer and Menon. Instead, the velocity distribution found by Rouyer and Menon may be just one out of the many distributions described here. The apparent universality may be due to the use of a narrow range of parameters (nearly elastic particles were used with $`\eta =0.93`$ and vary $`0.05\varphi 0.25`$). We find that varying $`\varphi `$ only has a minor effect on the shape of the distribution when $`\eta `$ is so close to the elastic limit. For low $`\eta `$, we instead find velocity distributions that look much like the distributions found by Blair and Kudrolli for gases with $`\eta 0.5`$.

Furthermore, we demonstrate that the distribution of velocities for dissipative gases, while not universal in form, depends only on two parameters: the coefficient of restitution $`\eta `$ (a material parameter) and $`q=N_H/N_C`$, the average ratio of heatings and collisions in the gas (a function of experimental conditions). We find that velocity distributions range from Gaussian for $`q1`$, where heating dominates dissipation, to strongly non-Gaussian for $`q1`$, where the dynamics of the gas is dominated by the dissipative collisions between particles. The differences in distributions we find between uniform and boundary heating can then be understood as different limits of $`q`$, for $`q1`$ and $`q1`$ respectively. We can control the parameter $`q`$ more easily using two-point heating and here we observe the transistion form Gaussian to strongly non-Gaussian behavior directly as we change $`q`$.

Furthermore, a simple model of a driven, inelastic gas without spatial degrees of freedom reproduces the entire family of velocity distributions we find in simulation, as we vary $`\eta `$ and $`q`$. This means that the velocity distributions are non-Gaussian not because of spatial correlations. Rather, it is the cascade of energy from a few high-energy particles to the slow-moving bulk of the gas that is the key determinant of the non-Gaussian velocity distributions. These observations should aid in the construction of a kinetic theory of dissipative gases.

We thank L. P. Kadanoff, N. Menon, A. Kudrolli, D. Blair and J. Bougie for useful conversations. This work is supported in part by the National Science Foundation under Grant PHY99-07949.


format: text

text:
# Manipulation of the Spin Memory of Electrons in n-GaAs

## Abstract

We report on the optical manipulation of the electron spin relaxation time in a GaAs based heterostructure. Experimental and theoretical study shows that the average electron spin relaxes through hyperfine interaction with the lattice nuclei, and that the rate can be controlled by the electron-electron interactions. This time has been changed from 300 ns down to 5 ns by variation of the laser frequency. This modification originates in the optically induced depletion of n-GaAs layer.

A long spin relaxation time for electrons, $`\tau _s`$,  could hold promise for magnetoelectronics and quantum information applications . This time can be controlled by an external magnetic field  , by choice of heterostructure geometry , by doping level , and by a gate voltage that changes the electron concentration in the two-dimensional channel .

This letter is devoted to the optical manipulation of the spin memory of electrons. We show, using optical orientation methods, that the spin relaxation rate in a GaAs-based heterostructure can be modified by a factor of 100 by variation of the laser frequency. The comparison between experiment and theory shows that the observed effect is the result of turning on or off of the hyperfine channel of spin relaxation of electrons localized on shallow donors. The giant changes of spin relaxation time allow one to extend the concept of spin dynamics engineering in quantum wells  to quasibulk semiconductors having much longer spin lifetimes.

The optical orientation experiments were carried out on a sample with a semi-insulating GaAs substrate, on which a 500 nm thick GaAs buffer and a 25 nm thick AlAs barrier separating the substrate from the main structure were MBE-grown. Following this is a 100 nm thick GaAs layer capped by a 25 nm thick Al0.3Ga0.7As barrier, above which a series of five quantum wells of different thickness were prepared (Fig.1b). The sample, though nominally undoped, has a background doping level of order $`10^{14}`$ cm-3 n-type from Hall measurements on test structures. The dominant photoluminescence (PL) at the wavelength close to the GaAs band edge is due to the 100 nm layer. The sample was placed in a liquid-helium cryostat and pumped by a tunable Ti-sapphire laser, with the circular polarization of light being alternated in sign at a frequency of 26.61 kHz with a photoelastic quartz modulator. The geomagnetic field was compensated to less than 0.1 G at the sample. The PL polarization was measured in the reflection geometry with a fixed circular-polarization analyzer by a double-grating spectrometer. The electronics provided measurement of the effective degree of circular polarization $`\rho =\left(I_+^+I_+^{}\right)/\left(I_+^++I_+^{}\right)`$, where $`I_+^+,I_+^{}`$ are the intensities of the $`\sigma ^+`$ PL component under $`\sigma ^+`$ and $`\sigma ^{}`$ pumping, respectively.

When the excitation energy, $`h\nu `$ , is below the absorption edge of the quantum wells, the PL spectrum of the 100 nm GaAs layer (curve $`1`$ in Fig.1a at $`h\nu =1.520`$ eV) is dominated by two lines, corresponding to the recombination of the free exciton (X) and of the exciton bound to a neutral donor (D0X) . The GaAs PL spectrum changes drastically when the laser frequency is in resonance with the heavy-hole exciton transition (HHX) of the quantum well (14 nm) nearest to the 100 nm layer (curve 2 at $`h\nu =1.534`$ eV, Fig.1a). The intensity of the exciton line becomes 3 times larger whereas the intensity of the D0X complex decreases. The same changes take place when one excites the light-hole (LHX) exciton ($`h\nu =1.540`$ eV). Curve 3 shows the PL excitation spectrum (PLE) of the GaAs layer exciton line when the laser frequency is tuned through the exciton resonances in the 14 nm QW.

One can see sharp enhancements in the PL intensity of the 100 nm GaAs exciton peak under resonant excitation of heavy-hole and light-hole excitons in the QW. The PLE spectrum of the D0X complex evolves with phase opposite to that of the exciton PLE. There is also a striking interplay between the exciton and trion lines in the spectrum of the QW PL. The solid line in Fig.3c shows that the PL trion peak (1.5326 eV) dominates over the HHX PL (1.534 eV) when one excites the LHX transition in the 14 nm QW. The full (open) circles show the PLE spectrum of the QW exciton (trion) peak. One can see that resonant excitation of excitonic transitions favors the trion line whereas the exciton line is stronger at excitation out of resonance. The quenching of D0X in the PL of bulk GaAs and the simultaneous sharp enhancement of the trion in the PL spectrum of the QW indicate that the GaAs layer is depleted of electrons whereas the 14 nm QW is enriched of them . This phenomenon disappears when the pump density is below 0.1 mW/cm2. In this case the PL of the GaAs layer consists of D0X and X lines and the QW PL reveals only the exciton line. These results can be understood qualitatively with the use of Fig.1b. Without excitation of the QW, the GaAs layer is charged negatively: besides donor-bound electrons there are extra free electrons escaping the donors from the surrounding barriers and wells. The excitation of the QW induces their redistribution, depleting the GaAs layer and filling the QW. One possibility is the tunneling of the QW photoexcited holes into the 100 nm GaAs layer and recombination with electrons. Another possibility is the annihilation of one QW exciton with the simultaneous ionization of another one through Auger-process. In this case the hole goes into the GaAs by jumping over the barrier while the electron goes back due to repulsion from the negatively charged layer.

Overall, the experimental data clearly show that the density of electrons in the 100 nm layer is reduced and that in the QW is increased independent of the specific mechanism responsible for the recharging effect. In this paper we will exploit this phenomenon to manipulate the 100 nm layer electron concentration optically in order to trace the electron spin relaxation in the 100 nm GaAs layer. We show that the depletion of the 100 nm GaAs switches on a fast electron spin relaxation mechanism (5 ns) involving lattice nuclei, whereas the presence of extra electrons in the layer strongly suppresses it (up to 300 ns).

The electronic spin relaxation in GaAs can be probed accurately with the use of the Hanle effect under optical pumping conditions . Circularly polarized light creates electronic spin polarization oriented along the exciting beam. The degree of circular polarization, $`\rho `$, is determined by the projection, $`S_z,`$ of the mean electron spin on the pump beam direction: $`\rho =S_z`$. An external magnetic field in Voigt geometry (i.e. a transverse field) decreases the steady state PL polarization due to Larmor precession with frequency $`\omega =\mu _Bg_eB/\mathrm{}`$ ($`\mu _B`$ is the Bohr magneton, $`g_e=0.44`$  is the electron $`g`$ -factor). The field dependence of $`\rho `$ is described by a Lorenzian curve with halfwidth, $`B_{1/2}=\mathrm{}/\mu _Bg_eT_s`$ (Hanle effect), where the electron spin lifetime $`T_s=\left(1/\tau +1/\tau _s\right)^1`$ is determined by recombination $`1/\tau `$ and spin relaxation $`1/\tau _s`$ rates. The electron lifetime in n-type semiconductors depends strongly on the pumping intensity (photo-hole concentration), and in the low pumping limit is very long $`\left(\tau >>\tau _s\right)`$. In this limit the spin lifetime coincides with the spin relaxation time, e.g. $`T_s\left(W=0\right)=\tau _s`$, and determines directly the linewidth of the Hanle curve.

Figure 2 shows the Hanle depolarization curves measured at different laser energies with fixed power density, $`W=40`$ mW/cm2. Curve 1 shows the data measured under quasiresonant creation of the exciton of the 100 nm layer. The halfwidth value, $`B_{1/2}=2.9`$ G, corresponds to $`T_s=90`$ ns. The spin lifetime becomes shorter when one excites above the GaAs energy gap (curve 2): $`T_s=40`$ ns ($`B_{1/2}=6.8`$ G). The insert in Fig.2 shows the dependence of $`1/T_s`$ on $`W`$ under quasiresonant excitation of the GaAs exciton. The low power limit of this dependence gives $`\tau _s=\left(290\pm 30\right)`$ ns.

Increasing the laser frequency, with power fixed, up to the heavy-hole exciton energy in the 14 nm QW brings no drastic changes to the halfwidth of the Hanle curve, $`B_{1/2}`$. The situation changes dramatically when either the QW heavy-hole or light-hole exciton is excited: $`B_{1/2}`$ broadens to 50 G. As an example we present the data obtained under resonant LHX excitation (curve 3 in Fig.2). When the exciting frequency is between the HHX and LHX, the QW is not excited and the polarization behavior of the $`100`$ nm layer luminescence is similar to curve 2. Further increasing the laser frequency above the QW resonances increases the QW excitation and $`B_{1/2}`$ again approaches to $`50`$ G. We conclude that there is a large quantitative difference in $`B_{1/2}`$ of the 100 nm layer when the QW is excited (the layer is depleted of electrons) or not (extra free electrons fill up the layer). The change of $`B_{1/2}`$ can not be connected with spin transfer from the QW across the 25 nm AlGaAs barrier during the recharging process. Direct evidence against this is that both the zero-field polarization values ($`+2\%`$) and the Hanle effects in the GaAs layer are identical under LHX (Fig.2, curve 3) and HHX (Fig.3, b) excitation of the 14 nm QW. However the QW exciton polarization has opposite signs in these two cases due to selection rules (experimental values are $`10\%`$ for LHX and $`+20\%`$ for HHX excitations for positive laser polarization).

Remarkably, qualitative differences between the two regimes appear when the magnetic field is applied in the Faraday geometry, i.e. under a longitudinal field. In Figure 3, polarization data taken in the Faraday and Voigt geometry are compared for the two regimes: when the QW is not excited and extra electrons remain in the 100 nm layer (Fig.3a); and when the QW is excited and the electron population is reduced (Fig.3b). With the extra electron population, there is practically no magnetic field dependence in the Faraday geometry (upper part of Fig.3a). In contrast, when the electron population is reduced, application of the magnetic field in the Faraday geometry increases the polarization by a factor of 2.5 (upper part of Fig.3b). Moreover, although there is no correlation in the change in polarization in the two geometries in Fig. 3a, the characteristic transverse and longitudinal magnetic field values that decrease and increase, respectively, the electron polarization are surprisingly similar in Fig. 3b. This difference in behavior reflects a qualitative difference in the spin relaxation process in the two experimental regimes, i.e. with high (Fig.3a) and low (Fig.3b) electron concentrations in the GaAs layer.

The experimental data shown in Fig.3 correspond to the two limiting cases of the electron spin relaxation in a rapidly fluctuating or static random magnetic field. The spin relaxation can be viewed as the result of the action of random local magnetic fields on electron spin  and depends both on the amplitude of the random field and on its fluctuation rate. In the static limit the characteristic time ($`\tau _c`$) of the local field fluctuation is much longer than the period of electron spin precession in this field ($`\omega _f^1`$) and the electron spins undergo precession in the static fields distributed randomly. The average spin of the ensemble relaxes to 1/3 of its initial value with a characteristic time $`\tau _s\omega _f^1`$ . The longitudinal external field eliminates spin relaxation, and the transverse field depolarizes electrons when it overwhelms the local field, i.e. $`\omega \omega _f`$. In the opposite limit: $`\tau _c<<\omega _f^1`$ (the so-called motional narrowing case ), dynamic averaging takes place and the spin relaxation slows down: $`\tau _s^1\omega _f^2\tau _c`$. In this case, the transverse magnetic field that depolarizes the electron spin (Hanle effect) satisfies the condition $`\omega \omega _t\tau _s^1\omega _f^2\tau _c`$. The longitudinal field required to suppress spin relaxation satisfies the condition: $`\omega \omega _l\tau _c^1`$,(see Ref.). It happens in much greater magnetic fields because $`\omega _t=\tau _s^1\omega _f^2\tau _c=\left(\omega _f\tau _c\right)^2\omega _l<<\omega _l`$. Therefore, the experimental result of Fig.3a, with additional electrons in the layer, corresponds to the motional narrowing case $`\left(\omega _l>>\omega _t\right)`$, and the result of Fig.3b, with electrons removed, corresponds to the static limit $`\left(\omega _l\omega _t\right)`$.

The experimental data obtained in the static case enables us to deduce that the main spin relaxation mechanism of electrons in this experiment is the hyperfine interaction of localized electrons with lattice nuclei. It follows from Fig.3b that the characteristic precession frequency of the electron spins in the random local field is $`\omega _f210^8s^1`$. Therefore $`\tau _c>\tau _s\omega _f^15`$ ns. None of the known spin relaxation mechanisms of free carriers and excitons (Chap.3 in Ref.) can have such long correlation times. However, bound electrons at low electron concentration and low temperature may spend a long time at a given donor. Moreover, the $`\omega _f`$ value agrees well with the calculated precession frequency of spins of the donor-bound electrons in the random hyperfine fields of surrounding nuclei as we now show.

To describe the average electron spin polarization in the low concentration limit we consider an evolution of the electron spin, $`𝐒_n`$, localized at the $`n`$-th donor. It is described by the Bloch equation  with the precession frequency, $`𝛀_n=\omega +\omega _{fn}`$:

$$\frac{d𝐒_n}{dt}=\frac{𝐒_i𝐒_n}{\tau _c}+𝛀_n\times 𝐒_n=0$$
(1)
where $`𝐒_i`$ is the initial electron spin value in the moment of trapping. During the correlation time, $`\tau _c`$, the electron spin dynamics is determined by the external and the local nuclear magnetic fields only. The static nuclear fields are described by a Gaussian distribution, $`\mathrm{\Phi }\left(B_N\right)=\pi ^{3/2}\mathrm{\Delta }^3\mathrm{exp}\left(B_N^2/\mathrm{\Delta }^2\right)`$, with a dispersion, $`\mathrm{\Delta }=54`$ G, given by the root-mean-square of the random nuclear field . The average nuclear polarization is zero because we are modulating the polarization of exciting light. Averaging the solution of Eq.1 over the nuclear magnetic field distribution  results in a theoretical dependence of circular polarization on the longitudinal and transverse magnetic field that are presented in Fig.3b by dashed and solid lines, respectively. The only fitting parameter for the curves in Fig.3b is $`\tau _c=17`$ ns, which satisfies the static condition, $`\tau _c>\omega _f^1`$. This value is determined by the processes of donor electron ionization, recombination, exchange scattering by free electrons, etc. The good agreement between experiment and theory confirms that the interaction with nuclei is the main mechanism of the electron spin relaxation for the localized carriers.

The much longer spin relaxation time in the motional narrowing case ($`\tau _s300`$ ns) arises from the increase in electron concentration. Recall that the presence of extra free electrons in addition to the “intrinsic” donor-bound ones, comes from donors in the barriers and quantum wells surrounding the GaAs layer. The suppression of spin relaxation is caused by the decrease of the correlation time (from $`\tau _c>`$5 ns to $`\tau _c0.1`$ ns ) and can be explained by scattering processes between the free and bound electrons discussed in Ref.. These processes average out the inhomogeneous nuclear field and suppress spin relaxation by two orders of magnitude. The exciton PL serves as a detector of spin polarization of donor-bound electrons .

In conclusion, we have shown that the spin relaxation time of donor-bound electrons in a thick GaAs layer depends strongly on the concentration of extra electrons. This enabled us to manipulate the spin relaxation rate by two orders of magnitude. The presence of extra free electrons in the conduction band suppresses strongly the spin relaxation of the donor bound electrons because the inhomogeneous nuclear field is averaged out, leading to an extremely long spin relaxation time of $`0.3\mu `$s. Decrease of the electron concentration breaks up the electron spin system into an ensemble of individual donor-bound electrons, and the electronic spin relaxation is the result of dephasing by the randomly distributed hyperfine fields. In our experiments, the average of the nuclear spin polarization is zero but its dispersion is approximately 54 G, corresponding to an inhomogeneous dephasing rate of 5 ns. The effect of inhomogeneous dephasing can be eliminated with the use of the spin echo technique Ref.. This case may be important in efforts to implement quantum computing using localized electronic spins, which require not only long coherence times but also must be localized in order to be individually addressed.

Partial support of the Russian Basic Research Foundation (grants 00-02-16991 and 01-02-17906), CRDF (grant RP1-2252), DARPA/SPINS, and ONR is gratefully acknowledged.


text:
# Exact solutions of time-dependent three-generator systems

## I Introduction

Exact solutions and geometric phase factor of time-dependent spin model have been extensively investigated by many authors. Bouchiat and Gibbons discussed the geometric phase for the spin-$`1`$ system. Datta et al found the exact solution for the spin- $`\frac{1}{2}`$ system  by means of the classical Lewis-Riesenfeld theory, and Mizrahi calculated A-A phase for the spin- $`\frac{1}{2}`$ system  in a time-dependent magnetic field. The more systematic approach to obtaining the formally exact solutions for the spin- $`j`$ system was proposed by Gao et al who made use of the Lewis-Riesenfeld quantum theory. In this spin- $`j`$ system, the three generators of the Hamiltonian satisfy the commutation relations of $`SU(2)`$ Lie algebra. In addition to the spin model, there exist many quantum systems whose Hamiltonian is also constructed in terms of three generators of various Lie algebras which we will illustrate in the following.

The invariant theory which can be better applied to solutions of the time-dependent Schrödinger equation was first proposed by Lewis and Riesenfeld in 1969. This theory is very appropriate for treating the geometric phase factor. In 1991, Gao et al generalized this theory and put forward the invariant-related unitary transformation method. Exact solutions for time-dependent systems obtained by using the generalized invariant theory contain the geometric and dynamical phase. This formulation was developed from the Lewis-Riesenfeld<sup>,</sup>s formal theory into a powerful tool for treating exact solutions of the time-dependent Schrödinger equation and geometric phase factor. In the present paper, we obtain exact solutions of various time-dependent three-generator systems based on these invariant theories.

This paper is organized as follows: In section 2, we set out several quantum systems or models which have three-generator algebraic structure. In section 3, use is made of the invariant theories and exact solutions of various time-dependent three-generator systems are obtained. In section 4, there are some discussions concerning the closure property of the Lie algebraic generators in the sub-Hilbert-space. In section 5, we conclude with some remarks.

## II The algebraic structures of various three-generator systems

In our previous work we showed that the time-dependent Schrödinger equation is solvable if its Hamiltonian is constructed in terms of the generators of Lie algebra. Then analyzing the algebraic structures of Hamiltonians plays significant role in obtaining exact solutions of the time -dependent systems. To our knowledge, a large number of quantum systems which have three-generator Hamiltonians have been investigated in the literature. In the following we set out these systems and discuss the algebraic structures of their Hamiltonians.

(1) Spin model. The time evolution of the wavefunction of a spinning particle in a magnetic field was studied by regarding it as a spin model whose Hamiltonian can be written

$$H(t)=c_0\{\frac{1}{2}\mathrm{sin}\theta \mathrm{exp}[i\phi ]J_++\frac{1}{2}\mathrm{sin}\theta \mathrm{exp}[i\phi ]J_{}+\mathrm{cos}\theta J_3\}$$
(1)
with $`J_\pm =J_1\pm iJ_2`$ satisfying the commutation relations $`[J_3,J_\pm ]=\pm J_\pm ,[J_+,J_{}]=2J_3.`$ Analogous to this case, in the gravitational theory of general relativity both spin-gravitomagnetic interaction and spin-rotation coupling are proved spin models. It can be verified that the investigation of the propagation of a photon inside the noncoplanarly curved optical fiber  is also equivalent to that of a spin model. The Hamiltonian of spin model is composed of three generators which constitute $`SU(2)`$ algebra.

(2) Two-coupled harmonic oscillator. The Hamiltonian of the two-coupled harmonic oscillator is of the form ( in the units $`\mathrm{}=1`$)

$$H=\omega _1a_1^{}a_1+\omega _2a_2^{}a_2+ga_1^{}a_2+g^{}a_2^{}a_1,$$
(2)
where $`a_1^{},a_2^{},a_1,a_2`$ are the creation and annihilation operators for these two harmonic oscillators, respectively; $`g`$ and $`g^{}`$ are the coupling coefficients. Set $`J_+=a_1^{}a_2,J_{}=a_2^{}a_1,J_3=\frac{1}{2}(a_1^{}a_1a_2^{}a_2),N=\frac{1}{2}(a_1^{}a_1+a_2^{}a_2),`$ then we can show that the generators of this Hamiltonian satisfy the commutation relations of $`SU(2)`$ algebra. Since $`N`$ commutes with $`H`$, i.e., $`[N,H]=0,`$we consequently say $`N`$ is an invariant (namely, it is a conserved generator whose eigenvalue is time-independent). In terms of $`J_\pm ,J_3`$ and $`N,`$ the Hamiltonian in the expression (2) can be rewritten as follows

$$H=\omega _1(N+J_3)+\omega _2(NJ_3)+gJ_++g^{}J_{}.$$
(3)
Another interesting Hamiltonian of the two-coupled harmonic oscillator is written in the form

$$H=\omega _1a_1^{}a_1+\omega _2a_2^{}a_2+ga_1a_2+g^{}a_1^{}a_2^{}.$$
(4)
If we take $`K_+=a_1^{}a_2^{},K_{}=a_1a_2,K_3=\frac{1}{2}(a_1a_1^{}+a_2^{}a_2),N=\frac{1}{2}(a_1^{}a_1a_2^{}a_2),`$ the generators of the $`SU(1,1)`$ group are thus realized. The commutation relations are immediately inferred as

$$[K_3,K_\pm ]=\pm K_\pm ,[K_+,K_{}]=2K_3.$$
(5)
(3) $`SU(1,1)_sh(4)`$ Lie-algebra system. A good number of quantum systems whose Hamiltonian is some combinations of the generators of a Lie algebra, e.g., $`SU(1,1)_sh(4)`$ ($`_s`$denotes a semidirect sum) which is used to discuss both the non-Poissonian effects in a laser-plasma scattering and the pulse propagation in a free-electron laser. The $`SU(1,1)_sh(4)`$ Hamiltonian is

$$H=AK_3+FK_++F^{}K_{}+Ba^{}+B^{}a+G,$$
(6)
where $`a`$ and $`a^{}`$ are harmonic-oscillator annihilation and creation operators, respectively.

(4) General harmonic oscillator. The Hamiltonian of the general harmonic oscillator is given by

$$H=\frac{1}{2}[Xq^2+Y(qp+pq)+Zp^2]+Fq,$$
(7)
where the canonical coordinate $`q`$ and the canonical momentum $`p`$ satisfy the commutation relation $`[q,p]=i.`$ The following three-generator Lie algebra is easily derived

$$[q^2,p^2]=2\{i(qp+pq)\},[i(qp+pq),q^2]=4q^2,[i(qp+pq),p^2]=4p^2.$$
(8)
(5) Charged particle moving in a magnetic field. The motion of a particle with mass $`\mu `$ and charge $`e`$ in a homogeneous magnetic field $`\stackrel{}{B}=(0,0,B)`$ is described by the following Hamiltonian in the spherical coordinates

$$H=\frac{1}{2\mu }(\frac{^2}{r^2}+\frac{2}{r}\frac{}{r}\frac{1}{r^2}L^2)+\frac{1}{8}\mu \omega ^2r^2\frac{\omega }{2}L_z$$
(9)
with $`\omega =\frac{B}{\mu }.`$ Since both $`L_z`$ and $`L^2`$ commute with $`H`$ and thus they are called invariants, only the operators associated with $`r`$ should be taken into consideration. We can show that if the following operators are defined

$$K_1=\mu r^2,K_2=\frac{1}{\mu }(\frac{^2}{r^2}+\frac{2}{r}\frac{}{r}\frac{1}{r^2}L^2),K_3=2i(\frac{3}{2}+r\frac{}{r}),$$
(10)
$`K_1,K_2`$ and $`K_3`$ form an algebra

$$[K_1,K_2]=2iK_3,[K_3,K_2]=4iK_2,[K_3,K_1]=4iK_1.$$
(11)
Apparently, $`H`$ can be rewritten in terms of the generators of this Lie algebra.

(6) Two-level atomic coupling. The model under consideration is consisted of two-level atom driven by the photons field. The interaction part of the Hamiltonian contains the transition operator $`|12|`$ and $`|21|`$, where $`|1`$ and $`|2`$ are the atomic operators of the two-level atom. Simple calculation yields

$`[|11||22|,|12|]`$ $`=`$ $`2|12|,[|11||22|,|21|]=2|21|,`$ (12)
$`[|12|,|21|]`$ $`=`$ $`|11||22|,`$ (13)
which unfolds that the Hamiltonian contains a $`SU(2)`$ algebraic structure.

(7) Supersymmetric Jaynes-Cummings model. In addition to the ordinary Jaynes-Cummings models, there exists a two-level multiphoton Jaynes-Cummings model which possesses supersymmetric structure. In this generalization of the Jaynes-Cummings model, the atomic transitions are mediated by $`k`$ photons . Singh has shown that this model can be used to study multiple atom scattering of radiation and multiphoton emission, absorption, and laser processes. The Hamiltonian of this model under the rotating wave approximation is given by

$$H(t)=\omega (t)a^{}a+\frac{\omega _0(t)}{2}\sigma _z+g(t)(a^{})^k\sigma _{}+g^{}(t)a^k\sigma _+,$$
(14)
where $`a^{}`$ and $`a`$ are the creation and annihilation operators for the electromagnetic field, and obey the commutation relation $`[a,a^{}]=1`$; $`\sigma _\pm `$ and $`\sigma _z`$ denote the two-level atom operators which satisfy the commutation relation $`[\sigma _z,\sigma _\pm ]=\pm 2\sigma _\pm `$ . We can verify that this model is solvable and the complete set of exact solutions can be found by working in a sub-Hilbert-space corresponding to a particular eigenvalue of the supersymmetric generator $`N^{^{}}`$

$$N^{^{}}=\left(\begin{array}{cc}a^k(a^{})^k& 0\\ 0& (a^{})^ka^k\end{array}\right).$$
(15)
It can be verified that $`N^{^{}}`$ commutes with the Hamiltonian in (14), $`N^{^{}}`$ is therefore called the time-independent invariant. Vogel and Welsch have studied the $`k`$-photon Jaynes-Cummings model with coherent atomic preparation which is time-independent. In the framework of the formulation presented in this paper, we can study the totally time-dependent cases of work done by Vogel and Welsch.

(8) Two-level atom interacting with a generalized cavity. Consider the following Hamiltonian

$$H=r(A_0)+s(A_0)\sigma _z+gA_{}\sigma _++g^{}A_+\sigma _{}$$
(16)
where $`r(A_0)`$ and $`s(A_0)`$ are well-defined real functions of $`A_0,`$ and $`A_0,A_\pm `$ satisfy the commutation relations $`[A_0,A_\pm ]=\pm mA_\pm \text{[29]}`$. One can show that this Hamiltonian possesses a three-generator algebraic structure.

(9) The interaction between a hydrogenlike atom and an external magnetic field. This model is described by

$`H`$ $`=`$ $`\alpha \stackrel{}{L}\stackrel{}{S}+\beta (L_z+2S_z)`$ (17)
$`=`$ $`\beta L_z+({\displaystyle \frac{1}{2}}\alpha L_z+\beta )\sigma _z+{\displaystyle \frac{1}{2}}\alpha (L_{}\sigma _++L_+\sigma _{})`$ (18)
with $`L_\pm `$ $`=L_x\pm iL_y.`$ It is evidently seen that this form of Hamiltonian is analogous to that in (16).

(10) Coupled two-photon lasers. The Hamiltonian of this model is in fact the combination of that of the two-coupled harmonic oscillator and general harmonic oscillator. One can show that there exists a $`SU(2)`$ algebraic structure in this model.

From what has been discussed above we can draw a conclusion that a number of typical and useful systems or models in laser physics, atomic physics and quantum optics can be attributed to various three-generator type. The analysis of these algebraic structures shows the solvability of these systems. It should be noted that, in the literature, most of above systems or models are only investigated in the stationary cases where the coefficients of the Hamiltonians are totally time-independent. Some systems are investigated for the Hamiltonian with partly time-dependent coefficients. In the present paper, we will give exact solutions of the time-dependent Schrödinger equation of all these systems or models where the coefficients of the Hamiltonians are totally time-dependent.

## III  Exact solutions of time-dependent Schrödinger equation

Time evolution of most above systems or models is governed by the Schrödinger equation

$$i\frac{|\mathrm{\Psi }(t)_s}{t}=H(t)|\mathrm{\Psi }(t)_s,$$
(19)
where the Hamiltonian is constructed by three generators $`A,B`$ and $`C`$ and is often given as follows

$$H(t)=\omega (t)\{\frac{1}{2}\mathrm{sin}\theta (t)\mathrm{exp}[i\varphi (t)]A+\frac{1}{2}\mathrm{sin}\theta (t)\mathrm{exp}[i\varphi (t)]B+\mathrm{cos}\theta (t)C\}$$
(20)
with $`A,B`$ and $`C`$ satisfying the general commutation relations of a Lie algebra

$$[A,B]=nC,[C,A]=mA,[C,B]=mB,$$
(21)
where $`m`$ and $`n`$ are the structure constants of this Lie algebra. Since all simple $`3`$-generator algebras are either isomorphic to the algebra $`sl(2,C)`$ or to one of its real forms, we treat these time-dependent quantum systems in a unified way. According to the Lewis-Riesenfeld invariant theory, an operator $`I(t)`$ that agrees with the following invariant equation

$$\frac{I(t)}{t}+\frac{1}{i}[I(t),H(t)]=0$$
(22)
is called an invariant whose eigenvalue is time-independent, i.e.,

$$I(t)|\lambda ,t_I=\lambda |\lambda ,t_I,\frac{\lambda }{t}=0.$$
(23)
It is seen from Eq. (22) that $`I(t)`$ is the linear combination of $`A,B`$ and $`C`$

and may be generally written

$$I(t)=y\{\frac{1}{2}\mathrm{sin}a(t)\mathrm{exp}[ib(t)]A+\frac{1}{2}\mathrm{sin}a(t)\mathrm{exp}[ib(t)]B\}+\mathrm{cos}a(t)C,$$
(24)
where the constant $`y`$ will be determined below. It should be pointed out that it is not the only way to construct the invariants. Since the product of two invariants also satisfies Eq. (22), there are infinite invariants of a time-dependent quantum system. But the form in Eq. (24) is the most convenient and useful one. Substitution of (24) into Eq.(22) yields

$`y\mathrm{exp}(ib)(\dot{a}\mathrm{cos}ai\dot{b}\mathrm{sin}a)im\omega [\mathrm{exp}(i\varphi )\mathrm{cos}a\mathrm{sin}\theta y\mathrm{exp}(ib)\mathrm{sin}a\mathrm{cos}\theta ]`$ $`=`$ $`0,`$ (25)
$`\dot{a}+{\displaystyle \frac{ny}{2}}\omega \mathrm{sin}\theta \mathrm{sin}(b\varphi )`$ $`=`$ $`0.`$ (26)
where dot denotes the time derivative. The time-dependent parameters $`a`$ and $`b`$ are determined by these two auxiliary equations.

It is easy to verify that the particular solution $`|\mathrm{\Psi }(t)_s`$ of the Schrödinger equation can be expressed in terms of the eigenstate $`|\lambda ,t_I`$ of the invariant $`I(t),`$ namely,

$$|\mathrm{\Psi }(t)_s=\mathrm{exp}[\frac{1}{i}\phi (t)]|\lambda ,t_I$$
(27)
with

$$\phi (t)=_{0I}^t\lambda ,t^{^{}}\left|[H(t^{^{}})i\frac{}{t^{^{}}}]\right|\lambda ,t^{^{}}_I𝑑t^{^{}}.$$
(28)
The physical meanings of $`_{0I}^t\lambda ,t^{^{}}\left|H(t^{^{}})\right|\lambda ,t^{^{}}_I𝑑t^{^{}}`$ and $`_{0I}^t\lambda ,t^{^{}}\left|i\frac{}{t^{^{}}}\right|\lambda ,t^{^{}}_I𝑑t^{^{}}`$ are dynamical and geometric phase, respectively.

Since the expression (27) is merely a formal solution of the Schrödinger equation, in order to get the explicit solutions we make use of the invariant-related unitary transformation formulation which enables one to obtain the complete set of exact solutions of the time-dependent Schrödinger equation (19). In accordance with the invariant-related unitary transformation method, the time-dependent unitary transformation operator is often of the form

$$V(t)=\mathrm{exp}[\beta (t)A\beta ^{}(t)B]$$
(29)
with $`\beta (t)=\frac{a(t)}{2}x\mathrm{exp}[ib(t)],\beta ^{}(t)=\frac{a(t)}{2}x\mathrm{exp}[ib(t)].`$ By making use of the Glauber formulae, lengthy calculation yields

$`I_V`$ $`=`$ $`V^{}(t)I(t)V(t)=\{{\displaystyle \frac{y}{2}}\mathrm{exp}(ib)\mathrm{sin}a\mathrm{cos}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}ax]{\displaystyle \frac{(\frac{mn}{2})^{\frac{1}{2}}}{n}}\mathrm{exp}(ib)\mathrm{cos}a\mathrm{sin}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}ax]\}A`$ (32)
$`+\{{\displaystyle \frac{y}{2}}\mathrm{exp}(ib)\mathrm{sin}a\mathrm{cos}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}ax]{\displaystyle \frac{(\frac{mn}{2})^{\frac{1}{2}}}{n}}\mathrm{exp}(ib)\mathrm{cos}a\mathrm{sin}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}ax]\}B`$
$`+\{\mathrm{cos}a\mathrm{cos}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}ax]+{\displaystyle \frac{(\frac{mn}{2})^{\frac{1}{2}}}{m}}y\mathrm{sin}a\mathrm{sin}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}ax]\}C.`$
It can be easily seen that when

$$y=\frac{m}{(\frac{mn}{2})^{\frac{1}{2}}},x=\frac{1}{(\frac{mn}{2})^{\frac{1}{2}}},$$
(33)
one may derive that $`I_V=C`$ which is time-independent. Thus the eigen-equation of the time-independent invariant $`I_V`$ may be written in the form

$$I_V|\lambda =\lambda |\lambda ,|\lambda =V^{}(t)|\lambda ,t_I.$$
(34)
Under the transformation $`V(t),`$ the Hamiltonian $`H(t)`$ can be changed into

$`H_V(t)`$ $`=`$ $`V^{}(t)H(t)V(t)V^{}(t)i{\displaystyle \frac{V(t)}{t}}`$ (35)
$`=`$ $`\{\omega [\mathrm{cos}a\mathrm{cos}\theta +{\displaystyle \frac{(\frac{mn}{2})^{\frac{1}{2}}}{m}}\mathrm{sin}a\mathrm{sin}\theta \mathrm{cos}(b\varphi )]+{\displaystyle \frac{\dot{b}}{m}}(1\mathrm{cos}a)\}C`$ (36)
by the aid of Baker-Campbell-Hausdorff formula

$$V^{}(t)\frac{}{t}V(t)=\frac{}{t}L+\frac{1}{2!}[\frac{}{t}L,L]+\frac{1}{3!}[[\frac{}{t}L,L],L]+\frac{1}{4!}[[[\frac{}{t}L,L],L],L]+\mathrm{}$$
(37)
with $`V(t)=\mathrm{exp}[L(t)].`$ Hence, with the help of Eq.(27) and Eq.(34), the particular solution of the Schrödinger equation is obtained

$$|\mathrm{\Psi }(t)_s=\mathrm{exp}[\frac{1}{i}\phi (t)]V(t)|\lambda $$
(38)
with the phase

$`\phi (t)`$ $`=`$ $`{\displaystyle _0^t}\lambda \left|[V^{}(t^{^{}})H(t^{^{}})V(t^{^{}})V^{}(t^{^{}})i{\displaystyle \frac{}{t^{^{}}}}V(t^{^{}})]\right|\lambda 𝑑t^{^{}}=\phi _d(t)+\phi _g(t)`$ (39)
$`=`$ $`\lambda {\displaystyle _0^t}\{\omega [\mathrm{cos}a\mathrm{cos}\theta +{\displaystyle \frac{(\frac{mn}{2})^{\frac{1}{2}}}{m}}\mathrm{sin}a\mathrm{sin}\theta \mathrm{cos}(b\varphi )]+{\displaystyle \frac{\dot{b}}{m}}(1\mathrm{cos}a)\}𝑑t^{^{}},`$ (40)
where the dynamical phase is $`\phi _d(t)=\lambda _0^t\omega [\mathrm{cos}a\mathrm{cos}\theta +\frac{(\frac{mn}{2})^{\frac{1}{2}}}{m}\mathrm{sin}a\mathrm{sin}\theta \mathrm{cos}(b\varphi )]𝑑t^{^{}}`$ and the geometric phase is $`\phi _g(t)=\lambda _0^t\frac{\dot{b}}{m}(1\mathrm{cos}a)𝑑t^{^{}}.`$ It is seen that the former phase is in connection with the coefficients of the Hamiltonian such as $`\omega ,\mathrm{cos}\theta ,\mathrm{sin}\theta ,`$ etc., whereas the latter is not immediately related to these coefficients. If the parameter $`a`$ is taken to be time-independent, $`\phi _g(T)=\lambda _0^T\frac{\dot{b}}{m}(1\mathrm{cos}a)𝑑t^{^{}}=\frac{\lambda }{m}[2\pi (1\mathrm{cos}a)]`$ where $`2\pi (1\mathrm{cos}a)`$ is an expression for the solid angle over the parameter space of the invariant. It is of interest that $`\frac{\lambda }{m}[2\pi (1\mathrm{cos}a)]`$ is equal to the magnetic flux produced by a monopole of strength $`\frac{\lambda }{4\pi m}`$ existing at the origin of the parameter space. This, therefore, implies that geometric phase differs from dynamical phase and it involves the global and topological properties of the time evolution of a quantum system. This fact indicates the geometric or topological meaning of $`\phi _g(t).`$

The expression (38) is a particular exact solution corresponding to $`\lambda `$ and the general solutions of the time-dependent Schrödinger equation are easily obtained by using the linear combinations of all these particular solutions. Generally speaking, in Quantum Mechanics, solution with chronological-product operator (time-order operator) $`P`$ is often called the formal solution. In the present paper, however, the solution of the Schrödinger equation governing a time-dependent system is sometimes termed the explicit solution, for reasons of the fact that it does not involve time-order operator. But, on the other hand, by using Lewis-Riesenfeld invariant theory, there always exist time-dependent parameters, for instance, $`a(t)`$ and $`b(t)`$ in this paper which are determined by the auxiliary equations (26). In the traditional practice, when employed in experimental analysis and compared with experimental results, these nonlinear auxiliary equations should be solved often by means of numerical calculation. From above viewpoints, the concept of explicit solution is understood in a relative sense, namely, it can be considered explicit solution when compared with the time-evolution operator $`U(t)=P\mathrm{exp}[\frac{1}{i}_0^tH(t^{^{}})𝑑t^{^{}}]`$ involving time-order operator, $`P`$; whereas, it cannot be considered completely explicit solution for it is expressed in terms of some time-dependent parameters which should be obtained via the auxiliary equations. Hence, conservatively speaking, we regard the solution of the time-dependent system presented in the paper as exact solution rather than explicit solution.

## IV Discussions

In the previous section, we obtain exact solutions of some time-dependent three-generator systems or models by using these invariant theories. In what follows there is the closure property of the Lie algebraic generators in the sub-Hilbert-space that should be discussed.

The generalized invariant theory can only be applied to the study of the system for which there exists the quasialgebra defined in Ref.. It is easily seen from (16) that there is no such quasialgebra for the Hamiltonian (16) of example (7)(supersymmetric Jaynes-Cummings model) and example (8)(two-level atom interacting with a generalized cavity). We generalize the method which has been used for finding the dynamical algebra $`O(4)`$ of the hydrogen atom to treat this type of time-dependent models. In the case of hydrogen, the dynamical algebra $`O(4)`$ was found by working in the sub-Hilbert-space corresponding to a particular eigenvalue of the Hamiltonian . In this paper, we will show that a generalized quasialgebra can also be found by working in a sub-Hilbert-space corresponding to a particular eigenvalue of $`\mathrm{\Delta }=A_0+m\frac{1+\sigma _z}{2}`$ in the time-dependent model of two-level atom interacting with a generalized cavity. This generalized quasialgebra enables one to obtain the complete set of exact solutions for the Schrödinger equation. It is easily verified that $`\mathrm{\Delta }`$ commutes with $`H(t)`$ and is a time-independent invariant according to Eq. (22). Then in order to unfold the algebraic structure of the Hamiltonian (16), the following three operators are defined

$$\mathrm{\Sigma }_1=\frac{1}{2[\chi (\mathrm{\Delta })]^{\frac{1}{2}}}(A_{}\sigma _++A_+\sigma _{}),\mathrm{\Sigma }_2=\frac{i}{2[\chi (\mathrm{\Delta })]^{\frac{1}{2}}}(A_+\sigma _{}A_{}\sigma _+),\mathrm{\Sigma }_3=\frac{1}{2}\sigma _z,$$
(41)
where $`\chi =n\left|A_+A_{}\right|n,|n`$ denotes the eigenstates of $`A_0.`$ It is easy to see all these operators commute with $`\mathrm{\Delta }`$ and the quasialgebra $`\{H,\mathrm{\Sigma }_1,\mathrm{\Sigma }_2,\mathrm{\Sigma }_3\}`$ is thus found. Therefore, this type of time-dependent models is proved solvable by working in a sub-Hilbert-space corresponding to the eigenstates of the time-independent invariant.

For the supersymmetric Jaynes-Cummings model, the commutation relations of its supersymmetric Lie- algebraic structure are

$`[Q^{},Q]`$ $`=`$ $`\lambda \sigma _z,[N,Q]=Q,[N,Q^{}]=Q^{},`$ (42)
$`[Q,\sigma _z]`$ $`=`$ $`2Q,[Q^{},\sigma _z]=2Q^{},`$ (43)
where

$$N=a^{}a+\frac{k1}{2}\sigma _z+\frac{1}{2},Q=(a^{})^k\sigma _{},Q^{}=a^k\sigma _+,$$
(44)
and $`\lambda `$ denotes the eigenvalue of the time-independent invariant $`N^{^{}}.`$ By the aid of (43) and (44), the Hamiltonian (14) of this supersymmetric Jaynes-Cummings model can be rewritten as

$$H(t)=\omega (t)N+\frac{\omega (t)\delta (t)}{2}\sigma _z+g(t)Q+g^{}(t)Q^{}\frac{\omega (t)}{2}.$$
(45)
Then this time-dependent model can be exactly solved by using the invariant-related unitary transformation formulation where the unitary transformation operator is of the form

$$V(t)=\mathrm{exp}[\beta (t)Q\beta ^{}(t)Q^{}].$$
(46)
It should be noted that the above approach to the time-dependent Jaynes-Cummings model is also appropriate for treating the periodic decay and revival of some multiphoton-transitions models which has been investigated by Sukumar and Buck.

It is readily verified that in the sub-Hilbert-space corresponding to a particular eigenvalue of the conserved generator (the time-independent invariant), the Hamiltonian of original two-level Jaynes-Cummings model (mono-photon case) possesses the $`SU(2)`$ Lie-algebraic structure, and three-level two-mode mono-photon model possesses the $`SU(3)`$ structure. The solution of the time-dependent case of $`SU(2)`$ Jaynes-Cummings model is easily obtained by taking the number of photons mediating in the process of atomic transitions $`k=1.`$ Since Shumovsky et al have considered the three-level two-mode multiphoton Jaynes-Cummings model  whose Hamiltonian is time-independent, it is also of interest to exactly solve the time-dependent supersymmetric three-level two-mode multiphoton Jaynes-Cummings model by means of invariant theories.

## V Summary

On the basis of the fact that all simple 3-generator algebras are either isomorphic to the algebra $`sl(2,C)`$ or to one of its real forms, exact solutions of the time-dependent Schrödinger equation of all three-generator systems or models in quantum optics, nuclear physics, solid state physics, molecular and atomic physics and laser physics are offered by making use of the Lewis-Riesenfeld invariant theory and the invariant-related unitary transformation formulation in the present paper. Since it appears only in systems with time-dependent Hamiltonian, the geometric phase factor would be easily investigated if the exact solutions of time-dependent systems had been obtained. In view of above discussions, the invariant-related unitary transformation formulation is a useful tool for treating the geometric phase factor and the time-dependent Schrödinger equation. This formulation replaces the eigenstates of the time-dependent invariants with those of the time-independent invariants through the unitary transformation. Apparently, it is also applicable to the algebraic structure whose number of generators is more than three. Additionally, it should be pointed out that the time-dependent Schrödinger equation is often investigated in the literature, whereas less attention is paid to the time-dependent Klein-Gordon equation. Work in this direction is under consideration and will be published elsewhere.

Acknowledgments This project is supported by the National Natural Science Foundation of China under the project No.$`19775040`$ and $`30000034`$. The authors thank Prof. Gao Xiao-Chun for offering the knowledge concerning the closure property of the Lie algebraic generators in sub-Hilbert-space.