Three-dimensional numerical model of the dynamics of photorefractive beam self-focusing in InP:Fe

We propose a time-dependent three-dimensional numerical model which successfully explains the complex spatiotemporal dynamic of photorefractive self-focusing of an infrared beam in InP:Fe. Characteristic behaviors are in good agreement with experiments previously reported. Intensity dependence and asymmetric photo-induced index proﬁle are explained. In addition self-trapping dynamics shows potential self-trapping response in the microsecond range.


Three-dimensional numerical model of the dynamics of photorefractive beam self-focusing in InP:Fe
We propose a time-dependent three-dimensional numerical model which successfully explains the complex spatiotemporal dynamic of photorefractive self-focusing of an infrared beam in InP:Fe. Characteristic behaviors are in good agreement with experiments previously reported. Intensity dependence and asymmetric photoinduced index profile are explained. In addition self-trapping dynamics shows potential self-trapping response in the microsecond range. DOI: 10.1103/PhysRevA.79.033823 PACS number͑s͒: 42.65.Jx, 42.65.Hw, 42.65.Tg

I. INTRODUCTION
Development of all optical routing and interconnections devices is one of the most challenging objectives in optical communications. Self-trapping of optical beams and soliton interactions in photorefractive ͑PR͒ semiconductor materials at telecommunication wavelengths offer solutions to design such devices. Indeed with proper doping, semiconductors such as GaAs, CdTe, or InP exhibit photorefractive properties in the near infrared with fast response at below milliwatt power level. However due to weak electro-optic coefficients enhanced PR properties are necessary to efficiently route a signal beam with self-induced waveguiding in semiconductors. Impressive experimental results have been reported in CdTe ͓1͔ and promising demonstrations have been reported in InP:Fe.
Since photorefractive effect has been observed in InP:Fe in the late 1980's ͓2͔, several attempts have been made to exploit this property. For instance, industrial inspection by ultrasonic motion detection was demonstrated based on In-P:Fe sensitivity and fast response ͓3͔ while tunable filters for telecommunications were developed thanks to its near infrared capability ͓4͔.
While the latter studies are based on two-wave mixing ͑TWM͒, beam self-trapping in InP:Fe also possesses attractive properties. Indeed, experiments exhibiting some interesting features such as anisotropic focusing of light, deviation of the focused beam, or intensity controlled switching from focusing to defocusing with low power continuous wave ͑cw͒ one-dimensional ͑1D͒ or two-dimensional ͑2D͒ light beams have been demonstrated ͓5,6͔. First characterization of buildup time on the order of microseconds of the selffocusing phenomenon has also been recently reported ͓7͔.
Few theoretical 1D models have been proposed to interpret both temporal and steady-state regime ͓8,9͔ but none give satisfying explanations of the richness of the phenomenon. Modeling complexity is linked to the influence of both electrons and holes and the problem is even more challenging when two transverse dimensions are considered.
In a recent paper we presented a time-dependent threedimensional ͑3D͒ PR numerical model that was successfully applied to describe large beam bending observed during formation of self-trapped beam in photonic grade LiNbO 3 ͓10͔. In addition this time-dependent 3D model was used to demonstrate the influence of the anisotropic tensorial electrooptical effect on propagation of a vortex beam in LiNbO 3 :Fe ͓11͔. In this work, the above mentioned one-carrier model has been refined in order to be applied to materials for which both electrons and holes compete to give PR effect. This two-carrier model successfully explains the complex dynamic of the PR effect in InP:Fe and gives results in good agreement with previously reported works.

II. THEORETICAL BACKGROUND AND NUMERICAL MODEL
From the Kukhtarev's band transport model ͓12͔, the usual system of equations to solve the space-charge field distribution E ជ for an arbitrary light intensity distribution I for a two-carriers single deep level photorefractive medium is given by the following equations ͓13,14͔: ‫ץ‬n T ‫ץ‬t = e p p T − e n n T − c p pn T + c n np T , ͑1g͒ * fabrice.devaux@univ-fcomte.fr; where n, p, n T , and p T represent, respectively, the densities of electrons, holes, ionized, and neutral deep iron traps. , J ជ n , and J ជ p define, respectively, the space charge, electrons, and holes current densities. These variables are functions of time and space. N D , N A are constants that represent, respectively, the shallow donor and acceptor densities. ͓͔ is the static dielectric tensor. n and p represent, respectively, the electron and hole mobilities. k B is the Boltzmann constant, T is the temperature, and e is the elementary charge of electron. c n and c p are the recombination rates for electrons and holes, respectively. e n and e p give the emission parameters for electrons and holes and are defined as e n = e n th + n ͑I + I d ͒, ͑2a͒ Where e n th and e p th are the thermal emission rates, n and p are the photoionization cross sections. I is the beam intensity distribution and I d is the background illumination intensity.
The peculiar photorefractive behavior of biased InP:Fe was revealed through TWM experiments ͓2͔ showing an enhanced space-charge field for a given intensity. The model developed later ͓14͔ showed the importance of electrons and holes competition in this resonant process. First physical explanations for beam self-trapping in InP have been naturally derived from TWM theory. However contrary to TWM, beam self-trapping implies localized illumination, large index modulation, and above all 2D intensity distribution so that many approximations valid for TWM are not suitable for illumination with narrow GAUSSIAN beams. 1D analytical models were proposed but under restrictive approximations ͓8,9͔ which describe partially the behaviors observed experimentally. We instead propose a numerical model that takes into account the two transverse dimensions for spatial derivation operators. Thermal diffusion is not neglected, thermal and optical excitations of electrons and holes are both considered. In addition the presented model includes the dynamic of the PR process.
To emphasize possible applications for telecommunications and to facilitate comparison with experiments described in Ref. ͓7͔, we consider propagation along the direction ͓110͔͑ i.e., z axis͒ of a 6 mm long InP:Fe crystal, of an cw infrared TEM 00 GAUSSIAN beam at 1.55 m focused to a 25 m beam waist at the input face of the crystal for peak intensities ranging from 5 to 10 4 W / m 2 . An electric field E 0 of 10 6 Vm −1 is applied in the ͓001͔ direction ͑i.e., y axis͒. Laser beam is polarized along the ͓110͔ direction ͑i.e., x axis͒. All calculations are performed with a 10 23 m −3 iron deep level concentration ͑N T ͒ among which 5% are ionized due to presence of shallow levels. Other parameters of the crystal taken in the literature are listed in Table I.
In order to obtain evolution in time and space of the charge density for a given intensity distribution, Eqs. ͑1e͒-͑1g͒ are iteratively solved starting from initial condi-tions with a time step ⌬t. Next step of the procedure is to deduce the electric field distribution produced by the newly calculated space-charge distribution.
Instead of using Poisson equation ͓Eq. ͑1͔͒, which is not straightforward to solve, we compute it from Eq. ͑3͒ by using three-dimensional discrete Fourier transform. Note that Eq. ͑1͒ is only solved in two transverse dimensions since variation of PR variables along the propagation axis are weak. However Eq. ͑3͒, which is a key feature in our model, is solved in three dimensions in order to take into account contribution of the whole space charge to the space-charge electric field formation. It thus gives all the components of the electrical field produced by charges distribution ͑r ជ͒ in the medium volume V ͑dV being an elementary volume͒. Then index perturbation induced by the electro-optical effect is calculated in the whole crystal volume. When transverse components of the space-charge electric field ͑E x , E y ͒ are both considered in electro-optic tensor, eigenvalue that gives index perturbation for light polarized along the ͓110͔ direction depends on E x and E y . Nevertheless, contribution of E x component to index perturbation is minor and can be neglected. Then, index perturbation for light polarized along the ͓110͔ direction is given by: ⌬n Ӎ − 1 2 n 0 3 r 41 E y , where n 0 is the refractive index, r 41 is the electro-optical coefficient, and E y is the component of the electric field along the ͓001͔ direction. While E x is neglected to calculate index perturbation, its contribution for charge transport ͓in Eqs. ͑1e͒ and ͑1f͔͒ is considered. Propagation of light in the perturbed medium obeys the following paraxial wave equation: where ⌬ Ќ = ‫ץ‬ 2 ‫ץ‬x 2 + ‫ץ‬ 2 ‫ץ‬y 2 , A͑r ជ͒ is the slowly varying amplitude of the beam and k the wave number. This equation is solved by a classical split-step-Fourier transform method that gives the new light intensity distribution. The whole procedure is repeated until steady-state regime is reached.

III. RESULTS
In order to understand the complex dynamic of self focusing in InP:Fe, we first launched a typical beam peak intensity of 200 W / m 2 with an intensity to dark irradiance ratio of 2. As later described in the paper, this intensity leads to efficient trapping and is consequently appropriate to illustrate both the dynamic and the distribution of the PR variables.

A. Early stage
First, let us examine the different distributions, normalized to initial ionized Fe traps density n T0 , of holes, elec-trons, ionized Fe traps, and space-charge densities in the transverse plane at the exit face of the crystal at an early stage of the induction process that is few microseconds after beam is switched on ͑Fig. 1͒. At the beginning of the process, because of electron and hole large mobilities, optically excited free charges drift in opposite direction along the applied field and move significantly far away from the beam ͓Figs. 1͑a͒ and 1͑b͔͒. Because of a larger optical cross section and for chosen n T0 p T0 ratio, holes are dominant and recombine on the upper side of the beam giving a lower n T concentration at that place at the expense of the central illuminated region ͓Fig. 1͑c͔͒. It thus leads to an asymmetric space-charge density ͓Fig. 1͑d͔͒. Electrons also contribute to the asymmetric space-charge distribution but to a less extent. Figure 2 gives the distribution of the two transverse space-charge field components E x and E y normalized to E 0 induced by the space-charge distribution depicted in Fig.  1͑d͒. Since refractive index change due to electro-optic effect is proportional to E y , it gives an asymmetric index modulation along y axis. We can note also that the regions where E y Ͻ E 0 ͑i.e., higher index area͒ and E y Ͼ E 0 ͑i.e., lower index area͒ are shifted with respect to the maximum intensity. These two observations are precursors of both the asymmetric focusing and bending of the beam as observed experimentally ͓5-7,9͔.

B. Dynamics of PR effect up to steady-state regime
After this initial stage the photorefractive process continues to rise. The temporal dynamics of the space-charge density, the electric field component E y , the index modulation, and the beam intensity along the y-axis are given respectively in Figs. 3͑a͒,3͑c͒,3͑e͒,and 3͑g͒. Figures 3͑b͒,3͑d͒,3͑f͒,and 3͑h͒ show the corresponding steady-state distributions in the transverse plane at the crystal output. The spacecharge density gradually increases, which induces a strong modulation depth of the electric field ͑i.e., difference between minimum and maximum values of the electric field͒. As a consequence, holes and electrons current densities are significantly perturbed by this inhomogeneous electric field. This perturbation leads to an accumulation of densely packed negative charges on the lower edge of the beam ͓Fig. 3͑b͔͒.It gives an electric field approaching zero in the upper region and an electric field increasing dramatically up to about three times the applied field E 0 on the lower edge of the beam ͓Fig. 3͑d͔͒. As a consequence, index modulation is strongly asymmetric and is composed of an abrupt index decay on one side of the beam and a larger and smoother region where the index is higher ͓Fig. 3͑f͔͒. Focusing occurs simultaneously with appearance of a beam displacement in the direction of E ជ 0 . It thus gives an asymmetric beam profile along y axis as it was clearly observed in experiments reported in ͓5,6͔. At last we would like to emphasize that the present model gives for the first time some insight on the 2D distribution of the index change. It notably reveals a beam-shaped low index region that surrounds one side of the beam ͓Fig. 3͑f͔͒. Figure 4 presents the beam barycenter displacement ͓15͔ along y axis and the focusing ratio defined as the ratio of the beam width along y axis in nonlinear regime over the width in linear regime at the output face. At the chosen intensity, an efficient self-trapping is thus observed with this short length crystal. The model reveals that a space-charge field considerably larger than the applied field is formed, which induces an index modulation strong enough to modify the beam propagation despite the low electro-optic coefficient of InP. Note that we consider the barycenter of the beam to evaluate the beam distortion and shift. It leads to a rather modest shift value compare to larger shifts experimentally observed ͓5,6͔. Nevertheless, we would like to emphasize that the model additionally shows that displacement of the beam peak intensity is very sensitive to parameters of the crystal such as doping concentration, ionization rate of iron traps, dark irradiance and, obviously, crystal length.

C. Influence of beam peak intensity
To complement the above development where a given beam intensity of 200 W / m 2 was considered, we performed additional calculations up to the steady-state regime for peak beam intensities ranging from 5 W to 10 kW/ m 2 . Figure 5 shows the focusing ratio and the beam displacement as a function of the input beam intensity. In the steady-state regime, an efficient focusing is obtained for intensities from 50 W to 1.5 kW/ m 2 , while a small defocusing effect for intensities lower than 30 W / m 2 and no focusing for intensities greater than 5 kW/ m 2 are observed. For reverse applied field a weak self-focusing is now predicted for weak intensities while defocusing is present for more intense beam intensities. This behavior showing that steady-state focusing occurs only for an intensity range is in accordance with reported experiments ͓5,6͔.
As depicted in Fig. 6͑a͒, the dynamic of self-focusing computed for three different typical intensities also reveals interesting features suggesting that fast self-focusing with a response time in the microsecond range effect should be observed in InP:Fe. In Fig. 6͑b͒ the corresponding normalized E y electric field profiles along y axis are plotted for steadystate regime. Even if the modulation of the electric field ͑i.e., index modulation͒ is small for low intensity, an area where E y Ͼ E 0 ͑i.e., low index͒ is located in the illuminated region which leads to small defocusing. Note that a higher index region is present on the right-hand side of the beam. When intensity raises, the modulation amplitude of the electric field increases and gives a strong asymmetry between the high and low index regions leading to self-focusing and displacement of the beam. While this index distribution gives focusing in steady-state regime for intensity ranging from 50 W to 1.5 kW/ m 2 , higher intensities tend to shift the low index region further from the beam while high index region widens. High intensity beams consequently do not exhibit selffocusing behavior in steady-state regime. However, a strong index change is expected to be present at some distance from the beam as depicted in Fig. 6͑b͒  where a local electric field approaching four times the applied field is predicted. It can be noted that this distribution has many similarities with the one calculated in the frame of a 1D model for periodic illumination with high contrast in Ref. ͓16͔ if we consider only one fringe.

IV. DISCUSSION AND CONCLUSION
This study enlightens the temporal and spatial properties of beam self-focusing in InP:Fe where both holes and electrons are involved in the PR effect. The 3D numerical model explains remarkably well several features previously ob-served experimentally but with an underlying physics unsuccessfully clarified until now. Furthermore, thanks to the timedependent model used, the whole complex spatiotemporal dynamic of the PR effect that leads to anisotropic focusing and bending of the beam is depicted up to the steady-state regime. Unexpected features of the transient regime are revealed especially at high intensity, which open up the possibility to induce fast self-trapping with microsecond response time. Experiments are in progress to visualize this fast dynamic. Because some crystal parameters of the crystal cannot be measured precisely, quantitative comparison between numerical results and experiments remains difficult. For in- stance, the intensity range for which focusing occurs depends on several parameters such as iron concentration or optical cross section. The presented numerical model can be usefully applied to extensively explore the influence of different parameters of the material. While in our model only the E y component of the electric field is considered for electro-optical effect, the magnitude of the E x component is sufficiently large to influence beam propagation in specific electro-optic configurations which should be studied. Finally this model can be easily adapted to other semiconductor materials either with one or two free carriers involved.