text:
# Studies of group velocity reduction and pulse regeneration with and without the adiabatic approximation

## I Introduction

Propagation of optical pulses in a resonant medium has long been an active field of research. The advances of laser technologies have stimulated further development and understanding of the propagation of optical pulses under various conditions, in particular in the field of ultra-short and extremely high intensity optical pulse generation and propagation under adiabatic conditions. Early works in this field include the interpretation of the nonlinear wave mixing in terms of an adiabatic following model , counter-intuitive pulse sequence  and induced transparency . Later, Eberly and co-workers  developed a theory for the propagation of optical field under the condition that atomic response can be treated adiabatically. Further development on the propagation properties of a pulse pair under the similar adiabatical conditions resulted in substantial understanding of the evolution of optical pulses in a resonant media. In particular, Grobe and co-workers first proposed the creation and recall of excitations distributed in media . Recently, the field has experience a renewed interest and resurgence following the demonstration of extremely slow propagation of optical pulse both in ultra-cold and high temperature atomic gases . Experimental studies on possible excitation of the atomic spin wave excitations with ultra-slow optical field, and the subsequent “revival” of an optical field with frequencies near that of the ultra-slow field have promised future applications in information technology \[7-9\]. Theoretical works on such spin wave excitations and the subsequent “revival” of an optical filed have also been formulated in a “dark-state polariton” theory , which basically is a field theoretical reformulation of the adiabatic approximation of pulse pair propagation in a resonant medium originally formulated before . These new treatments have provided some further understanding of the process, especially from the view point of atomic spin wave excitation and their relation to atomic coherence. Here, we present a study on the propagation and “retrieval” of ultra-slow optical fields. We investigate both analytically and numerically the propagation problem, with and without adiabatic approximation, for the cases of near resonance, on-resonance, and far detuned from resonance. In all cases, we show that the slow group velocity propagation and pulse “revival” is independent of the one-photon detuning by the probe laser, providing that the probe and coupling lasers are tuned to exact two-photon resonance and certain conditions are satisfied. That is, the results on off-resonance excitations are identical to the adiabatic theory formulation for the exact on-resonance case treated in . The main contribution of our work, however, is the detailed analysis, both analytical and numerical, on the extremely slow propagation of the optical field, the corresponding atomic response, and the “revival” of an optical field when a second time-delayed coupling pulse is injected into the medium. These results, to the best of our knowledge, are not contained in the original studies by Eberly and co-workers and all the subsequent studies, including all recent theoretical works, on the subject. In addition, we show that all the experimental observables can be well predicted with the usual treatment combining the classical electrodynamics and the three-state model, without invoking a full field theoretical methodology. This is not surprising, since when a detailed handling of spontaneous emission is not critical and the fields are not extremely weak, a mean-field approximation to the quantum treatment of the electromagnetic field is valid. In Section II, we first present an adiabatic theory for the case where non-vanishing one-photon detuning exists. This treatment is very similar to the original treatment given in  except that the latter deals exclusively with the case of on resonance excitation. In Section III, we examine the propagation of a pulse pair and give two examples on how to obtain the evolution of the optical fields. Section IV constitutes the main contribution of the present study. We discuss in detail the propagation of the pulse pair in the adiabatic limit and deduce analytically substantial insight and subtle understanding on the slow propagation of the probe field, the optical field-atomic spin cohenrence conversion, and the characteristics of the “retrieved” optical field. These results are not available in any of previous studies, and we believe they are very useful in providing a complete picture for the physics involved. In Section V, we compare predictions and estimates derived in Section IV with numerical simulations. Section VI is devoted to the case where the probe laser is tuned far from resonance, i.e. a Raman excitation scheme is employed. In this case, a different set of conditions are required for the adiabatic theory to be accurate. We show, however, that when certain conditions are satisfied, the solution is independent of the detuning, a surprising result that does not, to our best knowledge, exist in literature. Finally, the results are compared with numerical calculations for vigorous testing. In Section VII, we present a conclusion for our study.

## II  Preliminary

In this section, we examine the propagation of a pair of optical pulse in resonant media under the conditions where the adiabatic approximation can be made. We present here an expanded investigation based on the early treatment given in . The difference is the inclusion of the non-vanishing one-photon detuning since the original treatment assumed exact on one-photon resonance excitations. We will show that this non-vanishing one-photon detuning leads to correct predictions about pulse propagation, coherent population transfer and storage, and the regeneration of an optical field when a time delayed second pulse is provided from the coupling laser. It should be pointed out that all theoretical studies published recently on this subject also assume exact on resonance tuning by the probe laser as originally assumed in Ref.. Readers should consult a series pioneer studies on the on-resonance case given in . Consider a $`\mathrm{\Lambda }`$ system as depicted in Figure 1. We assume that a probe laser ($`E_p`$, frequency $`\omega _p`$) is tuned on or near resonance with the $`|1>|2>`$ transition, and a coupling laser ($`E_c`$, frequency $`\omega _c`$) is tuned so that exact two-photon resonance between states $`|1>`$ and $`|3>`$ is achieved. Assume atomic wave function of the form

$$|\mathrm{\Psi }(z,t_r)>=a_1e^{i\omega _1t}|1>+a_2e^{i\omega _2t}|2>+a_3e^{i\omega _3t}|3>,$$
(1)
we thus find the following atomic equations of motion

$`{\displaystyle \frac{a_1}{t_r}}`$ $`=`$ $`i\mathrm{\Omega }_pe^{ik_pz}a_2,`$ (2a)
$`{\displaystyle \frac{a_2}{t_r}}`$ $`=`$ $`i\mathrm{\Omega }_p^{}e^{ik_pz}a_1+i\mathrm{\Omega }_c^{}e^{ik_cz}a_3+i\left(\delta +i{\displaystyle \frac{\gamma _2}{2}}\right)a_2,`$ (2b)
$`{\displaystyle \frac{a_3}{t_r}}`$ $`=`$ $`i\mathrm{\Omega }_ce^{ik_cz}a_2,`$ (2c)
where $`a_j`$ and $`\gamma _j`$ are the $`jth`$ amplitude of the atomic wave function and decay rate, respectively. $`\mathrm{\Omega }_p`$ and $`\mathrm{\Omega }_c`$ are half-Rabi frequencies, given by $`\mathrm{\Omega }_p^{}=D_{21}E_p/(2\mathrm{})`$, $`\mathrm{\Omega }_c^{}=D_{23}E_c/(2\mathrm{})`$ with $`D_{ij}`$ being the dipole moment of the relevant transition. Also, $`\omega _p=\omega _2\omega _1+\delta `$$`\omega _c=\omega _2\omega _3`$, and the quantity $`t_r=tz/c`$ is the retarded time which is the combination of $`t`$ and $`z`$ that the laser field amplitudes depend on in vacuum. In order to correctly describe the propagation of the optical pulse, the atomic equations of motion (2) must be simultaneously solved self-consistently with Maxwell equations. In the limit of plane waves and slowly varying amplitudes, the positive frequency part of these fields satisfy

$`\left({\displaystyle \frac{\mathrm{\Omega }_p^{}}{z}}\right)_{t_r}`$ $`=`$ $`i\kappa _{12}\tau A_1^{}A_2,`$ (3a)
$`\left({\displaystyle \frac{\mathrm{\Omega }_c^{}}{z}}\right)_{t_r}`$ $`=`$ $`i\kappa _{32}\tau A_3^{}A_2,`$ (3b)
where the amplitudes $`A_1`$, $`A_2`$, and $`A_3`$ differ by position dependent phase factors from $`a_1`$, $`a_2`$, and $`a_3`$, and we have introduced notations $`\kappa _{12}=2\pi N\omega _p|D_{12}|^2/(\mathrm{}c)`$ and $`\kappa _{32}=2\pi N\omega _c|D_{32}|^2/(\mathrm{}c)`$ with $`N`$ being the concentration. In order to evaluate the polarization terms on the right hand sides of Eq.(3a,3b), it is necessary to obtain the amplitudes of the atomic wavefunction. The problem is intrinsically difficult because both Eq.(2) and (3) must be solved simultaneously. A great simplification can be achieved if one invokes the adiabatic approximation while evaluating the atomic wavefunction. The wavefunction obtained with the adiabatic limit will then be used to calculate the optical field in a self-consistent manner. This is the spirit of all adiabatic approximation based optical propagation problem and will be the method used in the present work.

In an adiabatic treatment we introduce $`A_1=a_1`$, $`A_2=e^{ik_pz}a_2`$, and $`A_3=e^{ik_cz}a_3`$ in order to remove the $`z`$ dependent phase factors. The adiabatic eigenvalues of the equations of motion can then be obtained from the characteristic equation given by

$$\lambda ^3\left(\delta +i\frac{\gamma _2}{2}\right)\lambda ^2(|\mathrm{\Omega }_p|^2+|\mathrm{\Omega }_c|^2)\lambda =0.$$

They are

$`\lambda _0`$ $`=0,`$ (4a)
$`\lambda _+`$ $`=`$ $`\left({\displaystyle \frac{\delta +i\frac{\gamma _2}{2}}{2}}\right)+\sqrt{|\mathrm{\Omega }_p|^2+|\mathrm{\Omega }_c|^2+\left({\displaystyle \frac{\delta +i\frac{\gamma _2}{2}}{2}}\right)^2},`$ (4b)
$`\lambda _{}`$ $`=`$ $`\left({\displaystyle \frac{\delta +i\frac{\gamma _2}{2}}{2}}\right)\sqrt{|\mathrm{\Omega }_p|^2+|\mathrm{\Omega }_c|^2+\left({\displaystyle \frac{\delta +i\frac{\gamma _2}{2}}{2}}\right)^2}.`$ (4c)
As can be seen, the inclusion of the non-vanishing one-photon detuning changes the eigenvalues. For sufficiently small $`\delta `$, this condition reduces correctly to that given in . The adiabatic condition requires that the eigenvalue $`\lambda _0\tau `$ always differ from $`\lambda _\pm \tau `$ by an amount that is large compared with unity. When $`|\delta \tau |>>|\mathrm{\Omega }_c|`$ and $`|\delta \tau |>>|\mathrm{\Omega }_p|`$, however, additional conditions are required. In fact, we require that the pulse length of the coupling laser be much longer than that of the probe laser, and $`|\mathrm{\Omega }_c\tau |^2/|\delta \tau |>>1`$ be satisfied. Therefore, during the entire pulse of the probe laser, the two-photon transition between $`|1>`$ and $`|3>`$ is shifted well outside the bandwidth of the probe laser. In another words, the overlap of the probe pulse band width with the resonance must be avoided. We now focus on the eigenvector corresponding to $`\lambda _0=0`$, since it is this adiabatic state that corresponds to the population all being in state $`|1>`$ before the laser pulses. In the lowest order of the adiabatic approximation Eqs.(2a) and (2c) imply that $`a_2=0`$. Eq.(2b) then yields

$$A_1=\frac{\mathrm{\Omega }_c^{}}{\mathrm{\Omega }_p^{}}A_3.$$
(5)
Since, for a closed system, $`|a_1|^2+|a_2|^2+|a_3|^2=1`$ we can write

$`A_1(z,t_r)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_c^{}(z,t_r)}{\sqrt{|\mathrm{\Omega }_c(z,t_r)|^2+|\mathrm{\Omega }_p(z,t_r)|^2}}},`$ (6a)
$`A_2(z,t_r)`$ $`=`$ $`0,`$ (6b)
$`A_3(z,t_r)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_p^{}(z,t_r)}{\sqrt{|\mathrm{\Omega }_c(z,t_r)|^2+|\mathrm{\Omega }_p(z,t_r)|^2}}}.`$ (6c)
It is important to realize that solutions (6a-6c) are not sufficiently accurate for many propagation problems. This is because that even though numerical solutions show that $`A_2`$ is much smaller than the amplitudes of the other states, it cannot be exactly zero without assuming the laser pulses are propagating at exactly the vacuum speed of light. In many situations, however, this assumption of the speed of laser pulses is not valid. This is particularly true in the case of sufficiently high optical density and very near a strong resonance excitation where significant modification of the propagation velocity is anticipated. In order to accurately account the propagation effect in the regime where significant dispersion effect is encountered, one must seek the next order correction to the lowest order approximation where the small quantities such as $`A_2`$ becomes critically important. To do this, one can use either Eq.(2a) or Eq.(2c) to derive a second approximation for $`A_2`$, thereby obtains the first order correction to the lowest order adiabatic approximation. Substitute Eq.(6a) into Eq.(2a) or Eq.(6c) into Eq.(2c) we obtain

$$A_2(z,t_r)=\frac{i}{\mathrm{\Omega }_p}\frac{A_1}{t_r}=\frac{i}{\mathrm{\Omega }_p}\frac{}{t_r}\left(\frac{\mathrm{\Omega }_c^{}}{\sqrt{|\mathrm{\Omega }_c|^2+|\mathrm{\Omega }_p|^2}}\right),$$
(7a)
and

$$A_2(z,t_r)=\frac{i}{\mathrm{\Omega }_c}\frac{A_3}{t_r}=\frac{i}{\mathrm{\Omega }_c}\frac{}{t_r}\left(\frac{\mathrm{\Omega }_p^{}}{\sqrt{|\mathrm{\Omega }_c|^2+|\mathrm{\Omega }_p|^2}}\right).$$
(7b)
Within the adiabatic approximation and the slowly varying amplitude approximation one can show that $`\mathrm{\Omega }_p(z,t_r)`$ and $`\mathrm{\Omega }_c(z,t_r)`$ can be chosen to be real. It is then trivial to show that two forms for $`A_2`$ as shown in Eqs.(7a,7b) are self consistent. It is interesting to note that even though having $`\delta `$ non-zero makes the criteria for the validity of the adiabatic approximation different, it does not change the solution for the state amplitudes. As usual with the adiabatic approximation used in a resonance situation, we must have $`\mathrm{\Omega }_c`$ already strong when $`\mathrm{\Omega }_p`$ starts to build up. This means that $`\mathrm{\Omega }_c\tau >>1`$, where $`\tau `$ is a measure of the time over which $`\mathrm{\Omega }_c`$, or $`\mathrm{\Omega }_p`$ change significantly. One scenario where the adiabatic approximation would hold through the whole laser pulse is if the two lasers peak at the same time, but the pulse length of the coupling laser is much longer. This feature of the pulse lengths and the requirement $`|\mathrm{\Omega }_c\tau |>>1`$ through the pulse length of the probe laser will make results based on the adiabatic approximation quite accurate. In essence, one must keep both $`|\lambda _+\tau |>>1`$ and $`|\lambda _{}\tau |>>1`$ throughout the pulse of the probe laser, e.g. there is a wide gap between $`\lambda _0=0`$ and the other eigenvalues so that curve crossing is avoided throughout the entire pulse length. We will now investigate the use of the adiabatic approximation in the propagation problem.

## III The propagation problem

In order to treat the propagation of the laser pulses through a resonant medium, the polarization of the medium as a function of position and time must be obtained. Using the atomic wave function obtained in Sec. II, we get

$`P(t)`$ $`=`$ $`N<\mathrm{\Psi }(z,t_r)|\widehat{D}|\mathrm{\Psi }(z,t)>`$ (8)
$`=`$ $`A_1^{}(z,t_r)A_2(z,t_r)e^{ik_pzi\omega _pt}ND_{12}+A_3^{}(z,t_r)A_2(z,t_r)e^{ik_czi\omega _ct}ND_{32}+c.c.`$
where $`N`$ is the concentration. For the positive frequency parts of the polarizations at the probe and coupling laser frequencies, we have

$`P_{\omega _p}^+`$ $`=`$ $`ND_{12}e^{ik_pi\omega _pt}A_1^{}A_2,`$ (9a)
$`P_{\omega _c}^+`$ $`=`$ $`ND_{32}e^{ik_ci\omega _ct}A_3^{}A_2.`$ (9b)
Using the first order correction to $`A_2`$, i.e. Eqs.(7a,7b), we find, for the amplitude of $`P_\omega ^+`$

$`P_{\omega _p,0}`$ $`=`$ $`iND_{12}{\displaystyle \frac{1}{\sqrt{|\mathrm{\Omega }_c|^2+|\mathrm{\Omega }_p|^2}}}{\displaystyle \frac{}{t_r}}\left({\displaystyle \frac{\mathrm{\Omega }_p^{}}{\sqrt{|\mathrm{\Omega }_c|^2+|\mathrm{\Omega }_p|^2}}}\right),`$ (10a)
$`P_{\omega _c,0}`$ $`=`$ $`iND_{32}{\displaystyle \frac{1}{\sqrt{|\mathrm{\Omega }_c|^2+|\mathrm{\Omega }_p|^2}}}{\displaystyle \frac{}{t_r}}\left({\displaystyle \frac{\mathrm{\Omega }_c^{}}{\sqrt{|\mathrm{\Omega }_c|^2+|\mathrm{\Omega }_p|^2}}}\right).`$ (10b)
Substitute these polarizations into Eq.(3a,3b), we obtain

$`\left({\displaystyle \frac{\mathrm{\Omega }_p^{}}{z}}\right)_{t_r}`$ $`=`$ $`{\displaystyle \frac{\kappa _{12}\tau }{\sqrt{|\tau \mathrm{\Omega }_c|^2+|\tau \mathrm{\Omega }_p|^2}}}{\displaystyle \frac{}{t_r/\tau }}\left({\displaystyle \frac{\tau \mathrm{\Omega }_p^{}}{\sqrt{|\tau \mathrm{\Omega }_c|^2+|\tau \mathrm{\Omega }_p|^2}}}\right),`$ (11a)
$`\left({\displaystyle \frac{\mathrm{\Omega }_c^{}}{z}}\right)_{t_r}`$ $`=`$ $`{\displaystyle \frac{\kappa _{32}\tau }{\sqrt{|\tau \mathrm{\Omega }_c|^2+|\tau \mathrm{\Omega }_p|^2}}}{\displaystyle \frac{}{t_r/\tau }}\left({\displaystyle \frac{\tau \mathrm{\Omega }_c^{}}{\sqrt{|\tau \mathrm{\Omega }_c|^2+|\tau \mathrm{\Omega }_p|^2}}}\right).`$ (11b)
Equations (11) are the key equations describing propagation of a pulse pair within the adiabatic approximation. By using the adiabatic approximation for the atomic response one needs to solve only two nonlinear partial differential equations, instead of solving five nonlinear partial differential equations Eq.(2a-2c) and Eq.(3a,3b). To further illustrate the effectiveness of the adiabatic approximation, we now examine two special cases. We proceed by first noticing that the quantity

$$F=\frac{|\mathrm{\Omega }_p|^2}{\kappa _{12}}+\frac{|\mathrm{\Omega }_c|^2}{\kappa _{32}},$$
(12)
represents the sum of the photon fluxes at $`\omega _p`$ and $`\omega _c`$ divided by the concentration of the medium through which the waves propagate. Differentiate $`F`$ with respect to $`z`$ while holding $`t_r`$ fixed and make use of Eq.(11), we obtain

$$\frac{}{z}\left(\frac{|\mathrm{\Omega }_p|^2}{\kappa _{12}}+\frac{|\mathrm{\Omega }_c|^2}{\kappa _{32}}\right)=0.$$
(13)
Thus, $`F`$ does not depends explicitly on $`z`$. This permits us to evaluate $`F`$ by evaluating it at $`z=0`$, the entrance to the atomic vapor cell, e.g.

$$F(z,t)=F(z=0,t)=\frac{|\mathrm{\Omega }_p(0,t)|^2}{\kappa _{12}}+\frac{|\mathrm{\Omega }_c(0,t)|^2}{\kappa _{32}}.$$
(14)
Thus, whenever $`|\mathrm{\Omega }_c(z,t_r)|^2/\kappa _{32}+|\mathrm{\Omega }_p(z,t_r)|^2/\kappa _{12}`$ occurs it can be replaced by $`F(t_r)`$, as determined in Eq.(14). An exact relation for this quantity can also be directly derived from Eq.(2a-2c) and Eq.(3a,3b). One obtains

$$\frac{}{z}\left(\frac{|\mathrm{\Omega }_p|^2}{\kappa _{12}}+\frac{|\mathrm{\Omega }_c|^2}{\kappa _{32}}\right)=\frac{|A_2|^2}{t_r}\gamma _2|A_2|^2,$$

This relation indicates that anytime $`|A_2|^2`$ is very small and slowly varying, and the decay rate of this state is not too fast, this total photon flux is close to depending only on $`t_r`$. The lack of dependency of this quantity on $`z`$ when the full set of equations is solved numerically is an important test for the validity of the adiabatic approximation.

### Case 1. $`\kappa _{12}=\kappa _{32}`$

In this situation the quantity $`|\mathrm{\Omega }_p(z,t_r)|^2+|\mathrm{\Omega }_c(z,t_r)|^2=\kappa _{12}F(t_r)`$. We now use this fact in Eqs.(11a,b) to obtain an analytical approximation that can be compared with numerical solutions to Eqs.(2a-2c) and (11a,b). Since $`\sqrt{|\mathrm{\Omega }_c|^2+|\mathrm{\Omega }_p|^2}`$ does not depend on $`z`$ when $`t_r`$ is held fixed, dividing both sides of Eqs.(11a,b) by this quantity yields

$`{\displaystyle \frac{}{z}}\left({\displaystyle \frac{\mathrm{\Omega }_p^{}\tau }{\sqrt{|\mathrm{\Omega }_c\tau |^2+|\mathrm{\Omega }_p\tau |^2}}}\right)`$ $`=`$ $`{\displaystyle \frac{\kappa _{12}\tau }{|\mathrm{\Omega }_c\tau |^2+|\mathrm{\Omega }_p\tau |^2}}{\displaystyle \frac{}{t_r/\tau }}\left({\displaystyle \frac{\mathrm{\Omega }_p^{}\tau }{\sqrt{|\mathrm{\Omega }_c\tau |^2+|\mathrm{\Omega }_p\tau |^2}}}\right),`$ (15a)
$`{\displaystyle \frac{}{z}}\left({\displaystyle \frac{\mathrm{\Omega }_c^{}\tau }{\sqrt{|\mathrm{\Omega }_c\tau |^2+|\mathrm{\Omega }_p\tau |^2}}}\right)`$ $`=`$ $`{\displaystyle \frac{\kappa _{12}\tau }{|\mathrm{\Omega }_c\tau |^2+|\mathrm{\Omega }_p\tau |^2}}{\displaystyle \frac{}{t_r/\tau }}\left({\displaystyle \frac{\mathrm{\Omega }_c^{}\tau }{\sqrt{|\mathrm{\Omega }_c\tau |^2+|\mathrm{\Omega }_p\tau |^2}}}\right).`$ (15b)
A few substitutions make an analytical solution to these equations relatively obvious. Let

$`W_p`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_p^{}\tau }{\sqrt{|\mathrm{\Omega }_c\tau |^2+|\mathrm{\Omega }_p\tau |^2}}},`$ (16a)
$`W_c`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_c^{}\tau }{\sqrt{|\mathrm{\Omega }_c\tau |^2+|\mathrm{\Omega }_p\tau |^2}}},`$ (16b)
$`v(t_r)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{t_r/\tau }}\left(|\mathrm{\Omega }_c(0,t_r^{^{}})\tau |^2+|\mathrm{\Omega }_p(0,t_r^{^{}})\tau |^2\right)d\left({\displaystyle \frac{t_r^{^{}}}{\tau }}\right),`$ (16c)
$`u(z)`$ $`=`$ $`{\displaystyle _0^z}\kappa _{12}\tau 𝑑z^{^{}}.`$ (16d)
In terms of these dimensionless quantities Eqs.(15) become

$`{\displaystyle \frac{W_p}{u}}+{\displaystyle \frac{W_p}{v}}`$ $`=`$ $`0,`$ (17a)
$`{\displaystyle \frac{W_c}{u}}+{\displaystyle \frac{W_c}{v}}`$ $`=`$ $`0,`$ (17b)
where general “travelling wave” type solutions are immediately obtained as

$`W_p`$ $`=`$ $`F_p(vu),`$ (18a)
$`W_c`$ $`=`$ $`F_c(vu).`$ (18b)
This special case provides a very useful test case for examining the validity of the adiabatic approximation for the “trapped” light problem. The functions $`F_p`$ and $`F_c`$ are easily determined by evaluating at $`z=0`$ (therefore, $`u=0`$) and using the fact that the laser fields are known as a function of $`t`$ and hence $`v`$. When the second coupling laser pulse is sent into the medium after a time delay, the predictions about the “revival” of the probe laser are all contained in this solution. The tabulation of $`F_p`$ remains the same as long as only a second coupling laser is sent into the medium after some delay time.

### Case 2. $`|\mathrm{\Omega }_p(0,t)\tau |<<|\mathrm{\Omega }_c(0,t)\tau |`$ with $`\mathrm{\Omega }_p\tau <<1`$

In this limit there can never be a significant population in state $`|3>`$. Therefore, to a good approximation we have $`\mathrm{\Omega }_c(z,t_r)=\mathrm{\Omega }_c(0,t_r)`$. That is, the coupling laser is not effected very much by the resonant medium and it propagates at speed $`c`$ without distortion. Thus, in this limit we only need to solve Eq.(11a) which becomes

$$\frac{}{z}\left(\frac{\mathrm{\Omega }_p^{}\tau }{|\mathrm{\Omega }_c(0,t_r)\tau |}\right)=\frac{\kappa _{12}\tau }{|\mathrm{\Omega }_c(0,t_r)\tau |^2}\frac{}{t_r/\tau }\left(\frac{\mathrm{\Omega }_p^{}\tau }{|\mathrm{\Omega }_c(0,t_r)\tau |}\right).$$
(19)
If we make the following dimensionless substitutions

$`W_p`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_p^{}}{|\mathrm{\Omega }_c(0,t_r)|}},`$ (20a)
$`v(t_r)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{t_r/\tau }}\left(|\mathrm{\Omega }_c(0,t_r^{^{}})\tau |^2\right)d\left({\displaystyle \frac{t_r^{^{}}}{\tau }}\right),`$ (20b)
$`u(z)`$ $`=`$ $`{\displaystyle _0^z}\kappa _{12}\tau 𝑑z^{^{}},`$ (20c)
we then find that the solution for $`W_p`$ is

$$W_p=F_p(vu),$$
(21)
with the function $`F_p`$ being determined by using $`z=0`$ at the entrance of the cell,

$$F_p(v(t))=\frac{\mathrm{\Omega }_p(0,t)}{|\mathrm{\Omega }_c(0,t)|},$$
(22)
where $`v(t)`$ is given in Eq.(20b), with $`t_r=t`$ at $`z=0`$. A table of $`F_p`$ as a function of its argument can be made, which will be used when $`z`$ is no longer zero and the argument is $`v(t_r)u(z)`$. Once we have tabulated $`F_p`$, the probe field at different $`z`$ can be trivially calculated using the known coupling field, e.g.

$$\mathrm{\Omega }_p^{}(z,t_r)=|\mathrm{\Omega }_c(0,t_r)|F_p(v(t_r)u(z)).$$
(23)
Again, the tabulation of $`F_p`$ is determined completely at all later times by Eq (22) - at least for the case where only a time delayed coupling laser pulse enters the material at $`z=0`$.

## IV Adiabatic approximation applied to optical field regeneration

In the precious sections, we have provided detailed descriptions and examples on how adiabatic approximation works. In this section, we will discuss various aspects of the adiabatic approximation, especially the validity in the context of probe pulse “revival”. These discussions provide many subtle understandings on the process which, in our view, are not discussed either in the recent works on “probe pulse revival” and the original studies of adiabatic approximation. This section, therefore, constitutes the main contributions of the present study. The validity of the adiabatic approximation depends on more than just the slow time dependence of the amplitudes of the probe and coupling lasers. For instance, the approximation will fail in dealing with the solution for the atomic response if the two pulses are sent simultaneously and the pulse widths and shapes are the same, even with both fields being turned on slowly. In this situation, at early times substantial populations will be promoted to state $`|2>`$, in contradiction to the adiabatic solution which, to the lowest order, predicts a zero population for this state. What is needed is to have the coupling laser built up and establish laser induced transparency before the probe laser starts to become intense. This could be done by either time delaying the probe laser relative to the coupling laser, or by making the coupling laser pulse length much longer than that of the probe laser and by having them peak in intensity at the same time. In the following treatment we will consider the latter situation. We therefore will choose the characteristics of $`\mathrm{\Omega }_p(0,t)`$ and $`\mathrm{\Omega }_c(0,t)`$ that would be conducive to making the adiabatic approximation work well, and demonstrate that this is indeed the case. We emphasize that our choice of the form of the laser fields is only one of many workable assumptions on the acceptable forms of the laser pulses that will make the adiabatic theory work well. Let us consider the functional forms for $`\mathrm{\Omega }_p(0,t)`$ and $`\mathrm{\Omega }_c(0,t)`$ as

$`\mathrm{\Omega }_p(0,t)`$ $`=`$ $`\mathrm{\Omega }_{p0}e^{(t/\tau )^2},`$ (24a)
$`\mathrm{\Omega }_c(0,t)`$ $`=`$ $`\mathrm{\Omega }_{c0}\left(e^{0.2(t/\tau )^2}+Re^{0.2(t/\tau x_0)^2}\right).`$ (24b)
With this functional form, the Rabi frequency due to the coupling laser rises up to produce laser induced transparency before the probe laser becomes large enough to have any effect, an important consequence of having a significantly longer pulse length for the coupling laser. Without this characteristic of the coupling laser, the probe laser would produce a population in state $`|2>`$ before the transparency was induced. The production of population in $`|2>`$ would mean that different adiabatic states have already been mixed and we have a failure of a strictly adiabatic picture in which only a single adiabatic state persists throughout the pulse. In Eq.(24), $`\mathrm{\Omega }_{p0}`$ and $`\mathrm{\Omega }_{c0}`$ are real constants characterising the peak amplitude of the two half-Rabi frequencies before the pulses enter the resonant medium. The parameter $`\tau `$ is a measure of the pulse length of the probe laser, $`R`$ is the ratio of the Rabi frequency at which the coupling laser recurs to its initial amplitude, $`x_0=t_d/\tau `$ is the value of $`t_r/\tau `$ at which the peak of the coupling laser recurs. We will begin by pointing out how the adiabatic approximation can be used to understand what happens for $`R=0`$. That is, for the case where there is no recurring coupling laser pulse. We take $`\mathrm{\Omega }_{p0}\tau =\mathrm{\Omega }_{c0}\tau =20`$, $`\kappa _{12}\tau =\kappa _{32}\tau =700cm^1`$, and $`\gamma _2\tau =0`$. We will soon see that in this numerical example the group velocity is sufficiently small so that the coupling laser dies away before the probe pulse can propagate through the vapor cell. When the coupling laser begins to die out, the intensity of the two lasers becomes proportional to each other. This is because when the coupling laser is beginning to become small, the following relation is appropriate (see Eq.(16c))

$$v(t_r)|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /2}+|\mathrm{\Omega }_{p0}\tau |^2\sqrt{\pi /2}.$$
(25)
That is, almost all of the area under the two Gaussian pulse shapes has already been included in the integration. Also, $`|\mathrm{\Omega }_p(0,t_r)|^2`$ has long been very small compared with $`|\mathrm{\Omega }_c(0,t_r)|^2`$. Thus, for $`t_r/\tau >3.0`$ we have, as a very good approximation

$$\mathrm{\Omega }_p^{}(z,t_r)\sqrt{|\mathrm{\Omega }_c(0,t_r)|^2+|\mathrm{\Omega }_p(0,t_r)|^2}F_p((|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /2}+|\mathrm{\Omega }_{p0}\tau |^2\sqrt{\pi /2})\kappa _{12}\tau z).$$
(26)
(In the case of an inhomogeneous medium $`\kappa _{12}\tau z`$ would be replaced by an integral over the $`z`$ dependent $`\kappa _{12}\tau `$.) At such a late time during the pulse, since the argument of $`F_p`$ depends only on the $`z`$ coordinate, therefore the time dependence of $`\mathrm{\Omega }_p^{}`$ is obviously exactly the same as that of $`\sqrt{|\mathrm{\Omega }_c(0,t_r)\tau |^2+|\mathrm{\Omega }_p(0,t_r)\tau |^2}`$. When $`|\mathrm{\Omega }_c(0,t_r)|`$ is several times larger than $`|\mathrm{\Omega }_p(0,t_r)|`$, as always is at such late times at $`z=0`$, this means that $`\mathrm{\Omega }_p^{}(z,t_r)`$ has the same time dependence as $`\mathrm{\Omega }_c^{}(0,t_r)`$ at such late times. The fact that the two pulses have the same time dependence as they begin to become very weak is the reason that populations are left in each of the states $`|1>`$ and $`|3>`$. If one of the laser fields died out much more rapidly than the other there would only be population left in only one state. ¿From Eq.(26) we see that if we choose $`z`$ such that

$$\kappa _{12}\tau z=\frac{1}{2}\left(|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /2}+\sqrt{\pi /2}|\mathrm{\Omega }_{p0}\tau |^2\right),$$

then the argument of $`F_p`$ will be the same as at $`z=0`$ and $`t_r=t=0`$. With this argument, the value of $`F_p`$ is $`F_p=\mathrm{\Omega }_{p0}/\sqrt{|\mathrm{\Omega }_{c0}|^2+|\mathrm{\Omega }_{p0}|^2}`$. This is the largest value that $`F_p`$ takes on. Thus, at this depth into the medium the value of $`|A_3(z,t_r)|`$ matches its largest value at $`z=0`$. This population persists at large $`t_r`$ until very slow collisional effects either destroy the coherence left behind in states $`|1>`$ and $`|3>`$, or until collisions mix the populations of the two states. This long persistence of a coherent mix of populations in states $`|1>`$ and $`|3>`$ is what leads to the “revival” of an optical field with frequency very close to that of the original probe laser when a coupling laser is re-injected into the medium after some time delay . The number of atoms left in $`|3>`$ is (within the adiabatic approximation) equal to the number of photons in the original probe pulse (recall that we have chosen parameters so that the probe pulse does not penetrate the entire medium). The creation of a coherent mixed stat is also what makes this situation interesting from the point of view of quantum computing. Next, we investigate the situation where $`R>0`$, i.e. a second coupling pulse enters the medium at a delayed time. Let $`\mathrm{\Omega }_{c0}\tau =20`$, $`\mathrm{\Omega }_{p0}\tau =5`$, $`\gamma _2\tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$. From the above consideration, we expect $`|A_3|`$ to have its largest value at large $`t_r`$ when $`z=2.86`$ cm. In Figure 2 we show a contour plot of $`|A_3(z,t_r)|`$ based on our adiabatic theory. The predicted $`|A_3|`$ is indicated by the line of constant color leading from $`t_r=0`$ and $`z=0`$ out to the horizontal path at $`z2.86`$ cm, as expected. In Figure 3 we show a surface plot of the probe field Rabi frequency $`\mathrm{\Omega }_p\tau `$ as functions of $`t_r/\tau `$ at $`z=3`$ cm. The long asymmetric tail of $`\mathrm{\Omega }_p`$ is a consequence of the behavior described in Eq.(26). In the region between $`2.5t_r/\tau 5`$ the ratio of the two half-Rabi frequencies is close to constant, averaging around 0.24. This ratio is also close to the adiabatic approximation for $`A_3`$ since $`|\mathrm{\Omega }_p|^2<<|\mathrm{\Omega }_c|^2`$. In a different set of parameters, we choose $`\mathrm{\Omega }_{c0}\tau =\mathrm{\Omega }_{p0}\tau =20`$, $`\gamma _2\tau =0`$, $`R=2.923`$, and $`\kappa _{12}\tau =\kappa _{32}\tau =700cm^1`$. In this example,our theory predicts that the peak value of $`|A_3(t_r,z)|`$ before the second coupling laser pulse enters the medium will occur at $`z=1.159`$ cm. Indeed, as can be seen clearly from Figure 4 that $`|A_3(z,t_r)|`$ takes on the value $`1/\sqrt{2}`$ at a depth of about 1.1-1.2 cm into the medium, as predicted. Figure 5 depicts the corresponding probe field predicted with adiabatic theory for this set of parameters. We now illustrate how one can take a set of initial conditions and very simply determine the characteristics of the “revived” probe pulse at the point it emerges from the atomic vapor.

First, note that with our choice of pulse characteristics and during the recurring coupling laser pulse

$`v(t_r)`$ $`=`$ $`S+{\displaystyle \frac{R^2}{2}}|\mathrm{\Omega }_{c0}\tau |^2\sqrt{{\displaystyle \frac{5\pi }{2}}}\left(1+erf\left(\sqrt{{\displaystyle \frac{2}{5}}}{\displaystyle \frac{(t_rt_d)}{\tau }}\right)\right),`$ (27a)
$`S`$ $`=`$ $`|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /2}+|\mathrm{\Omega }_{p0}\tau |^2\sqrt{\pi /2}.`$ (27b)
The reader will recall that $`F_p`$ was determined from the functional dependence of the probe and coupling laser pulses at $`z=0`$. In particular, $`W_p(v(t))=A_3(0,t)`$, where $`A_3(0,t)`$ is determined from Eq.(6c) and the pulse characteristics at $`z=0`$. In addition, the following properties of the function $`F_p`$ are very useful:

$`F_p(0)`$ $`=`$ $`0,`$ (28a)
$`F_p(S)`$ $`=`$ $`0,`$ (28b)
$`F_p(S/2)`$ $`=`$ $`F_p(v(0))={\displaystyle \frac{\mathrm{\Omega }_{p0}\tau }{\sqrt{|\mathrm{\Omega }_{c0}\tau |^2+|\mathrm{\Omega }_{p0}\tau |^2}}}.`$ (28c)
As will be seen later, we will invoke these properties when $`z=z_m`$ is reached, where $`z_m`$ is the value of $`z`$ at the end of the medium, to estimate the time interval over which the “revived” probe pulse exits, as well as the peak amplitude of $`\mathrm{\Omega }_p^{}(z_m,t_r)`$.

Before proceeding further, we will point out some obvious properties that must be satisfied by the parameter, $`R`$, that appears in Eq.(24b). We first note that $`R`$ must be large enough so that $`v(t_r)\kappa _{12}\tau z_m>0`$. That is, $`R`$ must be large enough so that the group velocity becomes large enough to allow most of the regenerated photons to exit the cell before the laser induced transparency ends. This group velocity depends on position and time as (see Eq.(15))

$$\frac{c}{v_g}=1+\frac{\kappa _{12}c\tau ^2}{|\mathrm{\Omega }_p(z,t_r)\tau |^2+|\mathrm{\Omega }_c(z,t_r)\tau |^2}.$$
(29)
We first write down the condition for a time $`t_{r1}`$ at which the “revived probe pulse” first reaches $`z_m`$. This is the earliest time at which the argument of $`F_p`$ goes from being negative and becomes zero. At this point $`F_p(0)=0`$, but at later times it will become non-zero. This time is determined by

$$v(t_{r1})\kappa _{12}\tau z_m=0.$$
(30)
At time $`t_{rm}`$ the value of $`F_p(S/2)`$ will be equal to $`A_3(0,0)`$. This is the maximum value $`F_p`$ takes on with the set of parameters chosen. This time is determined by

$$v(t_{rm})\kappa _{12}\tau z_m=S/2.$$
(31)
Finally, at the time $`t_{r2}`$ at which the “revived probe pulse” completes its exit from the medium. At this time, we have $`F_p(S)=0`$. Therefore,

$$v(t_{r2})\kappa _{12}\tau z_m=S.$$
(32)
Using Eqs.(27a,b) we can rewrite Eqs.(30-32) as

$`erf\left(\sqrt{2/5}{\displaystyle \frac{(t_{r1}t_d)}{\tau }}\right)`$ $`=`$ $`{\displaystyle \frac{2(\kappa _{12}\tau z_mS)}{R^2|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /2}}}1,`$ (33a)
$`erf\left(\sqrt{2/5}{\displaystyle \frac{(t_{rm}t_d)}{\tau }}\right)`$ $`=`$ $`{\displaystyle \frac{2(\kappa _{12}\tau z_mS/2)}{R^2|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /2}}}1.`$ (33b)
$`erf\left(\sqrt{2/5}{\displaystyle \frac{(t_{r2}t_d)}{\tau }}\right)`$ $`=`$ $`{\displaystyle \frac{2(\kappa _{12}\tau z_m)}{R^2|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /2}}}1,`$ (33c)
which can be further condensed into

$$erf\left(\sqrt{2/5}\frac{(t_rt_d)}{\tau }\right)=\frac{\alpha f\beta }{R^2}1,$$

where

$`\alpha `$ $`=`$ $`{\displaystyle \frac{2\kappa _{12}\tau z_m}{|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /2}}},`$
$`\beta `$ $`=`$ $`{\displaystyle \frac{2S}{|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /2}}},`$
with $`f=1`$ to determine $`t_{r1}`$, $`f=0.5`$ to determine $`t_{rm}`$, and $`f=0`$ to determine $`t_{r2}`$. Correspondingly, we obtain the peak value of $`\mathrm{\Omega }_p^{}(z_m,t_{rm})\tau `$ as

$$\mathrm{\Omega }_p^{}(z_m,t_{rm})=R\mathrm{\Omega }_{c0}e^{(t_{rm}t_d)^2/(5\tau ^2)}\frac{\mathrm{\Omega }_{p0}\tau }{\sqrt{|\mathrm{\Omega }_{c0}\tau |^2+|\mathrm{\Omega }_{p0}\tau |^2}}.$$
(34)
As an example, let $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$, $`z_m=8`$ cm, $`\mathrm{\Omega }_{c0}\tau =20`$, $`\mathrm{\Omega }_{p0}\tau =5`$, $`R=4`$, and $`\gamma _2\tau =0`$. With these parameters, we find $`\alpha =2.854598`$, $`\alpha /R^2=0.1784124`$, $`\beta =2.0559017`$, and $`\beta /R^2=0.128494`$. Thus, $`(t_{r1}t_d)/\tau =2.19`$, $`(t_{rm}t_d)/\tau =1.767`$, and $`(t_{r2}t_d)/\tau =1.50`$. The value of $`\mathrm{\Omega }_p^{}(z_m,t_{rm})`$ is thus estimated as

$$\mathrm{\Omega }_p^{}(z_m,t_{rm})=4(20)\mathrm{exp}(0.2(1.767)^2)\frac{5}{\sqrt{20^2+5^2}}=10.39.$$

The time $`t_{rm}`$ corresponds to the largest value of $`F_p`$, but not quite to the maximum of $`\mathrm{\Omega }_p^{}(z_m,t_r)`$. The actual maximum occurs closer to -1.74 with the maximum value of $`10.44`$. However, our simple estimate usually does very well if $`R`$ is large enough to easily let all of the regenerated probe laser photons escape. We next note that if $`R`$ is chosen to make the argument of $`F_p`$ exactly the same as its value at $`z=0`$ and $`t_r=t=0`$ at the time $`t_r=t_d`$, we then have

$$|\mathrm{\Omega }_{p0}\tau |^2\sqrt{\frac{\pi }{2}}+|\mathrm{\Omega }_{c0}\tau |^2\sqrt{\frac{5\pi }{2}}\left(1+\frac{R^2}{2}\right)\kappa _{12}\tau z_m=\frac{1}{2}\left(|\mathrm{\Omega }_{c0}\tau |^2\sqrt{\frac{5\pi }{2}}+|\mathrm{\Omega }_{p0}\tau |^2\sqrt{\frac{\pi }{2}}\right).$$
(35)
Thus,

$$|\mathrm{\Omega }_{p0}\tau |^2\sqrt{\pi /8}+|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /8}(1+R^2)=\kappa _{12}\tau z_m.$$
(36)
By choosing the value of $`R`$ such that this equation is satisfied, one achieves a situation where the peak value of $`\mathrm{\Omega }_p(z_m,t_d)`$ is given by (see Eq.(34))

$$\mathrm{\Omega }_{p0}^{}(z_m,t_d)=R|\mathrm{\Omega }_{c0}|\frac{\mathrm{\Omega }_{p0}\tau }{\sqrt{|\mathrm{\Omega }_{c0}\tau |^2+|\mathrm{\Omega }_{p0}\tau |^2}}.$$
(37)
The characteristics of this regenerated field can be estimated as follows. We first determine the value of $`t_r`$ such that the right hand side of Eq.(32) is zero or the full area under square of the incident probe and coupling lasers. This yields a range of time over which the probe rises from zero and returns to zero at the exit of the cell. One half of this time interval is close to the full-width at half-maximum of the exiting probe laser beam. We then note that our condition on $`R`$ yields

$$R^2=\frac{\kappa _{12}\tau z_m|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /8}|\mathrm{\Omega }_{p0}\tau |^2\sqrt{\pi /8}}{|\mathrm{\Omega }_{c0}\tau |^2\sqrt{5\pi /8}}.$$
(38)
Using these considerations we estimate the full-width at half-maximum pulse length, in the unit of the original probe pulse length $`\tau `$, to be

$$\mathrm{\Delta }_{1/2}=\frac{\sqrt{5\pi /2}}{R^2}\left(1+\frac{1}{\sqrt{5}}\left|\frac{\mathrm{\Omega }_{p0}}{\mathrm{\Omega }_{c0}}\right|^2\right).$$
(39)
If the original probe laser has only penetrated a small fraction of the thickness of the medium when the coupling laser pulse has passed by, then $`R`$ will turn out to be much larger than unity. This means that the width of the regenerated field is generally much smaller than the width of the initial probe pulse. Correspondingly, the bandwidth of the light will be larger. An immediate conclusion is thus available: the regenerated pulse is not the replica of the original probe pulse. It is coherent, to an extended degree, with the coupling laser because of the stimulated nature of the regeneration process. Therefore, the notion of “coherent storage” of the information carried in the original probe field is not accurate \[7-9\]. These features can be seen from Figures 3 and 5 which show $`\mathrm{\Omega }_p^{}`$ as a function of $`t_r/\tau `$ and $`z`$, as well as the time dependence of the regenerated field at the exit of the medium. In particular, as can be clearly seen from Figures 3 and 5, there is no optical field left in the region between the time when the original probe pulse and the companion coupling pulse enter the medium and the time when the second coupling pulse enters the medium. Therefore, photons are neither “stored” nor “stopped”. All photons from the original probe pulse are converted to the population of the state $`|3>`$. These figures also show that in “reviving the probe pulse” every regenerated photon comes at the expense of flipping population from state $`|3>`$ to state $`|1>`$. When the population of $`|3>`$ has been exhausted, there can be no further photon generated. Thus, in cases where the exiting regenerated field has a peak intensity much larger than the initial probe beam, the pulse width of this exiting pulse must be necessarily and correspondingly narrower. This is why the peak intensity of the regenerated field can go as $`R^2`$, while the width of the pulse tends to go as $`1/R^2`$.

## V Comparison between numerics and the adiabatic theory

We will now present some examples showing the degree of agreement between numerical calculations and the adiabatic theory developed in the previous sections. To obtain the numerical results, we solve numerically the coupled Maxwell and time dependent Schrödinger equations

$`\left({\displaystyle \frac{A_1}{t_r}}\right)_z`$ $`=`$ $`i\mathrm{\Omega }_pA_2,`$ (40a)
$`\left({\displaystyle \frac{A_2}{t_r}}\right)_z`$ $`=`$ $`i\mathrm{\Omega }_p^{}A_1+i\mathrm{\Omega }_c^{}A_3+i\left(\delta +i{\displaystyle \frac{\gamma _2}{2}}\right)A_2,`$ (40b)
$`\left({\displaystyle \frac{A_3}{t_r}}\right)_z`$ $`=`$ $`i\mathrm{\Omega }_cA_2,`$ (40c)
$`\left({\displaystyle \frac{\mathrm{\Omega }_p^{}}{z}}\right)_{t_r}`$ $`=`$ $`i\kappa _{12}A_1^{}A_2,`$ (41a)
$`\left({\displaystyle \frac{\mathrm{\Omega }_c^{}}{z}}\right)_{t_r}`$ $`=`$ $`i\kappa _{32}A_3^{}A_2.`$ (41b)
Here, $`A_1,A_2,`$ and $`A_3`$ are the same as in Eqs.(2a-2c) except for the removal of the $`z`$ dependent phase factors. First, recall the example that we have studied before where $`\mathrm{\Omega }_{c0}\tau =20`$, $`\mathrm{\Omega }_{p0}\tau =5`$, $`\gamma _2\tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$. The discussions based on the adiabatic approximation predicts that a population of $`|A_3|=1/\sqrt{17}=0.242536`$ at a depth of $`z=2.86`$ for $`t_r/\tau >6.0`$, as can be seen from Fig. 2. As a comparison, we show, in Figure 6, a contour plot of $`|A_3(z,t_r)|`$ for the same parameters given in Figure 2 but using numerical solution of Eqs.(40-41). The flat region in the middle, a region where both of the probe and coupling lasers have been extinguished, indicates that during this period of time the population in the two-photon state remains nearly constant, just as we have concluded before based on the adiabatic theory that a population is left behind in states $`|3>`$ and $`|1>`$ which persists for a long time. Correspondingly, we show, in Figure 7, a surface plot of the probe field resulted from the direct numerical solution of Eqs.(40-41) using the same parameters as in Figure 3. We note that there is no probe laser field between $`t_r/\tau =2`$ and the time of the injection of the second coupling laser pulse, the same period where the population in the two-photon state remains nearly constant. This indicates, as our adiabatic theory has shown, that photons are neither “stopped” or “stored” in the medium as claimed in Ref. \[7-9\]. As the original probe photons being slowed down, they are absorbed to produce the two-photon excitation with the above noted population in the state $`|3>`$. When all probe photons have been absorbed, there will be no further increase of the population in state $`3>`$. One may still argue that the probe photons are “stored” in the form of atomic coherence involving the two-photon excitation, and later when a second coupling pulse is injected into the medium, the probe pulse is “retrived” from this persist atomic coherence. This is, however, not correct. As we have shown the regenerated pulse has very different time-spectra characteristics from that of the original probe pulse. Therefore, the concept of “retriving the probe pulse” is incorrect. Indeed, the Fourier transform spectrum of all the regenerated fields have shown both different frequency content and noise characteristics that were not there in the original probe pulse. A simple argument based on our result (see Eqs.(37,39)) is sufficient to show that the regenerated pulse is very different from the original probe pulse, that is the regenerated pulse contains information in the second coupling pulse. Therefore, can not be the faithful replica of the original probe pulse alone.

Further comparison between the adiabatic theory and the numerical solution are presented in Figures 8 and 9 for the second set of parameters used to create Figures 4 and 5 which were obtained with adiabatic theory. Finally, we show a case where non-vanishing one-photon detuning is assumed. Figure 10 is surface plot of the probe field obtained from the adiabatic theory with $`\delta \tau =10`$, whereas Figure 11 is obtained by directly solving Eqs.(40-41) using the same parameters. All these comparisons show very good agreement between our theory and numerical solutions. The numerical calculations were also tested by making use of the fact that with $`\gamma _2\tau =0`$, the sum of the squares of the three state amplitudes should be unity. In all numerical examples described in this work, this condition of unity was preserved through at least seven significant figures if $`\gamma _2\tau =0`$. Also, in all cases when the number of $`t_r/\tau `$ and $`z`$ grid points were doubles for all figures shown, the values of the state amplitudes and half-Rabi frequencies repeated through five significant figures at all points where the functions were not near a zero.

## VI Case of far from resonance

In this Section, we examine the far off resonance excitation, a case has never been studied in literature with the same adiabatic treatment described in previous sections. Specifically, we investigate the condition for $`\lambda _0\tau `$ to stay well separated from $`\lambda _\pm \tau `$. We will show analytically that the formulism for far detuned excitation exhibits detuning independent feature when a specific condition is satisfied, a surprising result that to the best of our knowledge has never been reported in the literature. It is easy to see the condition for the eigenvalues to be well separated when $`|\delta |>>|\mathrm{\Omega }_c|`$ and $`|\delta |>>|\mathrm{\Omega }_p|`$. In this limit

$`\sqrt{|\mathrm{\Omega }_p|^2+|\mathrm{\Omega }_c|^2+\left({\displaystyle \frac{\delta +i\gamma _2/2}{2}}\right)^2}`$ $`{\displaystyle \frac{\delta +i\gamma _2/2}{2}}\left(1+2{\displaystyle \frac{|\mathrm{\Omega }_p|^2+|\mathrm{\Omega }_c|^2}{(\delta +i\gamma _2/2)^2}}\right)`$
$`=`$ $`\left({\displaystyle \frac{\delta +i\gamma _2/2}{2}}\right)+{\displaystyle \frac{|\mathrm{\Omega }_c|^2+|\mathrm{\Omega }_p|^2}{\delta +i\gamma _2/2}}.`$
The key element that could invalidate the adiabatic treatment is the eigenvalue

$$\lambda _{}=\frac{|\mathrm{\Omega }_c|^2+|\mathrm{\Omega }_p|^2}{\delta +i\gamma _2/2}\frac{|\mathrm{\Omega }_c|^2+|\mathrm{\Omega }_p|^2}{\delta },$$

which could become very close to the eigenvalue of the dark state $`lamda_0=0`$ at sufficiently large detuning. Since the length of the pulse due to the coupling field is chosen to be much longer than that due to the probe laser, this eigenvalue can only be well separated from the zero eigenvalue at times when the probe laser first starts to build up if $`|\lambda _{}\tau |>>1`$ at such times. Physically, what is required is for the ac Stark shift in level $`|3>`$ to satisfy

$$\frac{|\mathrm{\Omega }_c|^2}{\delta }\tau >>1.$$

This means that as the two-photon transition between $`|1>`$ and $`|3>`$ is driven, the final state, i.e. state $`|3>`$, must be shifted out of exact resonance by an amount that is much larger than the laser bandwidth. (We have assumed that the natural width of level $`|3>`$ is very small, being due to collisional relaxation in a very cold vapor. If the laser bandwidth is less than the width of the state $`|3>`$, the latter will replace $`\tau ^1`$ in the above equation.) In Figure 12, we first show a surface plot of $`\mathrm{\Omega }_p(z,t)`$ based on our adiabatic theory for zero detuning. Figure 13 is a corresponding plot where the full numerical evaluations of Eqs.(40-41) are carried out. As expected, the solution based on adiabatic theory agrees well with the numerical solution. We now compare Figures 12 and 13 with Figures 14 and 15 which are the surface plots of $`\mathrm{\Omega }_p(z,t)`$ with large one-photon detuning $`\delta \tau =120`$. Figure 14 is based on our adiabatic theory, therefore, should be compared with Figure 12, whereas Figure 15 is the result of full numerical calculation, thus should be compared with Figure 13. From these figures, we therefore conclude that when the appropriate adiabatic condition is met, the solution to the problem is insensitive to the detuning. Indeed, when the conditions given above are satisfied, the adiabatic approximation will be valid for all smaller $`\delta `$ as well, a surprising results that has not been reported in literature and could play an important role in the future experiments.

## VII Conclusion

In conclusion, we have shown in detail how the adiabatic approximation can be used to understand the propagation of a pair of optical pulses in a highly dispersive medium. Starting from the adiabatic solution of a three-level system interacts with two laser fields, we showed analytically that in the case where the one-photon detuning is small (i.e. $`0\delta <|\mathrm{\Omega }_c|`$), a coherent optical field with the frequency very close to that of the original probe can be regenerated when a coupling laser pulse in re-injected into the medium at a delayed time. We have provided detailed analysis on the conditions and characteristics of the regenerated field, including the estimate of the pulse width. Our analysis has shown that the regenerated pulse has a very different width therefore cannot be a perfect replica of the original probe pulse. Therefore, the notion of “coherent storage of the original probe pulse” is an incorrect description. As expected, the regenerated pulse is coherent, to an extended degree, with the coupling laser but not the original probe field. Subsequent numerical calculations have confirmed these analytical results. We have further developed the adiabatic theory to the cases where the one-photon detuning is larger enough so that $`|\mathrm{\Omega }_c|<\delta `$ and obtained a different type of adiabatic condition and solution. For time varying probe and coupling fields, the highly nonlinear wave equation for the recurring probe field is then investigated numerically. It is shown that even with this large one-photon detuning, as long as one chooses the pulse envelope and sequence properly so that there is a coherence left in the system, then there will be a coherent optical field generated when the coupling laser is turned on only.

We thank Professor R. Grobe for providing reprints.

## List of References

* D. Grischkowsky, M.M.T. Loy, and P.F. Liao, Phys. Rev. A12, 2514 (1975).
* J. Oreg, F.T. Hioe and J.H. Eberly, Phys. Rev A29,690 (1984); J.R. Kuklinski, U. Gaubatz, F.T. Hioe, and K. Bergmann, Phys. Rev. A 40 R6749 (1989).
* K.J. Boller, A. Imamoglu, and S.E. Harris, Phys. Rev. Lett. 66, 2593 (1991).
* R. Grobe and J.H. Eberly, Laser Physics 5, 542 (1995); J.H. Eberly, A. Rahman, and R. Grobe, Laser Physics 6, 69 (1996).
* J.R. Csesznegi and R. Grobe, Phys. Rev. Lett. 79, 3162 (1997); J.R. Csesznegi, B.K. Clark, and R. Grobe, Phys. Rev. A57, 4860 (1998); R. Grobe, Acta Physica Polica A93, 87 (1998).
* L.V. Hau et al., Nature (London) 397,594 (1999); M.M. Kash et al., Phys. Rev. Lett. 82, 5229 (1999); D.Budker et al., Phys. Rev. Lett. 83, 1767 (1999). For early works on the reduction of the group velocity of optical field in resonant media, see O.Schmidt, R. Wynands, Z. Hussein, and D. Meschede, Phys. Rev. A53, R27(1996); G. Muller, A. Wicht, R. Rinkleff, and K. Danzmann, Opt. Communi. 127, 37 (1996); A.M. Akulshin, S. Barreiro, and A. Lezama, Phys. Rev. Lett. 83, 4277 (1999).
* C. Liu et al., Nature (London) 409, 490 (2001).
* D.F. Phillips et al., Phys. Rev. Lett. 86, 783 (2001).
* E.A. Cornell, Nature (London), 409, 461 (2001).
* M. Fleischhauer, Opt. Express 4, 107 (1999); M. Fleischhauer, S.F. Yelin, and M.D. Lukin, Opt. Communi. 179, 395 (1999); M. Fleischhauer and M.D. Lukin, Phys. Rev. Lett. 84, 5094 (2000).
* Such a regeneration of an optical field is not unique to the process described in \[7-9\]. In fact, it is just a stimulated Raman generation in the context of electromagnetically induced transparency. It is not surprising at all that the frequency of the regenerated field is close to that of the original probe field because the second coupling pulse is identical to the first one (in frequency). The regenerated field, however, is not the faithful copy of the probe field as we have shown later. Therefore, the statement of “reviving the probe pulse” in Refs.\[7-9\] is inaccurate. In fact, one can mix the two lower states by pulsing a magnetic field, and then turning on a coupling laser to regenerate a very similar optical field.

## Figure Captions

* Energy level diagram showing relevant laser excitations. The decay rate of the state $`|3>`$ is assumed to be very small and hance neglected.
* Contour plot of $`A_3(z,t)`$ based on adiabatic theory. Parameters used: $`\mathrm{\Omega }_p\tau =5`$, $`\mathrm{\Omega }_c\tau =20`$, $`\gamma _2\tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$, $`R=4`$, $`t_d/\tau =11`$.
* Surface plot of $`\mathrm{\Omega }_p(z,t)`$ based on adiabatic theory. Parameters used: $`\mathrm{\Omega }_p\tau =5`$, $`\mathrm{\Omega }_c\tau =20`$, $`\gamma _2\tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$, $`R=4`$, $`t_d/\tau =11`$.
* Contour plot of $`A_3(z,t)`$ based on adiabatic theory. Parameters used: $`\mathrm{\Omega }_p\tau =\mathrm{\Omega }_c\tau =20`$, $`\gamma _2\tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =700cm^1`$, $`R=3`$, $`t_d/\tau =11`$.
* Surface plot of $`\mathrm{\Omega }_p(z,t)`$ based on adiabatic theory. Parameters used: $`\mathrm{\Omega }_p\tau =\mathrm{\Omega }_c\tau =20`$, $`\gamma _2\tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =700cm^1`$, $`R=3`$, $`t_d/\tau =11`$.
* Contour plot of $`A_3(z,t)`$ based on numerical solution of Eqs.(40-41). Parameters used: $`\mathrm{\Omega }_p\tau =5`$, $`\mathrm{\Omega }_c\tau =20`$, $`\gamma _2\tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$, $`R=4`$, $`t_d/\tau =11`$. This figure should be compared with the Figure 2.
* Surface plot of $`\mathrm{\Omega }_p(z,t)`$ based on numerical solution of Eqs.(40-41). Parameters used: $`\mathrm{\Omega }_p\tau =5`$, $`\mathrm{\Omega }_c\tau =20`$, $`\gamma _2\tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$, $`R=4`$, $`t_d/\tau =11`$. This figure should be compared to Figure 3.
* Contour plot of $`A_3(z,t)`$ based on numerical solution of Eqs.(40-41). Parameters used: $`\mathrm{\Omega }_p\tau =\mathrm{\Omega }_c\tau =20`$, $`\gamma _2\tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =700cm^1`$, $`R=3`$, $`t_d/\tau =11`$. This figure should be compared with Figure 4.
* Surface plot of $`\mathrm{\Omega }_p(z,t)`$ based on numerical solution of Eqs.(40-41). Parameters used: $`\mathrm{\Omega }_p\tau =\mathrm{\Omega }_c\tau =20`$, $`\gamma _2\tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =700cm^1`$, $`R=3`$, $`t_d/\tau =11`$. This figure should be compared with Figure 5.
* Surface plot of $`\mathrm{\Omega }_p(z,t)`$ based on adiabatic theory with non-vanishing one-photon detuning. Parameters used: $`\mathrm{\Omega }_p\tau =10`$, $`\mathrm{\Omega }_c\tau =20`$, $`\gamma _2\tau =0`$, $`\delta \tau =10`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$, $`R=4`$, $`t_d/\tau =11`$.
* Surface plot of $`\mathrm{\Omega }_p(z,t)`$ based on numerical solution of Eqs.(40-41) with non-vanishing one-photon detuning. Parameters used: $`\mathrm{\Omega }_p\tau =10`$, $`\mathrm{\Omega }_c\tau =20`$, $`\gamma _2\tau =0`$, $`\delta \tau =10`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$, $`R=4`$, $`t_d/\tau =11`$.
* Surface plot of $`\mathrm{\Omega }_p(z,t)`$ based on adiabatic theory with zero detuning. Parameters used: $`\mathrm{\Omega }_p\tau =10`$, $`\mathrm{\Omega }_c\tau =40`$, $`\gamma _2\tau =0`$, $`\delta \tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$, $`R=4`$, $`t_d/\tau =10`$.
* Surface plot of $`\mathrm{\Omega }_p(z,t)`$ based on numerical evaluation of Eq.s(40-41) with zero detuning. Parameters used: $`\mathrm{\Omega }_p\tau =10`$, $`\mathrm{\Omega }_c\tau =40`$, $`\gamma _2\tau =0`$, $`\delta \tau =0`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$, $`R=4`$, $`t_d/\tau =10`$. This figure should be compared with Figure 12.
* Surface plot of $`\mathrm{\Omega }_p(z,t)`$ based on adiabatic theory with large one-photon detuning. Parameters used: $`\mathrm{\Omega }_p\tau =10`$, $`\mathrm{\Omega }_c\tau =40`$, $`\gamma _2\tau =0`$, $`\delta \tau =1200`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$, $`R=4`$, $`t_d/\tau =10`$. This figure should be compared with Figure 12.
* Surface plot of $`\mathrm{\Omega }_p(z,t)`$ based on numerical evaluation of Eq.s(40-41) with large one-photon detuning. Parameters used: $`\mathrm{\Omega }_p\tau =10`$, $`\mathrm{\Omega }_c\tau =40`$, $`\gamma _2\tau =0`$, $`\delta \tau =120`$, $`\kappa _{12}\tau =\kappa _{32}\tau =200cm^1`$, $`R=4`$, $`t_d/\tau =10`$. This figure should be compared with Figure 13.


format: text

text:
# 1 Introduction

## 1 Introduction

If a quark-gluon system created in a heavy-ion collision at high energy undergoes a second-order phase transition to hadrons, one should expect large fluctuations in the hadron densities from region to region. Thus in any small interval of hadronization time the particles produced on the surface of the system, where the temperature $`T`$ is at the critical temperature $`T_c`$, should form patterns of clusters and voids of all sizes. The question is how such fluctuations can be detected. Since the detector integrates over the hadronization time, the hadron patterns formed at different times, when added, tend toward a uniform sum. In Refs.  we have studied the problem and found a way to overcome the difficulty. On the one hand, one must make severe cuts in the transverse momentum $`p_T`$ (accepting particles only in a narrow window of width $`\mathrm{\Delta }p_T`$), while on the other hand, appropriate measures that are insensitive to the randomization process following hadronization should be used to quantify the critical fluctuations. We proposed a void analysis that identifies the scaling properties of the fluctuations. The suggested measures in the analysis are in the two-dimensional space of rapidity $`(y)`$ and azimuthal angle $`(\varphi )`$. The analysis has been applied to the NA49 data of $`Pb`$-$`Pb`$ collisions at CERN-SPS , but remains to be considered for the data collected at BNL-RHIC.

In this paper we aim at simplifying the experimental problem of detecting the critical fluctuation, while making one of the theoretical parameters in the void analysis more experimentally relevant. They are accomplished by reducing the 2D space to 1D so that a void becomes a gap. The study of gaps in rapidity was undertaken by us earlier in connection with multiparticle production in hadronic collisions . The proposed gap analysis has been applied to the NA22 data  and provides a way of testing the predictions of the dynamical models on soft processes. We now apply the gap analysis to the phase transition problem in heavy-ion collisions. Besides showing what can be expected theoretically, we believe that the suggestions made here are more amenable to the experimental situation. An important outcome of this study is the identification of a numerical constant that characterizes the critical fluctuation. It can be checked by experiment.

Our theoretical tool, as in , is still the simulation of critical phenomenon on the 2D Ising lattice. Not only is the simulation simple to execute, it is also physically relevant, since the quark-hadron phase transition is very likely to be in the same universality class as the Ising system . We map the 2D Ising lattice to the surface of the quark-gluon plasma cylinder formed in a central collisions of heavy ions. We select a row in the lattice and map that to a strip on the plasma cylinder confined to a narrow interval in rapidity. Our proposed analysis is for the experimental examination of the fluctuation of gaps between particles in the azimuthal $`\varphi `$ variable.

## 2 Gap Analysis on the Lattice and for the Experiments

The relationship between the Ising problem and the hadron density problem has been discussed in detail in Refs. . We summarize here only the essentials for what is needed in this paper.

We start with a 2D lattice of size $`L^2`$, where $`L=288`$. On the lattice we define cells of size $`ϵ^2`$, where $`ϵ=4`$. The density $`\rho _i`$ at the ith cell is defined to be

$`\rho _i=\lambda c_i^2\theta (c_i),`$ (1)
where $`c_i`$ is the net spin at cell $`i`$ defined to be positive along the overall magnetization, i.e.,

$`c_i=\text{sgn}\left({\displaystyle \underset{jϵL^2}{}}\sigma _j\right){\displaystyle \underset{jϵi}{}}\sigma _j,`$ (2)
$`\sigma _j`$ being the lattice spin at site $`j`$, and $`\lambda `$ an unspecified factor relating $`c_i^2`$ to the particle density in $`y`$-$`\varphi `$ space. Near $`T_c`$, $`\rho _i`$ fluctuates from cell to cell; however, the existence of correlation of all length scales leads to clusters of all sizes formed out of neighboring cells having non-vanishing $`\rho _i`$.

In Refs.  we considered bins of size $`\delta ^2`$, with $`\delta `$ being integral multiples of $`ϵ`$. We regarded a bin to be empty if the average bin density, $`\overline{\rho }_b`$, is less than a threshold density $`\rho _0`$. We then defined a void to be a collection of contiguous empty bins, connected by at least one side between neighboring empty bins. The measure used to quantify the fluctuation of void sizes was shown to be insensitive to the value of $`\rho _0`$.

Without undercutting the importance of doing the void analysis on experimental data with high statistics, especially at high energies, we want now to consider in this paper the possibility of doing a similar analysis in 1D, since if clustering of hadrons occurs in 2D, they should also be visible in 1D. Not only will the statistics improve, there will also be no need to do binning to define voids, as we shall see.

Let us consider a row of $`C`$ cells, where we initially take $`C=L/ϵ=72`$. We map that row to a slice of the plasma cylinder in the central region with width $`\mathrm{\Delta }y`$ in rapidity. What $`\mathrm{\Delta }y`$ is will be discussed below. We also will limit $`p_T`$ to a narrow interval $`\mathrm{\Delta }p_T`$, as discussed extensively in Refs. . Thus we have only the azimuthal angle $`\varphi `$ as our space of variable. The row of $`C`$ cells is mapped to the interval $`0\varphi <2\pi `$.

We simulate $`N_e`$ configurations (where $`N_e=3\times 10^5`$) on the Ising lattice in 2D at $`T=T_c`$, and determine the density $`\rho _i`$ at each cell in each configuration. The simulation cannot be done in 1D because a 1D Ising system does not exhibit critical behavior. From a 2D configuration we choose an arbitrary row and examine the values of $`\rho _i`$ in each cell with $`1iC`$. The value of $`\lambda `$ in Eq. (1) is effectively set to 1 in what we do with $`\rho _i`$, since we always compare it to a threshold $`\rho _0`$ in units of $`\lambda `$. We define an occupied cell to be one in which

$`\rho _i>\rho _0,`$ (3)
and we place one (and only one) particle at cell $`i`$. If $`\rho _i\rho _0`$, then it is an empty cell. In mapping to the $`\varphi `$ variable, we place a particle at $`\varphi _i`$ randomly in the interval $`(i1)2\pi /C\varphi <i2\pi /C`$, when the ith cell satisfies Eq. (3). Thus after the whole row is mapped to the $`\varphi `$ variable, there are $`N`$ particles in the $`\varphi `$ space, where $`N`$ fluctuates from event to event. An event corresponds to a configuration on the lattice.

On the lattice the average number of occupied cells $`N`$ in a row, after averaging over all configurations, depends on the threshold $`\rho _0`$. In the experiment the average number of particles $`N`$ in the $`\varphi `$ space depends on $`\mathrm{\Delta }y\mathrm{\Delta }p_T`$. In our analysis we shall vary $`\rho _0`$. The predicted dependence of our measure on $`\rho _0`$ is to be checked by experiment by varying $`\mathrm{\Delta }y\mathrm{\Delta }p_T`$. In that sense the variable $`\rho _0`$ is not a theoretical parameter devoid of experimental relevance.

In an experiment of heavy-ion collisions the single-particle inclusive distribution in $`\varphi `$, $`dn/d\varphi `$, may not be uniform, because the impact parameter is not always exactly zero (even if a centrality cut is made), or because the detector efficiency is not symmetric in $`\varphi `$, mostly likely both. It is then better to use the cumulative variable $`X`$ , defined in our case by

$`X(\varphi )={\displaystyle _0^\varphi }{\displaystyle \frac{dn}{d\varphi ^{}}}𝑑\varphi ^{}/{\displaystyle _0^{2\pi }}{\displaystyle \frac{dn}{d\varphi ^{}}}𝑑\varphi ^{}.`$ (4)
The range $`0\varphi 2\pi `$ is thus mapped to $`0X1`$, and the density of particles in $`X`$, $`dn/dX`$, is uniform. The distribution $`dn/d\varphi `$ is determined after many events, but the variable $`X`$ is used for exclusive distributions event-by-event. In a simulation $`dn/d\varphi `$ is uniform, if $`N_e`$ is large; we nevertheless use the $`X`$ variable between 0 and 1 just so that our proposed gap analysis corresponds to what is to be done with experimental data.

Consider an event with $`N`$ particles in our window, labeled by $`i=1,\mathrm{},N`$, located in the $`X`$ space at $`X_i`$, ordered in accordance to $`X_i<X_{i+1}`$. Now define the distance between neighboring particles by

$`x_i=X_{i+1}X_i,i=0,\mathrm{},N,`$ (5)
with $`X_0=0`$ and $`X_{N+1}=1`$ being the boundaries of the $`X`$ space. We call these $`x_i`$’s gaps, which clearly satisfy

$`{\displaystyle \underset{i=0}{\overset{N}{}}}x_i=1.`$ (6)
Let us define for each event the moments

$`G_q={\displaystyle \frac{1}{N+1}}{\displaystyle \underset{i=0}{\overset{N}{}}}x_i^q,`$ (7)
which, for positive $`q`$, puts more emphasis on the large gaps than the small ones. The set $`\{G_q\}`$ with $`2qQ`$, $`Q`$ being some number less than 10, say, will be our quantification of the event structure. The set fluctuates from event to event, especially at the critical point. Our method is to use $`\{G_q\}`$ as a basis to construct a measure that can reveal the critical behavior.

Among the various erraticity measures considered in Ref.  we choose $`S_q`$ for our consideration here, since it has the simplest behavior. Starting with the definition

$`s_q=G_q\mathrm{ln}G_q,`$ (8)
where the average $`\mathrm{}`$ is performed over all events, $`S_q`$ is defined by

$`S_q=s_q/s_q^{st}`$ (9)
where $`s_q^{st}`$ is as defined in Eq. (8) but with the statistical contribution to $`G_q`$ only. Thus the deviation of $`S_q`$ from 1 is a measure of how strong the dynamical fluctuations are relative to the statistical fluctuations. In our simulation we generate $`N`$ random numbers between 0 and 1 to populate the $`X`$ space for the calculation of $`G_q^{st}`$. In an experiment one can do something similar to generate a random event.

We shall study the dependence of $`S_q`$ on $`q`$ for various values of $`\rho _0`$, using configurations simulated on the Ising lattice. The same can be done with the experimental data from heavy-ion collisions. If hadrons are formed in the latter due to a second-order phase transition, then the experimental $`S_q`$ should reveal some features that are similar to our theoretical $`S_q`$ to be presented below.

## 3 Results on Critical Fluctuations

Using $`N_e=3\times 10^5`$ configurations simulated on the Ising lattice at $`T=T_c=2.315`$ in units of $`J/k_B`$ , where $`J`$ is the near neighbor coupling and $`k_B`$ the Boltzmann constant, we have determined $`S_q`$ for a range of $`\rho _0`$. For $`\rho _0=20`$ the dependence of $`S_q`$ on $`q`$ is shown in a log-log plot in Fig. 1. The solid line is a straight-line fit of the calculated result, giving strong evidence for the power-law behavior

$`S_qq^\alpha .`$ (10)
The value of the exponent is $`\alpha =1.67\pm 0.02`$. The fact that $`\mathrm{ln}S_q`$ deviates unambiguously from 0 implies that $`S_q`$ is a statistically significant measure of nontrivial dynamical fluctuation. The power-law behavior is not necessarily a consequence of the dynamics of critical phenomenon, since similar behavior has been found before in the case of hadronic collisions . However, the exponent $`\alpha `$ is indicative of critical fluctuations; it is an order of magnitude larger than the value of $`\alpha =0.156`$ obtained for $`pp`$ collisions at $`\sqrt{s}=20`$ GeV.

To have one exponent $`\alpha `$ for all $`q`$ is a very economical description of the critical fluctuations. We can then investigate the dependence of $`\alpha `$ on $`\rho _0`$. Before so doing, let us first find a replacement of $`\rho _0`$ by some quantity that is directly measurable. In the preceding section we have related $`\rho _0`$ to the intervals $`\mathrm{\Delta }y\mathrm{\Delta }p_T`$ as the variable under experimental control that can be used to tune the average multiplicity $`N`$ accepted in the $`\varphi `$ window. Thus $`N`$ is a quantity that is both experimentally observable and theoretically computable on the lattice. It is therefore a suitable replacement for $`\rho _0`$. Nevertheless, we prefer to use an even better one that is the average number of gaps $`M`$, which is also observable. On the lattice in accordance to our convention of counting gaps in Eqs. (5) and (6), we have simply $`M=N+1`$. However, it should be recalled that in our simulation we have adopted the rule that, when $`\rho _i`$ exceeds $`\rho _0`$, only one particle is placed in a cell, not more, no matter how high $`\rho _i`$ is. That procedure makes sense in view of our measure $`G_q`$ of the event structure, since a tightly packed cell with many particles in it would have very small gaps that make negligible contribution in Eq. (7). Experimentally, if there are particles whose momenta are indistinguishable, or nearly so, whether they are separately counted or not also makes no significant difference in the calculation of $`G_q`$. Thus to allow for such possibilities it is better to count the number of distinguishable gaps, rather than the number of particles.

On the lattice we have simulated the configurations for a range of $`\rho _0`$ from 20 to 200. For every value of $`\rho _0`$ we can calculate the average number of gaps, $`M`$. Figure 2 shows the dependence of $`M`$ on $`\rho _0`$ at $`T_c`$. The straight-line segments are just the linear interpolations between neighboring points of $`\rho _0`$ where the calculation has been made. With such a definitive relationship at hand, any quantity that depends on $`\rho _0`$ will in the following be shown as a function of $`M`$, so as to render it amenable to experimental verification.

For each of the higher values of $`\rho _0`$ examined, we have found power-law behavior of $`S_q`$, as in Fig. 1. Thus the exponent $`\alpha `$ can be determined in each case. In Fig. 3 we show the dependence of $`\alpha `$ on $`M`$. Remarkably, the dependence is very linear. If we parameterize it as

$`\alpha =\alpha _0+\xi M,`$ (11)
we obtain

$`\alpha _0=0.258,\xi =0.055.`$ (12)
The nature of critical fluctuation is now seen to be reduced to a simple formula, Eq. (11), when the moments of gaps are used to describe the event structure. If we further put emphasis on the property that is independent of $`M`$, then the slope $`\xi `$ in Eq. (12) emerges as a numerical output of the theory that relies on no numerical input. This is perhaps the most succinct characterization of the critical phenomenon, beside the critical exponents. The latter depend on the temperature of a critical system near $`T_c`$. In heavy-ion collisions $`T`$ is not directly observable, so the corresponding critical exponents cannot be measured experimentally (see, however, Ref. ). Here we have an index $`\xi `$ that is eminently measurable and is the only numerical constant that can be meaningfully associated with critical fluctuation.

So far our study has been done only at $`T=T_c`$. To see how the results change when $`T`$ deviates from $`T_c`$, we repeat the above analysis for a range of $`T`$. Figure 4 shows $`S_q`$ vs $`q`$ in log-log plot for $`2.27<T<2.80`$ in units of $`J/k_B`$. Note that linearity is quickly lost when $`T`$ goes below $`T_c=2.315`$. For $`T>T_c`$ the linearity persists for a limited range $`2.315<T<3.1`$, but the slope is reduced. For $`T>3.2`$ the dependence bends over at high $`q`$ and the linearity is lost. For the range of $`T`$ where $`\alpha `$ can be determined from linear fits, we show in Fig. 5 $`\alpha (T)`$ up to $`T=2.8`$. The sharp peak that occurs at $`T=T_c`$ is a remarkable manifestation of the critical behavior. For $`T`$ less than $`T_c`$, so many hadrons are produced that the gap distribution rapidly tends toward the statistical. For $`T`$ larger than $`T_c`$, fewer hadrons are produced, and it takes more $`TT_c`$ difference for the statistical fluctuation to dominate. If $`T`$ were experimentally controllable, the measurement of $`\alpha (T)`$ would be an excellent way to determine the critical temperature of the quark-gluon system. That being not the case in reality, we can only learn from Figs. 4 and 5 that, if the system hadronizes in a range of $`T`$ around $`T_c`$, the most significant portion of the contribution to $`S_q`$ would come from the immediate neighborhood of $`T=T_c`$, and that only the $`\alpha `$ value around $`T_c`$ is experimentally relevant.

The above analysis is done for $`\rho _0=20`$. We can, of course, repeat the analysis for other values of $`\rho _0`$. In each case we find the peak value of $`\alpha `$ at $`T_c`$. At each $`T`$ where $`\alpha `$ can be meaningfully determined, we can investigate its dependence on $`M`$, just as we have done in Fig. 3 at $`T_c`$. In each case a linearity is found that allows us to determine the slope index $`\xi `$ in Eq. (11). The result is shown in Fig. 6, where $`\xi (T)`$ also exhibits a peak at $`T_c`$. We can now conclude that if hadronization is to occur around $`T_c`$, we expect the measurable value of $`\xi `$ to be around 0.05.

Having established the properties of the observables as functions of $`T`$ due to the dynamics of the critical system, we now investigate the question of stability with respect to changes in the size of the detector window, which is kinematical. We have so far considered a lattice of size 72x72 cells, from which we choose one row of $`C=72`$ cells from each configuration. We now want to vary the length and width of the row. We shall do so by setting $`T`$ at $`T_c`$. First, we consider a row of $`C=54`$ cells from each configuration, and later a row of $`C=36`$ cells. In each case it does not matter whether the row is mapped to a correspondingly shorter range in $`\varphi `$, since $`dn/d\varphi `$ is converted to $`dn/dX`$ in the $`X`$ space before the gap analysis is performed. However, the average number of gaps $`M`$ does change, so the simulation with shorter rows does correspond to a shorter ranges of $`\varphi `$ in the experiment. In Fig. 7 we show the result on $`\alpha `$ vs. $`M`$ for the $`C=72,54`$, and $`36`$ cells. The straight lines are linear fits, which are evidently very good. The values of the slope $`\xi `$ are :

$`\xi =0.055\pm 0.005,`$ $`C=72`$
$`0.057\pm 0.005,`$ $`C=54`$
$`0.048\pm 0.007,`$ $`C=36`$ (13)
The case of $`C=36`$ cells yields a slightly lower value and relatively larger errors, as is reasonable, since there are fewer points in a shorter range. Nevertheless, all three values of $`\xi `$ agree within errors. We therefore conclude that $`\xi `$ is stable against variations in the length of the $`\varphi `$ window.

Next, we consider variations in the width of the $`\varphi `$ window. In our simulation so far we have used only rows with 1-cell width. In an experiment the width $`\mathrm{\Delta }y`$ of a window in $`\varphi `$ will have to be adjusted in order to vary $`M`$. To check the stability of our result against variations in the width of our strip on the Ising lattice, we consider combinations of two and three rows. Let $`k`$ be the row index, and $`i`$ the cell index in the row, as before. Thus the location of a cell on the lattice is now denoted by $`ki`$. We take $`r`$ adjacent rows and combine them by performing vertical average of $`r`$ cells at each horizontal position $`i`$ , i. e., we define

$`\rho _i={\displaystyle \frac{1}{r}}{\displaystyle \underset{k=1}{\overset{r}{}}}\rho _{ki}`$ (14)
and then proceed in the gap analysis using $`\rho _i`$ as in the preceding section. It is important that the $`r`$ rows be adjacent so that the short-range vertical correlation of the cells can influence the value of $`\rho _i`$ so as to reflect its dynamical content; otherwise, if the $`r`$ rows are randomly chosen, the averaging in Eq. (14) tend to render $`\rho _i`$ more statistical and its distribution in $`i`$ more uniform. In short, taking $`r`$ adjacent rows corresponds to widening $`\mathrm{\Delta }y`$ in the experiment. In our calculation we consider only the case $`r=2`$ and $`3`$. Note that we take the average in Eq. (14) and then vary $`\rho _0`$ to change $`M`$. We could have added $`\rho _{ki}`$ (without dividing by $`r`$) and not vary $`\rho _0`$; that is equivalent to averaging $`\rho _{ki}`$ and dividing $`\rho _0`$ by $`r`$. Our chosen procedure allows more continous change in $`\rho _0`$.

The result of our analysis for $`r=1,2`$ and $`3`$, and $`C=72`$, are shown in Fig. 8. Roughly speaking, the slopes are essentially the same within errors, so $`\xi `$ may be regarded as independent of $`r`$. More precisely, we find

$`\xi =0.055\pm 0.005,r=1,`$
$`\xi =0.044\pm 0.003,r=2,`$
$`\xi =0.045\pm 0.003,r=3.`$ (15)
The slight decrease in $`\xi `$, as $`r`$ increases, is not unreasonable, since averaging over $`r`$ rows tends to suppress dynamical fluctuations. In the limit $`r`$ becoming very large, clearly only statistical fluctuations remain. We may summarize Eqs. (13), (15) and Fig. 6 by the grand average

$`\xi =`$ $`0.05\pm 0.01.`$ (16)
If this can be verified by experiments, then one should be able to claim that a signature of critical transition has been observed. If the value of $`\xi `$ is not confirmed, yet the power-law behavior of Eq. (10) is shown to exist in the data, with or without the linear dependence of $`\alpha `$ on $`M`$ in Eq. (11), that would still be an exciting experimental finding, suggestive of dynamical fluctuations.

There is the usual question about final-state interaction and the dilution of the dynamical signature due to randomization. The answer depends on the type of measure for that signature. The issue has been addressed in Ref. , where the dependence on the number of steps of final-state scatterings is examined, in Ref.  where configuration mixing is considered, and in Ref.  where different options in simulating time evolution are investigated. What one learns from all those studies is that the measures considered are not strongly affected by the final-state randomization. Since the gap analysis is a derivative of the void analysis , the same conclusion follows here. The one reminder that we should emphasize is that the window $`\mathrm{\Delta }p_T`$ in the transverse momentum should be kept as small as possible to minimize the overlap of particles emitted at different times.

## 4 Conclusion

We have studied the problem of even-to-event fluctuations of the hadronic patterns in phase space in heavy-ion collisions in search for detectable signatures of second-order quark-hadron phase transition. We have reduced the complication of voids in two dimensions to the simpler problem of gaps in one dimension. Using the moments of gaps to construct an entropy-like measure $`S_q`$, we have found a power-law dependence on $`q`$ with an exponent $`\alpha `$. It is the dependence of $`\alpha `$ on the average number of gaps that yields the index $`\xi `$, which serves to characterize critical fluctuation. We have found the stability of $`\xi =0.05\pm 0.01`$ against variations in the length and width of the detector window in $`\varphi `$. When the temperature of the system is moved away from $`T_c`$, the power-law behavior of $`S_q`$ on $`q`$ persists in a narrow range of $`T`$ around $`T_c`$, and the values of $`\alpha `$ and $`\xi `$ show strong peaks at $`T_c`$. Thus $`\xi `$ is a measure of the critical behavior and is a number that arises out of the study of fluctuations without any numerical input. It is highly significant that $`\xi `$ can be checked by experiments, since all measures leading to its determination are designed to be observable.

In heavy-ion experiments it is not difficult to make cuts in $`\mathrm{\Delta }p_T`$ and $`\mathrm{\Delta }y`$ to limit the average number of particles in a narrow strip in $`\varphi `$ to the range from 10 to 40. If those particles are found to be non-uniformly distributed for every event, the gap analysis proposed here is a way to quantify those fluctuations. When the exponent $`\alpha `$ and the index $`\xi `$ are found to exist, there are numerous variables under the control of the experiments to vary. We can mention, for example, the position of the $`y`$ slice, the value of $`p_T`$, the centrality of collisions, the nuclei sizes and the c.m. energy. It would be very interesting to see the dependence of $`\xi `$ on the total transverse energy $`E_T`$, since it can provide us with some idea of when the critical behavior is lost as the nuclei overlap gets to be too small to create a quark-gluon plasma. Also, for nonzero impact parameter, even if critical fluctuations exist, the index $`\xi `$ may depend on which sector of the $`\varphi `$ space the analysis is performed, since there is no $`\varphi `$ invariance in elliptic flow. One can envision a rich variety of phenomenological studies once the exponent $`\alpha `$ and index $`\xi `$ are found in the experimental data. They can provide valuable information about the quark-gluon system. The application of the gap analysis to the data is therefore strongly urged.

## Acknowledgment

This work was supported, in part, by the U. S. Department of Energy under Grant No. DE-FG03-96ER40972, and the Natural Science and Engineering Research Council of Canada and the Fonds FCAR of the Quebec Government.

## Figure Captions

$`S_q`$ vs $`q`$ in log-log plot for $`T=T_c`$ and $`\rho _0=20`$.

The average number of gaps $`M`$ vs $`\rho _0`$ at $`T_c`$.

The exponent $`\alpha `$ vs $`M`$ at $`T_c`$.

$`S_q`$ vs $`q`$ in log-log plot for a range of $`T`$ at $`\rho _0=20.`$

$`\alpha `$ vs $`T`$ at $`\rho _0=20`$.

$`\xi `$ vs $`T`$ at $`\rho _0=20`$.

$`\alpha `$ vs $`M`$ at $`T_c`$ for various number of cells in a row.

$`\alpha `$ vs $`M`$ at $`T_c`$ after averaging over various number of rows.