Numerical method for an analysis of nonlinear light propagation in photorefractive media - time nonlocal approach

Andrzej Ziółkowski

doi:10.1364/OE.22.0A1907

1. Introduction

Nonlinear phenomena are unique in many respects. Although not always noticeable, these phenomena are present in many physical systems of the surrounding world. Consequently, nonlinear phenomena play an essential role in all disciplines of modern science. One measureable outcome of research into these phenomena is a number of interesting and complex optical effects that occur during interactions of light with nonlinear crystals [1].

Photorefractive media [2–4] seem to be exceptionally interesting in the broad range of optical nonlinear materials. They have drawn researchers' attention from both fundamental and practical viewpoint. The former viewpoint relates to a complex nature of physical processes inherent in the mechanism of photorefraction and a variety of media that support this mechanism [5]. The latter is connected with a wide spectrum of potential applications of photorefractive media [6], from dynamic holography to medicine and telecommunications. Both theoretical and practical investigations run parallel to the need for developing effective numerical tools, offering possibly most complete analysis of the photorefraction mechanism.

To date, the numerical approach most frequently used for photorefractive phenomena has comprised two areas. One is related to analyses of physical processes induced by an interference pattern falling onto photorefractive material aimed at creating in it a grating of refractive index changes [7–10]. Numerical tools used in this area are mostly based on the solution of a system of materials equations, whose form and complexity depend on the examined material. This approach proves particularly valuable in cases going beyond the so called approximation of linearization and in testing materials with such complex photorefractive mechanism that an analytical approach cannot be used.

The other area where numerical approach is justified is the research related to nonlinear light propagation. In this case, methods of beam propagation serve two purposes: as a cognitive tool and for the verification of analytical approach in examination of, for instance, photorefractive solitons. Works on this problem are mostly based on a wave equation solution that takes into account available analytical relations describing light-dependent distribution of refractive index changes [11,12]. However, we can find many interesting research themes in areas where an analytical form of relations describing changes in the refractive index is not available or available only if necessary analytical approximations are made. One such example is the analysis of nonlinear light propagation in complex photorefractive materials, such as semiconducting compounds or media with complex structure of defects or dopants. In this case the nonlinear term of the wave equation should be determined by numerical techniques. If the method used for the purpose is sufficiently general, the numerical procedure will deliver information on light propagation as well as time evolution of essential physical quantities such as concentrations of free carriers or ionized dopants [13,14].

Interesting, relatively general numerical approach, designed to nonlinear light propagation analysis in photorefractive media can be based on algorithm discussed in [15]. Practical use of this method is presented in [16], containing results of the investigation of the influence of hot electrons on the propagation of solitons on the photorefractive gallium arsenide. The algorithm offers wide applicability due to its accuracy, no need to use analytical approximations and relatively high flexibility of implementation. However, it may be improved, as it does not consider time nonlocalities, essential in soliton propagation analysis. Latest research works on the dynamics of light self-trapping in photorefractive materials [17,18] has shown that in typical experimental conditions the time evolution of a self-trapping beam is more complex than was thought before. The process of forming a space-charge field in a given cross-section of a photorefractive crystal (along the z axis of propagation) depends on the dynamics of processes taking place in the crystal in previous places (for z^’ < z). This phenomenon, affecting the time of reaching a steady state and the process of forming spatial solitons, should be taken into account in the numerical algorithm used for an analysis of nonlinear light propagation in photorefractive materials.

This article presents a modified version of the computing procedure described in the work [15]. The new approach to this numerical problem has been designed taking into account nonlocal time phenomena. In particular, the article discusses and compares three numerical procedures. The first one is time local numerical approach presented in [15]. Here, a local algorithm is perceived only a basis for understanding the method which refers to time nonlocality. The concept of modifying the local algorithm is rather simple. It involves reversing the order of separated integration in domains of time and space, that has been used in paper [15]. However, this can be done in two different ways, resulting in subtly different numerical procedures. The possibility to address the issue of time nonlocality in two ways is attributed to the fact that, as it will be shown, it is possible (but not required) to completely separate carrier dynamics evaluation and wave equation calculation. This in turn, allows to choose a short integration time for carrier dynamics and a longer one to solve the wave equation. Such an approach if made with properly chosen parameters, can be considered as an improved version of the nonlocal numerical procedure as it leads to a shorter computation time.

The paper is organized as follows. Section 2 contains the theoretical background that describes band transport model and its most commonly used simplifications with particular emphasis on the approximation of time locality. For illustrative purposes, like in [15], two types of media are chosen: one with typical photorefractive properties (strontium barium niobate − SBN), and a formally more complex semi-insulating gallium arsenide which, demonstrate bipolar conductivity and hot electrons effect. Section 3 describes numerical procedures that can be used to simulate nonlinear light propagation in the aforementioned media. Section 4 presents results of calculations made for the comparison between time-local and nonlocal approaches. Section 5 compares the computational times for the presented algorithms.

2. Theoretical background

2.1 Model of the phenomenon

In the simplest approach used in descriptions of photorefraction, e.g. in ferroelectric SBN crystal, free electrons are induced by light from the deep level of donors, that fully compensate shallow acceptors. Excited carriers move due to electric field or diffusion until they permanently recombine with traps in dark areas. Since the distribution of light intensity is not homogeneous, the distributions of created charge and associated electric field in a photorefractive material are not homogeneous, either. Thus modulated electric field causes modulation of the refractive index through the electro-optic effect. In the formal approach, these processes are described by material equations, which for media such as an SBN crystal have this form:

\frac{\partial N_{D}^{+}}{\partial t} = [β_{n} + S_{n} (I + I_{B})] (N_{D} - N_{D}^{+}) - γ_{n} n N_{D}^{+},

J_{n} = q μ_{n} n E + μ_{n} k_{B} T grad n,

\frac{\partial}{\partial t} (N_{D}^{+} - N_{A}^{-} - n) = \frac{1}{q} div J_{n},

div E = - div (grad φ) = \frac{q}{ε ε_{o}} (N_{D}^{+} - N_{A}^{-} - n),

where n describes the concentration of free electrons,

N_{D}, N_{D}^{+}

and

N_{A}^{-}

are, respectively, total concentration of donors, concentration of ionized donors and concentration of ionized acceptors. It is assumed that at a room temperature all shallow acceptors are ionized. Besides: β is the coefficient of thermal generation of carriers from deep traps, S is a cross-section for photoionization of traps divided by photon energy, γ is the coefficient for recombination to traps. The other quantities in the equations are as follows: J – current density, E – intensity of electric field in the medium, φ - electric potential, I – light intensity, I_B – homogeneous background lighting allowing to control medium conductivity, μ - mobility of carriers, k_B –Boltzmann's constant, T- temperature, q – elementary charge, ε₀ – permittivity of free space, ε - relative permittivity of a medium.

The system of Eqs. (1a)-(1d) is appropriate for the description of typical photorefractive materials, but in many cases it is insufficient. Semiconductors, for example, require more sophisticated approach. Models of photorefractive materials such as indium phosphide, gallium arsenide or multiple quantum well structures based on such material are usually unique. Much attention has been recently paid to nonlinear light interaction with photorefractive indium phosphide [19–24]. This material allows to amplify light waves interacting in two-wave and four-wave mixing processes and to generate bright spatial solitons. Related experimental research is confirmed by a theory that unlike model (1a) – (1d), takes into account bipolar transport, where usually one photocarrier and one thermal carrier are considered.

Bipolar transport is one of the semiconductor properties that make models of these materials complex. Materials such as gallium arsenide are characterized additionally by the phenomenon of hot electrons [25,26]. Gallium arsenide belongs to multi-valley semi-conductors, whose conductivity band consists of a central valley Γ and higher side valleys L, characterized by large effective mass. In the thermodynamic equilibrium, most electrons are found in the central valley. If a semiconductor is placed in an electric field, at low values the energy is conveyed to the lattice by the mechanism of intervalley scattering on optical phonons. At higher values of electric field, the mean energy of electrons increases, and the energy relaxation becomes less effective. 'Hot' electrons are formed, with a corresponding effective temperature higher than the lattice temperature T_L [25–28]:

T_{n} (E) = T_{L} + \frac{2 q τ_{r} v_{n} (E)}{3 k_{B}} E,

where v_n(E) is drift velocity of electrons that depends on the electric field, τ_r is mean time of energy relaxation, ranging from 0.1 ps to 1 ps [26,28]. Further increase of the electric field is related to intervalley scattering on mostly optical phonons. The considerable majority of electrons reaches energy that enables them to pass to a side valley, where they gain a large effective mass. In this way in lower energy states of the central valley the electron occupation decreases, while the occupation of states with higher energies in side valleys rises. Because electrons in side valleys are less mobile, once a critical value of the electric field is exceeded, the mean rate of electron transfer begins to decrease. In the description of this phenomenon, most often the mobility of electrons is given as a mean weighted mobility:

μ_{n} (E) = μ_{n l} f (E) + μ_{n u} [1 - f (E)],

where μ_nl and μ_nu are, respectively, mobilities of carriers in the central and side valleys, f(E) is a distribution function describing the occupation of the central valley, given by this relation:

f (E) = {1 + R \exp [\frac{- Δ U}{k_{B} T_{n} (E)}]}^{- 1},

where ΔU is the energy difference between the central and side minimum, R is the ratio of the state density in the central to that in side valleys. For gallium arsenide R = 94, ΔU = 0.31 eV. If we take into account bipolar transport and hot electron effect, then the materials equations describing the photorefractive mechanism will take this form [3,15,16,27,29,30]:

\frac{\partial n}{\partial t} - \frac{1}{q} div J_{n} = [S_{n} (I + I_{B}) + β_{n}] (N_{D} - N_{D}^{+}) - γ_{n} n N_{D}^{+},

\frac{\partial p}{\partial t} + \frac{1}{q} div J_{p} = [S_{p} (I + I_{B}) + β_{p}] N_{D}^{+} - γ_{p} p (N_{D} - N_{D}^{+}),

J_{n} = q μ_{n} (E) n E + k_{B} grad [μ_{n} (E) T_{n} (E) n],

J_{p} = q μ_{p} p E - k_{B} μ_{p} T_{L} grad p,

\frac{\partial}{\partial t} (n + N_{A}^{-} - p - N_{D}^{+}) = \frac{1}{q} div (J_{n} + J_{p}),

div E = - div (grad φ) = \frac{q}{ε ε_{o}} (N_{D}^{+} + p - n - N_{A}^{-}),

where p stands for hole concentration. Indexes n and p refer to quantities characterizing electrons and holes, respectively.

Equations (1a)-(1d) and (5a)-(5f) describe physical processes coupled with a light beam propagating in a photorefractive medium, whose evolution is described by a nonlinear wave equation. In the configuration (1 + 1)D, for the light polarized in the direction x and propagating in the direction z, the electric component of light field: $E_{o p t} (x, z, t) = A (x, z) \exp (i ω t - i k z),$ fulfills the wave equation, which in paraxial approximation gets this form:

(\frac{\partial}{\partial z} + \frac{i}{2 k} \frac{\partial^{2}}{\partial x^{2}}) A (x, z) + \frac{i k}{n_{b}} Δ n (E) A (x, z) = 0,

where

k = 2 π n_{b} / λ

is wave number,

ω = 2 π c / λ

is circular frequency, λ is wavelength in vacuum, and

Δ n

is optically induced change of refractive index

n_{b}

. Analyses of spatial soliton propagation in photorefractive materials mainly make use of analytical approximations from the system of Eqs. (1a-1b) or (5a-5f), resulting in a relation between the space-charge field and distribution of light intensity. Thus a numerical analysis [12,31,32], and in some cases analytical analysis [33,34] of soliton states resulting from the propagation Eq. (6) is possible. Of standard analyses, particular attention can be paid to the approach presented in works [17,18,35] allowing to consider photorefractive time nonlocality, typical of solitons.

2.2 Analytical approximations

To justify the use of numerical methods in nonlinear propagation analysis, let us take a closer look at the generally made transition from materials equations to the relation between spatial-charge field and light intensity for bright spatial soliton propagation in a semiconducting material characterized by Eqs. (5a)-(5f). This will allow us to separate all standard analytical approximations and point out differences between the time local and nonlocal approach.

One approximation that considerably restricts an analysis of photorefractive solitons is an approximation of the linear transport of electrons. In this approach, both electron mobility and temperature do not depend on the electric field intensity:

μ_{n} = μ_{n l}, T_{n} = T_{L} .

This approximation is right only in the linear range of the relation between electron velocity v_n and electric field intensity E. The relation holds only for field intensity below a critical one, which for GaAs is about 3.5 kV/cm. This restriction is essential, because semiconductors are characterized by relatively weak electro-optic properties. This, in turn, necessitates the application of higher values of external fields. It is important, because we overcome formal difficulties by using a numerical approach. Nonlinear transport of electrons leads to strong effects of spatial nonlocality that play an important role in nonlinear light propagation [16] and photorefractive processes of wave mixing [3,27,28].

Another significant simplifying step refers to the dynamics of photorefraction. An analysis of this problem distinguishes two characteristic time quantities: mean lifetimes of carriers and time of dielectric relaxation describing the rate of space-charge field formation. These quantities usually differ substantially from each other, as the lifetime of carriers is negligibly short compared to the time of dielectric relaxation. Examining the macroscopic dynamics of the photorefractive process in the time scale comparable to dielectric relaxation time, we assume that processes of generation and recombination remain in equilibrium. This, for typical photorefractive materials given by the system of Eqs. (1a)-(1d), is expressed in relation resulting from Eq. (1a): $[β_{n} + S_{n} (I + I_{B})] (N_{D} - N_{D}^{+}) ≃ γ_{n} n N_{D}^{+}$ . If we use a similar approach to a medium described by Eqs. (5a-5f), assuming a priori that the left-hand sides of Eqs. (5a) and (5b) may be omitted, we obtain stationary concentrations of electrons and holes expressed by the relation:

n = \frac{[S_{n} (I + I_{B}) + β_{n}] (N_{D} - N_{D}^{+})}{γ_{n} N_{D}^{+}}, p = \frac{[S_{p} (I + I_{B}) + β_{p}] N_{D}^{+}}{γ_{p} (N_{D} - N_{D}^{+})} .

Although Eqs. (7) substantially simplify a formal analysis, their form is subject to further simplifications. In a low light intensity range, carrier concentrations reach small values compared to concentrations of traps, which allows to neglect them in the Gauss law (5f). This fact combined with an approximation of slow-varying electric field permits to assume that concentrations of ionized acceptors and donors are quasi-equal $N_{A}^{-} \approx N_{D}$ . As a result, we obtain carrier concentration distributions that depend linearly on the distribution of light intensity. Approximations leading to a simple linear relation between carrier concentrations and light intensity are rather weak, e.g. for media characterized by large drift lengths of carriers. Time evolution of free carriers in such cases obviously runs in two stages: first, comparable to lifetime of carriers, and second, longer stage comparable to the time of electric field formation [29]. For materials such as GaAs this assumption introduces a significant restriction, even if nonlocal phenomena are not considered.

Another frequently used simplification is when the process of diffusion is neglected. A formal analysis can be thus simplified when electric field applied to a photorefractive material has relatively high intensity, or when we consider the propagation of relatively broad beam of light. In both cases we deal with drift dominated carrier transport.

Based on the above simplifying assumptions, we can bring the time evolution of space-charge field to a problem of one differential equation. Combining the continuity Eq. (5e) with the Gauss's law (5f) we obtain:

\frac{\partial}{\partial x} (J_{n} + J_{p}) = - ε ε_{0} \frac{\partial^{2} E}{\partial t \partial x} .

Then, using the described simplifications, Eq. (8) can be reduced to this form:

\frac{\partial E}{\partial t} + \frac{1}{τ_{D}} (1 + \frac{I}{I_{B} + I_{d}}) E = \frac{E_{0}}{τ_{D}} .

In the above equation, the quantity describing dielectric relaxation time

τ_{D}

equals:

τ_{D} = \frac{ε ε_{0}}{q (μ_{n} a_{n} + μ_{p} a_{p}) I_{B} + q (μ_{n} b_{n} + μ_{p} b_{p})},

where

\begin{array}{l} a_{n} = \frac{S_{n} (N_{D} - N_{A})}{γ_{n} N_{A}}, a_{p} = \frac{S_{p} N_{A}}{γ_{p} (N_{D} - N_{A})}, \\ b_{n} = \frac{β_{n} (N_{D} - N_{A})}{γ_{n} N_{A}}, b_{p} = \frac{β_{p} N_{A}}{γ_{p} (N_{D} - N_{A})} . \end{array}

Besides, the quantity

I_{d}

is expressed by the relation:

I_{d} = \frac{μ_{n} b_{n} + μ_{p} b_{p}}{μ_{n} a_{n} + μ_{p} a_{p}} .

The solution of Eq. (9) describing the relation between electric field and light intensity can be presented in this form:

E (t) = E_{0} \exp [- \frac{1}{τ_{D}} \int_{0}^{t} (1 + \frac{I}{I_{B} + I_{D}}) d t^{'}] {1 + \frac{1}{τ_{D}} \int_{0}^{t} d t^{'} \exp [\frac{1}{τ_{D}} \int_{0}^{t^{'}} (1 + \frac{I}{I_{B} + I_{D}}) d t^{″}]} .

The expression (12) takes into account time nonlocality, characteristic of photorefractive response. However, if we assume that the light intensity distribution of an optical beam propagating in a medium is independent of time, Eq. (12) can be simplified to this form:

E (t) = E_{0} \exp [- \frac{t}{τ_{D}} (1 + \frac{I}{I_{B} + I_{D}})] + E_{0} {(1 + \frac{I}{I_{B} + I_{D}})}^{- 1} {1 - \exp [- \frac{t}{τ_{D}} (1 + \frac{I}{I_{B} + I_{D}})]},

which describes the relation between electric field and light intensity in time-local approximation. All above approximations can be avoided by using an appropriate numerical approach.

It is worth to add a comment here regarding similar works that have been made in the context of photorefractive indium phosphide [36,37]. For this semiconductor, authors refer in their papers to the theory based on similar approximations as presented in Eqs. (9)-(13), whereas they adopt the numerical method based on a technique different from the one discussed in this article.

3. Numerical method

The simplest numerical approach to nonlinear light propagation in PR media is based on the algorithm that has a two-block structure [15]. One block is a computing scheme for solving a system of differential equations, respectively (1a) – (1d) or (5a) – (5f). This part of the algorithm enables to determine the photorefractive response of a medium to interacting light with a preset intensity distribution. Formally determined in this computing procedure are the distributions of concentrations of free carriers, ionized traps and electric field. Information on the photorefractive response is passed to the other block, which supplements it by determining the distribution of optically induced changes of the photorefractive index: $Δ n (E) = - n_{b}^{3} r_{e f f} E / 2$ , where $r_{e f f}$ is the effective linear electro-optical coefficient. The principal part of this computing procedure is the solution of wave Eq. (6), in which a proper distribution of changes in the refractive index has been introduced.

It is worth emphasizing that such construction of the computing scheme executes separated integration in the domains of time and space. Two blocks of the computing procedure may co-operate in a variety of ways, yielding various results from the physical viewpoint. Let us take a closer look at the procedure presented in [15] and its possible modifications.

Following Eq. (6) we will consider the configuration (1 + 1) D for an optical beam polarized in the direction of axis x and propagating in the axis z direction. Let us assume that Q(x,z,t) represents all physical quantities describing the photorefractive response of a medium. We introduce a beam with intensity distribution $I (x, z = 0, t)$ at the medium input. We expect solutions which present a light intensity distribution and distributions of variables describing the photorefractive response for any time $τ$ in the whole plane of the computing field, i.e. $I (x, z, t = τ)$ and $Q (x, z, t = τ)$ . We also assume that spatial discretization is characterized by J computing nodes in the axis x direction (with index j) spaced at dx intervals from each other and K computing nodes in the axis z direction (with index k) spaced at dz intervals. Discretization in the time domain includes S nodes (with index s) spaced at dt time increments.

The computing procedure presented in [15] is so constructed that to obtain a light intensity distribution in a computing line (k + 1) for a set time $t = τ$ , first on the basis of intensity distribution in the previous line (k-th) we have to obtain the photorefractive response $Q (x, z = k d z, τ = S d t)$ by executing S time loops within the first block. Then, taking $Q (x, z = k d z, τ = S d t)$ and $I (x, z = k d z, τ)$ the wave Eq. (6) is solved, which allows to determine light intensity distribution in the line (k + 1) for a preset time $τ$ , that is $I [x, z = (k + 1) d z, τ]$ . The procedure is repeated for subsequent steps towards the propagation until the last Kth line of the computing field is reached. Importantly, in each line of the computing procedure, the determination of quantities characterizing the photorefractive response starts from preset initial values $Q (x, z, t = 0)$ . In the simplest case, quantities $Q (x, z, t = 0)$ may have uniform distribution in the whole computing field. The block diagram of the algorithm is shown in Fig. 1. Further in these considerations we will refer to this numerical approach as algorithm or procedure I. It is worth noting that a similar computational procedure has been used in the simulations describing the propagation of a short pulse in a photorefractive insulator [38]. However, in this work, authors were able (after applying appropriate simplifying assumptions) to reduce the problem (1a) –(1d) to a single differential equation. In more complex issues (as in the case of semiconductors), such an approach may not be feasible.

Fig. 1 Computing algorithm: time local version; numerical procedure I.

Download Full Size | PDF

Despite advantages of procedure I, such as overcoming the restrictions of typically used analytical approximations, procedure I has a disadvantage. The process aimed at determining the distribution of light intensity in k-th computing line (for z = z₀) through the wave equation depends on light distribution in previous lines, but it does not depend on the dynamics of processes occurring earlier (for z < z₀). The algorithm structure does not take account of time nonlocal phenomena. The results obtained by the approach presented in [15] correspond to those obtained using the model based on Eq. (13) rather than Eq. (12). The foregoing details of algorithm I describe the method of numerical integration in the domains of time and space as illustrated in Fig. 1. However, we can build an algorithm with the reverse order of integration. Having the initial conditions and light intensity distribution of the beam entering a medium, we can update the photorefractive response by making just one time step dt, obtaining $Q (x, z, t + d t)$ . Next, based on the response and using the beam propagation method, we can get an updated to the time t + dt distribution of light intensity in the whole computing plane $I (x, z, t + d t)$ .

To get a light intensity distribution for a preset time τ, the described procedure has to be repeated S times, updating each time the last obtained photorefractive response and light intensity distribution with the next time step. If the time increment dt is relatively small, the established computing procedure will have properties of a time nonlocal algorithm. The described algorithm is schematically shown in Fig. 2(a). For clarity, further in this article this numerical approach will be referred to as procedure II.

Fig. 2 Block diagram of a computing algorithm: time nonlocal version; a) numerical procedure II, b) numerical procedure III.

Download Full Size | PDF

It is essential from the physical and numerical viewpoints to determine the magnitude of time step dt. The computing procedure for determining the photorefractive response imposes certain restrictions in this area. For the procedure to be stable and yield correct results from the physical viewpoint, the value dt should not be greater than the lowest of values of mean lifetimes of electrons τ_n and holes τ_p. For a model with donors these values can be estimated from this relation:

τ_{n} = \frac{1}{γ_{n} N_{A}}, τ_{p} = \frac{1}{γ_{p} (N_{D} - N_{A})} .

For typical semiconductors, the estimates are one-digit nanoseconds. Such short times affect the relatively long time of computing of the numerical procedure as each update of the medium response requires the wave equation to be solved by the beam propagation method. For instance, if we simulate a time period τ = 1μs with

d t = 10^{- 10} s

steps, in each time loop of algorithm II, apart from updating the medium response the algorithm should make 10⁴ activations of the beam propagation method. If time domain integration did not require such a short step dt, the overall procedure could be much faster. The simplest method to obtain such effect is to build an algorithm in which the photorefractive response will be updated with a step dt required for stable operation of the program, while the wave equation will be solved after each step dτ, being its own multiplication

(d τ = H d t)

. Note that the step dτ should be short enough to properly take into account time-nonlocal effects. An algorithm of this type will be called herein the numerical procedure III. Its details are shown in Fig. 2(b), showing that the external loop of the procedure is executed R times. To reach the preset simulation time τ, conforming to procedures I and II, the relation

S = H \cdot R

must be satisfied. The next chapter describes properties of all presented numerical procedures and compares the results produced by them.

4. Results of simulations

4.1 One-carrier model

To test the proposed numerical procedures for their applicability in light self-trapping simulations with time nonlocality effects, let us distinguish, following [17], two conditions that reduce Eq. (9) to a solution of local character. One condition refers to an approximation of a small modulation, with leads to the linearization of Eq. (9). This condition is satisfied when the inequality $I ≪ I_{B} + I_{D}$ is satisfied, meaning weak nonlinear coupling between electric field and light intensity. The other condition refers to a situation in which the light intensity distribution does not show strong dependence on time. Both conditions are untypical for light self-trapping experiments. To go beyond the approximation defined by the former condition, in simulations (the results are presented below) the beams used had a maximum intensity comparable to background intensity. To assure relatively high dynamics of light intensity distribution changes, light propagation took place along a few diffraction lengths.

Let us first present the computation results concerning the simulation illustrating the process of spatial soliton formation in a material described by Eqs. (1a)-(1d). The parameters used in computations are characteristic of the ferroelectric SBN crystal. Their values are gathered in Table 1. Numerical tests were based on procedure II and the results compared to those obtained from the local approach (procedure I).

Table 1. Material parameters used in computations based on model (1a)-(1d).

View Table | View all tables in this article

To generate initial conditions for t = 0 s, the whole computing field was uniformly filled with the following values of variables included in materials Eqs. (1a)-(1d):

n_{0} = \frac{(S_{n} I_{B} + β_{n}) (N_{D} - N_{A})}{γ_{n} N_{A}}, N_{D}^{+}_{0} = N_{A} + n_{0}, E_{0} = \frac{φ_{0}}{L} .

where

φ_{0}

is the value of electric potential applied to one electrode of a sample (the other electrode maintained zero potential), L – sample width.

To establish the initial light intensity distribution $I (x, z, t = 0)$ , the wave Eq. (6) was solved by the BPM method, taking into account the initial uniform change of the refractive index: $Δ n_{0} = - n_{b}^{3} r_{e f f} E_{0} / 2$ . At the input of the computing plane a stationary soliton profile was introduced, numerically determined using the approach similar to that presented in works [31,32]. To ensure propagation in conditions typical of time nonlocalities, a normalized value of maximum light intensity $I_{M} / (I_{B} + I_{D})$ was assumed to be equal to unity. The soliton profile (for electric field E₀ = 700 V/cm) has a half-width Fwhm = 14.5 μm, which yields a diffraction length L_D = 2.18 mm. To ensure dynamic changes of light intensity during propagation, the length of computing field was established Z = 14 mm, giving the ratio Z/L_D = 6.4. The other parameters of numerical simulation are as follows: L = 1 mm, dx = 2.5 μm, dz = 2.0 μm, dt = 0.3 μs. The results for six selected transition states (τ = 0, 0.3, 0.9, 1.65, 2.55, 9 ms) are displayed in Fig. 3. For clarity and better illustration, the soliton formation process is recorded as movie (Media 1). A slight but visible in the steady state Fig. 3(f) deflection of soliton trajectory to the left is caused by a unsymmetrical diffusion field. It can be observed because the simulation was carried out for the parameters for which the effect of diffusion cannot be neglected.

Fig. 3 Time evolution of a soliton beam in a material described by the model (1a)-(1d), determined by procedure II. The results present transition states for a) 0 ms, b) 0.3 ms, c) 0.9 ms, d) 1.65 ms, e) 2.55 ms, f) 9 ms. For more effective illustration, the process has been recorded (Media 1).

Download Full Size | PDF

The operation of both algorithms was assessed in quantitative terms by an analysis of time dependence of the diffraction coefficient D(t), defined as the ratio of the light intensity half-width at the computing field output to the intensity distribution half-width at the input. Figure 4 (a) shows the time dependence of the diffraction coefficient D(t) for algorithms I and II. The changes of the diffraction coefficient in both cases is rather exponential, although in the nonlocal approach the relation is 'extended' in time. For the local approach, the stationary state of solitons is formed in a shorter time than in the nonlocal approach. The results seem to conform to experimental observations described in work [17]. As noted before, the effect of time nonlocality occurs when the changes of light intensity distribution are relatively dynamic, characteristic of phenomena related to optical beams self-trapping. In case of low dynamics, the nonlocality should be much less prominent during the beam evolution. To test the operation of algorithm II in computing conditions reflecting such case, simulation for shorter lengths of propagation were performed. The propagation length was shortened to Z = 3.5 mm, thus Z/L_D = 1.6. Figure 4(b) shows the evolution of diffraction coefficients computed by procedure II for both cases, and compares them to the procedure I results. The plot differences resulting from the local and nonlocal approach are slight and smaller than for propagation with a length Z = 14 mm.

Fig. 4 a) Comparison of diffraction coefficient changes during self-trapping of light beam, determined by algorithms I and II. b) Comparison of time evolution of the diffraction coefficient for low dynamics of self-trapping. Computed results produced by algorithms I and II for propagation lengths Z = 3.5 mm (Z/L_D = 1.6).

Download Full Size | PDF

In the computations with results presented above, the time step of numerical integration was determined using the stability criterion of the method for solving the system of Eqs. (1a-1d). The step dt = 0.3 μs and was shorter than the mean lifetime of electrons $τ_{n} ≃ 0.5 μs$ estimated from relation (14). The time step may be shorter than the used one, which leads to longer computing time. To check whether additional shortening of the step dt can affect the computation results, simulations were made with propagation length Z = 14 mm and the time step of 0.1 μs. Comparison of the computation results for both cases are shown in Fig. 5(a). It can be observed that a change of time step at this level may only lengthen the computing time, with no qualitative changes.

Fig. 5 a) Comparison of the diffraction coefficient obtained by algorithm II for various time steps dt = 1x10⁻⁷s and 3x10⁻⁷s. b) Comparison of the diffraction coefficient obtained by algorithm III for various values of parameter H to results obtained by numerical procedures I and II.

Download Full Size | PDF

It seems more interesting to consider a reverse situation, where the time step will be longer than that based on the stability criterion of the method solving the system of Eqs. (1a)–(1d). This approach may be implemented by procedure III. A relevant experiment is significant mainly because the procedure reduces computing time. To compare the performance of algorithm III to algorithms I and II, three series of computations were performed for various step lengths dτ. The first simulation had a loop updating the photorefractive response repeated H = 200 times with a time step dt = 0.3 μs. The solution of wave equation by the BPM method in this case updated light intensity distribution with the step dτ = Hxdt = 60 μs. The whole procedure was repeated R = 150 times, so the total simulation time was 9 ms.

In the second series of computations H equaled 1000, giving dτ = 300 μs. To obtain the same time of simulation, the procedure was repeated R = 30 times. In the last series the updating loop was repeated H = 2500 times, giving dτ = 750 μs. In this case the procedure had to be repeated R = 12 times. The obtained values of time dependent diffraction coefficient D(t) from the three series are shown in Fig. 5(b). As a reference, graphs of diffraction coefficients produced by algorithms I and II are drawn. The results are interesting because selecting a proper, sufficiently small number of loops H, we can obtain results convergent with the procedure that takes into account time nonlocalities. Establishing a properly large number of loops H, we get a reverse result, that is convergent with the results of procedure I. Further in this article the computing times of the three numerical algorithms will be compared.

4.2 Two-carrier model

To test the performance of the described methods for model (5a)-(5f), a series of simulations were performed, illustrating the propagation of light beam in the photorefractive gallium arsenide. One characteristic of this material is its small value of linear electro-optic coefficient. This means that to obtain non-diffractive light propagation in such material, experimental conditions should allow for propagation of relatively broad optical beams at the presence of relatively high intensities of electric field. The work [16] presents the results illustrating the self-trapping process of a Gaussian optical beam propagating in a semiconducting material described by equations similar to Eqs. (5a)-(5f). The only difference with the model used in this study is the type of traps introduced in the model. Deep traps in Eqs. (5a)-(5f) are donors compensating shallow acceptor traps (unlike the model used in [16] based on deep acceptors and shallow donors). In this way the model (5a) – (5f) is a natural expansion of the model (1a) – (1d), making this study more consistent. Note that both approaches allow to control the type of conductivity of the examined material by controlling the value of trap compensation coefficient: r = N_A/N_D. Consequently, both donor and acceptor models may lead to similar research situations. That is why both approaches are used in studies explaining nonlinear wave mixing in multiple quantum well structures created in the material system GaAs/AlGaAs [3,15,16,27,29,30]. Numerical experiments presented in [16] are characterized by strong spatial nonlocality occurring during propagation. The phenomenon manifests itself in effects of trajectory curvature of the propagating optical beam, which may be temporal or stationary. From the illustrations of numerical methods discussed herein the most interesting case is that of strong temporal self-bending of the optical beam. It depicts the high dynamics of changes in light intensity distribution. Such situation occurs in case of domination of holes, that within the model (5a) – (5f), taking into account material parameters collected in Table 2, may be attained for the trap compensation coefficient r = 0.5.

Table 2. Material parameters used in computations based on model (5a) – (5f).

View Table | View all tables in this article

In tests, whose results are presented above, a beam with Gaussian intensity distribution was introduced at the input of a computing field. The beam radius was w = 14.5 μm, equivalent to half-width Fwhm = 24 μm. To allow the beam to propagate in conditions in which a soliton of similar dimensions propagates, an electric field of 35 kV/cm should be applied to the material. In numerical experiments the whole computing area was illuminated with homogeneous background with a value selected so as to give the normalized maximum light intensity $I_{M} / (I_{B} + I_{D}) = 1$ . A relatively wide input beam requires a long computing field, in this case Z = 20 mm, giving the ratio Z/L_D = 9.2. To define the length of time step dt, mean carrier lifetimes, determined from material parameters, were calculated to have these values: $τ_{n} ≃ 0.213 ns$ and $τ_{p} ≃ 0.55 ns$ . Such values permit to obtain a stable operation of the method solving the system of Eqs. (5a)-(5f) with the time step dt = 0.2 ns. The other parameters of numerical simulation are as follows: L = 1 mm, dx = 2.5 μm, dz = 2.0 μm. As the simulations made for two-carrier model are demanding in terms of computing time, to reduce it we used algorithm III and set a low value of H = 20. Thus the light intensity distribution was updated with the time step dτ = Hxdt = 4 ns. The results for six transition states (τ = 0, 2, 4.8, 12.8, 25.2, 41 μs) are given in Fig. 6. Like for one-carrier model, the optical beam self-trapping process is demonstrated in a movie (Media 2).

Fig. 6 Time evolution of a Gaussian beam in a material described by the model (5a)- (5f), determined by algorithm III for H = 20. The results present transition states for a) 0 μs, b) 2 μs, c) 4.8 μs, d) 12.8 μs, e) 25.2 μs, f) 41 μs. For better illustration, the process is shown in a movie (Media 2).

Download Full Size | PDF

The case of temporal self-bending effect, can linked to experimental observations in a cubic photorefractive BTO crystal [39]. Although the model of the photorefractive effect in semiconductors such as GaAs, appears to be relatively more complex, the cause of demonstrated phenomena can be posited in a large photocarrier drift lengths, as suggested in [39].

Strong dynamics of light intensity distribution changes observed in transition states of the self-trapping process will entail the effect of strong time nonlocality. To compare the results to those obtained in the time locality approach, computations were made by numerical procedure I for selected transition states. Figure 7 displays profiles of light intensity distribution at the output of the computing field obtained by both methods. It is characteristic of both time evolutions that the distribution of light intensity goes through similar transition states, although the profile from algorithm III is delayed in comparison to profile from algorithm I. On this basis we can also state that the results presented in [16], even if obtained through the time local approach, should be qualitatively convergent with those taking into account time nonlocality. Nevertheless, an analysis based on time nonlocality seems to be essential from the fundamental and practical viewpoints. The use of such analysis is of major importance particularly when new devices with photorefractive materials are to be designed.

Fig. 7 Time evolution of light intensity distribution collected at the output of computing field for algorithm I and III (H = 20). The computation results present transition states for a) 0 μs, b) 2 μs, c) 4 μs, d) 8 μs, e) 16 μs, f) 20 μs, g) 24 μs, h) 32 μs, i) 48 μs.

Download Full Size | PDF

4.3 Computing efficiency

The described numerical algorithms are generally characterized by a relatively long computing time. Of many factors affecting simulation time, the complexity of a model describing a photorefractive medium is the major one. In this context materials described by model (5a) – (5f) will need much higher computing power than less complex media described by model (1a) – (1d). One method to reduce computing time of the former (as indicated in [15]) is to arrange parallel computations of hole and electron concentrations. In this study, tests were made to estimate the computing efficiency of the numerical procedures based on sequential algorithms. The tests mainly aimed at comparing the computing time of procedures taking into account time nonlocal effects.

The reference computing field used in the tests was 1mm x 1mm, and the discretization was the same as in simulations, whose results are given in sections 4.1 and 4.2: dx = 2.5 μm, dz = 2 μm, dt = 0.3 μs and dt = 0.2 ns, for models (1a) – (1d) and (5a) – (5f), respectively. The total simulation time τ = 1000dt. The tests, using procedure III, started for parameters H = 1 and R = 1000, which is equivalent to computations by procedure II. In subsequent tests the parameter H was gradually increased with such changes of parameter R that allowed to maintain the preset total simulation time τ. The last simulation for H = 1000 and R = 1 was equivalent to computations by procedure I. The computations were performed on an average class PC (AMD Athlon II X3 450 3.20 GHz) in the Windows environment. The results are collected in Table 3.

Table 3. Computing times determined by the numerical procedure III.

View Table | View all tables in this article

As expected, procedure II is the most and procedure I is the least demanding in terms of computing efficiency. Comparison of the presented computing times leads to a conclusion that the use of algorithm III is particularly attractive when the light intensity distribution is updated every H = 20 ÷ 100 time loops that update the photorefractive response. Such selection of algorithm parameters permits to reach computing time at a level attained by procedure I, with time nonlocalities taken into account relatively well.

5. Conclusions

The article presents a method of numerical analysis of nonlinear light propagation in photorefractive materials, taking account of time nonlocal phenomena. Two numerical approaches are discussed (algorithm II and III), modifications of the method presented in [15] (algorithm I), based on the time locality approach. Tests were made for a one-carrier model, characteristic of common photorefractive materials, such as ferroelectric strontium barium niobate and for a two-carrier model describing photorefraction in semiconductors with hot electron effect, e.g. gallium arsenide.

The tests have shown that by reversing the sequence of integration embedded in the structure of algorithm I, we obtain results taking into account time nonlocality that are congruent with results obtained experimentally. Besides, it has been proved that by using properly chosen parameters, algorithm III produces satisfactory results (taking into account time nonlocal phenomena), within computing time comparable to that reached by algorithm I. Tests of the process of self-trapping in semiconductors have also shown that results from local and nonlocal approaches are qualitatively congruent, and differences consist mainly in the time to reach the stationary state.

The presented methods offer an interesting tool for analysis of nonlinear light propagation in photorefractive media. The area of application of these methods comprises, first of all, analysis of phenomena occurring in materials with a relatively complex band transport model. It seems worthwhile working on the proposed numerical solutions towards the reduction of computing time, increasing the model dimensionality and incorporation of the associated anisotropy of a medium.

Appendix

Here we briefly describe the way of discretization used in the numerical procedure discussed in this paper. Formal structure of equations describing PR response is similar to the structure of basic semiconductor equations. Standard approach to solve the time-dependent semiconductor equations in the frame of finite difference method can be found in [40,41]. A modification of this approach to the needs of PR media described by Eqs. (1a)-(1d) was presented by Singh et al. [8]. The extension of this approach to the case of materials characterized by bipolar transport and hot electron effect is presented below. Firstly, combining the continuity Eq. (5e) with Eqs. (5a) and (5b) we can obtain:

\frac{\partial N_{D}^{+}}{\partial t} = (\frac{\partial n}{\partial t} - \frac{1}{q} \frac{\partial J_{n}}{\partial x}) - (\frac{\partial p}{\partial t} + \frac{1}{q} \frac{\partial J_{p}}{\partial x}) .

From Eq. (16) we can determine discrete value of ionized donor concentration in the next time-step:

N_{D}^{+}_{j}^{s + 1} = \frac{N_{D}^{+}_{j}^{s} + N_{D} [​ S_{n} (I_{j}^{s} + I_{B}) + β_{n} + γ_{p} p_{j}^{s}] d t ​ - ​ N_{D}^{+}_{j}^{s} [S_{p} (I_{j}^{s} + I_{B}) + β_{p} + γ_{n} n_{j}^{s}] d t}{1 + [S_{n} (I_{j}^{s} + I_{B}) + β_{n} + γ_{p} p_{j}^{s}] d t},

where s and j are, respectively, time and space (in the axis x direction) computing nodes. Values of free carrier concentrations, can be obtained from discrete form of continuity equations of electrons and holes. Using a similar approach as in [8,41] we can rewrite Eq. (5b) in form:

A^{+} p_{j + 1}^{s + 1} + B^{+} p_{j}^{s + 1} + C^{+} p_{j - 1}^{s + 1} = A p_{j + 1}^{s} + B p_{j}^{s} + C p_{j - 1}^{s} + D,

where:

\begin{array}{l} A^{+} = ξ (\frac{μ_{p}}{2 d x} E_{j + 1}^{s} - \frac{D_{p}}{d x^{2}}), B^{+} = \frac{1}{d t} + ξ \frac{2 D_{p}}{d x^{2}}, C^{+} = ξ (- \frac{μ_{p}}{2 d x} E_{j - 1}^{s} - \frac{D_{p}}{d x^{2}}), \\ A = (ξ - 1) (\frac{μ_{p}}{2 d x} E_{j + 1}^{s} - \frac{D_{p}}{d x^{2}}), B = \frac{1}{d t} + (ξ - 1) \frac{2 D_{p}}{d x^{2}} - γ_{p} (N_{D} - N_{D}^{+}_{j}^{s}), \\ C = (ξ - 1) (- \frac{μ_{p}}{2 d x} E_{j - 1}^{s} - \frac{D_{p}}{d x^{2}}), D = N_{D}^{+} [S_{p} (I_{j}^{s} + I_{B}) + β_{p}], \end{array}

In case of electrons, discrete form of continuity Eq. (5a) takes the form:

F^{+} n_{j + 1}^{s + 1} + G^{+} n_{j}^{s + 1} + H^{+} n_{j - 1}^{s + 1} = F n_{j + 1}^{s} + G n_{j}^{s} + H n_{j - 1}^{s} + M,

where:

\begin{array}{l} F^{+} = ξ (- \frac{μ_{n}_{j + 1}^{s}}{2 d x} E_{j + 1}^{s} + \frac{D_{n}_{(j + 1 / 2)}}{d x^{2}}), G^{+} = \frac{1}{d t} - ξ \frac{D_{n}_{(j + 1 / 2)} + D_{n}_{(j - 1 / 2)}}{d x^{2}}, \\ H^{+} = ξ (\frac{μ_{n}_{j - 1}^{s}}{2 d x} E_{j - 1}^{s} + \frac{D_{n}_{(j - 1 / 2)}}{d x^{2}}), F = (\frac{μ_{n}_{j + 1}^{s}}{2 d x} E_{j + 1}^{s} - \frac{D_{n}_{(j + 1 / 2)}}{d x^{2}}), \\ G = \frac{1}{d t} - (1 - ξ) \frac{D_{n}_{(j + 1 / 2)} + D_{n}_{(j - 1 / 2)}}{d x^{2}} - γ_{n} N_{D}^{+}_{j}^{s}, H = (ξ - 1) (\frac{μ_{n}_{j - 1}^{s}}{2 d x} E_{j - 1}^{s} + \frac{D_{n}_{(j - 1 / 2)}}{d x^{2}}), \\ M = [S_{n} (I_{j}^{s} + I_{B}) + β_{n}] (N_{D} - N_{D}^{+}_{j}^{s}), D_{n}_{(j + 1 / 2)} \approx \frac{D_{n}_{j + 1} + D_{n}_{j}}{2}, D_{n}_{(j - 1 / 2)} \approx \frac{D_{n}_{j} + D_{n}_{j - 1}}{2} . \end{array}

In the above equations

D_{p}

and

D_{n}

are the diffusion coefficients,

ξ

is an auxiliary parameter (for

ξ = 0

the scheme is explicit, while for

0 < ξ < 1

the scheme is implicit). To determine value of electron mobility

μ_{n}

and electron diffusion coefficient

D_{n}

we used Eqs. (2)-(4) with approximate form of the distribution function [15,30]:

f (E) = {[1 + {(\frac{E}{E_{C}})}^{5.6}]}^{- \frac{1}{4}},

where

E_{C} = 2.8 k V / c m .

Furthermore, we have assumed that the electron diffusion coefficient is slowly varying in space. Using known values of electron, hole and ionized donor concentrations, we can determine electric potential distribution (and in the next step electric field distribution) from the discrete form of Eq. (5f):

φ_{j + 1}^{s + 1} - 2 φ_{j}^{s + 1} + φ_{j - 1}^{s + 1} = - \frac{q d x^{2}}{ε ε_{0}} (N_{D}_{j}^{s + 1} + p_{j}^{s + 1} - n_{j}^{s + 1} - N_{A}),

The difference Eqs. (18)-(20) are solved using tridiagonal matrix solver. All Eqs. (17)-(20) are the results of the block I. In block II, these equations are supplemented by the discrete form of the wave Eq. (6), which can be solved by Crank-Nicolson method.

Acknowledgments

This work has been partially supported by the National Center for Science under the project awarded by decision number DEC-2011/01/B/ST7/06234

References and links

1. R. W. Boyd, Nonlinear Optics, Electronics & Electrical (Elsevier Science, 2003).

2. P. Yeh, Introduction to Photorefractive Nonlinear Optics (John Wiley & Sons, 1993).

3. D. D. Nolte, Photorectractive Effects and Materials (Kluwer, 1995).

4. L. Solymar, D. J. Webb, and A. Grunnet-Jepsen, The Physics and Applications of Photorefractive Materials (Clarendon, 1996).

5. P. Günter and J. P. Huignard, Photorefractive Materials and Their Applications 2: Materials (Springer-Verlag, 2006).

6. P. Günter and J. P. Huignard, Photorefractive Materials and Their Applications 3: Applications (Springer-Verlag, 2007).

7. G. A. Brost, “Numerical analysis of photorefractive grating formation dynamics at large modulation in BSO,” Opt. Commun. 96(1-3), 113–116 (1993). [CrossRef]

8. N. Singh, S. P. Nadar, and P. P. Banerjee, “Time-dependent nonlinear photorefractive response to sinusoidal intensity gratings,” Opt. Commun. 136(5-6), 487–495 (1997). [CrossRef]

9. J. G. Murillo, “Photorefractive grating dynamics in Bi12SiO20 using optical pulses,” Opt. Commun. 159(4-6), 293–300 (1999). [CrossRef]

10. T. Gatlin and N. Singh, “Nonlinear frequency response of a moving grating with an applied field in bismuth silicon oxide,” Opt. Lett. 24(22), 1593–1595 (1999). [CrossRef] [PubMed]

11. N. Fressengeas, D. Wolfersberger, J. Maufoy, and G. Kugel, “Build up mechanisms of (1+1)-dimensional photorefractive bright spatial quasi-steady-state and screening solitons,” Opt. Commun. 145(1-6), 393–400 (1998). [CrossRef]

12. J. Maufoy, N. Fressengeas, D. Wolfersberger, and G. Kugel, “Simulation of the temporal behavior of soliton propagation in photorefractive media,” Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 59(55 Pt B), 6116–6121 (1999). [CrossRef] [PubMed]

13. F. Devaux, V. Coda, M. Chauvet, and R. Passier, “New time-dependent photorefractive three-dimensional model: application to self-trapped beam with large bending,” J. Opt. Soc. Am. B 25(6), 1081–1086 (2008). [CrossRef]

14. F. Devaux and M. Chauvet, “Three-dimensional numerical model of the dynamics of photorefractive beam self-focusing in InP:Fe,” Phys. Rev. A 79(3), 033823 (2009). [CrossRef]

15. A. Ziółkowski, “A numerical approach to nonlinear propagation of light in photorefractive media,” Comput. Phys. Commun. 185(2), 504–511 (2014). [CrossRef]

16. A. Ziółkowski, “Self-bending of light in photorefractive semiconductors with hot-electron effect,” Opt. Express 22(4), 4599–4605 (2014). [CrossRef] [PubMed]

17. C. Dari-Salisburgo, E. DelRe, and E. Palange, “Molding and stretched evolution of optical solitons in cumulative nonlinearities,” Phys. Rev. Lett. 91(26), 263903 (2003). [CrossRef] [PubMed]

18. E. DelRe and E. Palange, “Optical nonlinearity and existence conditions for quasi-steady-state photorefractive solitons,” J. Opt. Soc. Am. B 23(11), 2323–2327 (2006). [CrossRef]

19. G. Picoli, P. Gravey, and C. Ozkul, “Model for resonant intensity dependence of photorefractive two-wave mixing in InP:Fe,” Opt. Lett. 14(24), 1362–1364 (1989). [CrossRef] [PubMed]

20. M. Chauvet, S. A. Hawkins, G. J. Salamo, M. Segev, D. F. Bliss, and G. Bryant, “Self-trapping of planar optical beams by use of the photorefractive effect in InP:Fe,” Opt. Lett. 21(17), 1333–1335 (1996). [CrossRef] [PubMed]

21. M. Chauvet, S. A. Hawkins, G. J. Salamo, M. Segev, D. F. Bliss, and G. Bryant, “Self-trapping of two-dimensional optical beams and light-induced waveguiding in photorefractive InP at telecommunication wavelengths,” Appl. Phys. Lett. 70(19), 2499–2501 (1997). [CrossRef]

22. H. Leblond and N. Fressengeas, “Theory of photorefractive resonance for localized beams in two-carrier photorefractive systems,” Phys. Rev. A 80(3), 033837 (2009). [CrossRef]

23. C. Dan, D. Wolfersberger, and N. Fressengeas, “Experimental control of steady state photorefractive self-focusing in InP:Fe at infrared wavelengths,” Appl. Phys. B 104(4), 887–895 (2011). [CrossRef]

24. N. Fressengeas, C. Dan, and D. Wolfersberger, “Microsecond infrared beam bending in photorefractive iron doped indium phosphide,” Opt. Laser Technol. 48, 96–101 (2013). [CrossRef]

25. S. M. Sze and K. N. Kwok, “Physics and Properties of Semiconductors - A Review,” in Physics of Semiconductor Devices (John Wiley & Sons, 2002).

26. K. Seeger, Semiconductor Physics - An Introduction (Springer-Verlag, 2004).

27. Q. N. Wang, R. M. Brubaker, and D. D. Nolte, “Photorefractive phase shift induced by hot-electron transport: multiple-quantum-well structures,” J. Opt. Soc. Am. B 11(9), 1773–1779 (1994). [CrossRef]

28. D. D. Nolte, S. Balasubramanian, and M. R. Melloch, “Nonlinear charge transport in photorefractive semiconductor quantum wells,” Opt. Mater. (Amst) 18(1), 199–203 (2001). [CrossRef]

29. A. Ziółkowski, “Temporal analysis of solitons in photorefractive semiconductors,” J. Opt. 14(3), 035202 (2012). [CrossRef]

30. M. Wichtowski and A. Ziółkowski, “Interband photorefractive effect in semiconductors with hot-electron transport at arbitrary modulation depth,” Opt. Commun. 300, 257–264 (2013). [CrossRef]

31. M. Segev, G. C. Valley, B. Crosignani, P. DiPorto, and A. Yariv, “Steady-state spatial screening solitons in photorefractive materials with external applied field,” Phys. Rev. Lett. 73(24), 3211–3214 (1994). [CrossRef] [PubMed]

32. D. N. Christodoulides and M. I. Carvalho, “Bright, dark, and gray spatial soliton states in photorefractive media,” J. Opt. Soc. Am. B 12(9), 1628–1633 (1995). [CrossRef]

33. A. Ziółkowski and E. Weinert-Raczka, “Spatial solitons in biased photorefractive media with quadratic electro-optic effect,” Opto-electronics Rev. 13, 135–140 (2005).

34. A. Ziółkowski and E. Weinert-Raczka, “Dark screening solitons in multiple quantum well planar waveguide,” J. Opt. A, Pure Appl. Opt. 9(7), 688–698 (2007). [CrossRef]

35. E. DelRe, B. Crosignani, and P. Di Porto, “Photorefractive solitons and their underlying nonlocal physics,” Prog. Opt. 53, 153–200 (2009). [CrossRef]

36. N. Fressengeas, N. Khelfaoui, C. Dan, D. Wolfersberger, G. Montemezzani, H. Leblond, and M. Chauvet, “Roles of resonance and dark irradiance for infrared photorefractive self-focusing and solitons in bipolar InP:Fe,” Phys. Rev. A 75(6), 063834 (2007). [CrossRef]

37. D. Wolfersberger, N. Khelfaoui, C. Dan, N. Fressengeas, and H. Leblond, “Fast photorefractive self-focusing in InP:Fe semiconductor at infrared wavelengths,” Appl. Phys. Lett. 92(2), 021106 (2008). [CrossRef]

38. D. Wolfersberger, F. Lhomme, N. Fressengeas, and G. Kugel, “Simulation of the temporal behavior of one single laser pulse in a photorefractive medium,” Opt. Commun. 222(1-6), 383–391 (2003). [CrossRef]

39. P. A. M. Aguilar, J. J. S. Mondragon, S. Stepanov, and V. Vysloukh, “Transient self-bending of laser beams in photorefractive crystals with drift nonlinearity,” Phys. Rev. A 54(4), R2563–R2566 (1996).

40. S. Selberherr, Analysis and Simulation of Semiconductor Devices (Springer-Verlag, 1984).

41. M. Reiser, “Large-scale numerical simulation in semiconductor device modelling,” Comput. Methods Appl. Mech. Eng. 1(1), 17–38 (1972). [CrossRef]

Parameter	Value
Cross section for electron photogeneration divided by photon energy	S_n = 4x10⁻⁶ m²/J
Cross section for recombination	γ_n = 1x10⁻¹⁶ m³/s
Electron mobility	μ_n = 2.7x10⁻⁶ m²/Vs
Donor concentration	N_D = 4.7x10²⁴ m⁻³
Acceptor concentration	N_A = 1.7x10²² m⁻³
Dielectric constant	ε = 850
Wavelength	λ = 514 nm
Refractive index	n_b = 2.35
Electro-optic coefficient	r₃₃ = 224 pm/V

Parameter	Value
Cross section for photogeneration from traps	S_n = 1⋅10⁻¹⁷ cm²	S_p = 1⋅10⁻¹⁶ cm²
Cross section for recombination to traps	γ_n = 4.7⋅10⁻⁷ cm³/s	γ_p = 1.8⋅10⁻⁷ cm³/s
Linear carrier mobility	μ_nl = 5000 cm²/Vs	μ_p = 300 cm²/Vs
Electron mobility in upper valley	μ_nu = 300 cm²/Vs
Carrier thermal emission rate	β_n ≈0.5 s⁻¹	β_p ≈1 s⁻¹
Concentration of donor traps	N_D = 2x10¹⁶ cm⁻³
Concentration of acceptor traps	N_A = 1x10¹⁶ cm⁻³
Dielectric constant	ε = 12.58
Wavelength	λ = 1.06 μm
Refractive index	n_b = 3.48
Electro-optic coefficient	r₄₁ = 1.43 pm/V

	H = 1 R = 1000	H = 5 R = 200	H = 20 R = 50	H = 100 R = 10	H = 200 R = 5	H = 1000 R = 1
Model (1a - 1d) τ = 1000dt = 300 μs	429 s	158 s	108 s	95 s	93 s	91 s
Model (5a - 5f) τ = 1000dt = 0.2 μs	856 s	574 s	519 s	503 s	502 s	499 s

Parameter	Value
Cross section for electron photogeneration divided by photon energy	S_n = 4x10⁻⁶ m²/J
Cross section for recombination	γ_n = 1x10⁻¹⁶ m³/s
Electron mobility	μ_n = 2.7x10⁻⁶ m²/Vs
Donor concentration	N_D = 4.7x10²⁴ m⁻³
Acceptor concentration	N_A = 1.7x10²² m⁻³
Dielectric constant	ε = 850
Wavelength	λ = 514 nm
Refractive index	n_b = 2.35
Electro-optic coefficient	r₃₃ = 224 pm/V

Parameter	Value
Cross section for photogeneration from traps	S_n = 1⋅10⁻¹⁷ cm²	S_p = 1⋅10⁻¹⁶ cm²
Cross section for recombination to traps	γ_n = 4.7⋅10⁻⁷ cm³/s	γ_p = 1.8⋅10⁻⁷ cm³/s
Linear carrier mobility	μ_nl = 5000 cm²/Vs	μ_p = 300 cm²/Vs
Electron mobility in upper valley	μ_nu = 300 cm²/Vs
Carrier thermal emission rate	β_n ≈0.5 s⁻¹	β_p ≈1 s⁻¹
Concentration of donor traps	N_D = 2x10¹⁶ cm⁻³
Concentration of acceptor traps	N_A = 1x10¹⁶ cm⁻³
Dielectric constant	ε = 12.58
Wavelength	λ = 1.06 μm
Refractive index	n_b = 3.48
Electro-optic coefficient	r₄₁ = 1.43 pm/V

Numerical method for an analysis of nonlinear light propagation in photorefractive media - time nonlocal approach

Abstract

1. Introduction

2. Theoretical background

2.1 Model of the phenomenon

2.2 Analytical approximations

3. Numerical method

4. Results of simulations

4.1 One-carrier model

4.2 Two-carrier model

4.3 Computing efficiency

5. Conclusions

Appendix

Acknowledgments

References and links

Supplementary Material (2)

Cited By

Figures (7)

Tables (3)

Equations (33)

Optics Express