## Abstract

We present a one-dimensional iterative predictor-corrector finite-difference time-domain method for modeling of broadband optical pulse propagation and interaction with inhomogeneously broadened materials. The simulator is used to demonstrate two- and three-pulse photon echoes resulting from bandwidth limited pulse and matched chirp interactions with a material modeled with hundreds of equally spaced, discrete spectral lines of detuning. The results are illustrated as Bloch-sphere evolution movies.

©2005 Optical Society of America

## 1. Introduction

Nonlinear optical processes have potential applications in lasers, optical memory, fiber optics, and optical signal processing (OSP). Optical coherent transients (OCT) nonlinear processes involve pulses interacting with two-level systems (TLS), that is, media whose atomic state can be varied between two distinct energy levels by the application of an optical field at the appropriate resonance frequency. In solid or crystalline hosts, due to fluctuations of the local environment surrounding each TLS, the homogeneous resonance line of each individual two-level system can shift, with the effect that the absorbtion lineshape is said to be inhomogeneously broadened, since different subsets of atoms will have different resonance frequencies. Another mechanism that can cause inhomogeneous broadening is the Doppler shift due to atomic motion in gasses.

In order to understand the behaviour of inhomogeneously broadened absorbers (IBA) better, several simulation schemes of varying complexity have been developed, with the newest generation of simulations allowing for phase modulated signals interacting at multiple different angles of propagation [1]. One possible shortcoming of most simulation methods is the use of the rotating wave approximation (RWA), in which terms oscillating at twice the optical frequency are omitted. We have developed a simulator for IBA materials which avoids using the RWA and could thus be used to compare high bandwidth simulations based on other techniques to verify their validity. The simulator is based on an extension of Ziolkowski’s iterative predictor-corrector finite difference time-domain (IPC-FDTD) scheme for simulating a single TLS interaction with a short pulse optical field [2], but has been extended to calculate the Maxwell-Bloch system for hundreds of lines of detuning simultaneously in order to approximate an inhomogeneously broadened lineshape continuum. This approach does not make the RWA, and thus incorporates the physical effects of the higher frequency terms. Of particular interest is the interaction of a three-pulse sequence with the IBAs, since this results in a three pulse photon echo, whose behaviour is governed (for small pulses) by a causal correlation-convolution integral, but can be evaluated with the techniques presented here for arbitrarily large pulses [3]. The three-pulse echo can be used to directly correlate high bandwidth optical signals, or can be used to perform a number of OSP tasks [4, 5, 6, 7, 11]. In this paper we will outline the Bloch system for an inhomogeneously broadened absorber, the iterative predictor corrector FDTD simulation and visualization techniques, and present simulation results for a standard two- and three-pulse photon echo, as well as chirped two- and three-pulse echos.

## 2. Maxwell-Bloch system

The Maxwell-Bloch system for IBAs consists of Maxwell’s equations for the magnetic and electric field propagation, with the appropriate terms for coupling to the Bloch equations, which govern the two-level system. The analytic form of the 1-D Maxwell’s equations is given below, where we have assumed a z-propagating electric field polarized along the x-dimension (*E*_{x}
), along with a y-polarized magnetic field (*H*_{y}
).

$$=-\frac{1}{{\epsilon}_{0}}{\partial}_{z}{H}_{y}-\sum _{j}\frac{{N}_{j}\gamma}{{\epsilon}_{0}{T}_{2}}{\rho}_{1}^{j}\left(z\right)+\frac{{N}_{j}\gamma {\omega}_{j}}{{\epsilon}_{0}}{\rho}_{2}^{j}\left(z\right)$$

In these equations, *∂* represents a partial derivative with respect to the subscripted variable, *µ0* and *ε*_{0}
are the permeability and permittivity of the host material, repectively, and the *P*_{x}
term indicates that a polarization in the material is coupling back into the electric field. *N*_{j}
is the density of TLS at the *j*th line center (so that *N*_{atoms}
=∑
_{j}*N*_{j}
), *γ* is the dipole matrix element coupling coefficient, *ω*_{j}
is the resonance frequency of the *j*th transition and *T*_{2}
is the coherence lifetime. The polarization term is due to the TLS that is provided by the IBA, and is governed by the Bloch equations for each resonance frequency *ω*_{j}
:

The Bloch equations consist of three coupled evolution equations for *ρ*
_{1,j},*ρ*
_{2,j}, and *ρ*
_{3,j}, which are also functions of space. *ρ*
_{1,j} and *ρ*
_{2, j} represent the in-phase (dispersion) and quadrature (absorbtion) components of the coherence between the E-field and the TLS wave function, while *ρ*
_{3,j} represents the TLS inversion of a subset of atoms at a particular resonance frequency. In other words, *ρ*
_{1,j} and *ρ*
_{2,j} provide a measure of the phase relationship between the upper and lower atomic state coherence of the *j*th detuned absorbtion line and the driving optical field. Also, *ρ*
_{1,j} and *ρ*
_{2,j} are the terms which couples energy back and forth between the two-level system and the electric field. *ρ*
_{3,j} provides a measure of what fraction of the atomic ensemble is in the upper/lower level. Without decoherence and the limited excited state lifetime, ${\rho}_{1}^{2}$+${\rho}_{2}^{2}$+${\rho}_{3}^{2}$=1 for each *j*, resulting in the famous Bloch sphere, in which the atomic TLS dynamics are governed by the evolution of the Bloch vector on the surface of a unit sphere. The north pole (*ρ*
_{3}=1) represents full inversion, while the south pole (*ρ*
_{3}=-1) represents the ground state. The rotating dynamics in the equitorial plane represent the relative phase of the ground and excited state, and causes the Bloch vectors to spin around the polar axis at an angular rate of *ω*_{j}
. Decoherence of the atomic system can be envisioned in the Bloch sphere picture as an exponential shrinking of the *ρ*
_{1,2} plane with a time constant of *T*
_{2} (coherence time), while the excited state lifetime manifests itself as an exponential term which drives *ρ*
_{3} towards -1 (the ground state, or -*ρ*30 for a non-zero equillibrium inversion) with a time constant of *T*
_{1}, the excited state lifetime.

In previous work, Ziolkowski thoroughly outlined a derivation of a numerical Maxwell-Bloch solver for a 1D, single homogeneous resonance line using an IPC-FDTD code, which can be duplicated with our simulation by letting *ω*_{j}*=ω*_{0}
[2]. We will start with Ziolkowski’s numerical equations, and outline how the system can be modified to incorporate multiple inhomogeneous lines, as indicated in the analytic equations above. It should be noted that a few typographical errors in the original numerical equations were corrected in the equations below. The discretized numerical equations for the typical FDTD leap-frog algorithm, in which the z-propagating magnetic field (*H*_{y}
) and electric field (*E*_{x}
) are spatially separated by Δ*z*/2 and temporally separated by Δ*t*/2, are shown in Eq. (6–7), in which *m* and *n* denote the spatial and temporal discretization steps, respectively. The Bloch-vector values describing the in-phase and quadrature components of the polarization, as well as the population of a particular resonance line, are given by Eq. (8–10) as *u*
_{1}, *u*
_{2}, and *u*
_{3}, respectively. These values are numerically constrained to ${u}_{1}^{2}$+${u}_{2}^{2}$+${u}_{3}^{2}$=1, and will later be multiplied by the corresponding exponentials to transform them to *ρ*1,*ρ*2, and *ρ*3, thus simulating decoherence and excited state lifetime decay.

*Maxwell equations*

*Bloch equations*

In these equations, Δ*z* is the spatial increment, and should typically not exceed $\frac{\lambda}{100}$ for stable, dispersion-free numerics of a nonlinear problem, thus also locking the value of Δ*t*, since the Courant condition requires that Δ*t*<Δ*z/c*. For more robust FDTD simulations, Δ*t* is commonly made smaller, and in our case Δ*t*=Δ*z*/2*c* was chosen. Other constants are the permeability and permittivity of the material host, *µ*
_{0} and *ε*
_{0}, while coefficients used in the equations above are given by

In our simulation, we chose some illustrative values for the parameters which simplify and mitigate the enormous computational burden, and help illustrate the physics as efficiently as we could, even though the chosen values do not precisely model any particular physical system. In particular, *N*_{atom}
=10^{24}m^{-3} is the number of two-level atoms involved in the interaction at a particular frequency bin, *γ*=10^{-29}
*Cm* is the dipole coupling coefficient, *ω*
_{0}=2*π f*_{0}*, f*_{0}
=2.0×10^{14}s^{-1} is the resonance frequency of a two-level system, *h̄* is Planck’s constant, and *T*
_{1}=10^{-10}s and *T*
_{2}=3×10^{-12}
*s* are the excited state lifetime and coherence times, respectively. Finally, *ρ*30 gives the initial ground state population, and is usually set to -1, indicating that the system is in the ground state before the pulse interaction occurs. It should be noted that the $n+\frac{1}{2}$ superscript indicating a time dependence of the variables A, B, *C*+, *C*-, and D has been omitted for brevity. Values used in the simulations are identical to those in Ziolkowski’s paper, unless stated otherwise, and are appropriate for low coherence time systems such as organics or semiconductors. Since the electric field equations (Eq. (7)) and the Bloch equations (Eq. (8–10)) are coupled, they can not be integrated in the typical FDTD leap-frog scheme. These equations do, however, lend themselves to be integrated via a predictor-corrector (PC) scheme [8]. In this method of numerical integration, the unknown future value of a certain variable on the RHS of the equation is predicted, and the LHS is calculated. The calculated value is then compared to the predicted value, and if the two differ by a more than a maximum allowable error, the process is repeated using the calculated value from the previous PC iteration as the predicted value on the RHS. This is represented in Eq. (16),

where *X* represents the variables of interest (*E*_{x}
, *u*
_{1}, *u*
_{2} and *u*
_{3}), and *f* is the final bracketed term in Eq. (8–10) and Eq. (7) or (17) (depending on whether the simulated system is homogeneously or inhomogeneously broadened). This process is repeated for all four variables at each spatial and temporal step until the coupled equations yield acceptable results. It was found that 3–4 PC steps result in differences of predicted and calculated values that were on the order of 0.001%. Once a predictor-corrector step is complete, the dephasing effect due to the limited coherence time is implemented by multiplying *u*
_{1} and *u*
_{2} by a decaying exponential so that ${\rho}_{1,2}={u}_{1,2}{e}^{-n\Delta t\u2044{T}_{2}}$. Similarly, a slowly decaying inversion of the form ${\rho}_{3}=\left({u}_{3}+1\right){e}^{-n\Delta t\u2044{T}_{1}}-1$ is generated, which exhibits the effect of slow decay into the ground state. In the simulations shown later, *ρ*
_{1,2,3} are plotted for each *j*, and the coherence-decay is easily visible. The lifetime effect is less pronounced, since the lifetime of the material is hundereds of times longer than the simulation time. The code was tested with only one material resonance, and pulses of areas *π* and 2*π* verified proper operation for an inversting pulse and self-induced transparency, reproducing Ziolkowski’s results and the well-known analytic solutions [9].

## 3. Inhomogeneous media

After verifying the model’s response with a single homogeneous absorbtion line, and comparing the results with Ziolkowski’s, the model was expanded to handle an inhomogeneously broadened absorber. Typically, inhomogeneous broadening is a result of Doppler shift in gaseous media (0.1–1 GHz), or the result of variations in the local crystal environment in rare-earth ion-doped crystal lattices (10–100 GHz), or the result of amorphous environment fluctuations in glasses, polymers, or liquids (1–10 THz). The result is that the homogeneous absorption lines are detuned from some central wavelength, as shown in Fig. 1. Since the simulation will be used with either Fourier-limited pulses with only a few optical cycles (corresponding to a pulse duration of as little as 35 fsec) or chirped pulses spanning the same spectral range (but lasting up to 450 fsec), the spectral bandwidth (or inhomogeneous linewidth) of the material should be at least Γ_{IN}=30 THz. This is the absolute lowest limit of the spectral bandwidth of the material, and we found that extending the spectral domain by a factor of 3–4 gives superior results, since the tailing ends of the pulse’s spectra are thus also adequately sampled. For this reason, we chose a material bandwidth of 120 THz. This choice is not physically realistic, but is used to better illustrate the dynamics on the Bloch sphere and shorten the simulation time, since higher bandwidth can accomodate shorter pulses.

At the optical carrier frequency of 200 THz, the 120 THz bandwidth corresponds to a fractional bandwidth of 0.6, and in both the short- and chirped-pulse case, 301 equally spaced spectral lines were used to cover the inhomogeneous bandwidth, corresponding to about one modeled line per homogeneous linewidth (1/*T*
_{2}=0.3 THz). The modeled material has a rectangular inhomogeneous absorption profile, which was achieved numerically by giving all the *N*_{j}
terms in Eq. (11 & 12) the same value, even though a different lineshape (such as a Gaussian) could easily be created by using the appropriate *N*_{j}
weighting terms. The modified rate equation for the electric field is given below, while the rate equations for *u*
_{1}, *u*
_{2} and *u*
_{3} are practically unchanged, except for the resonance frequency *ω*_{j}
(previously *ω*_{0}
), which now varies for different values of detuning.

$$+{B}_{j}\frac{1}{2}\left[{u}_{2,j}{\mid}_{m}^{n+1}+{u}_{2,j}{\mid}_{m}^{n}\right]\left\}\right]$$

Different inhomogeneous absorbtion profiles could easily be implemented by applying the appropriate weighting to each of the homogeneous, detuned lines via *N*_{j}
. In the simulations results shown below, further items of interest are the material lifetime *T*
_{1}=100 ps, the coherence time *T*
_{2}=3 ps, the length of the material L=20*µ*m, and *αL*=0.015 (by comparing the electric field before and after it had propagated through the imulated material). The electrical field of the driving term was injected at one end of the simulation grid, which in turn contained 1,200 spatial points (with 301 spectral detunings each), and was terminated with a perfectly matched layer absorbing boundary condition on the other. The total simulated times were 1.2 psec for the brief-pulse simulation and 2.5 psec for the chirped pulse case, and required about a day on a workstation. The spatial and temporal lengths of the simulations had to be kept small due to the computational load that the IPC-FDTD presents for this inhomogeneous case with 301 absorbtion lines. The sampling requirements for FDTD techniques at optical frequencies drove even the small 1-D simulation space to more than 100,000 time-steps and more than 1,000 spatial steps, at each of which 301 lines of detuning have to be calculated approximately 4 times to obtain good IPC-convergence. The computational load is probably this method’s biggest shortcomming, and it will most likely prohibit extension of this technique to two or three spatial dimensions. However, in the 1-D case, the simulations are still manageable, and could offer a reference to other simulation techniques which perform the rotating wave approximation, and operate only on the electric field envelope. A second concern that arises in the use of this IPC-FDTD method, which we have not yet investigated, is the propagation of error. Since the stringent temporal sampling requirement draws with it a need to perform several orders of magnitudes more calculations than would be necessary to simply calculate the envelope of the electric field, the integrated numerical error may in fact prove to outweigh any small effects that the non-RWA simulation technique could provide. The final problem with this technique is that the coherence decay starts at the beginning of the simulated time-window, not at the occurance of the applied pulses. This happens because Eq. (11–15) start decaying immediately, and will result in slightly different solutions for echo-response depending on when the pulse sequence occurs. This error should be minimal as long as the pulse sequence occurs soon after the simulation starts. Since our first pulses start within 0.1 ps of the start of the simulation, we expect the resulting error to be a shortening of the *u*
_{1} and *u*
_{2} vectors to ${e}^{-\frac{0.1}{3}}\approx 97\phantom{\rule{.2em}{0ex}}\%$ of their intended length. This effect can be minimized by either using longer coherence times (*T*
_{2}) or short programming pulses. Additionally, since the phase of the vectors is the governing aspect as to the timing of the echo pulses, a slight length error can easily be tolerated.

Our primary goal in performing these simulations was to produce two-pulse and three-pulse photon echoes using bandwidth limited pulses. The second goal was to visually illustrate the simultaneous Bloch-vector evolution for a slice of the simulated material. Both of these items will be discussed in more detail.

## 4. Photon echo using brief pulses

In this section we present the results for both the two-pulse (2P) and three-pulse (3P) photon echo simulations using 35 fs bandwidth-limited hyperbolic-secant pulses. The interaction of two optical pulses which are seperated by a time *τ*
_{12} will write a spectral grating in an inhomogeneously broadened material [7], assuming *τ*
_{12} is less than the coherence time of the material. These two pulses also result in the emission of an optical pulse by the material, delayed after the second pulse by *τ*
_{12}, known as the 2-pulse photon echo. This process is commonly explained by two different methods.

The first method evolves all of the inhomogeneously detuned lines on the same Bloch sphere at their particular phase evolution rates. After the first pulse hits the material, the Bloch vectors of all the detunings are excited, and thus rotate around the *ρ*
_{2} axis by an amount that depends both on the spectral content of the pulse, and the area or the pulse. For the hyperbolic secant pulse with a hyperbolic secant spectrum, the spectral excitation follows the spectral hyperbolic secant, and thus the Bloch vectors are somewhat distributed in the *ρ*
_{1,3} plane. Since the different detunings have linearly varying phase evolution rates, the corresponding Bloch vectors quickly fan out and spread apart as they rotate around the *ρ*
_{3} axis at a rate proportional to their resonance frequency, undergoing what is called inhomogeneous dephasing. The effect of this dephasing is that the vector sum of all the Bloch vectors, which is responsible for the polarization induced in the material, sums to approximately zero once the Bloch vectors have evenly spread throughout the *ρ*
_{1,2} plane. However, since this inhomogeneous dephasing process is deterministic and not irreversible like homogeneous dephasing, and the ensemble of Bloch vectors is controllable by subsequent applications of optical pulses, one can rephase the dephasing Bloch vectors with the following approach. By applying a second pulse to the material, (easiest to visualize is a *π* pulse, although other pulse areas have a similar though less efficient result), all the Bloch vectors with rotate around the *ρ*
_{2} axis by 180° and then continue to process around the *ρ*
_{3} axis in the same direction at their individual rates. This has the effect of setting up the individual inhomogeneous lines to fan back together so that they all re-phase at a time *τ*
_{12} after the second pulse on the opposite side of the Bloch sphere than they were initially excited to. This rephasing of the Bloch vectors will then induce a polarization in the material, and cause the two-pulse photon echo to be emitted coherently [9]. The second method is based on Fourier analysis. Because the two pulses act upon the material within its coherence time, a spectral grating with period $\frac{1}{{\tau}_{12}}$ is recorded. The second half of pulse two can read out this grating, causing the 2P echo at a time *τ*
_{12} after the second pulse [7]. This second approach also lends itself to explaining the 3P photon echo, since any pulse applied within the lifetime of the excited state will diffract off the spectral grating, causing the three pulse photon echo. The Bloch sphere picture of the 3P echo offers less insight, but one can regard the second pulse as a mechanism that converts coherence information in the equitorial plane into much longer-lived excited state information along the polar axis. The third pulse acts as a rephasing event, turning the spectral grating (excited state information) back into coherence information, allowing an echo to form a time *τ*
_{12} later. This process is best illustrated by observing the simultaneous evolution of the Bloch vectors for all detunings, since the vector sum of these is directly related to the strength of the polarization induced in the material. Figure 2 shows all the pertinent information from a short-pulse PE simulation using the IPC-FDTD code. The electric field shown in Fig. 2(a) is launched into the left end of the simulation grid, where at every location in the grid (at successively delayed propagation times) it creates a time-varying excited state inversion as shown in Fig. 2(b). This excited state inversion for the different detunings numerically ranges from -1 to 1, and corresponds directly to the *ρ*
_{3} value of the Bloch vector. Once the atomic ensemble has been lifted out of the ground state, the inhomogeneous lines can dephase at their particular rate, according to the value of inhomogeneous detuning, while slowly decaying due to the decoherence goverened by *T*
_{2}. This is illustrated in Fig. 2(c), in which the data has been down-converted into a rotating coordinate frame corresponding to the central detuning to remove the confusion which a 200 THz temporal carrier would cause. In essence, this down-conversion puts the viewer into a rotating coordinate frame, so that only the relative motion between the Bloch vectors is visible, but the overall continuous spinning around the Bloch-sphere at 200 THz has been eliminated. Finally, Fig. 2(d) shows the electric field detected at the output of the material, which demonstrates that the inhomogeneous lines have collectively generated a polarization that is momentum matched to efficiently radiate into the propagating electric field. This process is the underlying mechanism for both the 2PE and 3PE. The cause of the echos is easiest to see in Fig. 2(c) - the reversal of the dephasing of the different lines causes them to undo the initial dephasing, so that they come back into coincidence when the echo occurs, which is seen as a Moiré-pattern at the time of the echo. An artificial decoherence was implemented to shorten the simulation time. In an experiment, one would normally wait several *T*
_{2}’s before applying the third pulse, since this assures that the third pulse does not interact coherently with the first two. Since this idle waiting makes the simulations unbearably long, we multiplied the values for *u*
_{1,2} by a smoothly varying decay term 0.6 ps before applying the third pulse.

Another way of viewing the same information is with an animated Bloch sphere evolution, such as that depicted in the movie referenced by Fig. 3. The data displayed in this movie is the same as that which is displayed in Fig. 2, with the exception that the electric field has also been placed in a rotating coordinate frame (mixed to DC with the central frequency carrier) for easier visualization of the data. The left part of the movie shows a perspective view of the Bloch vectors evolving, where each of the different homogeneously broadened lines is represented by a colored vector whose detuning frequency is represented by a rainbow coloring scheme. This perspective view shows the state of the inversion relatively well, while the right part of the movie shows a projection of the *ρ*
_{1,2}, or in-phase and quadrature coherence plane. This is the plane in which the polarization in the material (which is responsible for the photon echo) is generated by the vector sum of all the Bloch vectors. Below the two sets of evolving Bloch spheres is a plot of the rotating coordinate frame electric field that is measured at the same location in the material where the shown Bloch vectors are valid. This measurement location is at the end of the material, since the echo requires some propagation distance to form, but other such phase delayed Bloch spheres are evolving at each spatial sampling point in the simulation grid and contributing to the final radiated output. The movie starts with all the Bloch vectors in the ground-state. When the first pulse hits the material, the Bloch vectors rotate up towards the equitorial plane, with the central absorbers rotating into the equatorial plane, indicating a $\frac{\pi}{2}$ pulse. The Bloch vectors with resonance frequencies at the spectral wings are less efficiently excited due to the hyperbolic secant nature of the spectrum. The phase of all 301 Bloch vectors progresses as a spinning around the vertical axis at a rate proportional to their frequency until the second pulse hits. Next, the second pulse (also $\frac{\pi}{2}$) drives the Bloch vectors with a positive *ρ*
_{1} upwards, while others which have processed around to have negative *ρ*
_{1} are driven downwards. Since the phase shifts in processing around the Bloch sphere are linearly proportional to detuning, this is the process by which the spectral grating in Fig. 2(b) is generated. The phase evolution after pulse two causes the two-pulse photon echo to occur. The coherence plane does not show this coincidence of the Bloch vectors very clearly, since we used a $\frac{\pi}{2},\frac{\pi}{2},\frac{\pi}{2}$ pulse sequence, and the Bloch vectors that are mostly responsible for the 2P echo are hidden below other Bloch vectors in the *ρ*
_{1,2} plane view. This set of pulses was chosen to preferentially illustrate the three-pulse echo response, at the cost of making the two-pulse response less visible. When the simulation is run using a $\frac{\pi}{2}$, *π* pulse sequence, the 2P echo is much more apparent. However, if one carefully watches the movie, one can see the Bloch vectors come into coincidence in the lower half of the right plot in the movie. After the 2P echo, time evolves, and the coherence plane once again looks like a random progression of vectors whose sum is close to zero. Note that the values of *ρ*
_{1} and *ρ*
_{2} slowly shrink, in accordance with the coherence decay based on *T*
_{2}. In order to keep the simulation time shorter, we added an additional coherence decay to the simulation, as indicated in the movie, which collapses the *ρ*
_{1,2} plane completely, leaving behind only the spectral grating in *ρ*
_{3}(*ω*). The same effect could be achieved by simulating a much longer delay between pulses two and three, but we found that this increases the simulation run time to an unacceptable duration. Fortunately, we have not seen any adverse effects due to the artificial de-cohering as long as it is done adiabatically. Artificially decohering too quickly causes an echo-like numerical artifact to occur. After the Bloch vectors have lost all coherence information, and only the inversion information (spectral grating) remains, the third pulse illuminates the material. The third pulse acts as another $\frac{\pi}{2}$ pulse, and folds the spectrally varying inversion back into the coherence plane where they renew their phase progression by spinning around the *ρ*
_{3} axis at a rate proportional to their frequency/detuning (depending on whether the information is viewed in the stationary/rotating coordinate frame, respectively). At the moment when the third pulse is acting on the Bloch vectors one can see a type of polar representation of the positive and negative components of the spectral grating folding back into the *ρ*
_{1,2} plane. The most interesting moment in the simulation occurs at *τ*
_{12} after the third pulse (and is shown in Fig. 3), when most of the Bloch vectors come into coincidence in one quadrant of the coherence plane, causing the three-pulse photon echo. The 3-pulse echo magnitude is proportional to the vector sum of all of the detuned lines, which on average cancel out at all times, except during the time of the photon echo, where all of the rotating frame coherence vectors in the plane of coherence point (on average) into the same quadrant.

## 5. Photon echo using chirped pulses

In this section, we will investigate the behaviour of the IPC-FDTD simulator when the IBA is illuminated with chirped optical pulses. Chirped pulses are technologically important, since they permit large bandwidth programming or readout using modulated low-power CW lasers, while potentially maintaining the power-spectral density of high-power bandwidth-limited pulses [10]. This is illustrated by our second set of simulations, in which the time-bandwidth product of the pulses is on the order of TB=10 (see Fig. 4(a)), and one can see that the peak electric field strength is one third of that used in the brief-pulse case. Hence the optical power is about nine times smaller than in the brief pulse case, but after the first pulse the Bloch vectors for most detunings undergo a $\frac{\pi}{2}$ rotation, indicating that the peak power spectral density is the same as it was in the brief-pulse case (this can be seen in Fig. 4(b) and the movie linked by Fig. 5). In fact, since the chirp’s power spectrum is flatter than the sech pulse, it does a better job of transferring the various Bloch vectors of different detunings into the equitorial plane than did the sech pulse. The particular pulse sequence chosen here consists of three chirps, where the chirp rate of the first pulse is half that of the following two chirps (200 THz/ps and 400 THz/ps, respectively). These chirp rates can be directly seen as the sloped features in the excitation sequence of the detunings vs. time images (Fig. 4(b) & (c)). This sequence was chosen because it results in efficient, bandwidth limited 2P and 3P echoes as brief pulse outputs.[10] As is visible in Fig. 4(b), the first chirped pulse causes the two level transitions of the various detunings to excite sequentially, so that the second chirp interferometrically generates a chirped excited state grating as each corresponding spectral component arrives. Figure 4(c) shows the evolution of the in-phase components of all the Bloch-vectors in the coherence plane, and one can again see the re-phasing as a hyperbolic Moiré pattern in the in-phase component of the coherence, which results in a pulse-compressed 2P echo output. After 1.6 ps, the artificial decoherence is applied, leaving behind the chirped spectral grating along *ρ*
_{3}, and followed by the matched chirped readout pulse. As before, the third chirp reads out the spectral grating, and since it is a $\frac{\pi}{2}$ pulse across its bandwidth, it efficiently transforms the spectral grating into a uniform population at a time that linearly increases with detuning. As the corresponding spectral components arrive, the chirp rates are matched as $\frac{1}{{b}_{1}}=\frac{1}{{b}_{2}}+\frac{1}{{b}_{3}}$ like a lens focusing relation, which thus results in the pulse compressed 3P photon echo. This pulse compression process is particularly well illustrated using the spectral domain analytic representation valid in the small pulse regime, [3, 7] giving electric field output

$$=\int {e}^{-i\frac{2\pi}{{b}_{1}}{\omega}^{2}}{e}^{i\left[\frac{2\pi}{{b}_{2}}{\omega}^{2}+\omega {\tau}_{12}\right]}{e}^{i\left[\frac{2\pi}{{b}_{3}}{\omega}^{2}+\omega {\tau}_{3}\right]}{e}^{-i\omega t}d\omega \approx \delta \left(t-{\tau}_{12}\right),$$

where *S*
_{1} through *S*
_{3} are the chirp Fourier spectra, which, when chirp rate and phase matched, yield an impulsive photon echo. For the chirped pulse interaction, the animated Bloch sphere evolution (linked by Fig. 5) shows the linearly progressing excitation process for the various detunings very well. Initially, the lower detunings (red) are excited to the equatiorial plane, with the higher detunings (purple) following several hundred femtoseconds later. The two-pulse echo is clearly visible in the movie, with most Bloch-vectors re-phasing in one quadrant of the coherence-plane. Next, the artificial decoherence collapses the Bloch-vectors onto the vertical, *ρ*
_{3} axis, leaving behind a chirped spectral grating sampled by the 301 modeled detuned absorbers. This chirped grating is difficult to see in the movies due to the overlap on the *ρ*
_{3} axis, but can easily be seen in Fig. 4(b). Finally, the third chirp reads out the spectral grating sequentially, transfers the spectral components to the equitorial plane, and the three-pulse echo occurs when most of the Bloch-vectors re-phase in one quadrant of the sphere. This modeling and illustration technique could be used for a variety of other processing techniques, but the computational inefficiency of the FDTD technique precludes modeling of multidimensional spatial-spectral holography or the precise modeling of GHz bandwidth operations currently of great experimental interest for RF signal processing applications [4, 5, 6, 7, 11].

## 6. Conclusion

We have successfully implemented an IPC-FDTD simulator designed to predict the interaction of short optical pulses with inhomogeneously broadened media. Two- and three-pulse photon echoe simulations were demonstrated from brief- and chirped-pulse sequences inside a material with 301 lines of inhomogeneous detuning. The complicated dynamics of these 301 spectral components were illustrated with plots as well as color-coded movies of multiple simultaneous Bloch vectors on the Bloch sphere as a powerful visualization technique. The formation of 2P and 3P echoes was illustrated on the Bloch-sphere, and the detailed rephasing of all the inhomogeneous lines that leads to the formation of the photon echoes was graphically illustrated.

## Acknowledgments

We wish to acknowledge support from W. Schneider of DARPA (MSU subcontract MDA972-03-1-0002), as well as valuable discussions withW.R. Babbitt, K. Anderson and M. Mohamed.

## References and links

**1. **T. Chang, M. Tian, Z.W. Barber, and W. R. Babbitt, “Numerical modeling of optical coherent transient processes with complex configurations - II. Angled beams with arbitrary phase modulations,” J. Luminescence **107**, 138–145 (2004). [CrossRef]

**2. **R. Ziolkowski, J. Arnold, and D. Gogny, “Ultrafast pulse interaction with two-level atoms,” Phys. Rev. **52**, 3082–3094 (1995). [CrossRef]

**3. **M. Mitsunaga and R. G. Brewer, “Generalized perturbation theory of coherent optical emission,” Phys. Rev. A **32**, 1605–1613 (1985). [CrossRef] [PubMed]

**4. **K. Merkel, R. K. Mohan, Z. Cole, T. Chang, A. Olson, and W. Babbitt, “Multi-Gigahertz radar range processing of baseband and RF carrier modulated signals in Tm:YAG,” J. Luminescence **107**, 62–74 (2004). [CrossRef]

**5. **V. Lavielle, F. D. Seze, I. Lorgere, and J.-L. L. Gouet, “Wideband radio frequency spectrum analyzer: Improved design and experimental results,” J. of Luminescence **107**, 75–89 (2004). [CrossRef]

**6. **Y. S. Bai, W. R. Babbitt, N. W. Carlson, and T. W. Mossberg, “Real-Time Optical Waveform Convolver Cross Correlator,” Appl. Phys. Lett. **45**, 714–716 (1984). [CrossRef]

**7. **T. W. Mossberg, “Time-Domain frequency selective optical data storage,” Opt. Lett. **7**, 77–79 (1982). [CrossRef] [PubMed]

**8. **S. Nakamura, *Applied Numerical Methods in C*, 1st ed., Prentice-Hall Inc, Englewood Cliffs, New Jersey (1993).

**9. **L. Allen and J. Eberly, *Optical Resonance and Two-Level Atoms*, 2nd ed., General Publishing Ltd., Don Mills, Toronto (1987).

**10. **K. D. Merkel and W. R. Babbitt, “Coherent transient optical signal processing without brief pulses,” Appl. Opt. **35**, 278–285 (1996). [CrossRef] [PubMed]

**11. **F. Schlottau and K. H. Wagner, “Demonstration of a continuous scanner and time-integrating correlator using spatial-spectral holography,” J. Luminescence **107**, 90–102 (2004). [CrossRef]