Optics Communications 281 (2008) 5884–5888
Contents lists available at ScienceDirect
Optics Communications journal homepage: www.elsevier.com/locate/optcom
Photorefractive hologram writing with high modulation depth in photovoltaic media under different boundary conditions W.K. Lee *, T.S. Chan Department of Physics, The Chinese University of Hong Kong, Tai Po Road, Hong Kong, PR China
a r t i c l e
i n f o
Article history: Received 10 February 2008 Received in revised form 26 July 2008 Accepted 6 August 2008
PACS: 42.65.Hw 42.40.Eq
a b s t r a c t Photorefractive recording response to sinusoidal intensity pattern in photovoltaic media under different boundary conditions, namely, open-circuit, crystal-resistor-circuit, short-circuit and crystal-voltage supply-circuit, are investigated by solving steady-state Kukhtarev equations numerically. In high modulation depth m limit and when photovoltaic field becomes sufficiently large, we find that both the magnitude and imaginary part of fundamental space-charge field E1 as a function of m may transit from superlinear increment to sublinear increment for short-circuit condition. The role of photovoltaic field in holographic storage and two-wave mixing is also discussed in comparison with an earlier experimental literature [G. Cook et al., Opt. Commun. 192 (2001) 393]. We further generalize Gu et al.’s remark on photovoltaic field and dc space-charge field, which is for small m only [C. Gu et al., J. Appl. Phys. 69 (1991) 1167], to whole range of m in the limit donor density much larger than ionized donor density that the effect of photovoltaic field cannot be cancelled or considered as an externally applied field for an arbitrary non-zero m. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction Buildup of space-charge field in photorefractive media under illumination of light intensity pattern has been studied both experimentally and theoretically by many investigators. The key for analyzing this phenomenon is a set of coupled differential equations known as Kukhtarev material equations obtained in 1979 [1]. Early works solve the Kukhtarev material equations analytically by linearization technique in the limit of small modulation of light intensity m [1–3]. Consequence of this limit is that the analyses are not applicable to a number of hologram recording experiments that have been done in large m for higher diffraction efficiency. Several investigators have extended the theoretical work of small m limit under a number of other conditions [4–6], or by introducing an empirical function called modulation-depth-dependent-correction functions f(m) in studies of two-wave mixing [7,8]. For review of works on photorefractive response at large m prior to 1987, we recommend Section 1 of Ref. [6]. Review of works on various forms of f(m) can be found in Section 1 of Ref. [8]. Regrettably, owing to nonlinearity of Kukhtarev material equations, exact analytical solution of the space-charge field in large-m limit has not been reported. Investigators then seek numerical method that can solve Kukhtarev material equations efficiently for large m including m = 1. A commonly used semi-analytical method, known as Powell method, ex* Corresponding author. Tel.: +86 852 2609 6395; fax: +86 852 2603 5204. E-mail address:
[email protected] (W.K. Lee). 0030-4018/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2008.08.015
pands the space-charge field into truncated Fourier series, and then substituting the series into the Kukhtarev equations to find Fourier components from a set of resultant nonlinear algebraic equations [9] or similarly, expands the charge–carrier density into truncated Fourier series, and then substituting the series into the Kukhtarev equations with algebraic manipulation afterwards [16]. Several investigators have reported their studies on photorefractive response for large m using these calculations. These reports include photorefractive response in BSO crystal upon incidence of two frequencydetuned beams studied by Au and Solymar [9], photorefractive response in BSO crystal upon incidence of two beams of the same frequency studied by Bledowski et al. [10] and by Kukhtarev et al. [16], and photorefractive response in LiNbO3 crystal upon illumination of two laser beams of same frequency studied by Rupp et al. [11] and Serrano et al. [12]. The photovoltaic field Ep in BSO crystal was not taken into account in studies of Au and Solymar, Kukhtarev et al., and Bledowski et al. due to its insignificance. On the other hand, typical photovoltaic field of LiNbO3 ranges from several to more than 100 kV/cm and its effect on photorefractive response is not negligible even in small m limit, which has been studied theoretically and experimentally by Gu et al. Nonetheless, the value of photovoltaic field in Ref. [13] is much smaller than typical values. Also, dc space-charge field E0, which is related to boundary condition of the crystal, is set zero throughout Ref. [13]. We note that Kukhtarev et al. also reported a compact closed-form solution with photovoltaic field considered in 1997 [17]. However, the solution was found under the assumption that
W.K. Lee, T.S. Chan / Optics Communications 281 (2008) 5884–5888
photovoltaic current density is directly proportional to light intensity, which is true only when ionized donor density is much smaller than donor density or ionized donor density is approximately a constant. These insufficiencies of previous works motivate us to study the photorefractive response at large m in photovoltaic media under different boundary conditions, which remains essentially uninvestigated analytically and numerically with realistic parameters, by using a simple but powerful numerical method that is applicable to solve Kukhtarev equations with the least assumptions. By comparing our results with the case in which glass constant j = 0, the role of photovoltaic field in holographic storage and two-wave mixing is discussed.
2. Theory and method The photorefractive response to illumination of light in photovoltaic media is described by the Kukhtarev material equations [1]:
oNiD ¼ SðI þ Id ÞðN D NiD Þ cnNiD ; ot on oN iD 1 ¼ r j; ot e ot j ¼ enlE þ kB T lrn þ jsðND NiD ÞI^c;
r eE ¼ eðn þ NA
NiD Þ;
ð1Þ ð2Þ ð3Þ ð4Þ
where s is the photoionization cross section, c the charge carrier recombination rate, e the magnitude of electron charge and l the mobility of the charge carrier, ND the donor density, N iD the ionized donor density, NA the acceptor density, n the charge-carrier density, b the dark generation rate, j the electric current density, j the glass constant, cˆ the +c-axis unit vector, E the space-charge field, kB the Boltzmann’s constant, T the temperature, t the time, e the dielectric tensor, I the light intensity, Id = b/s the dark irradiance. For steady state o/ot|t = 1 = 0, total current density J is equal to j and in the limit that all variables changes predominantly along caxis, with the assumptions n NA and sðI þ Id ÞN D jcðN iD Þ2 1, one may simplify Eqs. (1)–(4) into the equation below by following key steps in deriving Eq. (3) in Ref. [14]:
# " c e oE 1 jc J þE ðI þ Id Þ NA þ ND NA þI e oy esl e oy el " # kB T o lnðI þ Id Þ e oE 1 e oE 1 e o2 E ND NA þ ND NA þ e oy e oy e oy e2 oy2 1
e oE
¼ 0;
ð5Þ
where I and E are functions of y that is along c-axis, and J is a constant (by equation of continuity) to be determined by boundary condition, which is given by Ohm’s law:
Va
Z
0
Edy ¼ JSR;
ð6Þ
L
where Va is the externally applied voltage across the crystalline faces perpendicular to c-axis, R is the resistance of external resistor, L is the transverse length along c-axis, and S is the surface area of crystalline faces perpendicular to c-axis. For a given set of crystal parameters and normalized intensity profile, the solution of E(y) depends on two parameters, namely, J/Id and Imax/Id. When the value of J is given, E(y) becomes the only function to be solved in Jacobi relaxation method. For open-circuit condition, J = 0 everywhere, which leads to dc space-charge field E0 EP when Id is negligible and N D N iD everywhere. For crystal-resistor circuit condition, Va = 0. Short-circuit condition RL (R = 0), leads to 0 Edy ¼ 0 and hence E0 = 0. The value of J is ad-
5885
justed such that the solution of E satisfies this condition. Similar procedures are used for solving non-zero R and/or Va – 0. The main advantage of relaxation method over Powell method is its simplicity and flexibility. Similar to Powell method, a good initial guess can be sometimes important to obtain a solution. However, relaxation method is not limited to solving periodic solution but still; periodic solution can be simply obtained by applying periodic boundary condition at computational windows. Also, the incident light can have intensity of arbitrary form. Accompanied with Helmholtz equation, simulation of two-wave mixing of Gaussian beams and even more complicated beam profiles becomes possible. In this communication, the following parameters of LiNbO3 are used in the calculations [15]: l = 0.01 cm2 V1 s1, Id = 1 mW cm2, s = 3.6 cm2 s1 W1, ND = 1018 cm3, NA = 5 1016 cm3, sR = 27.6 ps, c = 1/(sRNF) = 7.25 107 cm3 s1, T = 300 K, j = 12 A ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cm s, e = 30 e0. We obtain Debye screening wave number 1028p ffi K D ¼ e2 N A ðN D N A ÞjekB TN D ¼ 3:34 107 m1 , Debye screening distance LD ¼ 2p=kD ¼ 0:188 lm, and Ep ¼ kcN A =el ¼ 27:2 kV=cm. In darkness, n sId ðN D N A Þ=cN A 106 cm3 N A 1016 cm3 and sId N D =cN 2A 1012 1. The two assumptions are checked and ensured case by case in this communication. We shall solve E(y) for the following intensity profile I(y):
IðyÞ ¼ I0 ½1 þ mcosðKyÞ ;
ð7Þ
where K is the grating spacing, and I0 is the dc light intensity, K ¼ 2p=K. We choose K = 0.188 lm, which gives K ¼ kD , the diffusion field Ed ¼ kB TK=e = 8.58 kV/cm, and the saturation field Eq ¼ eeN A ðN D NA Þ=KND ¼ 8:58 kV=cm. Moreover, we choose I0 = 500 Id, which is large enough such that space-charge field saturates. 3. Results and discussion Earlier works expand space-charge field into truncated Fourier series for solving Kukhtarev equations and the first Fourier component as functions of some variables are usually plotted. For better comparison with earlier works, we shall decompose space-charge field solutions, which are obtained by Relaxation method and originally presented in a list of field values of different positions, into Fourier series:
E0 þ
n 1X ½Er expðjrKyÞ þ Er expðjrKYÞ ; 2 r¼1
ð8Þ
where Er of the rth Fourier component of the space-charge field, Er ¼ E r . And we define ur to be the phase angle (relative to the light intensity) of the rth Fourier component of the space-charge field. Results for m = 0.1, 0.3, 0.5, 0.7, 0.9 are plotted in Fig. 1 for a range of dc space-charge field E0. Instead of plotting Im(E1), which contributes to energy transfer in two-wave mixing, and field amplitude |E1|, we plot |E1|/m and Im(E1)/m, which are independent of m when m is sufficiently small, for better scale arrangement. The asymmetrical nature of |E1| and Im(E1) shown in Fig. 1 is the characteristic of the presence of photovoltaic field. For m = 0.1, results are very close to those of small modulation, which are given by the theory of Gu et al. [3]. In small m limit, ionized donor density N iD is usually assumed equal to NA in by Eq. (4). As m becomes large, intensity changes more rapidly along y. Larger intensity gradient leads more photogenerated free charge carriers (electrons) to diffuse from bright regions to dark regions. Hence, for dark regions, more ionized donors recombine with electrons and ‘depletion’ of NiD is more significant when m is close or equal to 1. The break down of the assumption N iD N A leads to the development of asymmetry of space-charge field (Fig. 2) and |E1|, Im(E1) gradually deviate from the theory of Gu et al. as expected. Bledowski et al. have found that field amplitude of fundamental
5886
W.K. Lee, T.S. Chan / Optics Communications 281 (2008) 5884–5888
a
18
a
E0=-27.2kV/cm
16
m=0.3
6
14
|E1|/m (kV/cm)
Eac (kv/cm) 1
m=0.1 and Gu's result
12
4
10
m=0.5 m=0.7
m=0.9
2
8 6
0.2
4 E0=0
2
0.6
10
0.8
1
yΛ
-2
0 -150
-100
-50
0
50
100
b
E0(kV/cm)
b 100 m=0.1,0.3,0.5,0.7,0.9 and Gu's result
80
1
1.5
10
1
60
φ1 (deg)
N D NA
1.25
0.75
40
0.5 20
0.25
0
0.2
-20 -150
-100
-50
0
50
100
E0(kV/cm)
c
14
0.4
0.6
0.8
1
yΛ
Fig. 2. J = 0. (a) AC component of space-charge field profile. (b) N iD =N A profile. m values of curves 1 to 10: 1, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1. The curves are well ordered, and for clarity, labels of curves 2 to 9 are suppressed.
E0=-27.2kV/cm E0=0
m=0.3
12
m=0.1 and Gu's result
10
m=0.9 m=0.5 m=0.7
8 1
Im(E )/m (kV/cm)
i
1.75
6 4 2 Ep=-27.2kV/cm
0 -100
0
100
E0(kV/cm) Fig. 1. Ep = 27.2 kV/cm. (a) |E1|/m versus E0. (b) u1 versus E0. (c) Im(E1)/m versus E0.
Fourier component of space-charge field as a function of m increases superlinearly in non-photovoltaic media with lots of combinations of K and E0 including short-circuit condition E0 = 0 [10]. Serrano et al. have found similar results for E0 = 0 in photovoltaic media, though the value of Ep they used is much smaller than typical values [13]. In contrary to previous reports, with typical photovoltaic field value of LiNbO3, we have found sublinear relation for |E1| and Im(E1) versus m for E0 = 0 and depending on the value of E0, resultant |E1| or Im(E1) can be larger or smaller than results obtained in small m limit, which is clearly depicted in Fig. 1a and c. Further increment of E0 up to 80 kV/cm does not change the sublinear relation. However, decrement of E0 down to about 10 kV/cm and 5 kV/cm reverse the sublinear relation to superlinear relation for |E1| and Im(E1), respectively. There exists such a range of E0 that at large m, beam coupling due to fundamental harmonic is weaker than that predicted by small m theory while field amplitude of fundamental harmonics is larger than that predicted by small m theory. Further decrement of E0 down to 140 kV/cm does not change the superlinear relation and results of open-circuit condition (J = 0, E0 Ep = 27.2 kV/cm) is in the superlinear regime. For the parameters we used, the theory of Gu
et al. predicts that there exist some values of E0 at which |E1|, Im(E1) maximized or minimized. We find that for large m, |E1|, Im(E1) also maximize or minimize at approximate the same values. We repeat the same calculations except setting j = 0 for studying the photovoltaic field free case (Fig. 3). We have found superlinear relation for |E1| and Im(E1) versus m for E0 up to 120 kV/cm. No similar sublinear relation is found for a wide range of E0. Photovoltaic field plays an important role in the superlinear-sublinear relation transition reported above. We remark that the phase angle u1, under all conditions, is insensitive to m. We are particularly interested in open-circuit and short-circuit boundary conditions. These two conditions are the most commonly used configurations in experiments. We vary j(Ep) so as to further investigate its role in recording process under the two experimental conditions. For open-circuit condition (Fig. 4), for same m no matter how large it is, larger Ep leads to smaller |E1| and Im(E1), resulting smaller field amplitude and weaker beam coupling. At the same time, we have superlinear relation for |E1|, Im(E1) versus m for Ep up to 27.2 kV/cm. The short-circuit case is more interesting (Fig. 5). We have found superlinear relation for |E1|, Im(E1) versus m when Ep is small. It transits from superlinear relation to sublinear relation at Ep 16 kV/cm. For same m no matter how large it is, larger Ep leads to larger |E1|. However, whether larger Ep will lead to stronger beam coupling depends on the value of m. For m = 0.1, 0.3, 0.5, larger Ep will lead to larger Im(E1) and hence stronger beam coupling while for m = 0.7, 0.9, larger Ep will lead to smaller Im(E1) and hence weaker beam coupling. It is interesting to note the phase angle u1 (20° at Ep = 27.2 kV/cm) obtained by us is close that obtained experimentally by Rupp et al. [11] in the low doping limit. To summarize, for open-circuit condition, existence of photovoltaic field leads to smaller field amplitude and tends to weaken the energy transfer due to fundamental harmonics. For short-circuit condition, existence of photovoltaic field leads to larger field amplitude. However, as field amplitude increases more slowly at large m but phase difference between intensity and fundamental
5887
W.K. Lee, T.S. Chan / Optics Communications 281 (2008) 5884–5888
a
a
9
J=0, E0 ≅ -Ep
5.0
m=0.1 and Gu's result
|E 1|/m (kV/cm)
8
|E1|/m (kV/cm)
5.5
m=0.9
m=0.7
7
m=0.5
6 5
Ep=0
4.5 4.0 m=0.3
3.5 m=0.9
3.0
2.0 0
20
40
60 70 |E0| (kV/cm)
100
m=0.1 and Gu's result
2.5
4
0
5
10
120
15
20
25
30
Ep (kV/cm)
b
b
m=0.7 m=0.5
J=0, E0 ≅ -Ep
90
90
m=0.1,0.3,0.5,0.7,0.9 and Gu's result
80
70 1
φ1 (deg)
85
φ (deg)
80
60 50
75 40
Ep=0 70
m=0.1,0.3,0.5,0.7,0.9 and Gu's result
30 0 0
20
40
60
70
100
E0 (kV/cm) 9
c
m=0.7
7
m=0.5
6
25
30
J=0, E0 ≅ -Ep m=0.3
4.5 4.0 m=0.5
3.5 3.0 2.5 2.0
5
20
5.0
1
m=0.1 and Gu's result
15
5.5
m=0.9
8
10
Ep (kV/cm)
Im(E )/m (kV/cm)
Im(E1)/m (kV/cm)
c
5
120
m=0.1 and Gu's result
m=0.9 m=0.7
1.5
Ep=0 4
0 0
20
40
60
70
100
120
E0 (kV/cm) Fig. 3. Ep = 0. (a) |E1|/m versus E0. (b) u1 versus E0. (c) Im(E1)/m versus E0.
harmonic of space-charge field decreases dramatically as photovoltaic field increases, it leads to stronger beam coupling for small m and weaker beam coupling for large m. For the parameters we used and the two boundary conditions we considered, photovoltaic field only enhance beam coupling for short-circuit condition at relatively small m. Cook et al. believe that large power transfer in two-wave mixing in LiNbO3 is due to photovoltaic effect and they have reported experimental evidence to support their belief [13]. They wrote, ‘‘From Fig. 3, at about 140 kV/cm, the DOD reaches a minimum where the applied field just cancels out the generated photovoltaic field. At this point, the beam coupling is dominated by the diffusion mechanism”. To summarize, the main idea of their experiment was to cancel the effect of photovoltaic field by an externally applied field and to compare it with those cases with the presence of photovoltaic field. Before agreeing on their conclusion, it is essential to think carefully whether the effect of photovoltaic field can be cancelled out by an externally applied field. Gu et al., who first wrote down the equa-
5
10
15
20
25
30
Ep (kV/cm) Fig. 4. J = 0. E0 Ep. (a) |E1|/m versus Ep. (b) u1 versus Ep. (c) Im(E1)/m versus Ep.
tions (Eqs. (1) and (2) of Ref. [13]) that Cook et al. made use of, which is valid in the small m limit, remarked that the solution of spacecharge field cannot be expressed as a linear combination of photovoltaic field and E0. The consequence is that photovoltaic field cannot be thought as an equivalent applied field. This remark actually contradicts with what is claimed to be done, namely, cancellation of the effect of photovoltaic field by an externally applied field, in the experiment of Cook et al. Though the experimental results are in good agreement with fitted theoretical results, as long as m at the input face of crystal is not small and the photovoltaic effect cannot be convincingly proved to have been cancelled, the conclusion of Cook et al. should be considered with some caution. It is interesting to ask whether the remark of Gu et al. can be extended to large m limit. In the limit N D N iD , Eq. (5) with the use of Eq. (4), reduces to the following equation:
e oE c jc ðI þ Id Þ1 NA þ þI J þE e oy eslND el " # 1 kB T olnðI þ Id Þ e oE e o2 E NA þ ¼ 0: þ e oy e oy e2 oy2
ð9Þ
5888
W.K. Lee, T.S. Chan / Optics Communications 281 (2008) 5884–5888
a
E0=0
14
m=0.3 m=0.5 m=0.7 m=0.9
m=0.1 and Gu's result
|E1 |/m (kV/cm)
12 10 8 6 4 0
5
15
10
20
30
25
Ep (kV/cm)
b
E0=0
90 80
φ 1 (deg)
70 60 50 40
m=0.1,0.3,0.5,0.7,0.9 and Gu's result
30 20 0
5
10
15
20
25
30
Ep (kV/cm)
Im(E1)/m (kV/cm)
c
5.4
m=0.1 and Gu's
E0=0
5.2
m=0.3
5.0 m=0.5
4.8 4.6
Kukhtarev equations numerically with realistic parameters. Our analyses allow estimations of diffraction efficiency in experiments carried out in photorefractive media with strong photovoltaic field under different boundary conditions. In high modulation depth m limit and when photovoltaic field becomes sufficiently large, we find that both the magnitude and imaginary part of fundamental space-charge field E1 as a function of m may transit from superlinear increment to sublinear increment for short-circuit condition. This finding suggests a sublinear modulation-depth-dependentcorrection function f(m) for gain amplification in two-wave mixing may probably be needed when photovoltaic field is sufficiently large for some E0. We further investigate the role of photovoltaic field for two commonly used experimental boundary conditions, open-circuit and close-circuit conditions. What is not expected is that for short-circuit condition, while the theory of Gu et al. predicts that increment of Ep will lead to larger Im(E1), actual Im(E1) becomes smaller at large m because the increment of |E1| becomes more slowly than the prediction of Gu et al. and phase difference between intensity and fundamental harmonics of space-charge field decreases rapidly. So far, all of our findings support that the effect of photovoltaic field cannot be cancelled by externally applied field, which contradicts with what is claimed to be done experimentally, cancellation of photovoltaic effect by an externally applied field, in the report of Cook et al. We have discussed the arguments of Cook et al. We also generalize the remark of Gu et al. on photovoltaic field and dc space-charge field, which is for small m only to whole range of m in the limit donor density much larger than ionized donor density that the effect of photovoltaic field cannot be cancelled or considered as an externally applied field for arbitrary non-zero m. Acknowledgements This work was supported, in part, by a Hong Kong RGC (Grant No. CUHK 4035/02P). We would like to thank Dr. X.S. Wang for part of his contribution to derivation of Eq. (5). We are also grateful to Prof. W.L. She for his helpful comments. We also thank M.W. Chan for allocation of computer resources for us.
m=0.7
4.4
References
4.2
m=0.9
0
5
10
15
20
25
30
Ep (kV/cm) Fig. 5. E0 = 0. (a) |E1|/m versus Ep. (b) u1 versus Ep. (c) Im(E1)/m versus Ep.
For non-zero arbitrary m, as J is a constant, J eslcND þ I jelc, a term which is related to the contribution of photovoltaic field, cannot always be zero or any constant. In this limit, the effect of photovoltaic field cannot be cancelled or considered as an externally applied field for an arbitrary non-zero m. 4. Conclusion We have studied photorefractive response to sinusoidal intensity pattern in photovoltaic media under different boundary conditions for modulation depth m, from close to 0 to 1 by solving
[1] N.V. Kukhtarev, V.B. Markov, S.G. Odulov, M.S. Soskin, V.L. Vinetskii, Ferroelectrics 22 (1979) 949. [2] G.C. Valley, J. Opt. Soc. Am. B 1 (1984) 868. [3] C. Gu, J. Hong, H. Li, D. Psaltis, P. Yeh, J. Appl. Phys. 69 (1991) 1167. [4] M.G. Moharam, T.K. Gaylord, R. Magnusson, J. Appl. Phys. 50 (1979) 5642. [5] E. Ochoa, F. Vachss, L. Hesselink, J. Opt. Soc. Am. A 3 (1986) 181. [6] F. Vachss, L. Hesselink, J. Opt. Soc. Am. A 5 (1988) 690. [7] Ph. Refregier, L. Solymar, H. Rajbenbach, J.P. Huignard, J. Appl. Phys. 58 (1985) 45. [8] F. Pauly, O. Sandfuchs, F. Kaiser, M.R. Belic, Opt. Commun. 218 (2003) 385. [9] L.B. Au, L. Solymar, Opt. Lett. 13 (1988) 660. [10] A. Bledowski, J. Otten, K.H. Ringhofer, Opt. Lett. 16 (1991) 672. [11] R.P. Rupp, R. Sommerfeldt, K.H. Ringhofer, E. Kratzig, Appl. Phys. B 51 (1990) 364. [12] E. Serrano, V. Lopez, M. Carrascosa, F. Agullo-Lopez, IEEE J. Quantum Electron. 30 (1994) 875. [13] G. Cook, J.P. Duignan, D.C. Jones, Opt. Commun. 192 (2001) 393. [14] A.A. Zozulya, G. Montemezzani, D.Z. Anderson, Phy. Rev. A 52 (1995) 4167. [15] M. Segev, G.C. Valley, M.C. Bashaw, M. Taya, M.M. Fejer, J. Opt. Soc. Am. B 14 (1997) 1772. [16] N.V. Kukhtarev, P. Buchhave, S.F. Lyuksyutov, Phy. Rev. A 55 (1997) 3133. [17] N.V. Kukhtarev, S.F. Lyuksyutov, P. Buchhave, T. Kukhtareva, K. Sayano, P.P. Banerjee, Phy. Rev. A 58 (1998) 4051.