Introduction
Nuclear-pumped lasers (NPL) are systems where population inversion in a lasing medium is achieved through ionizing radiation from nuclear reactions. The ideal NPL system would be incorporated in a self-critical reactor configuration, hence products from fission reactions are the most commonly studied. Virtually all NPL systems studied are gas based due to the insensitivity of the lasing medium to radiation damage.
The most common NPLs studied were gaseous where the laser medium was contained within a metal tube. The fuels usually consisted of uranium (Helmick et al., Reference Helmick, Fuller and Schneider1975; McArthur and Tollefsrud, Reference McArthur and Tollefsrud1975; Voinov et al., Reference Voinov, Dovbyshov, Krivonosov, Mel'nikov, Podmoshenskii and S1979; Alford and Hays, Reference Alford and Hays1989; Hebner and Hays, Reference Hebner and Hays1993) and B-10 (DeYoung et al., Reference DeYoung, Wells, Miley and Verdeyen1976; Akerman et al., Reference Akerman, Miley and McArthur1977; Prelas et al., Reference Prelas, Akerman, Boody and Miley1977) film coatings or He-3 (Jalufka et al., Reference Jalufka, DeYoung, Hohl and Williams1976; De Young et al., Reference DeYoung, Jalufka and Hohl1977; DeYoung et al., Reference DeYoung, Jalufka and Hohl1978; Shaban and Miley, Reference Shaban and Miley1993; Williams and Miley, Reference Williams and Miley1993), BF3 (Prelas and Loyalka, Reference Prelas and Loyalka1981), and UF6 (Hassan, Reference Hassan1980) gasses mixed with the lasing medium. The gas is kept at around atmospheric pressure giving the nuclear products a range of a few centimeters. The tubes were designed such that the energy deposited in the gas is nearly uniform down the central axis. The gasses used in the previous studies had high-power thresholds thus requiring a reactor to be pulsed to produce a large enough neutron flux to reach the lasing threshold. In the end, such lasers had efficiencies barely above 1% (Prelas et al., Reference Prelas, Watermann, Wisniewski, Neher and Weaver Iii2014). The goal of these NPL systems is to group enough tubes together to create a self-critical reactor. To do so with this sort of design leads to large reactor systems (Prelas, Reference Prelas2016).
A semiconductor-based NPL is a solid-state system which takes thin layers of a semiconductor material and directly interface them with a thin layer of fuel. Such layers could be stacked to create hundreds or thousands of micro-cells. This NPL would require a lower power threshold to induce lasing and would be compact. An NPL of this sort could potentially create a self-critical system with a volume on the order of a cubic meter or less.
An NPL of this size would be a candidate for a space-based laser. The potential high-power CW beam produced could be used for missile defense (Carter, Reference Carter1984), space mining (Gibbings et al., Reference Gibbings, Vasile, Hopkins, Burns and Watson2012), space propulsion (Phipps and Luke, Reference Phipps, Luke and Phipps2007), low orbit debris elimination (Phipps et al., Reference Phipps, Albrecht, Friedman, Gavel, George, Murray, Ho, Priedhorsky, Michaelis and Reilly1996), and even meteor deflection (Prelas et al., Reference Prelas, Watermann, Wisniewski, Neher and Weaver Iii2014).
This paper introduces a model for a transverse electron beam-pumped semiconductor laser (see Fig. 1) as an analog to a fission-based system. Our goal is to show the validity of such a model. This is achieved by benchmarking a transport model with the proposed experimental studies using an e-beam. The next goal is to extend the model to fission fragments. The similarity between pumping with energetic electrons and fission fragments is that the key species in the laser, electron–hole pairs, are produced mainly by secondary and higher order elections thus leading to similar W values (the energy required to produce an e–h pair through particle collision, defined as eV/ion pair).
Transverse electron beam laser diagram. Electrons impinge upon a semiconductor crystal parallel to the x-axis. Electrons deposit their energy into the shaded region creating the population inversion. The crystal is cleaved on two parallel planes one of which is mirrored to 100% reflectivity. Laser output is along the z-axis. The electron beam is assumed uniform such that the electron–hole population does not vary along the y and z axes.
To model volumetric energy deposition, transport codes such as GEANT, PENELOPE, and MCNP can calculate them with a 3D mesh. The results of these simulations predict where the inversion is achieved. The volumes of these regions are irregularly shaped, depending both on the particle type and energy, and can simply not be represented generally by empirical equations.
The model for the electron beam provides data which can be used in benchmarking experiments. Part of the benchmarking process is to use the model to predict laser intensity and far-field distribution. By comparing the results from computational EBL models and physical experiments, the properties of a fission-based system can be extrapolated from a theoretical model.
The use of an electron beam is attractive because NPL experiments are complex and costly. A nuclear reactor is required and the experiments with fissile material are nearly impossible at the non-governmental level due to the cost and regulatory burden. The construction of the model incorporates several modules as described below.
Power density module
Electron beam-pumped lasers are able to reach oscillation with pure semiconductor materials. However, a model based on the population difference between electrons and holes is not sufficient. Additionally, due to the nature of free carrier generation through ionization, the distribution of e–h pairs is not uniform across the volume. By considering a temporally and spatially dependent complex permittivity coefficient, these issues can be accounted for. This is the approach Bogdankevich et al. took to solve these problems (Bogdankevich et al., Reference Bogdankevich, Goncharov, Lavrushin, Letokhov and Suchkov1967; Bogdankevich et al., Reference Bogdankevich, Letokhov and Suchkov1969). The relative permittivity is defined as
$${\rm \varepsilon} \,{\rm (}x,t{\rm )} = {\rm \varepsilon} _{\rm r} + {\rm \delta \varepsilon} \,{\rm (}x,t{\rm )} + i({\varepsilon} ^{\prime \prime}_0-{\varepsilon} ^{\prime \prime}(x,t))$$where ε″ relates the electron population which contributes to gain
$${\rm {\varepsilon}^{\prime \prime}}\,{\rm (}x,t{\rm )} = \displaystyle{\sigma \over k}N(x,t)$$where k is the wavenumber (k 2 = ω2ε r/c 2) and σ is, as stated by the cited authors, “the cross-section for radiative recombination averaged over the linewidth.” To keep with this definition, it is defined here as
$${\rm \sigma} = \langle{{\rm \sigma} \lpar {\rm \nu} \rpar } \rangle = \displaystyle{1 \over {\Delta {\rm \nu}}} \int\limits_{{\rm \nu} _0-\Delta {\rm \nu} /2}^{{\rm \nu} _0 + \Delta {\rm \nu} /2} {\displaystyle{{c^2} \over {8{\rm \pi} {\rm \nu} ^2 {\tau}_1}}g_{\rm \nu} ({\rm \nu} )d{\rm \nu}} $$where ν0 is the center line of the laser frequency, Δν is the linewidth due to broadening, and g ν(ν) is the lineshape. For the purposes of this study, natural broadening was considered with a Lorentzian lineshape.
$$\Delta {\rm \nu} = \displaystyle{1 \over {2{\rm \pi}}} \left( {\displaystyle{1 \over {{\tau}_1}} + \displaystyle{1 \over {{\tau}_2}}} \right)$$
$$g_{\rm \nu} ({\rm \nu} ) = \displaystyle{{\lpar {\Delta {\rm \nu} /2{\rm \pi}} \rpar } \over {{\lpar {{\rm \nu} -{\rm \nu}_0} \rpar }^2 + {\lpar {\Delta {\rm \nu} /2} \rpar }^2}}$$where τ1 and τ2 are the lifetimes of upper and lower lasing states for an arbitrary laser system. In this case, the lower state is stable hence τ2 = ∞. Inserting Eq. (5) into Eq. (3) and making the change of variables ν = ν0f inside the integral
$${\rm \sigma} = \displaystyle{1 \over { {\tau}_1{\rm \nu} _0^3}} \left( {\displaystyle{c \over {4{\rm \pi}}}} \right)^2\int\limits_{1-\displaystyle{{\Delta {\rm \nu}} \over {{\rm \nu} _0}}\displaystyle{1 \over 2}}^{1 + \displaystyle{{\Delta {\rm \nu}} \over {{\rm \nu} _0}}\displaystyle{1 \over 2}} {\,f^{-2}{\left( {{\lpar {\,f-1} \rpar }^2 + {\left( {\displaystyle{{\Delta {\rm \nu}} \over {{\rm \nu}_0}}\displaystyle{1 \over 2}} \right)}^2} \right)}^{-1}df} $$The term ε0″ relates the permittivity to the photon lifetime in the cavity τp.
$${\varepsilon} ^{\prime \prime}_0 = \displaystyle{1 \over {{\rm \omega} _0 {\tau}_{\rm p}}},\,\, {\tau}_{\rm p} = \displaystyle{{2L_{\rm c}} \over c}\left( {\ln \displaystyle{1 \over {R_1R_2}} + 2{\rm \alpha} L_{\rm c}} \right)^{-1}$$where R 1 and R 2 are the reflective coefficients for the front and back mirrors where the backing mirror has been assumed to be perfect, L c is the length of the resonator cavity, ω0 is the angular frequency of the principle laser line, and α is the linear loss coefficient. The δε term accounts for the change in refractive index from the inhomogeneous electron–hole plasma, given by
$$\eqalign{{\rm \delta \varepsilon} \,{\rm (}x,t{\rm )} & = -\displaystyle{{4{\rm \pi} e^2} \over {m{^\ast}{\rm \varepsilon} _0{\rm \omega} _0^2}} N\,(x,t) = -\displaystyle{{4{\rm \pi} e^2} \over {m{^\ast}{\rm \varepsilon} _0{\rm \omega} _0^2}} \displaystyle{k \over {\rm \sigma}} {\rm {\varepsilon}^{\prime\prime}}\,{\rm (}x,t{\rm )} \cr & = -{\rm \chi} _1{\rm {\varepsilon}^{\prime \prime}}\,{\rm (}x,t{\rm )}}$$where m* is the effective mass of electrons in the cavity and e is the electron charge. The equations for the field and complex permittivity in the cavity areFootnote a
$$\displaystyle{{\partial E} \over {\partial t}} = \displaystyle{{{\rm \omega} _0} \over 2}\left[ {\displaystyle{i \over {k^2}}\displaystyle{{\partial^2E} \over {\partial x^2}} + \displaystyle{{{\rm {\varepsilon}^{\prime \prime}}\,{\rm (}x,t{\rm )} + i{\rm \delta \varepsilon} \,{\rm (}x,t{\rm )}-{{\rm {\varepsilon}^{\prime \prime}}}_0} \over {{\rm \varepsilon}_{\rm r}}}E} \right]$$
$$\displaystyle{{\partial {\rm {\varepsilon}^{\prime \prime}}} \over {\partial t}} + \displaystyle{1 \over { {\tau}_1}}{\rm {\varepsilon}^{\prime \prime}} = \displaystyle{{\rm \sigma} \over k}g\,(x)-2{\rm \sigma} I{\rm {\varepsilon}^{\prime \prime}}$$where τ1 is the free carrier lifetime and I is the field intensity at the resonator face related to the field by I = (ħω0)−1(c/2)ε0|E|2 = χ0|E|2. The function g(x) is the e–h population distribution curve. The system can be reduced to a non-dimensionalized form by creating new dimensionless variables: x = ξ/k, t = τ1s
$$\displaystyle{2 \over {{\rm \omega} _0 {\tau}_1}}\displaystyle{{\partial U} \over {\partial s}} = i\displaystyle{{\partial ^2U} \over {\partial {\rm \xi} ^2}} + \displaystyle{{(1-i{\rm \chi} _1){\rm {\varepsilon}^{\prime \prime}}-{{\rm {\varepsilon}^{\prime \prime}}}_0} \over {{\rm \varepsilon} _{\rm r}}}U$$
$$\displaystyle{{\partial {\rm {\varepsilon}^{\prime \prime}}} \over {\partial s}} = Q\,({\rm \xi} )-(1 + {\rm \delta} ^2\vert U \vert ^2){\rm {\varepsilon}^{\prime \prime}}$$where several new terms were defined
$$E\,(x,t)\to u_0U\,({\rm \xi}, s),\,\,\,Q\,({\rm \xi} ) = {\tau}_1\displaystyle{{\rm \sigma} \over k}g({\rm \xi} ),\,\,\,{\rm \delta} ^2 = 2 {\tau}_1{\rm \sigma} {\rm \chi} _0u_0^2 $$where u 0 is a scaling constant which carries the dimensions of the field. Its value is arbitrary and is set to 1 V/m. The initial and boundary conditions in Eqs. (11) and (12) are
$$\eqalign{& U\,(0,t) = U\,(Lk,t) = 0 \cr & U({\rm \xi}, 0)\approx 0 \cr & {\rm {\varepsilon}^{\prime \prime}}({\rm \xi}, 0) = 0} $$where L is the maximum distance into the crystal. Equation (11) is homogeneous and as such cannot be set identically to zero for its initial condition. Instead all the interior points of the initial guess can just be set arbitrarily small.
Next the source distribution term must be defined. The work by several previous authors (Bogdankevich et al., Reference Bogdankevich, Letokhov and Suchkov1969; Johnston, Reference Johnston1971) simplified the model by assuming ε″(x)~Sech2(x) because this allowed for exact analytical steady-state solutions in the form of hypergeometric functions.
To find a true distribution curve, MCNP was employed. A 1 cm diameter beam of 100 keV electrons was transported through the slab of GaAs with dimensions 1 × 1 × 0.006 cm which was divided into a 100 × 100 × 100 3D mesh. The full simulation transported 5 × 108 electrons with an approximate run time of 3 weeks.
Because MCNP accounts for reflections and electrons scattering out of the medium, it was necessary to calculate what fraction of the beam energy was deposited into the volume. The energy deposition was calculated using the TMESH tally in MCNP which gives results in the units of MeV/(cm3 source particle). The energy fraction deposited into the volume is then
$$\eqalign{e_{{\rm frac}} &= \left[ {\sum\limits_{i,j,k} {e_{ijk}\left[ {\displaystyle{{{\rm MeV}} \over {{\rm c}{\rm m}^{\rm 3}\cdot {\rm source}\,{\rm particle}}}} \right]} \times {\rm Cell\_Volume}\lsqb {{\rm c}{\rm m}^3} \rsqb } \right] \cr & \quad / {\rm Electron\_Energy}\left[ {\displaystyle{{{\rm keV}} \over {{\rm particle}}}} \right]}$$where e ijk is the energy tally in the ijk-th cell. The result is still in the units of energy and is then divided by the electron energy to find the total fraction of the beam energy which is deposited into the volume. This factor was calculated to be 51.8%. The simulation data are then normalized
$$f_{ijk} = e_{ijk} \times \left[ {\sum\limits_{i,j,k} {e_{ijk}\left[ {\displaystyle{{{\rm MeV}} \over {{\rm c}{\rm m}^{\rm 3}\cdot {\rm source}\,{\rm particle}}}} \right]}} \right]^{-1}$$This allows the power deposition in any single cell to be related to the total electron beam power: Total Power × f ijk × e frac. The carrier generation rate in a single cell is then easily calculated. A 1 Amp electron beam of 100 keV electrons has a power output of 105 W. To convert this into e–h pairs per second, divide the beam power by the material's W value. The W value can be approximated for a wide variety of wide bandgap semiconductor materials by the Klein formula (Revankar and Adams, Reference Revankar and Adams2014)
$$W = 2.8 \times E_{\rm g} + 0.5\,\lsqb {{\rm eV}} \rsqb $$where E g is the material bandgap. In the case of GaAs, E g = 1.42 eV yielding W = 4.4 eV/pair. Thus the final form of the data is
$$\eqalign{g_{ijk} &= f_{ijk}/v_{\rm cell}\lsqb {\rm cm}^3 \rsqb \times \lpar {10}^5\,{\rm W/A} \rpar /{\rm W}\lsqb {4.4\,{\rm eV/pair}} \rsqb \times e_{\rm frac} \cr & = f_{ijk} \times \left[ {9.64597 \times {10}^{30}{\rm pairs} \over {{\rm cm}^3\cdot {\rm s}\cdot {\rm A}}} \right] = f_{ijk}g_0}$$The model can now be easily manipulated for any electron beam current.
While in principle this model could be applied to two or three dimensions, for the sake of simplicity, the data were averaged into a single vector. Specifically, because the data are circularly symmetric and the width of the electron beam is much greater than the straggling distance of electrons in GaAs, the data in any direction perpendicular to the beam axis are approximately constant and was summed out into a 2D matrix. What this matrix now represents is a 100 × 100 grid of columns whose data represent the fraction of beam energy deposited into a specific column and thus you must divide g 0 in Eq. (18) by 100. The centerline of this 2D matrix was then chosen for having the maximum energy deposition and furthest penetration depth. With the data in proper form, a curve was fitted to the normalized data with the following form
$$f\,(x) = a\exp \lsqb {-b^2\,{(x-\mu )}^2} \rsqb \times Sech \,(cx)$$The fitting parameters are (a, b, μ, c) = [465.076, 824.326 cm−1, 5.629393 × 10−3 cm, −5633.573 cm−1]. Plot x compares the fitted curve to the MCNP data as shown in Figure 2.
Electron–hole generation distribution curve. The step plot is taken directly from MCNP data, the smooth plot is the fitted curve.
Finally the pump can be scaled to any electron beam power by multiplying by a factor P in the units of Amps.
$${\rm PUMP}\,{\rm (}x,P{\rm )} = Pg\,(x) = Pf\,(x)g_0/10^2$$ Equations (11) and (12) can then be solved for any arbitrary beam current. The intensity of the beam will not change the shape of the distribution, only its magnitude. Because the entire beam spot size is no longer being considered, it will be more useful to discuss the results in terms of electron beam current density, 
$\tilde{P}$ = P/[π·(0.5 cm)2].
To reiterate the equations to be now solved are
$$\displaystyle{2 \over {{\rm \omega} {\rm \tau} _1}}\displaystyle{{\partial U} \over {\partial s}} = i\displaystyle{{\partial ^2U} \over {\partial {\rm \xi} ^2}} + \left( {\displaystyle{{\lpar {1-i{\rm \chi}_1} \rpar {\rm {\varepsilon}^{\prime \prime}}-{{\rm {\varepsilon}^{\prime \prime}}}_0} \over {{\rm \varepsilon}_{\rm r}}}} \right)U$$
$$\displaystyle{{\partial {\rm {\varepsilon}^{\prime \prime}}} \over {\partial {\rm s}}} = PQ-(1 + {\rm \delta} ^2\vert U \vert ^2){\rm {\varepsilon}^{\prime \prime}}$$The full electron beam current, P, is still used in Eq. (22) because the choice for Q accounts for the limited area which the beam interacts.
A table of pertinent properties and evaluations of the constants is given in Table 1.
Defined and evaluated constants
Solutions of steady-state equation
The system of Eqs. (21) and(22) was solved semi-implicitly by a method used by Kubíček and Hlaváček (Reference Kubíček and Hlaváček1983). In this method, the permittivity can be exactly solved for implicitly via backwards Euler method
$${\rm \varepsilon} _l^{n + 1} = \displaystyle{{{\rm \varepsilon} _l^n + \Delta sPQ_l} \over {1 + \Delta s\lpar {1 + {\rm \delta}^2{\vert {U_l^{n + 1}} \vert }^2} \rpar }}$$where the double primes on the permittivity have been dropped for brevity. The field equation is solved semi-implicitly where the differential operator is evaluated at the “new” time step and the non-linear terms are calculated at the “old” time step.
$$\eqalign{& \displaystyle{2 \over {{\rm \omega} {\tau}_1}}\displaystyle{{U_l^{n + 1}} \over {\Delta s}}-i\displaystyle{{U_{l + 1}^{n + 1} -2U_l^{n + 1} + U_{l-1}^{n + 1}} \over {\Delta {\rm \xi} ^2}} \cr & = \left( {\displaystyle{{\lpar {1-i{\rm \chi}_1} \rpar {\rm \varepsilon}_l^n -{{\rm {\varepsilon}^{\prime \prime}}}_0} \over {{\rm \varepsilon}_{\rm r}}} + \displaystyle{2 \over {{\rm \omega} {\tau}_1\Delta s}}} \right)U_l^n}$$Here the basic second order finite difference was applied to the derivative. In this scheme, the system of equations generated by Eq. (24) is solved and the resulting solutions are inserted into Eq. (23). The advantage with this scheme is no iterative submethod is necessary for each time step as opposed to a fully implicit scheme. Solving the linear system was carried out using Mathematica.
The maximum distance into the crystal considered was L = 60 µm. The solution of interest is the field intensity immediately outside the resonator restated here as
$$I\,(x) = \lpar {1-R} \rpar \varepsilon _0\displaystyle{c \over 2}\vert {u_0U} \vert ^2\,\,\left[ {\displaystyle{{\rm W} \over {{\rm c}{\rm m}^{\rm 2}}}} \right]$$where the factor (1–R) is the transmission coefficient for the front mirror. In the neighborhood of threshold, the average laser intensity has the plot seen in Figure 3 as a function of the electron beam current density.
Laser threshold curve. Each dot represents the average integral of the resultant steady-state fields squared for the given pumping power.
According to the model, it appears that the threshold is approximately 1.8 nA/cm2. Near the threshold, the full-time-dependent field square has the form seen in Figure 4.
Full plot of field intensity with respect to space and time.
Along the peak, the time-dependent curve has the form seen in Figure 5.
Curve along the point of maximum intensity from Figure 4.
For pumping intensities well beyond the threshold, the solution curves take the form seen in Figure 6.
Electric field squared at various pumping powers. Maximum field intensity is at x = 16 µm.
Despite the spatial inhomogeneity, the average laser intensity increases linearly (Fig. 7), which can be approximated by
$$I_{{\rm avg}}\lpar {\tilde{P}} \rpar = \displaystyle{1 \over L}\int_0^L {I\lpar {\tilde{P},x} \rpar dx} \approx 36080.3\tilde{P} + 0.046955\,\,\,\left[ {\displaystyle{{\rm W} \over {{\rm c}{\rm m}^{\rm 2}}}} \right]$$
where 
$\tilde{P}$ is the electron beam current density in A/cm2.
Demonstration of the linear relationship between electron beam current and laser intensity for large currents.
The overall efficiency versus pumping power is derived as such
$$\eqalign{\eta \lpar {\tilde{P}} \rpar &= \displaystyle{{L \times \lpar {{\rm cell}\,{\rm width}} \rpar \times I_{{\rm avg}}} \over {\tilde{P} \times 2R \times \lpar {{\rm cell}\,{\rm width}} \rpar \times \lpar {{10}^5\,{\rm W/A}} \rpar }} \cr & = \displaystyle{L \over {2R}}\displaystyle{{36080.3\tilde{P} + 0.046955} \over {\tilde{P} \times \lpar {{10}^5\,{\rm W/A}} \rpar }} \approx \displaystyle{L \over {2R}}\displaystyle{{36080.3} \over {{10}^5}}\approx 0.0022} $$where this applies when the width of the active region is much smaller than the beam diameter.
Conclusions
A model for an electron beam-pumped laser, which accounts for spatial inhomogeneities, has been presented and solved for a variety of pumping strengths. The forcing term which accounts for these inhomogeneities was obtained directly from particle transport codes. It was shown that, despite the inhomogeneities, the power output from the laser increases linearly for a wide range of pumping powers. At high pumping strengths, the maximum efficiency was calculated to be η = 0.22%.
This model, in which the change in permittivity from free electrons is used to calculate laser properties, can be applied to any similar system regardless of how the electron–hole plasma is created. Therefore, a system where the electron–hole plasma is created by fission products is applicable. Additionally, the effective W values for fission products will be similar to those of electron beams because the majority of ionizations from fission products come from secondary and higher order ionizations. An aspect in which an electron beam will not be comparable is that of damage to the lattice due to displacements. Displacements result from non-ionizing energy loss (NIEL). Fission products will produce significantly more NIEL events over the same exposure time than an electron beam. A more comprehensive model can be derived from the approach used in this paper to simulate the physics of a solid-state nuclear-pumped laser.
Appendices
A. Field equation derivation
To derive Eq. (9), begin with the typical electric field wave equation in a non-magnetic material (μ = 1).
$$\nabla ^2{\rm E} = \displaystyle{1 \over {c^2}}\displaystyle{{\partial ^2} \over {\partial t^2}}\lpar {\varepsilon {\rm E}} \rpar $$It is assumed that the time derivatives of ε are small compared to those of ξ, thus ε can be taken out of the derivative (Suchkov, Reference Suchkov1966). We are only interested in the transverse wave and thus assume ξ = E(x,y,t)exp[i(ωt-kz)] (Johnston, Reference Johnston1971) where E is the transverse electric field, furthermore we are not interested in the case where the field could vary in the y direction and hence it is ignored. The field equation now becomes
$$\displaystyle{{\partial ^2E} \over {\partial x^2}}-k^2E = \displaystyle{{\rm \varepsilon} \over {c^2}}\left( {\displaystyle{{\partial^2E} \over {\partial t^2}}-2i{\rm \omega} \displaystyle{{\partial E} \over {\partial t}}-{\rm \omega}^2E} \right)$$It is assumed that the time derivatives of E are slower than the inverse frequency of the laser. Laser transients are typically on the order of nanoseconds while ω = 2π(343 THz), making this assumption completely justifiable. Taking this into account, the second order time derivative can be removed. Next it is a matter of applying the definitions of k and ε and a few algebraic manipulations
$$\eqalign{& \displaystyle{{\partial ^2E} \over {\partial x^2}}-\displaystyle{{{\rm \varepsilon} _{\rm r}{\rm \omega} ^2} \over {c^2}}E = \displaystyle{{{\rm \varepsilon} _{\rm r} + \delta {\rm \varepsilon} + i\lpar {{{\rm {\varepsilon}^{\prime \prime}}}_0-{\rm {\varepsilon}^{\prime \prime}}} \rpar } \over {c^2}}\left( {-2i{\rm \omega} \displaystyle{{\partial E} \over {\partial t}}-{\rm \omega}^2E} \right) \cr & \displaystyle{{c^2} \over {{\rm \varepsilon} _{\rm r}}}\displaystyle{{\partial ^2E} \over {\partial x^2}}-{\rm \omega} ^2E = \displaystyle{{{\rm \varepsilon} _{\rm r} + \delta {\rm \varepsilon} + i\lpar {{{\rm {\varepsilon}^{\prime \prime}}}_0-{\rm {\varepsilon}^{\prime \prime}}} \rpar } \over {{\rm \varepsilon} _r}}\left( {-2i{\rm \omega} \displaystyle{{\partial E} \over {\partial t}}} \right)-{\rm \omega} ^2E-\displaystyle{{{\rm \delta \varepsilon} + i\lpar {{{\rm {\varepsilon}^{\prime \prime}}}_0-{\rm {\varepsilon}^{\prime \prime}}} \rpar } \over {{\rm \varepsilon} _r}}{\rm \omega} ^2E \cr & {\rm \omega} ^2\left( {\displaystyle{{c^2} \over {{\rm \varepsilon}_{\rm r}{\rm \omega}^2}}\displaystyle{{\partial^2E} \over {\partial x^2}} + \displaystyle{{{\rm \delta \varepsilon} + i\lpar {{{\rm {\varepsilon}^{\prime \prime}}}_0-{\rm {\varepsilon}^{\prime \prime}}} \rpar } \over {{\rm \varepsilon}_{\rm r}}}E} \right) = -2i{\rm \omega} \displaystyle{{{\rm \varepsilon} _{\rm r} + {\rm \delta \varepsilon} + i\lpar {{{\rm {\varepsilon}^{\prime \prime}}}_0-{\rm {\varepsilon}^{\prime \prime}}} \rpar } \over {{\rm \varepsilon} _{\rm r}}}\displaystyle{{\partial E} \over {\partial t}} \cr & \displaystyle{{\rm \omega} \over 2}\left( {\displaystyle{i \over {k^2}}\displaystyle{{\partial^2E} \over {\partial x^2}} + \displaystyle{{{\rm {\varepsilon}^{\prime \prime}} + i{\rm \delta \varepsilon} -{{\rm {\varepsilon}^{\prime \prime}}}_0} \over {{\rm \varepsilon}_{\rm r}}}E} \right) = \displaystyle{{\partial E} \over {\partial t}}}$$where ε/εr~1 because εr>>δε, ε0″, ε″.
It should also be noted that the equations stated here are not exactly the same as those in references Bogdankevich et al. (Reference Bogdankevich, Goncharov, Lavrushin, Letokhov and Suchkov1967), Bogdankevich et al. (Reference Bogdankevich, Letokhov and Suchkov1969) and Johnston (Reference Johnston1971). The signs on ε″ and 
${\rm {\varepsilon}^{\prime \prime}}_0$ are in fact reversed. It can be quickly justified that the form given here is correct. First, assume there is no spatial dependence, Eq. (9) becomes
$$\displaystyle{{\partial E} \over {\partial t}} = \displaystyle{{{\rm \omega} _0} \over 2}\left[ {\displaystyle{{{\rm {\varepsilon}^{\prime \prime}(}t{\rm )} + i{\rm \delta \varepsilon} \,{\rm (}t{\rm )}-{{\rm {\varepsilon}^{\prime \prime}}}_0} \over {{\rm \varepsilon}_{\rm r}}}E} \right]$$Next, multiply Eq. (31) by E* and add the complex conjugate of Eq. (31) multiplied by E
$$\eqalign{& \displaystyle{{\partial E} \over {\partial t}}E{^\ast} + E\displaystyle{{\partial E{^\ast}} \over {\partial t}} = \displaystyle{{{\rm \omega} _0} \over 2}\left[ {\displaystyle{{{\rm {\varepsilon}^{\prime \prime}}\,{\rm (}t{\rm )} + i{\rm \delta \varepsilon} \,{\rm (}t{\rm )}-{{\rm {\varepsilon}^{\prime \prime}}}_0} \over {{\rm \varepsilon}_{\rm r}}}EE{^\ast} + \displaystyle{{{\rm {\varepsilon}^{\prime \prime}}\,{\rm (}t{\rm )}-i{\rm \delta \varepsilon} \,{\rm (}t{\rm )}-{{\rm {\varepsilon}^{\prime \prime}}}_0} \over {{\rm \varepsilon}_{\rm r}}}EE{^\ast}} \right] \cr & \displaystyle{\partial \over {\partial t}}\vert E \vert ^2 = {\rm \omega} _0\left[ {\displaystyle{{{\rm {\varepsilon}^{\prime \prime}(}t{\rm )}-{{\rm {\varepsilon}^{\prime \prime}}}_0} \over {{\rm \varepsilon}_{\rm r}}}{\vert E \vert }^2} \right]}$$ Notice that at t = 0, ε″ = 0, and therefore the derivative is negative and there continues to be no gain until 
${\rm {\varepsilon}^{\prime \prime}} \gt {\rm {\varepsilon}^{\prime \prime}}_0$. If the signs are reversed, this would mean the system begins inverted and as ε″increases the system dies out. Clearly the latter is not physically accurate.
