Time-dependent nonlinear photorefractive response to sinusoidal intensity gratings

.--_ bf!!!i

cQ&% ELSEVIER

1 April 1997

OPTICS COMMUNICATIONS Optics Communications

136 (1997) 487-495

Full length article

Time-dependent nonlinear photorefractive response to sinusoidal intensity gratings Nagendra Singh, S.P. Nadar, Partha P. Banerjee Depurfmmt r$Elecfricul

und Compurrr Engineering,

Uniurrsiry of Alubumu in Hunrsoille, Huntsville. AL 35899, USA

Received 28 May 1996; revised 18 September

1996; accepted 6 November 1996

Abstract The time-dependent, nonlinear Kukhtarev equations are numerically solved to determine the temporal and spatial dependence of the space-charge field for a sinusoidal intensity profile with an arbitrary modulation depth M. The intensity profile is spatially bounded simulating two interfering laser beams of finite spot size. The numerical solution allows us to critically examine the validity of analytical theories for linear and nonlinear photorefractive effects. The parametric dependence of the characteristic time constant r for writing a grating is studied including the effects of quiescent intensity I,,, modulation index m and the thermal ionization rate p. We find that the Kukhtarev equations inherently yield a sublinear dependence of the characteristic time constant T on I, even for a single trap (acceptor) species when I,, is above a critical value I,. Expression for I, is given in terms of material properties, and the sublinearity index is determined for BSO. The exact dependence of T on the modulation index (0 < M 5 1) is determined for BSO, and the nonsinusoidal waveform of the space-charge field for large m is fully characterized both analytically and numerically.

1. Introduction Photorefractive materials are generally analyzed using a set of material equations first put together by Kukhtarev and now known as Kukhtarev equations [l]. These equations describe the temporal and spatial evolution of the space charges and associated space-charge electric fields in a photorefractive crystal when it is illuminated by an inhomogeneous light intensity. The set includes equations for the generation of free charge carriers, their transport and the Gauss law for the space-charge electric field. In early theoretical studies [l] on PR materials, electrons were assumed to be the only free charge carriers. However, in later studies both electrons and holes have been included

M. The set of Kukhtarev equations and its further modifications to include additional charge carriers [2] are basically a set of coupled partial differential equations involving time and space. The equations are intrinsically nonlinear. In order to solve such equations a variety of approximations have been made [3]. A commonly used approxi-

mation is to linearize the equations by assuming that the relevant physical quantities can be separated into zero- and first-order quantities. The first-order quantities depend on the strength of the spatial modulation in the light intensity. For example, for the light intensity pattern I(x) = 1,[1 + m cos(kRx)] with m as the modulation index and k, as the grating wave number, ml, determines the first-order quantities, while the zero-order quantities depend on the average intensity I,. Under the linear approximation for m -+c I, the temporal evolution of the zero and first-order quantities are calculated. However, the temporal evolution of the first-order quantities are inseparably coupled with the zero-order quantities. Since the zero-order quantities also depend on time, the analytical solutions of the firstorder quantities critically depends on the time scale (r,,) at which the zero-order quantities evolve. Typically 7, is given by the recombination time 7, y (y,N,)- ‘, where ‘y, is the recombination constant and N, the acceptor concentration. If this time scale T, -=c pi, the time scale for the evolution of the first-order quantities determined by the dielectric relaxation time, it is possible to obtain an analyt-

0030-4018/97/$12.00 Copyright 0 1997 Elsevier Science B.V. All rights reserved. PII SOO30-40 18(96)007 14-6

488

N. Singh et al./ Optics Communications

136 (1997) 487-495

ical solution for the temporal evolution of the first-order quantities. This regime (7, GC 7,) has been studied by a number of authors in the field of photorefractive materials [l-5]. On the other hand, when r0 > T, for large intensities I,, time-dependent analytical solutions are not possible even for the linear case (m -=x I> and numerical solutions of the Kukhtarev equations are needed. Another motivation for using numerical solution of Kukhtarev equations is that the linearization process described above is valid only for a relatively weak inhomogeneity (m < I) in the light intensity pattern. When the intensity pattern has relatively strong modulation, for example, a pattern produced by two equally strong laser beams, it is not possible to obtain analytical solutions for the temporal evolution of the physical quantities. It is surprising to find that in the past there has been little effort on numerical solutions of the Kukhtarev equations, which are a set of coupled nonlinear partial differential equations. The only reference we can find on such an effort is that by Moharam et al. [4]. The purpose of this paper is to present a systematic numerical study of the original set of Kukhtarev equations [I]. The numerical method presented here can be easily extended to the modifications in the original set, extending the model to include electron-hole competition and different levels of traps for the charge carriers [2,5]. We present here results from the numerical solution and compare them against available analytical results. Time-dependent analytical results on the evolution of the space-charge electric field are generally available under the linear approximation for relatively small modulation index m -GC1 [l-3]. We find that the analytical results on the temporal evolution even for small m -C I are generally valid only when the average intensity I, < I,, where I, is a critical intensity depending on the material properties. For I, < I,, the time (7) taken to reach a steady state for the space-charge electric field is found to behave as ~a I;‘. When I, > I,, 7 decreases much slower (subinversely) with increasing

electric field for arbitrarily large m, including the effects of thermal excitation, are given and compared against available analytical theory. Section 3.3 deals with the time evolution of the electric field for arbitrary m and thermal excitation rate p. Main conclusions are given in Section 4.

1,.

p;+’

For relatively large m, for which the linearization breaks down, the waveform of the space-charge electric field becomes highly nonsinusoidal in contrast to the sinusoidal behavior expected from steady-state solutions for small m [3]. The temporal evolution to the steady state for relatively large values of m --) I shows that the different parts of the spatial waveform of the field evolves quite differently; the maxima in the waveform evolve the slowest. The numerical and analytical [3,6] results on the steady-state waveform of the space-charge electric fields are found to agree extremely well for all values of I, and m. The rest of the paper is organized as follows. The numerical method is described in Section 2. Section 3 deals with the numerical results. In Section 3.1, numerical results are compared against analytical results for the parameters for which latter results are available for the linear case m -=z 1. In Section 3.2 results on the steady-state

2. Numerical model Our numerical model is based on the following Kukhtarev equations [I] 3Pp: ~=(SI+p)(p,-p;)-(y,/e)p:p,, 3P, -=at

3Pd at

d*@

dE

dX2

ax

-=__=__

(1) a%

+l&p,E)+DS,

:(

set of

d-p;

-Pi).

(2) (3)

where TV= eNd, pi = eNd+, p; = eN‘,, pi = eN; = eN,; Nd, Nd+, N, are the donor, acceptor, photo-ionized donor, and electron concentrations, respectively; @ is the electrostatic potential and electrostatic (space-charge) field E = -G/dx; I is the light intensity; p is the thermal excitation rate, S is photoionization cross-section, -y, is the recombination rate, D and p are the effective electronic diffusion coefficient and effective mobility, respectively, and E is the effective electrostatic permittivity. In the above equations t is time and x is the distance along the PR crystal as shown in Fig. 1. The diffusion constant D and the mobility p are related by D = pk,T/e, where k, is the Boltzmann constant, T is the crystal temperature, and e is the magnitude of the electronic charge. In order to seek time dependent solutions of the above set of couple equations, we convert them into the following difference equations [7,8] N,,

=

(sI+p)eN,Ar-(r,/e)p,kAt+p~ 1 +At(sI+p)

-~~jp,“,:,’+bjP,k,+’

- f_Yjp,“,,‘,’= djy

(4)

(5)

Eqs. (4), (5), and (6) are the difference equations for Eqs. (1). (2), and (3). respectively. The superscript k represents quantities at a time step kAt and the subscript j represents the position jAx, where At is the time step used in advancing solutions starting with a set of initial conditions at time r = 0, and Ax is the grid spacing shown in Fig. 1. The coefficients appearing in Eq. (5) are defined by

(7) b,=

-2Dp/Ax’-

I/At,

(8)

N. Singh et al./Optics

Communications I36 (1997) 487-495

489

(9) d,=

p,“,+l(l -p)A

- 24

(1 -p)D/Ax*

+ I/At

( -3p:)

+(I-p)(D/Ax’-&Ej+,)

(‘0) In writing the difference equation for pe, we have differenced the term containing a*pJdx* using a mix of k and k+ 1 time levels. The mixing is determined by p; when p = 0 the scheme is explicit, while for p = 1 it is implicit. In our solutions here we have fixed p = 0.5, which gives good numerical stability [7,8]. The difference equation (4) yields pj” from known values of p,k and pj at the kth time level. The difference equations (5) and (6) are solved using a tridiagonal matrix solver [8]. The boundary and initial conditions in solving (1) and (3) over 0 I x I L are as follows. The initial conditions inside the crystal are obtained from (1) and (2); assuming I = 0 for t < 0, Eq. (1) gives p( pd - p:) = (y,./e)p:p, and pJ= pi + pi. The solutions of these two equations give the initial values of pi and p, for a fixed value of pi. The boundary values of p:, p; and pi at x = 0 and x = L remain the same as the initial values because the ends of the crystal are kept in the dark by properly shaping the light beams as explained shortly. For the electrostatic potential @, we assume that E(x = 0) - &D/axl,= o = 0 and @(x = L) = 0. The geometrical configuration of our calculation is well suited for applying an external electric field to the crystal simply by changing the boundary condition @(x = ,C) to a non-zero value if desired. However, in this paper we limit our discussion to the case for a zero applied electric field. The numerical results presented in this Letter are mainly for BSO with the following parameters [9]: Nd= lO25 mm’, N, =0.95 X IO** m-s, p = 10e5 m2/V s, y,= 1.65 X lo-l7 m3 SC’, S= 1.06 X 10e5 m*/J, and E= 9.76,. We present here results for a prescribed light intensity pattern I( X> = I,[ 1 + m cos(k, x)]G( x), where G(x) de-

Fig. I. Geometry of problem. Kukhtarev (x = 0) to j= M(x Ax = A/60, where

our calculation. We solve a one-dimensional equations are solved over a grid from j= I = x,,,). Typically we used a grid spacing A is the grating wavelength.

A o.o& 5

-5 0

2 4 X (x10e4m)

6

Fig. 2. (a) Interference pattern at the center of the PR crystal for two super-Gaussian beams for I, = lo4 W/m*, p = 0.2 s- ‘, and m = 0.1. (b) Space charge field E(x) corresponding to the interference pattern of Fig. 2a. Note that the amplitude of E(x) is nearly uniform inside the intensity pattern, and near the edges the electric fields are relatively large.

scribes the finite size of the interfering beams producing the intensity variation and k, = 2n/A, with A as the grating wavelength. This intensity pattern is realized by the superposition of two obliquely incident light beams. The factor g(x) in this paper is modeled by a super-Gaussian G(x) =exp[-(x-x,,)*/(+*] where x,,,~~ is the midpoint of the crystal of length x,,. Such a super-Gaussian effectively approximates a rect function. We have used x,, = 600 km, A = 10 p,m, and CT= 120 km. As an example, Fig. 2a shows the light intensity pattern for I, = lo4 W/m* and m = 0.1. The corresponding waveform of the space-charge electric field for p = 0 is shown in Fig. 2b. Note that in the central region of the interfering light beams the amplitudes of both the light intensity and space-charge electric field patterns are constants over 250 p,rn < x < 350 km. This is the spatial region from which we present results in the following discussion. We also note from Fig. 2b that the large space-charge electric fields develop near the edges of the intensity pattern. These large fields are essentially determined by the relatively large gradients in the intensity pattern, as explained later in

490

N. Singh et al. / Optics Communications 136 f 1997) 487-495

Section 3.2. However. we note that fairly large electric fields are created near the edges of the beams. The choice of the super-Gaussian is motivated by the following two reasons. The first is that the laser beams are of submillimeter thickness, and the above choice of CTgives a beam size of about 0.2 mm. The second factor deals with numerical implementation of the boundary conditions at the ends of the crystal. When the ends of the crystal are kept in the dark (see Fig. 2a), all the relevant quantities remain constant near the ends allowing for an easy implementation of the boundary conditions on p,. On the other hand, when the ends are also illuminated, the perturbation in ~~ leads to spurious oscillations near the boundaries, and they are troublesome numerically.

o.oA 0.0

0.2

0.4

0.6

0.8

1.0

t(ms) 2.0 L

1

0

20

40

60

80

100

%a) 3. Numerical results 3.1. Relatively low thermal excitation (p << SIO) and m << I

We first establish the validity of our numerical method discussed above by comparing numerical results against accurate and reliable analytical results. As pointed out in the introduction, most of the analytical studies dealing with the temporal evolution of the refractive index gratings in a PR material are based on perturbation methods in which it is assumed that the light intensity inhomogeneity is small allowing for the linearization of the Kukhtarev equations. This leads to time-dependent equations for the zero- and first-order quantities, the former equations depend on the uniform part of the intensity while the latter ones depend on the amplitude of the inhomogeneity in the intensity pattern. It turns out that the time-dependent zero-order quantities enter into the equation for the first-order ones. This leads to difficulty in analytically solving even the linearized equations. However, it turns out that for a set of parameters of the problem the zero-order quantities reach steady-state at a much faster time scale than the first-order ones. In this case, it is possible to assume that at the time-scale of the evolution of the first-order quantities, the zero-order quantities are constant. Such linearized solutions are given by Yeh [3] for /3 +C SI,, +C y, N, and E/.L< N/E -=++ y,.N, where N, = (NJ - NJ SIJY,N,>

(1’)

and it is the steady-state, zero-order electron concentration. N, evolves with a time constant T,.E (y,N,)- ‘, which is the recombination time of electrons with the ionized donors. On the other hand, the first-order electron density evolves with a time constant rdi G (y/e~N,)y,N,/SI,,

(‘2)

where 7di is the dielectric relaxation time with an electron concentration N,. When r, < rdi, N, becomes nearly constant for the purpose of evolution of the first-order

Kx( rad) Fig. 3. Comparison of numerical and Yeh’s analytical transient growth of the electric field amplitude for BSO (a) I, = IO3 W/m’, and (b) I, = lo4 W/m2 when p = 0 and m = 0.1. The numerical and analytical solutions are shown by solid and broken line curves, respectively. (c) Comparison of numerical and Yeh’s analytical steady-state waveforms of E(x) for the parameters of (b). The waveforms are shown over two grating periods in the central region of the crystal.

electron density and the first-order space-charge electric field, allowing for an analytical solution for the temporal evolution of the refractive index waveform in a PR crystal. Although not explicitly discussed, the above assumptions are the basis of all perturbational analytical studies on the coupling of two laser beams in a PR material. The condition T, +X rdi can be cast in terms of intensity I, as ‘yrN,

1

IO< I,, = s

(

Na

_

N,)

( dek) W/m*.

(13)

This implies that the linear analytical results on the time evolution are valid only for sufficiently low intensities I,, < I,,, where the critical intensity I,, depends on the material properties. Therefore, in this section we present results for low intensities I < I,, and at the same time we assume that thermal excitation rate j3 = 0. As mentioned earlier, the comparison of such results against the corresponding analytical ones will demonstrate the accuracy of the numerical scheme. For BSO, rr G 6 ps, rdi = l/SIC, and I,, = IO4 W/m* or SIC, = 0. I. In order to check the validity of our numerical results we compare them against those from Yeh [3] for m=O.l, SI”
solid line shows the numerical results while the dashed curve is the corresponding analytical result from Yeh [3]. The latter result is given by E,(t) = E,(l - ee’/‘y) where E, is the amplitude of the electric field waveform and TV is the time constant given by 131

(1 +E&) ( 1+ %/Eq>’

(‘4)

ry = TLir

E, g eNa( Nd - NJ/( k, l Nd), and p. For our parameters E,,/E, < 1 and T? = rdi(l + E,/E,) = (1.4/Q,) KS. Fig. 3a shows that for SI, = 0.01 the time evolution of the space-charge electric field agrees well with the analytical prediction. This agreement, along with some others described later, shows that our numerical method is sufficiently accurate. When the light intensity I, is increased to SI, = 0.1, the agreement with the time evolution for E, deteriorates. This is shown in Fig. 3b. Numerically we see a somewhat slower evolution; we find that TV= 14 ps while the numerically obtained time constant r,, = 18 p.s. As I, increases above the critical intensity I,,, T,, and ry increasingly deviate from each other and we find that T, > r,, > rr. We performed a series of calculations by varying I, to study the dependence of r on I, for p < Sl,. Fig. 4 shows the comparison between the numerically determined time constant T, against ry for BSO, and also for KNbO,, which is discussed later. The numerical results are shown by the symbol ’ + ‘. We find numerically that for BSO there is an excellent agreement for I, < I,, = lo4 W/m2 or SI, < 0. I, i.e., ~a 1; ’ in this range. When SI, > 0.1, the numerically obtained time constant is generally larger than ry. The dependence of T on I, (> I,,) can be approximately described by r a 1; y, with y being numerically estimated to be 0.7. We point out that the results shown here for p = 0 are likely to be modified when /3 is non-zero, as discussed in a later section. Fig. 3c shows the comparison of the analytically and

where

Ed = k, k,T/e,

E, = ‘y, N,A/k,

numerically obtained waveforms of the space-charge electric fields; as in Figs. 3a and 3b, the solid and broken lines represent numerical and analytical results, respectively. It is seen that the waveform agrees quite well. This agreement between steady-state waveforms is found to be true generally even when the time evolution obtained numerically is generally slower than the analytically predicted evolution for I > I,,. We have considered some other PR materials as well. For BaTiO,, S = lo-’ m2/J, y, = 5 X lO-‘4 m2 s, NA = 2 X 10z2 m-j, Nd= 100 NA, p = IO-’ m2/V s, and E = 2OOe,, for which I,, = lOI W/m2. This intensity is too high to be realized experimentally. Therefore, the analytical results [I] based on Kukhtarev equations should be applicable to BaTiO, over a wide range of intensity. However, it has been found that BaTiO, is not adequately described by the original set of Kukhtarev equations as given by Eqs. (l)-(3); these need to be supplemented by equations for holes as charge carriers and different levels of traps [5]. For KNbO,, S = 0.03 m2/J, yr = IO- I5 rn’ s, N, = 5 X 102’ rne3, Nd = 6N,, /L = 8 X lop6 m2/V s, and E= 55c,, for which I,, E IO* W/m2. As noted earlier, a comparison of the numerically and analytically obtained time constants for KNbO, is shown in Fig. 4, which reveals that the numerical and analytical time constants are in excellent agreement for II lo7 W/m2. For I,, = IO8 W/m2 the numerically obtained time constant is 0.55 us, while its analytical value r, = 0.2 JLS. Calculations were not performed for I > IO*- W/m2 because for such large intensities the time steps for obtaining the solution become excessively small. This is particularly true for KNbO, for which S = 0.03 m2/J. This large value of S in combination with I, > IO* W/m2 makes the numerical solution of Eq. (1) a difficult task. For a PR material with relatively low values of S and I,,, the numerical method described here can be an invaluable tool. 3.2. Effects of thermal excitation

(non-zero

/3) and arbi-

trary m

123456789 log,o(Io)

(w/m”)

Fig. 4. Comparison of analytical CT,.) (broken lines) and numerical (T,,) time constants as functions of the intensity I, for m = 0.1 and /3 = 0, for BSO and KNbO,. The dotted and dashed lines are the analytical results, while the symbol ’ + ’ shows numerical values. Note T,, becomes systematically larger than 7,” above the critical for intensity I, - lo4 W/m’ for BSO and I,- 10’ W/m* KNbO,.

The results described above are for small values of m because of the common practice of linearizing the Kukhtarev equations for analytical studies. Here we present results for arbitrary values of m including the effects of thermal excitation which implies p f 0. We first performed calculations for several values of m by still assuming that /3 < SI,. These calculations are intended to determine the range of values of m for which the linear approximation can be useful. Fig. 5 shows the comparison between the analytical (solid line) and numerical (*) results for the amplitude of the space-charge electric field. The analytical result is from Yeh [3], and it is the amplitude of the electric field waveform given by E, = mE,,( 1 + Ey/Ed)

_

’,

(‘5)

492

N. Singh et al./ Optics Communicutions

136 (1997) 487-495

where E, and Ed are previously defined. As expected, the analytical result predicts a linear dependence of E, on m. We see from Fig. 5 that for m = 0.1, the analytical and numerical results are nearly identical with an error less than 1%. However, the error steadily increases with increasing m; when m = 0.5 the analytical theory underestimates the electric field amplitude by about 20%. Our numerical study therefore shows that the analytical theory can only be useful in predicting the amplitude of the space-charge electric field for m < 0.4 with an error of less than 10%. It turns out that for p # 0 and for an inhomogeneity of arbitrarily large strength in 1(x), a steady-state expression for the space-charge electric field can be derived for relatively large grating periods compared to the Debye length (u 0.25 km for I, = IO4 W/m* for BSO). This was first done by Banerjee et al. [lo] and E(x) is given by

[p+sqx)]-

I

a1 s-g.

The foregoing expression shows that E(x) is proportional to the gradient in the intensity pattern; the large electric fields (Fig. 2b) near the edges of the intensity pattern shown in Fig. 2a are the result of this behavior of E(x). For the sinusoidal intensity waveform I(x) = I,( 1 + m x cos kx), mSI, k, sink, x /3+S1,(1

(‘7)

+mcosk,x)

This expression

reduces to that of Yeh [3] with an amplitude given by (15) when p +z Sl,, m -=z 1 and E, -=KE,. The latter condition implies that the grating wave number k -C k,, where k, is the Debye wave number for the material. When p -=c SI, and m + 1, E(x) -+ 3~ when k, x --f (2n + l)a, with n as an integer. As remarked earlier about the agreement between analytical and numerical results on the steady-state waveform of E(x), we further demonstrate here that the steady-state

2.0x104

*

analytmal * numerica,

*

0.4 0.6 Modulation 1ndex.m

Fig. 5. Comparison of numerical and Yeh’s analytical steady-state electric field amplitude E, versus m for BSO when p = 0 s-’ and I,= lo4 W/m’.

-0 -10

-il

m=l

0

I x Kx (rad)

2x

37T

Fig. 6. (a) Interference pattern shown over two grating periods in the central region of a BSO crystal for m = I, 6 = 1 s- ’ , I, = IO6 W/m’. (b) Comparison of analytic and numerical steady-state waveforms of E(x) for the interference pattern in (a). (c) Numerically obtained steady-state waveforms of E(x) for p = I sm’ , and I, = IO6 W/m2 when m = 0.1, 0.5 and 1. In order to show the curves clearly for m = 0.5, and 0.1, the values are amplified by 2 and 5, respectively. (d) Space charge distribution for the electric field shown in (b).

expression (17) for arbitrary roles of m and /3 agrees remarkably well with the corresponding numerical results. Fig. 6a shows the intensity pattern over two grating periods for I, = IO6 W/m’, p = 1, and m = 1. Fig. 6b shows the comparison between the analytical (Eq. (17)) and numerical waveforms of the steady-state electric field for m= 1, p= 1 s-‘, and I,= lo6 W/m*. The two waveforms, shown by the solid and broken liens for the numerical and analytical results, respectively, match almost perfectly. We note from (17) that the waveforms for ECx) and the intensity I(x) are out of phase by 90” only for m + 0. As m increases, the electric field waveform becomes increasingly nonsinusoidal and the phase difference between the maxima of 1E(x)1 and I(x) depends on the relative magnitudes of /?, Sl,, and the modulation index m. It is easy to show from (17) that the maxima and minima in the waveform of E(x) occur when k,x=

(2/+

l)rrfcos-’

(‘8)

where / is an integer. The minus sign in (18) corresponds to the maxima in E(x), while the plus sign is for the

minima. When tn -+ 0, (18) yields the well known result that I( x> and I!?(x> differ in phase by 90”. The locations of the maxima and minima of E(x) lie below and above the locations of the minima of I(x), respectively. The location of the adjacent pairs of maximum and minimum tend to approach the minimum of I(x) as m + 1. Fig. 6c shows the waveforms for m = 0.1, 0.5, and I, when /I= 1 SC’, and I, = lo6 W/m2. In order to show the curves clearly for m = 0.1 and 0.5, the values are amplified by factor of 5 and 2, respectively. Note that the waveform distorts from being nearly sinusoidal to highly nonsinusoidal as m increases from 0.1 to I. Fig. 6b and 6c show that as m + 1 the electric field waveform steepens greatly near b = (2n + 1)~. This steepening suggests that the grating created by the sinusoidal intensity pattern with a periodicity of A has a rich variety of wave numbers, which play crucial roles in coupling of waves affected by the grating. One obvious application could be readout of strongly modulated gratings stored in the PR material with the reading beam at non-Bragg (as defined by the fundamental grating period) incidence. The generation of the higher wave numbers (> k, = 2n/A) is the consequence of the nonlinearity of the system. The steepening of the electric field waveform near kx = (2n + 1)TI indicates that a large space charge builds up near such locations inside the crystal. This space charge is shown in Fig. 6d. Note that a large negative space charge is created near the nulls in the intensity pattern I(x). The physical reason for this charge build up is as follows. Since near the nulls photoionization does not take place, the ionized donor concentration remains at its initial value Nl(r = 0). Therefore, the electrons reaching there by the process of diffusion are not affected by the recombination process. Even if the recombination occurs, the negative space charge builds up there due to the unneutralized electrons captured by the acceptor sites (N,). The saturation in the space charge near the nulls occur only when the electron diffusion is balanced by the opposing process associated with the electronic drift in the evolving self-consistent electric field, that is, pN,E=

-DaN,/ax,

or (19) The maximum electric field strengths seen from Fig. 6b and 6c are well borne out by the estimates given by (19). For small m, N,(x) is sinusoidal and (19) gives (15) for Ed<< E,. 3.3. Temporal evolution for arbitrary p and m For the case of arbitrary p and m, it is difficult to derive an analytical expression for the temporal evolution

1.00

r ,’ :.

.

0.75 i::’ ,. 0.50 ;: ”

----_

:____

*1=10.0.f9=1.0.m=1.0

,:

::IIL--_-l 0

10

20 30 t(w)

40

50

Fig. 7. Time evolution of the electric field at three locations. Solid line shows the evolution of the maximum electric field at k, x = 2~/ +0.86~. Dotted and dashed lines show the evolution of the electric field at k,Vx = 27rL’ + 7r/2 and k,yx = 2~!’ + r/8, respectively. of the space-charge electric field. However, Jarem and Banerjee [6] have given an approximate estimate for the time constant for the evolution of the electric field by making a variety of approximations. As usual, one of the major approximations is that the zero-order charge carrier density reaches steady state at a time scale much faster than the space-charge electric field. Their approximate time constant is given by

The above expression shows that rjb is a function of I(x), that is, the different parts of the electric field waveform evolve differently and the slowest evolution occurs near the minima of I(x). This is qualitatively borne out by our numerical calculations. However, we will shortly show that TV,, greatly underestimates the actual time constant. Fig. 7 shows the temporal evolution of the space-charge electric field at three locations for m = 1, I, = lo6 and p = 1; the solid curve shows the temporal evolution of the maximum electric field at k, x = 2/‘7r + 0.86~. The dotted curve is for the electric field at k, x = 2~TTef 7r/2 and the dashed curve is for k,x= 2~!+ 7r/8. The electric fields shown are normalized with respect to their respective maximum values during the course of their evolution to the steady state. As pointed out earlier, Fig. 7 shows that the closer the points are to a maximum of I(x) the faster the fields reach the steady state. Fig. 7 also shows that the electric fields at k,x= 2aL’+ 7r/8 and 27r/+ 7r/2 first overshoot above their respective steady-state values before settling down at slightly smaller values. This is caused by a relatively slow redistribution of the initially created space charges as the spatial distribution of electrons evolves self-consistently through drift and diffusion and the associated changes in the production rate through photoionization and recombination. This redistribution of space charge is associated with the weakening and strengthening of the electric field in the neighborhoods of the maxima and minima of I(X), respectively, leading to the formation of the sawtooth waveform.

494

N. Singh et ul./

Oprics Communications 136 (1997) 487-495

Table 2 Effects of p and m on the time constant 7 for I, = lo6 W/m’

BSO

x

x

0

p @‘I

0. I

1

10

0.1 m = 0.5 T??= 1.0

1.00 1.18 ?

I .25 I .45

1.56 I .83 4.10

m=

0.0

0.2

0.4

0.6 n-l

0.8

I.0

Fig. 8. Comparison of numerical (X) and analytical (0) time constants including the effect of p as a function of the modulation index m for BSO when I, = lo6 W/m’. Note that the numerically obtained time constants are about an order of magnitude larger than the analytical values for all m.

Fig. 8 shows the time constants as a function of m for the evolution of the maximum electric fields at the locations given by (18) for p = 1 SC’ and I, = IO6 W/m*. The data points shown by crosses are the results from our numerical calculations, while the squares show the time constants from the analytical formula given by (20). Fig. 8 clearly reveals that the analytical values of the time constant are nearly by an order of magnitude smaller than the numerical values. As a matter of fact the results in Fig. 8 are extensions of the results shown in Fig. 4 to larger values of M. Since for /? = 1 SC’ and I, = lo6 W/m’. the approximation p < SI, is well met, the point in Fig. 4 for I, = IO6 W/m* is the same as that for m = 0.1 in Fig. 8. As discussed earlier, it is apparent that the analytical result for the time constant given by (20) is a good approximation only for a sufficiently low ionization rate depending on both p and I,. An estimate of relatively low light intensities for which the ionization rate is sufficiently low can be obtained from the requirement that the value of T,,, > T,, which yields

(21) This defines a new critical intensity including the effects of the thermal excitation /3 and the modulation index nr, that is,

(22) where I0 is given by (13). Normally p lies in the range 0.1 < p < 2. For BSO with S = IO-j, I,( p,m) = 0, (20) is not valid even for relatively small intensities. For KNbO, it turns out that I,, z=+/3/S and I,( /3,m) = I,,,/(1 + m).

4.93

Table 1 shows the numerically determined time constants for BSO when p = I s - ‘, and M and I, are varied. It is seen that for I, = IO4 W/m2, 7, is relatively constant at about 7.3 ps for all values of m, and it is slightly larger than the recombination time 7, = 6 KS in contrast TV,,z 0.82 ps. The independence of T,, on m is simply the consequence of /3 >> H(x) for I, = lo4 W/m’. As the intensity increases the time constant decreases for a given m, but for a given intensity the time constant increases with increasing m. This latter effect takes place only when St,, 2 j3. We point out that the numerical results discussed here are for BSO for which p 2 St,,,, and (20) is not applicable even for small intensities because I,( /3,m) z 0. On the other hand, for example for KNbO,, I,( p,nz) = I,,/(1 + m) and (20) can be useful for intensities I,, < I,,,/(1 + m), as shown in Fig. 4. Table 2 shows the variation of the time constant with p and m when the intensity I,, is fixed at IO6 W/m*, for which SI,, = 10. It is seen that for /I < SI,, (also see Table I> the time constant sharply increases as m + 1. Furthermore, it turns out that for /3 < Sl,,, and m G 1, the evolution of the maximum space-charge electric field is too slow to reach the steady state within a reasonable computational time. This is indicated by a question mark for ,L3= 0.1 and m = I in Table 2. We notice a rather curious effect of increasing time constant with increasing /3 as long as m is not close to unity; although the increase in the time constant is only slight for 0.1 < /3 I IO. As /3 increases, so does the zero-order electron density making the dielectric relaxation time smaller than the recombination time. In such a case the time evolution is controlled by the combined action of recombination and the charge relaxation. For m G 1, the time constant decreases with increasing /3 because for small /3 < SZ,, the time constant tends to be large near I(X) = 0, where the electric field maximum occurs.

4. Conclusion and discussion Table I Effects of I, and m on the time constant

T in ps for p of I SC’

I, (W/m* )

IO4

105

106

m= 0.1 m = 0.5 m= I.0

7.30 7.30 7.30

4.58 4.83 5.63

I .25 1.45 4.93

The main purpose of this paper is to report numerical solutions of Kukhtarev equations for photorefractive materials. We have considered the case of two interfering light beams for producing a refractive index grating. For the interference pattern I( x1 = I,[ 1 + m cos(k, xl], it is possible to linearize the Kukhtarev equations for m -=x 1, for

N. Singh et al./ Optics Communications

which several authors have analytically studied the temporal evolution of the space-charge electric field. We have shown here that even for this linear regime, the analytical results have limited validity. Specifically, we show that analytical results are valid over a limited range of the light intensity I,. When Z, is sufficiently large, the analytical theory fails because the dielectric relaxation time becomes smaller than the electron recombination time. In such a case the time evolutions of zero- and first-order quantities become coupled, making it difficult to obtain time-dependent analytical solutions. We have given an expression for the critical intensity I, above which analytical theory is not valid. I, depends on material properties including the thermal ionization rate and the modulation index m. We have further shown that the linear approximation is good for m < 0.3 within a tolerance of about 10% error in the estimates of the space-charge electric field. For large m, the linear approximation fails, and there is a lack of studies giving rigorous solutions for the temporal evolution of the space-charge electric field. Jarem and Banerjee [6] have given an estimate of the time constant for the temporal evolution. This estimate is based on the approximation that the zero-order electron density evolves at a much faster rate than the higher-order perturbation in the electric field. We have shown here that for BSO this condition is not met even at low light intensities. Consequently, numerical results on the time evolution yield time constants by an order of magnitude larger than that given by the analytical result 161. Our study suggests that for some other materials such as BaTiO, and KNbO,, it is possible that at sufficiently low intensities I < I,. where I, is given by Rq. (221, the analytical result in (20) may give a useful estimate for the time constant. We have shown that when m is not sufficiently small, different parts of the electric field waveform evolve with different time scales, with slowest evolution occurring for the maximum electric fields occurring near the minima of I(x). The fastest evolution occurs near the maxima of I(X), where electric fields are relatively small. We found that the steady-state analytical solutions for E(x) agree well with the corresponding numerical solutions. It is found that as the modulation index m increases, the steady-state waveform of E(x) increasingly deviates from the sinusoidal waveform given by the light intensity pattern I(X). In the limit m + I, the waveform becomes a sawtooth with maximum electric fields occurring near the minima of I(x); the location of such maximum electric

136 (19971487-495

495

fields are given by Eq. (18). The strength of the maximum electric field is determined by the balance between the diffusion and drift of electrons. It is worth mentioning that the numerical method is applied here to the original set of Kukhtarev equations with electrons as the only charge carrier, but it can be easily extended to the modified form of Kukhtarev equations including both electrons and holes as charge carriers and different levels of shallow and deep traps for them [2,5]. In the latter cases, the analytical solutions are available only for the linear case of m < 1, and even in that case only for situations where the zero-order charge carrier concentrations evolve at a rate much faster than that for the first-order space-charge electric field. The modified form of the Kukhtarev equations await a numerical treatment to examine both linear and nonlinear regimes for the useful range of light intensities. In the present study we have treated a one-dimensional problem assuming a thin crystal. The numerical methods can be generalized to three dimensions, and when augmented by an appropriate equation for the modification in the light intensity, the numerical method can be a powerful tool in rigorous studies of PR devices.

References

[II M.V. Kukhtarev,

Pis’ma Zh. Tekh. Fiz. 2 (1976) I1 14; J. Tech Phys. Len. 2 (1976) 438; M.V. Kukhtarev. V.B. Markov, SC. Odulov, MS. Sorkin and V.L. Vinetskii, Ferroelectrics 22 (I 979) 949. Dl G.C. Valley, J. Appl. Phys. 59 (1986) 3363. Nonlinear Optics 131P. Yeh, Introduction to Photorefractive (Wiley, New York, 1993) ch. 3; G.C. Valley, IEEE Trans. Quantum Electron. QE- 19 ( 1983) 1637. 141 M.G. Moharam, T.K. Gaylord and R. Magnuson, J. Appl. Phys. 50 (1979) 5642. El D. Mahgerefteh and J. Feinberg, Phys. Rev. Len. 64 (1990) 2195; P. Tayebati, J. Opt. Sot. Am. B 9 (1992) 415. (61J.M. Jarem and P.P. Banerjee, Opt. Eng. 34 (1995) 2254. 171 M. Reiser, Comp. Methods App. Mech. Eng. I (1972) 17. 181R.D. Richtmyer and K.W. Morton, Difference Methods for Initial Value Problems (Wiley, New York, 1957). 191P.P. Banerjee, Q.W. Song and J.J. Liu, Optics Comm. 83 (1991) 115. 1101 P.P. Banerjee and R.M. Misra, Optics Comm. 100 (1993) 166.