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

F. Devaux, M. Chauvet

► To cite this version:

F. Devaux, M. Chauvet. Three-dimensional numerical model of the dynamics of photorefractive beam self-focusing in InP:Fe. Physical Review A, 2009, 79 (3), pp.033823. 10.1103/PhysRevA.79.033823. hal-00441882

HAL Id: hal-00441882

https://hal.science/hal-00441882v1

Submitted on 23 Apr 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Three-dimensional numerical model of the dynamics of photorefractive beam
self-focusing in InP:Fe

Fabrice Devaux* and Mathieu Chauvet

Département d'Optique P.M. Duffieux, Institut FEMTO-ST, UMR CNRS 6174, Université de Franche-Comté,

16 route de Gray, F-25030 Besançon Cedex, France

(Received 2 October 2008; published 17 March 2009)

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 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 proper-
ties 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 semiconduc-
tors. 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 infra-
red 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 self-focusing 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 three-
dimensional (3D) PR numerical model that was successfully
applied to describe large beam bending observed during for-
mation of self-trapped beam in photonic grade LiNbO₃ [10].
In addition this time-dependent 3D model was used to dem-
onstrate the influence of the anisotropic tensorial electro-
optical effect on propagation of a vortex beam in LiNbO₃: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 dy-
namic 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]:

$$
\vec{\nabla} \cdot \{[\varepsilon]\vec{E}\} = \rho, \qquad (1a)
$$

$$
\rho = e(N_D - N_A + p - n - n_T), \quad (1b)
$$

$$
\vec{J}_n = e \mu_n n \vec{E} + \mu_n k_B T \vec{\nabla} n, \quad (1c)
$$

$$
\vec{J}_p = e \mu_p p \vec{E} + \mu_p k_B T \vec{\nabla} p, \quad (1d)
$$

$$
\frac{\partial n}{\partial t} = e_n n_T - c_n n p_T + \frac{1}{e} \vec{\nabla} \cdot \vec{J}_n, \quad (1e)
$$

$$
\frac{\partial p}{\partial t} = e_p p_T - c_p p n_T - \frac{1}{e} \vec{\nabla} \cdot \vec{J}_p, \quad (1f)
$$

$$
\frac{\partial n_T}{\partial t} = e_p p_T - e_n n_T - c_p p n_T + c_n n p_T, \quad (1g)
$$

<footnote>*fabrice.devaux@univ-fcomte.fr; http://www.femto-st.fr/fr/
Departements-de-recherche/OPTIQUE/</footnote>
$$N_T = n_T + p_T, \qquad (1h)$$

where $n$, $p$, $n_T$, and $p_T$ represent, respectively, the densities of electrons, holes, ionized, and neutral deep iron traps. $\rho$, $\tilde{J}_n$, and $\tilde{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. $[\varepsilon]$ is the static dielectric tensor. $\mu_n$ and $\mu_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} + \sigma_n(I + I_d), \qquad (2a)$$

$$e_p = e_p^{th} + \sigma_p(I + I_d). \qquad (2b)$$

Where $e_n^{th}$ and $e_p^{th}$ are the thermal emission rates, $\sigma_n$ and $\sigma_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₀₀ 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⁴ W/m². An electric field $E_0$ of 10⁶ V m⁻¹ 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²³ m⁻³ 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-

<table><caption>TABLE I. Crystal parameters.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>μ<sub>n</sub>(m²/V s)</td><td>0.4</td></tr><tr><td>μ<sub>p</sub></td><td>0.015</td></tr><tr><td>e<sup>th</sup><sub>n</sub>(s⁻¹)</td><td>16.32</td></tr><tr><td>e<sup>th</sup><sub>p</sub></td><td>10<sup>-4</sup></td></tr><tr><td>σ<sub>n</sub>(m²/J) @ 1.55 μm</td><td>3.1 × 10<sup>-4</sup></td></tr><tr><td>σ<sub>p</sub></td><td>7.8 × 10<sup>-4</sup></td></tr><tr><td>c<sub>n</sub>(m³/s)</td><td>4.1 × 10<sup>-14</sup></td></tr><tr><td>c<sub>p</sub></td><td>1.6 × 10<sup>-14</sup></td></tr><tr><td>r<sub>41</sub>(pm/V)</td><td>1.68</td></tr><tr><td>ε<sub>r</sub></td><td>12.5</td></tr></tbody></table>

tions with a time step $\Delta t$. Next step of the procedure is to deduce the electric field distribution produced by the newly calculated space-charge distribution.

$$\vec{E}(\vec{r}) = \vec{E}_0 + \frac{1}{4\pi[\varepsilon]} \iint_V \rho(\vec{r}') \frac{\vec{r} - \vec{r}'}{|\vec{r} - \vec{r}'|^3} dV. \qquad (3)$$

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 $\rho(\vec{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: $\Delta n = -\frac{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:

$$\left\{ \frac{\partial}{\partial z} - \frac{i}{2k} \Delta_{\perp} \right\} A(\vec{r}) = i \frac{k}{n_0} \Delta n(\vec{r}) A(\vec{r}), \qquad (4)$$

where $\Delta_{\perp} = \frac{\partial^2}{\partial z^2} + \frac{\partial^2}{\partial y^2}$, $A(\vec{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.
<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color) Distributions of (a) normalized holes, (b) electrons, (c) relative ionized Fe traps, and (d) space-charge densities in the transverse plane at the exit face of the crystal at early stage of the PR process. Black circles: beam location.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. (Color) (a) $E_y$ and (b) $E_x$ space-charge field components normalized to $E_0$ induced by the space-charge distribution depicted in Fig. 1(d). Black circles: beam location.</figcaption></figure>
### 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² 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 $\frac{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

<figure><img src="image_4.png" /><figcaption>FIG. 3. (Color) (a) Temporal dynamics of the space-charge density, (c) $E_y$ component of the electric field, (e) index modulation, and (g) beam intensity along the y axis and respective distributions at steady-state regime in the transverse plane at the crystal output.</figcaption></figure>
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 space-charge 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 $\vec{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² was considered, we performed additional calculations up to the steady-state regime for peak beam intensities ranging from 5 W to 10 kW/m². 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², while a small defocusing effect for intensities lower than 30 W/m² and no focusing for intensities greater than 5 kW/m² 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 steady-state 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², 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 self-focusing 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) for 10 kW/m² intensity

<figure><img src="image_5.png" /><figcaption>FIG. 4. Displacement of the beam barycenter along y axis (dashed line) and focusing ratio (full line) as a function of time for an launched intensity of 200 W/m².</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 5. Focusing ratio and beam barycenter displacement as a function of the input beam intensity. Dotted and dotted-dashed horizontal lines represent, respectively, focusing ratios for a soliton regime and free diffraction.</figcaption></figure>

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 time-dependent 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-

<figure><img src="image_7.png" /><figcaption>FIG. 6. (a) Temporal dynamic of the focusing ratio for three different typical intensities. (b) Corresponding variation of the normalized electric field component $E_y$ along the y axis at steady-state regime.</figcaption></figure>
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.

[1] T. Schwartz, Y. Ganor, T. Carmon, R. Uzdin, S. Shwartz, and M. Segev, Opt. Lett. **27**, 1229 (2002).

[2] B. Mainguet, Opt. Lett. **13**, 657 (1988).

[3] P. Delaye, A. Blouin, D. Drolet, L. de Montmorillon, G. Roosen, and J.-P. Monchalin, J. Opt. Soc. Am. B **14**, 1723 (1997).

[4] E. Vourc'h and D. Hervé, Opt. Lett. **28**, 1105 (2003).

[5] M. Chauvet, S. A. Hawkins, G. J. Salamo, M. Segev, D. F. Bliss, and G. Bryant, Opt. Lett. **21**, 1333 (1996).

[6] M. Chauvet, S. A. Hawkins, G. J. Salamo, M. Segev, D. F. Bliss, and G. Bryant, Appl. Phys. Lett. **70**, 2499 (1997).

[7] D. Wolfersberger, N. Khelfaoui, C. Dan, and N. Fressengeas, Appl. Phys. Lett. **92**, 021106 (2008).

[8] R. Uzdin, M. Segev, and G. Salamo, Opt. Lett. **26**, 1547 (2001).

[9] N. Fressengeas, N. Khelfaoui, C. Dan, and D. Wolfersberger, G.

Montemezzani, H. Leblond, and M. Chauvet, Phys. Rev. A **75**, 063834 (2007).

[10] F. Devaux, V. Coda, M. Chauvet, and R. Passier, J. Opt. Soc. Am. B **25**, 1081 (2008).

[11] R. Passier, F. Devaux, and M. Chauvet, Opt. Express **16**, 7134 (2008).

[12] N. Kukhtarev, V. Markov, S. Odulov, M. Soskin, and V. V networks, Ferroelectrics **22**, 949 (1979).

[13] J. Strohkendl, J. Jonathan, and R. Hellwarth, Opt. Lett. **11**, 312 (1986).

[14] G. Picoli, P. Gravey, C. Ozkul, and V. Vieux, J. Appl. Phys. **66**, 3798 (1989).

[15] M. R. Belić, A. Stepken, and F. Kaiser, Phys. Rev. Lett. **82**, 544 (1999).

[16] N. Wolffer, I. Gravey, and R. Coquillé, J. Appl. Phys. **78**, 6375 (1995).