Temporal analysis of optical beams in biased photorefractive materials in the context of solitonic solutions: microscopic and macroscopic approach

Fres-sengeas et al. [6] presented the explicit wave propagation equation including temporal behavior of spatial solitons. In all mentioned papers [2–6] the starting point was the K–V band transport model, where analogous simplifying assumptions were formulated to determine the refractive index variation. The solution given in [6] was regarded as able to describe screening and photovoltaic bright solitons Abstract The validity of the commonly used time-depend-ent wave equation describing the propagation of screening one-dimensional solitons in photorefractive materials is discussed. Concentrating attention on temporal development of the space-charge field, we show that the widely used standard solution follows from a phenomenological description, which is consistent with the band transport model equations only in specific cases. The exact analytical solution for the localized optical beam is derived within the microscopic model under a low contrast approximation. The numerical modeling of photorefractive response to an arbitrary contrast is performed and compared with standard solutions. The range of applicability of the macroscopic approach for three basic classes of photorefractive crystals is discussed.


Introduction
It is generally accepted that the band transport model of Kukhtarev-Vinetskii (K-V) [1] gives the correct theoretical description of the photorefractive (PR) phenomena. The rate equations which make up the K-V model take into account the basic physical processes as photoexcitation, transport and trapping of free carriers. The K-V model has been employed primarily to describe PR response for two kinds of light intensity distributions: the sinusoidal light pattern and the single localized beam. Research connected in generic PR media. To obtain a searched for solution, the time-dependent relation between the space-charge field and the intensity was first established. Inclusion of this relationship in the paraxial wave equation allowed to carry out in a relatively simple manner the time-dependent analysis of behavior of an optical beam converging into a soliton state. This standard solution is rather controversial, but offering the neat and simple formula, the presented approach was widely accepted by other researchers and up to now has been the basis of one-dimensional PR solitons analysis. Although the original solution had been developed to describe bright solitons, later it was generally agreed that it could also be applicable to analyze dark solitons. Approximations, which underlie the considered solution, make up the basis of numerous theoretical analyses of solitonic beams behavior including: analytical and numerical modeling of evolution of bright and dark screening solitons [7][8][9][10][11][12][13][14], temporal formation of photovoltaic solitons [15,16], stability of dark screening solitons [17,18], analytical approximation for the soliton profile [19,20], concept of photovoltaic-screening solitons and dissipative solitons [21,22]. Also, we can find these assumptions in chapters in books and general reviews [23][24][25][26].
The purpose of the present work is to indicate that the standard time-dependent wave propagation equation for bright and dark optical beams is highly problematic in some cases. Particular attention is paid to the correctness of describing the temporal development of the space-charge field in PR materials. It is shown that the solution obtained from the K-V equations under two widely applied approximations leads in fact to the macroscopic solution with limited range of applicability. It is notable that for the PR grating analysis an approach based on the microscopic model is used almost without exception, while for localized beams the phenomenological model is exploited. Two different models used for two various light patterns provide an incoherent description of the photorefractive material response.
Taking into account three main families of PR materials, it turns out that the phenomenological approach yields a correct approximation only for bright beams in ferroelectric crystals. It can also be applied conditionally for bright beams in sillenites and semiconductors. On the other hand, investigating dark beams we find significant discrepancies in comparison to the results from the microscopic model equations. In this case, the macroscopic solution fails for all PR crystals herein considered.
The paper is organized as follows. In the first part, after introducing the K-V equations and formulas for bright and dark beams, we repeat the basic steps leading to the standard solution of the temporal space-charge field formation, paying special attention to assumptions made to simplify the analysis. The characteristic response time for optical beam is compared with the well-known solution for the PR sinusoidal grating, and the exact analytical solution for localized optical beam under the low contrast approximation is presented. The numerical solutions of time evolution of the space-charge field are shown for bright and dark beams at arbitrary contrast value. Finally, in the light of experimental results given in the literature we try to explain why, up to now, the inconsistency between theory and data has not been reported. The standard spatial-temporal wave equation describing the propagation of screening onedimensional solitons in PR materials is discussed and final conclusions are formulated.

Kukhtarev-Vinetskii model: bright and dark beams
For many PR materials, we can use the simple version of the K-V model with the single trap level and one type of charge carrier. Such model constitutes a compromise between the complexity of the PR effect and the simplicity of mathematical description. In this case one assumes two interband levels of impurities: one of photoactive donors with concentration N D , and one of shallow acceptors with concentration N A . In absence of light it is assumed that all acceptor sites are ionized, partially compensating donors, that is N D . Under the influence of illumination with the light distribution I(x) electrons are excited from nonionized donors N D 0 = N D − N D + to the conduction band, where carriers are transported due to diffusion and drift and eventually are trapped mainly in darker regions. As a result, the space-charge distribution of donors and electrons is created in a crystal which produces the nonuniform space-charge field E sc (x). For one-dimensional beams and neglected photovoltaic effect, these processes are represented by the set of material equations [6,27,28]: In above equations symbols ∂ t and ∇ stand for ∂ ∂t and ∂ ∂x, respectively. I s is the power intensity of the optical beam (the signal beam), I b = I d + I B , where I d is the so-called equivalent dark irradiance corresponding to the rate of electron thermal generation, I B is the background uniform illumination used to artificially increase I d , and S = s/hν is the photoionization cross-section s per photon energy hν. N D , N D of donors, ionized donors, acceptors and free electrons, respectively; γ is the recombination coefficient, μ-electron mobility, q-elementary charge, E is the total electric field inside a crystal, U T = k B T/q is the thermal potential (U T ~ 26 mV at room temperature), where k B is the Boltzmann's constant, and T-the absolute temperature; j is the current density; ε = ε 0 ε r , where ε 0 and ε r are the vacuum and relative low frequency dielectric constants, respectively. If a constant voltage V is applied to the crystal of length L, Eqs. (1a)-(1e) are supplemented by the bias condition The total electric field E(x) in a crystal is the sum of the external electric field E a and the spatially modulated spacecharge field E sc (x): E(x) = E a + E sc (x). The change in the refractive index n m induced by the linear electro-optic effect is proportional to the electric field: where r eff denotes the effective electro-optic coefficient. Further, we consider the one-dimensional light beam with the intensity distribution I(x) in the x-direction. In a general case of soliton beams investigations a total power intensity of light falling into a crystal is the sum of the signal beam intensity and the dark intensity I d together with the background intensity I B . In this work, the signal beam is assumed to be a Gaussian beam described by the distribution of light intensity in the form: where w is the half-width at a half maximum (HWHM).
For the bright beam, the signal intensity vanishes at infinity, thus I s (x → ∞) = I ∞ = 0, and total intensity is: The dark beam is a dark dip in the background of constant intensity. For this case, we have I s (x → ∞) = I ∞ = I m and the total intensity can be written as: In order to uniformly describe bright and dark beams, it is convenient to introduce a parameter m, herein called the contrast beam parameter or alternatively, brightness/darkness coefficient. The parameter m is given by m = I m /I b (m > 0) for the bright beam and m = I ∞ /(I b + I ∞ ), (0 < m ≤ 1) for the dark beam. Note that in the latter case, the coefficient m does not correspond to the grayness parameter defined in [14] as I(0)/I ∞ . In what follows, most calculations are carried out to simulate the PR behavior of Bi 12 SiO 20 (BSO)a non-photovoltaic sillenite crystal with fast PR effect and for strontium-barium niobate (SBN)-ferroelectric crystal often used in PR soliton experiments. Material parameters of BSO are quite well known, and for sillenites it is common to apply a high external voltage which allows investigating the PR response in relatively wide range of applied fields. On the other hand, SBN has a very large dielectric constant; hence, screening solitons can be formed in a rather small field of the order of 1 kV/cm. The considered single defect level model agrees quite well with experiments in BSO and SBN crystals. The values of relevant material parameters for BSO and SBN are gathered in Table 1.

The space-charge field equation
For typical PR materials the inequality N D ≫ N + D , N A is fulfilled and if the light intensity is not too high the concentration of photoexcited electrons obeys the relation n ≪ N + D , N 0 D , N A . As a result, N + D can be neglected with respect to N D in Eq. (1a) and electron density can be omitted in comparison to (N + D − N A ) in Eq. (1e). To make equations more legible, we introduce the characteristic length X N = εE a /(qN A ) and two normalized quantities: the normalized electric field and the function u(x) of relative light intensity, given, respectively, as: Moreover, for the description of the temporal behavior of PR material two characteristic times are defined: which stand for, respectively: the carrier recombination time and the dielectric relaxation time corresponding to the background illumination. In Eq. (5d) n b denotes the free-electron density n 0 in regions with I total (x → ±∞) = I ∞ + I b .
First, we will derive the equation for the temporal spacecharge evolution without applying any additional approximations. On the basis of Eqs. (1d) and (1e) under the above assumptions regarding impurities and carrier densities we obtain relations: Substituting Eqs. (6a) into (1a) leads to the expression for the electron density: By introducing Eqs. (1b) and (6b) into the continuity Eq. (1c) we arrive at the following equation: Integration of Eq. (8) over the spatial variable x leads to the Ampere's law: The total current density J total can be determined by exploiting the appropriate boundary conditions for the light intensity together with the condition given by Eq. (1f). The spatial x-width of the typically focused optical beam (~10 µm) in slab-soliton experiments is much less than the width L of the crystal (a few mm). Hence, at a sufficient distance from the beam center (formally at x → ±∞) the electric field E ∞ ≈ E a = V/L, that is e ∞ ≈ 1, and I ∞ → constant. Together with the background illumination I b , the resulting light intensity is thus given by I total (x → ±∞) = I ∞ + I b .
Taking these assumptions into account and averaging spatially both sides of Eq. (9) over the crystal length L we obtain: �ne� ≈ n b e ∞ ≈ n b , �∂ t e� = 0, �∇n� = 0. Hence, J total = qµn b and including (5b) we finally get: where the ratio n/n b according to Eqs. (7a) and (7b) is equal: Equation (10) involving Eq. (11) yields the looked for time-dependent equation for the total electric field E(x, t), hence the space-charge field E sc (x, t). Equation (10) does not have in general an analytical solution and needs to be solved numerically. In the stationary state on the basis (10) and (11), we get equation [23][24][25][26]32]: Experiments with screening solitons require the application of external field of the order of kV/cm, so that the second term in Eq. (12) associated with diffusion transport, for typical light beams with FWHM ~10 µm can be usually neglected. As a result, Eq. (12) reduces to: .

Approximations leading to the standard space-charge field expressions
In works [4][5][6] two strong assumptions were adopted that allow to significant simplification of Eq. (10) and obtain the analytical solution for the internal electric field E(x, t). These two assumptions: slowly varying of the space-charge profile and so-called adiabatic approximation will be discussed in more detail below.

Slowly varying of the space-charge field distribution
For typical optical beams with transverse size of the order of 10 µm, the power intensity distribution I(x) can be regarded as a slowly varying function; hence, the hypothesis was derived [2][3][4][5][6] that the induced electric field profile in a crystal should be also described by a slowly varying function of x. This suggests that inequality is satisfied and the term X N ∇e can be disregarded in Eqs. (11)- (13). In that case, Eq. (11) simplifies to the form: and the electric field distribution from Eq. (13) can be calculated according to the relation: The assumption given by Eq. (14) and resulting Eq. (16) is commonly applied in the theoretical analysis concerning screening solitons . However, it should be stressed that this approximation has a severely limited range of (14) applicability as regards dark beams. To prove this, note that Eq. (6a) reads: As seen, the considered approximation (14) is equivalent to the assumption that depletion of ionized donor traps is always negligible over all regions of the light intensity distribution-in other words N D + ≈ N A . Now, let us take into account the black beam described by an arbitrary function u(x), where the term "black" indicates that its minimum intensity is zero at any point x p . According to Eq. (13) we obtain: From (18) it can be seen immediately that for a beam with u(x p ) = 0 (for the symmetric beam x p = 0) one gets N D + (x p ) = 0. Thus, in the center of a black beam we always have 100 % depletion of donor traps regardless of the shape of the function describing the profile u(x). When background illumination is added, the depletion will be smaller than 100 %, although the above analysis reveals that for dark beams the analysis should be carried out on the basis of Eq. (13). Only for bright beams the approximation (16) is correct. The use of the discussed approximation to determine the stationary space-charge field distribution in the frame of K-V model can lead to large errors regarding dark beams. To illustrate this conclusion, Fig. 1a, b present profiles of the internal electric fields generated by dark and bright optical Gaussian beams with HWHM = 5 µm. Material parameters for cerium-doped strontium-bariumniobate (SBN) crystal (Table 1) has been adopted under the assumption of an applied field E a = 1 kV/cm. The difference between shapes of electric field profiles obtained from Eqs. (13) and (16) is apparent. Approximated conformity is obtained only for small values of the darkness coefficient, whereas for values of m close to unity we find the stark discrepancy between the considered solutions.
We can set a different argument, that for dark beams the assumption (14) cannot be accepted as correct. Let us consider the function u(x) given by Eq. (5b) defined for Gaussian bright and dark beams with HWHM = w described, respectively, by Eqs. (4a) and (4b) and, according to (16) let us determine HWHM = W E of the electric field distribution e(x). The normalized half-width is equal to: where the signs +, − refer to, respectively, bright beam (m > 0) and dark beam (0 < m ≤ 1). The graph of the function W E (m)/w is shown in Fig. 2.
For the bright beam, the width HWHM of the field distribution E(x) is greater than the HWHM of the optical beam and increases along with the coefficient m; hence, the slowly varying assumption of the function E(x) is justified. However, for the dark beam an increase in contrast coefficient leads in the limit m → 1, to the formation of field distribution E(x) with the width HWHM → 0. In such case, the donor density N D + calculated from the Gauss law (1d) takes high negative values. This result as a non-physical should be rejected.

"Adiabatic" approximation
In the analysis of photorefractive process dynamics we can distinguish two fundamental time scales determined by the carrier recombination time (τ r ) and the dielectric relaxation time (τ die )-see Eqs. 5c and 5d. For PR materials in typical CW-experiments the inequality holds: τ r ≪ τ die . Referring to this property, in many works [6, 9-11, 15, 34] authors assume that time evolution of carriers can be completely neglected, thus putting ∂N D + /∂t = 0 in Eq. (1a) is justified. Following [23] we call this substitution the "adiabatic" approximation, although as will be shown later it is an incorrect name. Insertion of this assumption into Eq. (11) in connection with the condition (14) yields Considering a case when the light is abruptly switched on at time t = 0, i.e., u(x, t) = u(x)θ(t), where θ(t) is the Heaviside's step function, Eq. (20) means that the carrier distribution attains the steady state immediately and during the space-charge field formation the carrier distribution does not change in time. In this event Eq. (10) takes the form: This is a simple differential equation which has the following solution [6,9,16,34,35] According to Eq. (22), the space-charge distribution varies monotonously to the stationary state without any oscillations with the response time: For bright beams, the response time decreases with the increasing brightness parameter, contrary to dark beams for which the response time grows with the grayness parameter. Thus, Eq. (23) yields the prediction that for bright beams the time formation of the space-charge field is below the dielectric relaxation time, whereas for dark beams the time τ sc can considerably exceed the time τ die . In the second case, the problem appears in the point x p in which u(x p ) goes to zero. Then, according to Eq. (23) τ sc (x p ) → ∞ and from Eq. (22), we obtain in the limit t → ∞ the nonphysical singular solution e(x p ,∞) → ∞. Note, that the solution (22) was originally formulated for an analysis of bright PR solitons. However, later this formalism was adapted by some researches to study dark solitons [12-15, 17, 18]. In the following, the expressions (16) and (22) will be referred to as the standard solutions. It is notable that none of the above considered assumptions is used in the theoretical analysis of photorefractive grating formation. Thus, two different ways of describing PR response are applied. To show main differences in these two approaches, we will recall briefly the derivation of the standard temporal grating solution.

Linearized solutions for harmonic grating
Consider a sinusoidal light pattern which is switched on at time t = 0: where I 0 is the average intensity, θ(t) is the temporal Heaviside's step function, K = 2π � is a grating constant for the spatial period Λ, m is the modulation index (fringe contrast). The solutions of the K-V equations in the linear regime assuming the light-driven function in the form (24) is well known. For a small modulation depth m the photorefractive rate Eqs. (1a)-(1f) are solved by applying the perturbation approach in the parameter m to the first order. The zero-order solutions, i.e., under homogenous illumination, accordingly for the ionized donor and carrier density, are where the steady state amplitude of photoexcited electrons is n 0 (∞) = SI 0 N D τ r . The first-order equations can be written in the matrix form: where n 1 , N 1 + designate amplitudes of electron and ion harmonic distributions, whereas elements of the square matrix represent characteristic transition rates (s −1 ): Γ Re = τ r −1 = γN A -the carrier recombination rate, Γ Ie = γn 0 (∞)-the ion recombination rate, Γ die = τ −1 die = qμn 0 (∞)/ε-the dielectric relaxation rate, Γ tote = Γ Re + µKE d − i·µKE a and Γ Ph = mSI 0 N D -the photogeneration rate. The set of Eq. (26) can be easily solved by Laplace transformation. Having solutions of (26), we can calculate the space-charge field amplitude from Gauss's law: E 1 ≈ (q/iεK)N 1 + . Eigenvalues s 1 , s 2 of the transition rates matrix determine time constants describing the temporal behavior by the relation τ j = −s j −1 , j = 1, 2. Taking into account, the amplitude of electron density under the initial condition n 1 (0) = 0, we achieve the complete solution: For PR materials, the inequality Re(1/τ 1 ) ≫ Re(1/τ 2 ) is generally valid; therefore, the second term in Eq. (28) can be neglected in comparison with the third one, so that the equation can be cast in the simplified form: Omission in the solution of the root s 1 is equivalent to putting dn 1 /dt = 0 in Eq. (26). This assumption, called the adiabatic approximation, allows to transform Eq. (26) into a single differential equation: where the steady state fundamental amplitude of the spacecharge field E 1 (∞) and the grating formation time τ N = τ 2 are given by expressions: The transient solution for the grating amplitude of the space-charge field describes the exponential function: The formulas (31a) and (31b) are written in terms of characteristic electric fields defined as follows: The real part of Eq. (32) yields the full solution of the space-charge field distribution [28]: Equation (34) represents a superposition of the stationary harmonic distribution and the exponentially damped distribution which moves at a speed ω sc /K. Unlike the local solution (22), Eq. (34) shows that the evolution of the distribution of E sc is dynamically non-local in respect to the intensity pattern during the grating formation process. The expression (34) is fully in accordance with the rigorous numerical solution. 6 Macroscopic approach and its applicability According to Eq. (31b) the grating build-up time constant is a function of the applied field, the gradient of light intensity distribution (through the grating constant K) and microscopic parameters as the recombination time τ r and the acceptor density N A . In contrast, the characteristic time given by Eq. (23) does not exhibit any of these dependences. This feature can be explained by noting that the formula (22) comes in fact from the macroscopic (phenomenological) approach. In other words, at any point we do not need to refer to the microscopic Eqs. (1a) and (1b) of the full K-V model. In the frame of the macroscopic scheme, the key assumption is setting that photoconductivity is locally proportional to the light intensity pattern , where C is the parameter of proportionality. The above equality can be expressed as is the photoconductivity sufficiently far from the beam.
Taking into account macroscopic Eqs. (1c) and (1d), we obtain the equation: Writing the dependence for the current density in the form: j = σ E + U T ∇σ and taking advantage of Eq. (36) after integration over the spatial variable x we arrive at the equation: By defining the dielectric relaxation time as τ die = ε/σ b Eq. (37) takes the form of Eq. (21). Note, that in the above scheme only two macroscopic parameters are exploited: the dielectric constant ε and the diffusion voltage U T . If we apply the analogous macroscopic approach to find the linear solution for a PR grating one obtains the equation which is the counterpart to Eqs. (21) and (30), and has the solution: Comparison of (39)-(30) together with (31) enables us to establish conditions for the correctness of the phenomenological description [36]. As regards, the amplitude E 1 (∞), Eq. (39) gives the same result as Eq. (30), if inequalities E d E q ≪ 1 and E a E q ≪ 1 are fulfilled. Considering the characteristic response time τ N , the phenomenological approach can be applied if the conditions E d E q ≪ 1 and E a E q ≪ 1 hold true, as well as E d E µe ≪ 1 and E a E µe ≪ 1. In Table 1 are listed the material parameter values for crystals which are representative of the three most common classes of PR material family: strontiumbarium niobate (SBN) for ferroelectrics, Bi 12 SiO 20 (BSO) for cubic crystals (sillenites) and GaAs:Cr for compound semiconductors. As follows from the calculated ratios of characteristic electric fields given in Table 1, we can draw the conclusion that within the small signal approximation the macroscopic approach may be successfully applied only to ferroelectric crystals.
It should be stressed, that limitations of validity of the macroscopic model presented here as regards the linear harmonic grating solution have a more general character. In fact, considering a large contrast of the light distribution, disagreement between the microscopic and macroscopic model is even higher particularly in relation to the description of PR response dynamics.

Adiabatic approximation in physics of photorefractive effect
The assumption that carrier distribution does not undergo time evolution by setting ∂ t N D + = 0 in Eq. (1a) leads to radical simplification of the K-V model equations. Thus, the free carrier density is determined from the relation (20). This assumption needs to be considered in more detail. From the standpoint of mathematics, this statement is of course not consistent with the continuity Eq. (1c). In this event, Eq. (6a) implies ∂ t E = 0. Note, that such approximation should not be called "adiabatic." Generally, in physics the adiabatic process conveys two meanings. The term adiabatic process mostly occurs in the context of thermodynamics. The other meaning of such process refers to sufficiently slow changes of the external conditions in a dynamic system [37,38]. In this case there are two characteristic times: the "external" time (T ext ) corresponding to rate of external parameter changing and the "internal" time (T int ) describing the speed of the process considered. An adiabatic process is one for which T ext ≫ T int . In photorefractive phenomena, photogeneration and recombination proceed in a much shorter time scale than the process of building spatial charge in traps; hence, this process can be treated as adiabatic. As a consequence, after appearing of an external light pattern, the photoexcited carrier distribution attains in the recombination time scale the quasi-equilibrium state with respect to the light intensity and space-charge in traps and follows temporal changes in these distributions with a very small delay. The adiabatic approximation in the analysis of the PR effect relies on completely neglecting that delay. This assumption breaks down for short light pulses comparable or shorter than the carrier recombination time.
It should be stressed that the adiabatic approximation does not mean that the distribution of free carriers instantly reaches a steady state but evolves from an initial state n(x, t = 0) to a steady state, remaining in equilibrium with the ionized donor distribution. It is a direct consequence of the fact that the transport of carriers causes a spatial charge separation. To illustrate this, let us recall the analysis of elementary hologram record presented in Sect. 5, assuming BSO parameters values as per Table 1. According to adiabatic approximation, instead of the complete solution (28), we apply Eq. (29) written in the form: where E 1 (t) is given by Eq. (22). Figure 3 compares the solutions obtained from the complete Eq. (28) and simplified Eq. (40). Figure 3 presents the time evolution of the free carrier density amplitude according to Eq. (28), involving two time constants. It is notable that contrary to the prediction of the expression (22) for the localized beam, the time in which the steady state is attained is of the order of hundreds of dielectric times. The use of adiabatic approximation neglects a short period of temporal evolution immediately after the light is switched on, which formally changes the initial condition from n 1 (t = 0) = 0 to n 1 (t = 0) ≠ 0, so that the plot n 1 (t) in Fig. 4 begins from the initial value Note that the value n 1 (0) may be much lower than the steady state value n 1 (∞) (in Fig. 4 the ratio n 1 (∞)/n 1 (0) ≈ 25 times). From that moment on the solutions of (28) and (40) become practically indistinguishable. Considering the intensity distribution of light in the form of localized optical beam, adopting of the adiabatic approximation in the K-V model implies that the carrier density distribution is determined by relation (7a). As a result, the derivative dn/dt does not exist in the continuity Eq. (1c), which is equivalent to the assumption that dn/ dt = 0. For this reason (7a) is an equivalent of Eq. (40) for the grating. Following from Eq. (7a) one can see that if X N ∇e ≪ 1, as for bright beams, then n(x, t) evolves in time in the same way as the field E sc time derivative. For dark beams with high grayness coefficient, for which we cannot neglect the term X N ∇e in the denominator, the change in time n(x, t) takes place differently from that for the bright beam. The related conclusion is that the assumption in (1a): ∂ t N D + = 0 leading to the solution (22) does not fit the adiabatic approximation definition, and should not be termed as such. This approximation is for macroscopic description. The solution of equation system (10) and (11) in the general case can only be achieved numerically. However, we can obtain the exact linear solution, for small values of beam contrast coefficient. Optical beams used for creating soliton states are usually characterized by large contrast, i.e., for dark beams m → 1, for bright beams m ≫ 1. However, the linear solution is useful for at least two reasons. Firstly, like for the harmonic grating, it allows comparing the exact solution with the approximated macroscopic solution. Secondly, it estimates the time of space-charge formation time in the material, which is particularly useful for dark beams, allowing determining the minimum response time.

Linearized solution for the localized optical beam
The linear solution of the K-V model describing the PR response dynamics for localized optical beam can be obtained as a superposition of solutions from the theory of harmonic grating given in Sect. 5. We commence by multiplying the function I(x) describing the optical beam power distribution, creating a periodic function in the form: where p = 0, ±1, ±2, … and X denotes the space period of the function Ĩ (x), where X should be significantly greater than the full beam width (X ≫ FWHM). This condition follows from the fact that during the temporal evolution Fig. 3 Exact solution (28) starts from the initial condition n 1 (0) = 0, whereas the solution (40) in the adiabatic approximation-from the value indicated by the dashed line; from this moment both plots become practically indistinguishable. Material parameters assumed for BSO changes in space-charge distribution in the vicinity of the function I(x + pX) should not affect charge formation at the next function I(x + pX ± X). The periodic function Ĩ(x) can be expanded in a Fourier series with a period of X: It is assumed that the beam profile I(x) is symmetrical relative to point x = 0; therefore, we use the expansion only according to the cosine function. For the sake of calculations, the truncated Fourier series is used, where the minimum number of components should be chosen so as to ensure convergence of the series.
In the linear approximation, each harmonic light pattern with a modulation index m n = m·Ĩ n leads to the formation of independent grating of the space-charge field E sc (x, t) are given by Eq. (34), where the substitution K → K n = n(2π/X) is made. Taking the superposition of harmonic components described by Eqs. (31) and (32), we obtain the PR material response to optical beam illumination as the expression for time-dependent complex amplitude of the field E sc (t): where τ (n) N is estimated from Eq. (31b). Calculating the module of the amplitude (44): we find out temporal changes of a minimum (bright beam) or maximum (dark beam) value of the field E sc in the spatial distribution. In the linear approximation the temporal evolution runs in the same way as in the case of illumination by bright or dark beam. The complete solution, including the variable x, showing how the space-charge field distribution changes in time: where τ (n) sc and ω (n) are calculated according to the formulas (33a), (33b) with the replacement of K by nK = n·2π/X. Writing the relations (33) in the explicit form we get: where R (n) ab denote ratios of the respective characteristic fields defined in Sect. 5: Let us underline that the solution (46) is fully consistent with the numerical solution of the K-V model equations. Fig. 4 presents a temporal change in the field amplitude while a steady state is being reached, calculated by Eq. (45), and a profile of space-charge field distribution at three selected instants according to Eq. (46). Parameters of two different PR materials were taken for calculations (Table 1): SBN and BSO. For comparison, Fig. 4a, c also include the approximate solution from Eq. (22). For the BSO transient state changes occur in the field E sc (x, t) profile, along with a space non-locality and oscillations. For SBN, the field changes monotonically without any oscillations and locally, relative to the light intensity distribution. In this case, the macroscopic approach gives a very good approximation of the microscopic solution. If in Eq. (45) the terms cos(ω (n) t) describing oscillations is omitted, a monotonic change of field amplitude is obtained. For thus defined function, we can estimate the characteristic time τ avg of reaching the steady state by defining this time as a weighted average of harmonic gratings response times: Thus we can write: The relation (49) is represented in Fig. 5a as a dotted line. The time τ avg is for sillenites and semiconductors as GaAs much longer than the dielectric relaxation time. From the practical viewpoint, it is useful to determine a different characteristic time τ1 avg , which defines the initial rate of (47a) τ (n) . an electric field rise to the value corresponding to a steady state. Postulating the exponential relation: to determine the time τ1 avg we adopted the criterion of mean square root error minimization between the function (50) and relation (45), i.e., [|E beam sc (t)| − F sc (t)] 2 time = min. Seeking the minimum relative to the time τ1 avg we calculate: After rejecting the small terms one finds that the looked for time is given approximately by the relation: It should be noted that the times τ avg and τ1 avg are inversely proportional to the mean light intensity, similarly to the relation (23). The function (50) is plotted in Fig. 4a, using the formula (52). In the considered case for BSO we get τ1 avg ≈ 4τ die .
For bright and dark beams with low contrast, the macroscopic approach generally is not correct for sillenites.

Numerical solutions for bright and dark beams
To model the temporal evolution of the space-charge field under the arbitrary beam contrast coefficient, we performed the numerical calculations according to the full K-V equations. The numerical algorithm was based on the implicit, finite difference method which follows partially the scheme proposed in [39]. The outcome was next compared to macroscopic solutions. Figure 5a, b show the time-dependent development of the internal electric field for a dark beam with a high darkness coefficient, adopting material parameter values for BSO and SBN crystals. The corresponding plots are displayed in Fig. 5c, d, but these reflect the case of a bright beam.  Fig. 5a (BSO) and Fig. 5c (SBN)-on the basis of Eq. (46). Transition states of the internal field profile for three different times- Fig. 5b, d. The external fields E a = 10 and 3 kV/cm were adopted, respectively, for BSO and SBN crystal. The phenomenological solution given by Eq. (22) and the microscopic solution obtained from (45) are compared. Time is normalized to the dielectric relaxation time Figure 5a, b reveal very large discrepancies between microscopic and macroscopic solutions in both sillenite and ferroelectric material. Two features are worth noting. First, for dark beams the electric field reaches the steady state without oscillations. For BSO the microscopic solution predicts much greater time to reach the steady state than it stems from the macroscopic model. The opposite behavior takes place in the case of SBN, where the divergence between microscopic and macroscopic solutions is also evident. Moreover, in both cases the phenomenological description predicts much greater magnitude of the internal electric field amplitude. Figure 5c, d show, for the BSO crystal, the time-dependent changes of electric field amplitude in the center of the light distribution with different values of the contrast coefficient. It is conspicuous that the photorefractive response for a bright beam is clearly different from the response for a dark beam. The electric field evolves toward the stationary state in an oscillatory way, however for higher values of the coefficient m oscillations disappear and the duration of the transition state is shortened. Furthermore, the time needed to attain the maximal screening of the external field, denoted in Fig. 5c, d by T m , tends to the value given by formula (23). Concluding, at a sufficiently high value of m, the analytical expression (22) can be regarded as a plausible approximation as far as bright beams are concerned.
In reference to the K-V model equations, the observed in Fig. 5 large differences in the macroscopic and numerical description of material response dynamics result from the fact that in Eq. (7a) and the second term in the numerator may become comparable to the first term; hence, it cannot be omitted. We can derive the approximate criterion of applicability of the macroscopic approximation for bright beams. To estimate the order of magnitude of ∂(∇e) ∂t in Eq. 7a, it can be assumed that the field in a crystal with spatial variation of FWHM magnitude reaches a value E max ~ E a in the time of τ die order. With such assumptions, the second term in the numerator in Eq. (7a) can be written as follows: Taking into account Eqs. (5d) and (7b) the dielectric relaxation time may be expressed as: Putting (54) into (53) and comparing with the first term of the numerator in (7b), and assuming it is much greater than the second term, we obtain: where the function u(x) is given by the relation (5b). This is the looked for criterion of macroscopic model applicability for the bright beam. For instance, assuming for the BSO crystal the field E a = 10 kV/cm, FWHM = 10 µm, and as per Table 1 the product of mobility and recombination time µτ r = 4 × 10 −11 (m 2 /V) we find that the expression on the right-hand side of Eq. (55) is of the order of unity. Therefore, for the bright beam, if the optical power of the beam is much higher than the background intensity, we have u(x) ≫ 1 and the phenomenological model yields a good approximation. In case of ferroelectric materials, a small value of the product µτ r makes Eq. (22) remain valid, even for a small coefficient of beam contrast.

The soliton wave propagation equation and experiments
In theoretical works dealing with an analysis of PR screening solitons 1D, the description of beam propagation is based on the paraxial envelope wave equation in the form [2,3,[6][7][8][9][10][11][12][13][14][15][16][17][18]: where φ = φ(x, z, t) designates the normalized slowly varying complex amplitude of the beam intensity I s = I b ·|φ| 2 which is allowed to diffract along the x-direction and propagates along the z-axis. The plus sign corresponds to a selffocusing nonlinearity and bright beam trapping, while the minus sign corresponds to a self-defocusing nonlinearity and possibility to obtain a dark soliton. In linear electrooptic crystals the variation of Δn m induced by propagating optical beam is proportional to the electric field E(x, z, t) according to Eq. (2). In the standard approach this internal field is calculated on the basis of Eq. (22). Taking into account that I s = I b ·|φ| 2 , neglecting diffusion due to strongly biased PR crystals, and inserting expression (22) into Eq. (56), the general form of the dynamic wave equation is derived [6, 8-12, 15, 34]: where A = (1/2)k 0 n 3 m r eff E a and u(φ) = (1 + |φ| 2 )/ (1 + I ∞ /I b ). Equation (57) is taken in the literature to analyze bright solitons and dark solitons and constitutes the basis of many numerical investigations [6, 8-11, 12, 15, 34]. However, in the context of herein presented arguments employing Eq. (57) we should bear in mind the limited range of applicability of the phenomenological approximation as regards the determination of the nonlinear term Δn m (E). As shown, the macroscopic approach yields a correct quantitative description in the case of bright beams with high contrast coefficients. However, to study dark solitons  [11,16,[40][41][42][43][44][45][46]. As it was indicated, for dark solitons, the macroscopic approach gives the space-charge field solution to a high degree inconsistent with the K-V model. However, one can show that solutions for dark beams obtained from the wave Eq. (57) (in a stationary state) reveal an property which can be called the "quasi-degeneration" of solitonic solutions. Although the nonlinear term Δn m (E) in Eq. (56) determined from the macroscopic model is considerably different from the microscopic solution, it turns out that intensity profiles of dark solitons exhibit in both cases a high similarity. In particular, half-widths of profiles are nearly identical [33].
Here, two points are important to note. The complete beam shape is not determined in the experiments, but the most relevant measured parameter is the intensity FWHM of a beam, which serves to plot the soliton existence curve-a basic tool for comparing experimental results with predictions of theoretical models. Furthermore, the characteristic feature of photorefractive experiments is the presence of relatively large temporal fluctuations in the output optical signal [47]. For these reasons, small differences in soliton profiles are difficult to notice. Taking this into account, the more reliable test of correctness of the theoretical model should be the investigation of the temporal development of PR screening solitons. Such experiments carried out for bright solitons with a large contrast coefficient confirmed the validity of the approximated description based on Eq. (22). As it was remarked above, the disagreements between experiment and theory can be expected only for bright beams with a low value of the brightness parameter.
According to conclusions in Sect. 6, the conclusive experiment verifying the correctness of Eq. (57) is the examination of the transition state for dark beams by measuring the time of reaching the steady solitonic state. According to our knowledge, results of such experiment have not been published to date. In experiments with dark photovoltaic and screening solitons, temporal changes of a beam width at the crystal output were studied, but the direct measurements of soliton formation time have not been reported. In the light of the above considerations, the result of such experiment should be negative for the phenomenological model.
In closing, the most important conclusion is that the case of dark screening solitons should not be analyzed in an analogous manner as bright solitons, particularly in relation to the dynamical behavior analysis. The dark beam is a more complex issue than the bright one, and the analysis of the dark solitons propagation in PR media should be carried out by combining BPM methods with a numerical technique used for determining the light-induced space-charge field [48].

Summary
It is generally adopted that the theory for photorefractive (1 + 1)D screening solitons, known for over two decades, is well established and satisfactorily confirmed by numerous experiments. In the present article, we dispute partially the above statement taking a critical look at commonly used simplifying assumptions underlying the description of the nonlinear term in the wave propagation equation.
Using the K-V model, we took a closer look at the dynamics of PR material response to the illumination of the crystal by spatially localized optical beams, bright and dark. Two main assumptions adopted in theoretical works: slowly varying of beam profile and "adiabatic" approximation were discussed. These assumptions lead to an approximate solution describing the dynamics of space-charge field formation. To assess the applicability of these approximations we have presented an exact analytical time-dependent solution for the optical beam in the linearized approximation and a numerical solution for high contrast beams. In the context of these considerations, we may formulate the following conclusions.
(1) The commonly used analytical solution derived from the K-V model for the determination of refractive index change distribution is a solution based on the macroscopic approach that has a limited applicability. It constitutes a satisfactory approximation for the description of bright soliton formation in sillenites and semiconductors, and very good approximation for ferroelectrics. (2) The macroscopic solution reveals high disagreement compared with the microscopic solution for dark beams. Therefore, examining dark solitons in each of the three types of PR crystals herein considered the microscopic model should be used, especially in the transition state analysis. (3) The standard form of envelope wave equation is formally incorrect for dark solitons because the slowly varying assumption of the optical beam profile does not imply a similar property for the spacecharge field profile. However, the use of such form of equation does not lead to major errors in an analysis of soliton beams in a stationary state. (4) Numerical simulations of the time-dependent state of dark soliton beams should not be based on Eq. (57) which employs the "adiabatic" approximation. (5) As the conclusive test verifying the validity of the macroscopic model in dark beam analysis it is suggested the direct measurement of the black soliton formation time.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.