As originally predicted by Letokhov^[1] in 1968, random lasers have attracted much attention in both theoretical and experimental work^[2–10]. Lasing action arises in disordered media with the combination of multiple scattering and light amplification. In a random laser, the multiple scattering process plays an important role in the lasing process. The lasing modes are determined by multiple scattering and not by a conventional laser cavity. Light scattering for an ultrashort pulse in disordered media has been studied in recent years^[11–14]. These works focus on the role of localized states of ultrashort laser pulse propagation in random media, eventually including nonlinear effects. However, multiple light scattering only provides the feedback mechanism. In order to achieve lasing action, the light amplification process should be considered.

Random laser under ultrashort pulse excitation has been realized in different types of disordered media^[15,16]; however, there is as yet no theory to describe lasing action in the ultrashort regime. A complete model must include the dynamics of the gain mechanism and the interference effects. By combining Maxwell’s equations with rate equations of the electronic population, Jiang^[4] developed the time-dependent theory for the one-dimensional (1D) case, and Sebbah^[17] extended it to the two-dimensional (2D) case. Many characteristics of random lasing have been successfully explored via this model. In early experiments, nanosecond or picosecond laser pulses were usually used as the pumping source. For numerical simulation, only a few picoseconds are required to reach the stable state in the system, so the pumping rate was always taken as a fixed value within this model. However, in the ultrashort pulse regime, the pumping process only last tens to hundreds of femtoseconds, and therefore a fixed pumping rate is inadequate. This model has been modified in our previous work and many new phenomena are explored, such as the fact that the random lasing action depends strongly on the pumping process^[18]. Despite the fact that most of the results are in agreement with previous experiments pumped by femtosecond lasers, we only used a phenomenological method to simulate the ultrashort pumping process in our model. Here we use two individual electrical fields to describe the pumping pulse and the emission pulse. By using the polarization equation, the pumping pulse was coupled into the rate equations. The pumping light will be scattered in the random media and form some spatial distribution modes, which leads to the spatially nonuniform gain in the sample. This is more accurate than using the special uniform gain in the entire sample. Our method not only allows us to investigate the spatial distribution of the pumping pulse and the random lasing emission independently, but also the relationship (or overlapping) of the spatial modes between them.

In this Letter, we report the computational results of random lasing action excited by an ultrashort pulse by changing some configuration parameters of the disordered media, including the surface-filling fraction, scatterer radius, and refractive index of the scatterers.

The 2D square random system is considered here with size $L^2$ in the $xy$–plane. It consists of circular particles with a radius $r$ and refractive index $n_2$, which are randomly distributed in a background medium with a refractive index $n_1$. The scattering strength is varied by adjusting the index contrast $\mathrm{\Delta }n=n_2-n_1$ and the surface-filling fraction $\mathrm{\Phi }(\mathrm{\Phi }=N·\mathrm{\pi }{r}^2/L^2)$, where $N$ is the amount of the particles. This system is essentially a 2D simplification of real experiments including a random collection of cylinders oriented along the $z$-axis. This model can also be considered as the multiple scattering light is confined in the $xy$–plane among a three-dimensional (3D) sample, resulting in a quasi-2D type of light transport.

The background medium is chosen as the active part and modeled as a four-level atomic system. The electrons in the ground Level 0 are transferred to the upper Level 3 by an external laser pulse excitation. Electrons in Level 3 flow downward to Level 2 by means of nonradiative decay process with a very short time. The intermediate Levels 2 and 1 are the upper and lower levels of the laser transition, respectively. The time evolution of the four-level atomic system is described by rate equations, as follows $$\frac{d{N}_3}{dt}=-\frac{N_3}{{\mathrm{\tau }}_{32}}-\frac{E_1}{�{\omega }_1}·\frac{d{P}_1}{dt},$$(1)$$\frac{d{N}_2}{dt}=\frac{N_3}{{\mathrm{\tau }}_{32}}+\frac{E_2}{�{\omega }_2}·\frac{d{P}_2}{dt}-\frac{N_2}{{\mathrm{\tau }}_{21}},$$(2)$$\frac{d{N}_1}{dt}=\frac{N_2}{{\mathrm{\tau }}_{21}}-\frac{E_2}{�{\omega }_2}·\frac{d{P}_2}{dt}-\frac{N_1}{{\mathrm{\tau }}_{10}},$$(3)$$\frac{d{N}_0}{dt}=\frac{N_1}{{\mathrm{\tau }}_{10}}+\frac{E_1}{\hslash {\omega }_1}·\frac{d{P}_1}{dt},$$(4)where $N_i$ is the population density in Level $i$, $i=0–3$, $P_i$ ($i=1–2$) is the polarization density, and $\hslash =h/2\pi$. The lifetime of the energy levels are chosen as ${\mathrm{\tau }}_{10}=5\times {10}^{-12}\mathrm{s}$, ${\mathrm{\tau }}_{21}={10}^{-10}\mathrm{s}$, and ${\mathrm{\tau }}_{32}={10}^{-13}\mathrm{s}$. The stimulated transition rate is given by the term $(E_2/\hslash {\omega }_2)·(d{P}_2/dt)$, where ${\omega }_2$ is the transition frequency between Levels 1 and 2 and is chosen as ${\omega }_2=2\mathrm{\pi }\times 5.13\times {10}^{14}\mathrm{Hz}$ (${\lambda }_2=585\mathrm{nm}$). Distinct from previous work, we use the term $(E_1/\hslash {\omega }_1)·(d{P}_1/dt)$ in the rate equations to describe the excitation process instead of the pumping rate $W_p$, where ${\omega }_1$ is the center frequency of the excitation light and chosen as ${\omega }_1=2\mathrm{\pi }\times 5.64\times {10}^{14}\mathrm{Hz}$ (${\lambda }_1=532\mathrm{nm}$). By coupling the polarization equation and rate equations, the excitation pulse can be described with an individual electric field intensity. Therefore, the electrical field of the pumping pulse $E_1$ and that of the emission pulse $E_2$ were represented separately in the random medium. We considered a 2D transverse magnetic (TM) field in the $xy$–plane; thus the Maxwell’s equations read as ($i=\mathrm{1,2}$) $${\mu }_0\frac{\partial H_{ix}}{\partial t}=-\frac{\partial E_{iz}}{\partial y},$$(5)$${\mu }_0\frac{\partial H_{iy}}{\partial t}=\frac{\partial E_{iz}}{\partial x},$$(6)$${\epsilon }_j{\epsilon }_0\frac{\partial E_{iz}}{\partial t}+\frac{\partial P_{iz}}{\partial t}=\frac{\partial H_{iy}}{\partial x}-\frac{\partial H_{ix}}{\partial y}.$$(7)The polarization density $P_i(i=\mathrm{1,2})$ obeys the following $$\frac{d^2{P}_1}{d{t}^2}+\mathrm{\Delta }{\omega }_1\frac{d{P}_1}{dt}+{\omega }_1^2{P}_1={\mathrm{\kappa }}_1(N_3-N_0)E_1,$$(8)$$\frac{d^2{P}_2}{d{t}^2}+\mathrm{\Delta }{\omega }_2\frac{d{P}_2}{dt}+{\omega }_2^2{P}_2={\mathrm{\kappa }}_2(N_1-N_2)E_2,$$(9)where $\mathrm{\Delta }{\omega }_1=6.67\times {10}^{12}\mathrm{Hz}$ is the linewidth of the excitation light, $\mathrm{\Delta }{\omega }_2=1/{\mathrm{\tau }}_{21}+2/T_2=5\times {10}^{13}\mathrm{Hz}$ is the linewidth of the atomic transition, $T_2=2\times {10}^{-14}\mathrm{s}$ is the collision time, the constant ${\mathrm{\kappa }}_2$ is given by ${\mathrm{\kappa }}_2=6\mathrm{\pi }{\epsilon }_0{c}^3/{\omega }_2^2{\mathrm{\tau }}_{21}$. The constant ${\mathrm{\kappa }}_i(i=1,2)$ is given by ${\mathrm{\kappa }}_1=6\mathrm{\pi }{\epsilon }_0{c}^3/{\omega }_1^2{\mathrm{\tau }}_{32}$ and ${\mathrm{\kappa }}_2=6\mathrm{\pi }{\epsilon }_0{c}^3/{\omega }_2^2{\mathrm{\tau }}_{21}$, respectively. By coupling the electric field of the excitation pulse and population density in atomic levels through the polarization density equation, we can describe the pumping process by the external excitation pulse, once the excitation ended, and the pumping process stopped. Amplification occurs when the external pumping mechanism produces population inversion.

The time function of the excitation pulse is introduced as a Gaussian $$E_1(t)=E_{\text{peak}}\exp (-\frac{4\mathrm{\pi }(t-t_0{)}^2}{{\mathrm{\tau }}^2}),$$(10)where $\mathrm{\tau }$ is the width of the Gaussian pulse, $E_{\text{peak}}$ is the peak value of the pumping, and $t_0$ is the time corresponding to the peak value.

The electromagnetic fields in the 2D active random medium can be obtained by finite-difference time domain (FDTD) methods. A perfectly matched layer (PML) is used to model an open system as absorbing conditions.

In what follows we will analyze the temporal properties of the emission light under different sample parameters, including the surface-filling fraction, scatterer size, and refractive index of the scatterer. Before this, we would like to distinguish two concepts; ultrashort pulse propagation and excitation. In the early work by Conti and coworkers^[12], they studies the sample with a distribution of spherical scatterers obtained by molecular dynamics (MD) simulations. A Gaussian transverse electromagnetic linearly $y$-polarized input pulse was impinged on the $xy$–face at normal incidence. By solving the Maxwell equations, the transmission behavior of the output pulse can be investigated, such as the spectrum, decay time, and energy velocity. In our work, a ultrashort pulse with TM polarization was introduced onto the 2D sample as the excitation source. The input light was scattered by the particles and form a spatially nonuniform gain. We focus on the temporal properties of the stimulated radiation process induced by this external pumping pulse. For each set of FDTD simulations, the electrical field intensity $E_2$ was integrated with respect to the entire $xy$–plane, and then the temporal evolution curve of the emission pulse can be obtained. We will discuss the laser action in the time domain with different system parameters, including the surface-filling fraction, refractive index, and scatterer size.

In a 2D random system, the density of the scatterers can be described by the surface-filling fraction $\mathrm{\Phi }$, which plays an important role in light localization; once $\mathrm{\Phi }$ exceed a critical value, the transition from diffusion to a localization state occurs. The parameters of the excitation pulse are chosen as follows: $\mathrm{\tau }=200\mathrm{fs}$, $E_{\text{peak}}=1\times {10}^9\mathrm{V}/\mathrm{m}$, and $t_0=200\mathrm{fs}$. This pumping pulse will be used in the following simulation work. Here we select a 2D sample with $S=5\times 5{\mathrm{\mu m}}^2$, $r=60\mathrm{nm}$, $n_1=1.4$, and $n_2=2.7$. For the pulse pumping case, the stimulated emission exhibits a pulse envelope with peak intensity $I_{\text{peak}}$ and width ${\mathrm{\tau }}_E$ following full-width at half-maximum (FWHM). The delay time $\mathrm{\Delta }T$ is defined as the interval between the excitation pulse and the emission pulse (Fig. 1, inset).

We start from the surface-filling fraction with ${\mathrm{\Phi }}_1=35\%$ and the result is shown in Fig. 1. As can be seen clearly, the peak intensity of the emission pulse is quite low and the pulse envelope exhibits a relaxation oscillation behavior at the trailing edge. As $\mathrm{\Phi }$ increases gradually, the peak intensity becomes higher and the pulse width becomes narrower, and at the same time the relaxation oscillation behavior disappears. The delay time $\mathrm{\Delta }T$ decreases as $\mathrm{\Phi }$ increases. The corresponding dependence of the peak intensity, pulse width, and delay time of the emission pulse on the surface-filling fraction is present in Fig. 2. As aforementioned, the system will transform from diffusive to a localization state as the surface-filling fraction exceeds a critical value. When ${\mathrm{\Phi }}_1=35\%$, the system is still in a diffusive case and the multiple scattering process is not quite intense; this leads to a stimulated radiation along with relaxation oscillation behavior. If the surface-filling fraction increases further, a more intensive multiple scattering occurs among the system and results in a strong stimulated radiation. The transport mean free path $l_t$ can be estimated by the formula $l_t=1/\rho ·{\sigma }_s$, where $\rho$ is the scatterer concentration and ${\sigma }_s$ is the effective scattering cross sectional area. For our 2D sample, the scatterer concentration can be represented by the surface-filling fraction $\mathrm{\Phi }$. The transport mean free path $l_t$ decreases with increasing surface-filling fraction $\mathrm{\Phi }$, which will lead to a more intensive light scattering process and form a localized state. Therefore, as $\mathrm{\Phi }$ increases, the stimulated emission occurs more quickly and lasts a shorter time, and simultaneously the residual oscillations disappear.

Light localization comes from the interacting between the photons and scattering particles. In previous work, Wiersma^[19] demonstrated that the emission wavelength can be controlled by changing the particle size. This indicates that scatterer radius may change the multiple scattering processes. In what follows, we will investigate how the scatterer size influences the light localization level, which is represented by the temporal properties of the emission pulse. As discussed previously, the surface-filling fraction plays an important role in light localization. In order to eliminate this influence, we choose three different scatterer radii ($r_1=50\mathrm{nm}$, $r_2=60\mathrm{nm}$, and $r_3=70\mathrm{nm}$) and keep the surface-filling fraction constant ($\mathrm{\Phi }=40\%$). A 2D sample with $S=5\times 5{\mathrm{\mu m}}^2$, $n_1=1.4$, and $n_2=2.7$ is used here.

As can be seen in Fig. 3, when the scatterer radius $r_1=50\mathrm{nm}$, the peak intensity and width of the emission pulse are about $I_{\text{peak}1}=1.17\times {10}^{10}(\text{a.u.})$ and ${\mathrm{\tau }}_{E1}=1.83\mathrm{ps}$, respectively, with a delay time $\mathrm{\Delta }{T}_1=8.33\mathrm{ps}$. When the radius is increased to $r_2=60\mathrm{nm}$, the peak intensity increases to $I_{\text{peak}2}=1.92\times {10}^{10}(\text{a.u.})$ and the pulse width and delay time decrease to ${\mathrm{\tau }}_{E2}=1.55\mathrm{ps}$ and $\mathrm{\Delta }{T}_2=5.54\mathrm{ps}$, respectively. As the radius increases further to $r_3=70\mathrm{nm}$, a more intense and shorter pulse emerges. Light localization can only take place in optical materials that are extremely strongly scattering. It requires the transport mean free path $l_t$ and wavevector $\mathbf{k}$ to satisfy $\mathbf{k}·l_t\le 1$, which is known as the Ioffe–Regel criterion^[7]. According to the Mie scattering theory, the transport mean free path $l_t$ is inversely proportional to the effective scattering cross sectional area ${\sigma }_s$, whereas ${\sigma }_s$ increases as the particle size becomes larger. This means that the transport mean free path $l_t$ decreases as the scatterer radius increases. As aformentioned, a shorter $l_t$ will lead to a more strongly multiple scattering and higher localized state, which will produce a sharper emission pulse.

Light scattering comes from the contrast of the refractive index between the particles and background. With the FDTD method in this work, we will study the multiple scattering of light and the stimulated emission behavior in 2D disordered system for a varying scattering strength, as obtained by changing the scatterer refractive index $n_2$. We then characterize the stimulated emission behavior by the pulse parameters including the peak intensity, pulse width, and delay time. The refractive index of the scatterers $n_2$ is chosen between 1.5 and 2.2 in accordance with some previous works^[12].

Figure 4 displays the emission light intensity versus time with the scatterer refractive index $n_2$ as a parameter. It shows that the refractive index has an effect on the temporal properties of the emission pulse. We start from the weak scattering regime for $n_2=1.5$; a weak pulse appears along with some residual oscillations, which indicates that the stimulated radiation process is not quite strong. As the refractive index increases to $n_2=1.8$, the trailing edge of the emission pulse becomes smooth and a Gaussian shape envelope of the pulse emerges. As the refractive index increase further, the peak intensity becomes stronger and the pulse width becomes narrower, along with a shorter delay time. Results demonstrate that the stimulated radiation becomes stronger as the multiple scattering strength increases, which comes from the increasing of the scatterer refractive index. This has been discussed by Vanneste^[20] in previous work. He investigated the lasing action in active random media in the weak and strong scattering regime by changing the refractive index of the scatterers. When the refractive index is low, lasing modes primarily consist of traveling waves and are spatially extended over the entire system. If the refractive index increases further, lasing modes transit to the localized regime. When the localized state forms, the stimulated emission process becomes stronger and leads to a more intense emission pulse.

In conclusion, the previous time-dependent theory of random lasers is extended to the ultrashort pulse pumping regime by introducing an independent electric field to describe the excitation pulse. By use of the FDTD method, a numerical study on the pulse shape of laser emission in 2D active random media is performed. The results show that the temporal properties of the emission pulse depend strongly on the sample parameters, including the refractive index of the scatterer, surface-filling fraction, and particle size. As the surface-filling fraction and refractive index of the particles increase, the pulse width and delay time decrease together, whereas the peak intensity of the emission pulse increases. However, the variation trend of the emission pulse behaves in the opposite manner as the scatterer radius increases. When the system is in the weak scattering regime, the lasing pulse exhibits some residual oscillations at the trailing edge. This phenomenon will disappear while the system transits to the strong scattering regime by changing some sample parameters. All these results give us good guidance in relevant experiments, especially in controlling the lasing emission in the temporal regime under ultrashort pulse excitation.