Photorefractive（PR）spatial solitons have been extensively investigated in the past two decades，in light of their unique features of formation on the order of a few $\mathrm{m}\mathrm{W}$ and their important potential applications^［1-4］. Various types of PR spatial solitons that arise from the change in refractive index due to only the linear electro-optic effect（Pockels effect）in noncentrosymmetric PR crystals or the quadratic electro-optic effect（DC Kerr effect）in centrosymmetric PR crystals have been investigated both theoretically and experimentally^［5-8］. Moreover，under proper conditions，PR spatial solitons and soliton pairs（or soliton families）governed by both the linear and quadratic EO effects have been proven to exist in many noncentrosymmetric PR crystals when the crystal temperature is close to its phase-transition temperature^［9-12］.

Modulational instability（MI）refers to the interplay between nonlinearity and diffraction or dispersion in the spatial or temporal domain，which occur in most nonlinear optical wave systems^［13-16］. For a plane wave or broad optical beam that propagates in a nonlinear optical medium，spatial MI causes the incident beam to disintegrate during propagation，which in turn leads to the formation of multiple wave filaments. The filaments resulting from the MI process can be considered as ideal soliton trains. In other words，solitons are tightly connected to the MI. Because MI typically occurs in the same parameter region in which bright solitons are observed，it is considered a precursor to soliton formation. To date，MI has been systematically investigated in biased photorefractive crystals owing to single- or two-photon photorefractive effect^［17-20］. Previous studies of MI in the context of photorefractive processes have been limited to photorefractive materials，in which the change in the refractive index is governed solely by the linear or quadratic electro-optic effect. In fact，incident beams propagating in photorefractive media with both linear and quadratic electro-optic effects experience refractive index modulation，which can strongly influence MI.

In this paper，we present a theoretical study of the MI of broad beams in biased PR crystals with both linear and quadratic electro-optic effects by treating the space-charge field globally and locally，whereby both the one-dimensional global and local MI growth rates are obtained. The properties of these MI growth rates that differ from previous results are discussed，and relevant examples are provided.

In order to investigate the modulational instability of a broad optical beam in biased PR crystals with both the linear and quadratic electro-optic effects，we consider a broad optical beam that propagates along the $z$-axis，where the PR crystal is placed with its principal axes aligned with the $x$，$y$ and $z$ directions of the system. The polarization of the broad beam and the external bias electric field are both assumed to be parallel to the $x$-axis. For simplicity，only $x$-axis diffraction will be considered，and any loss effects are neglected in our analysis.

The optical field of the incident beam is expressed as the slowly varying envelopes $\stackrel{⃑}{E}\left(x,z\right)=\widehat{x}\phi \left(x,z\right)$$\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}kz\right)$，where $k=k_{\mathrm{0}}{n}_{\mathrm{e}}=\left(\mathrm{2}\mathrm{\pi }/{\lambda }_{\mathrm{0}}\right)n_{\mathrm{e}}$ with $n_{\mathrm{e}}$ being the unperturbed index of refraction and ${\lambda }_{\mathrm{0}}$ the free-space wavelength，and $\widehat{x}$ is the unit vector along $x$. Under the above conditions，the optical beam^［2，21］ satisfies the following equation

The change of nonlinear refractive index $\mathrm{\Delta }n$ is governed by

where $\phi \left(x,z\right)$ is the slowly varying envelope of the optical beam and $r_{\mathrm{33}}$ and $g_{\mathrm{e}\mathrm{f}\mathrm{f}}$ are the linear and effective quadratic electro-optic coefficients of the PR crystal，respectively. Further，$k=k_{\mathrm{0}}{n}_{\mathrm{e}}=\left(\mathrm{2}\mathrm{\pi }/{\lambda }_{\mathrm{0}}\right)n_{\mathrm{e}}$ with $n_{\mathrm{e}}$ the unperturbed index of refraction and ${\lambda }_{\mathrm{0}}$ the free-space wavelength，and ${\epsilon }_{\mathrm{0}}$ and ${\epsilon }_{\mathrm{r}}$ denote vacuum and relative dielectric constants，respectively. By substituting Eq.（2） into Eq.（1），we obtain the following paraxial equation：

where ${\beta }_{\mathrm{1}}=k_{\mathrm{0}}{n}_{\mathrm{e}}^{\mathrm{3}}{r}_{\mathrm{33}}/\mathrm{2}$，${\beta }_{\mathrm{2}}=k_{\mathrm{0}}{n}_{\mathrm{e}}^{\mathrm{3}}{g}_{\mathrm{e}\mathrm{f}\mathrm{f}}{\epsilon }_{\mathrm{0}}^{\mathrm{2}}{\left({\epsilon }_{\mathrm{r}}-\mathrm{1}\right)}^{\mathrm{2}}/\mathrm{2}$，$U={\left(\mathrm{2}{\eta }_{\mathrm{0}}{I}_{\mathrm{d}}/n_{\mathrm{e}}\right)}^{-\mathrm{1}/\mathrm{2}}\phi \left(x,z\right)$，and the power density of the broad beam has been normalized to the so-called “dark-irradiance” $I_{\mathrm{d}}$，i.e.，$I=\left(n_{\mathrm{e}}/\mathrm{2}{\eta }_{\mathrm{0}}\right){\left|\phi \right|}^{\mathrm{2}}$ with ${\eta }_{\mathrm{0}}={\left({\mu }_{\mathrm{0}}/{\epsilon }_{\mathrm{0}}\right)}^{\mathrm{1}/\mathrm{2}}$. Based on the transport model of Kukhtarev et al.，the space-charge field $E_{\mathrm{s}\mathrm{c}}$^［3］ in the material is approximately given by

For the bright type，${\left|U\right|}^{\mathrm{2}}=\mathrm{0}$ at $x\to \pm \mathrm{\infty }$，where $E_{\mathrm{0}}$ represents the value of the space-charge field at $x\to \pm \mathrm{\infty }$，i.e.，$E_{\mathrm{0}}=E_{\mathrm{s}\mathrm{c}}\left(x\to \pm \mathrm{\infty },z\right)$，$e$ is the charge，$K_{\mathrm{B}}$ is Boltzmann's constant，$T$ is the absolute temperature，and $N_{\mathrm{A}}$ is the acceptor or trap density. It is worth noting that the diffusion term ${\left(\mathrm{2}k\right)}^{-\mathrm{1}}\left({\partial }^{\mathrm{2}}U/\partial x^{\mathrm{2}}\right)$ and the spatial derivatives of $E_{\mathrm{s}\mathrm{c}}$ and $U$ can be omitted in Eqs.（3）and（4），because $U$ remains relatively constant over a large range of $x$ for broad beams. Under the above conditions，the space-charge field $E_{\mathrm{s}\mathrm{c}}$ is given by

We begin our analysis by treating the space-charge field equation in Eq.（4） globally. In this case，

where $r$ is defined as $r=I_{\mathrm{m}\mathrm{a}\mathrm{x}}/I_{\mathrm{d}}=I\left(\mathrm{0}\right)/I_{\mathrm{d}}$. In what follows，we discuss the stability of the above solution by making the following ansatz：

here $\epsilon \left(x,z\right)$ represents an added weak spatial perturbation term and satisfies the condition $\left|\epsilon \left(x,z\right)\right|\ll$$r^{\mathrm{1}/\mathrm{2}}$. Substitution of Eq.（7） into Eqs.（3）and（4）yields:

where $E=E_{\mathrm{s}\mathrm{c}}-E_{\mathrm{01}}$，$E_{\mathrm{01}}=E_{\mathrm{0}}/\left(\mathrm{1}+r\right)$，$\upsilon =E_{\mathrm{01}}$$\left({\epsilon }_{\mathrm{0}}{\epsilon }_{\mathrm{r}}/e{N}_{\mathrm{A}}\right)$，$\mathrm{\Delta }=\left(K_{\mathrm{B}}T/e\right)\left({\epsilon }_{\mathrm{0}}{\epsilon }_{\mathrm{r}}/e{N}_{\mathrm{A}}\right)$. The space-charge field in the spatial-frequency space $\widehat{E}$ can then be obtained by employing the Fourier transform，satisfying

where $\widehat{\epsilon }$ is the Fourier transform of the spatial perturbation in the spatial-frequency space. The spatial perturbation $\epsilon \left(x,z\right)$ can be expressed as the sum of the following two terms：

It is easy to show that

where $\delta \left(k_x\right)$ is the delta function. Substituting this form of $\left(\widehat{\epsilon }+{\widehat{\epsilon }}^{•}\right)$ back into Eq.（10） and taking an inverse Fourier transform allows the space-charge field in real space to be obtained from the following equation：

In order to simplify calculations，$G\left(p\right)$ is defined as

Eq.（13） then reduces to

By substituting Eqs.（11）and（13）into Eq.（8） and keeping only the synchronous terms，we obtain to the following coupled differential equations：

We then decouple Eqs.（16） into an equivalent set of ordinary differential equations as follows：

From these two equations，we can directly obtain the global modulational instability gain $g_{\mathrm{g}\mathrm{l}}$ as

where $\mathrm{R}\mathrm{e}\left\{\cdot \right\}$ denotes the real part of a complex variable. From Eq.（13），it is clear that the MI gain is an even function of $p$，and that its value reaches zero at $p=\mathrm{0}$.

The local MI process deserves special consideration. In the next section，the local MI process is investigated theoretically by neglecting higher-order effects in the space-charge field. Under strong bias conditions，for a broad incident optical beam the diffusion effect can be neglected；that is，all terms associated with the diffusion process（i.e.，$K_{\mathrm{B}}T$ terms）may be omitted in Eq.（4）. Additionally，the dimensionless term $\left({\epsilon }_{\mathrm{0}}{\epsilon }_{\mathrm{r}}/e{N}_{\mathrm{A}}\right)\left(\partial E_{\mathrm{s}\mathrm{c}}/\partial x\right)$ is typically significantly less than unity^［17］. Accordingly，the paraxial Eq.（3） reduces to

Eq.（19） takes the form of a nonlinear Schrödinger equation with a saturable nonlinearity. In our preliminary work，we obtained the solitary wave solutions of dark，bright，and grey solitons in a steady-state regime，various of whose characteristics and properties have been discussed in detail elsewhere^［9-10］. In what follows，by introducing Eq.（7） into Eq.（19） and linearizing it in $\epsilon$，the local MI properties of Eq.（19） can be investigated using the following evolution equation.

Substituting Eq.（11） into Eq.（20），we obtain

The above equations can be decoupled into an equivalent set of ordinary differential equations：

and the local MI gain can be obtained directly from Eq.（23）：

Moreover，the maximum MI gain and its associated spatial frequency can be readily determined as follows：

To illustrate our results，we consider a single PMN-0.33PT crystal that exhibits maximal transparency，very good optical clarity，and low propagation loss. The parameters of PMN-0.33PT are $n_{\mathrm{e}}=\mathrm{2.562}$，$r_{\mathrm{33}}=\mathrm{182}\mathrm{p}\mathrm{m}/\mathrm{V}$，$g_{\mathrm{e}\mathrm{f}\mathrm{f}}=\mathrm{0.06}{\mathrm{m}}^{\mathrm{4}}\cdot {\mathrm{C}}^{-\mathrm{2}}$，${\epsilon }_{\mathrm{r}}=\mathrm{5378}$，$N_{\mathrm{A}}=\mathrm{3.7}\times {\mathrm{10}}^{\mathrm{22}}{\mathrm{m}}^{-\mathrm{3}}$^［22-25］. The other parameters are set as ${\lambda }_{\mathrm{0}}=\mathrm{632.8}\mathrm{n}\mathrm{m}$ and $x_{\mathrm{0}}=\mathrm{40}\mathrm{\mu }\mathrm{m}$. Based on these parameters，${\beta }_{\mathrm{1}}=\mathrm{0.0152}$，${\beta }_{\mathrm{2}}=\mathrm{1.1354}\times {\mathrm{10}}^{-\mathrm{8}}$，and $k=\mathrm{2.5439}\times {\mathrm{10}}^{\mathrm{7}}$.

Figs. 1 and 2 show the dependence of the global and local MI gains（i.e.，$g_{\mathrm{g}\mathrm{l}}$ and $g_{\mathrm{l}\mathrm{c}}$）on $p/k$ with three different values of $E_{\mathrm{0}}$ for the same $r$，where the dimensionless ratio $p/k$ represents the angle（in radians）at which the plane-wave components of the $\epsilon \left(x,z\right)$ perturbation propagate with respect to the quasi-plane-wave optical beam. From these two figures，it is evident that the two MI gains increase with increasing $E_{\mathrm{0}}$. Given that both gains are symmetrical about $p$，only the positive branch will be considered in our subsequent analyses. From Fig. 1，we can observe that the global MI gain curve reaches two different peaks：one appears in the low spatial-frequency domain，defined as $g_{\mathrm{p}\mathrm{1}}$，and the other occurs in the high-frequency region，defined as $g_{\mathrm{p}\mathrm{2}}$. Moreover，it is evident that both $E_{\mathrm{0}}$ and $r$ affect these two different peaks. Next，we study the effects of $E_{\mathrm{0}}$ and $r$ in isolation on MI gains using a variable-controlling approach.

We begin with the effect arising from the biased field $E_{\mathrm{0}}$. Fig. 3 depicts the curve of the global MI gains versus $E_{\mathrm{0}}$ associated with the two global peaks $p_{\mathrm{1}}/k$ and $p_{\mathrm{2}}/k$. Unlike the cases investigated previously，we find that the global MI gains are asymmetric with respect to the polarity of $E_{\mathrm{0}}$，that is，the global MI gains depend not only on the absolute strength of $E_{\mathrm{0}}$，but also on the polarity of the external bias electric field $E_{\mathrm{0}}$. This is because the global MI gains are governed by both linear and quadratic electrooptic effects. By altering the polarity of $E_{\mathrm{0}}$，the sign of the linear electro-optic term $\left({\beta }_{\mathrm{1}}{E}_{\mathrm{0}}\propto E_{\mathrm{0}}\right)$ changes as well. However，the sign of the quadratic electro-optic term $\left({\beta }_{\mathrm{2}}{E}_{\mathrm{0}}^{\mathrm{2}}\propto E_{\mathrm{0}}^{\mathrm{2}}\right)$ is not influenced by the polarity change of $E_{\mathrm{0}}$，so in the case of an externally biased field of equal magnitude but opposite polarity，the photorefractive effect is weakened and even counteracted by the interaction between the linear and quadratic electric-optic effects. Global MI gains can also be adjusted by altering the polarity of $E_{\mathrm{0}}$ in addition to changing its strength. Furthermore，Fig. 3 shows that $g_{\mathrm{p}\mathrm{1}}$ exceeds $g_{\mathrm{p}\mathrm{2}}$ when $E_{\mathrm{0}}>\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{V}/\mathrm{m}$ with the positive bias field.

In addition，Fig. 4 illustrates the dependence of the global MI gain peak $g_{\mathrm{p}\mathrm{1}}$ and $g_{\mathrm{p}\mathrm{2}}$ on $r$ for $E_{\mathrm{0}}=\mathrm{1.8}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{V}/\mathrm{m}$. As shown in Fig. 3，$g_{\mathrm{p}\mathrm{1}}$ and $g_{\mathrm{p}\mathrm{2}}$ attain their maxima at $r\approx \mathrm{1}$. Moreover，the peak $g_{\mathrm{p}\mathrm{1}}$ decreases rapidly when $r\ll \mathrm{0.1}$ and $r\gg \mathrm{10}$，and all $g_{\mathrm{p}\mathrm{1}}$ and $g_{\mathrm{p}\mathrm{2}}$ tend to be stable for $r\gg \mathrm{10}$.

Fig. 5 shows the variation of the two gain peaks $g_{\mathrm{p}\mathrm{1}}$ and $g_{\mathrm{p}\mathrm{2}}$ as a function of both $E_{\mathrm{0}}$ and $p/k$ when $r=\mathrm{1}$，from which we can see that in the low bias voltage region $g_{\mathrm{p}\mathrm{1}}<g_{\mathrm{p}\mathrm{2}}$，and $g_{\mathrm{p}\mathrm{1}}$ will exceed $g_{\mathrm{p}\mathrm{2}}$ when $E_{\mathrm{0}}$ is higher than a certain value. Moreover，in the low spatial-frequency domain $g_{\mathrm{p}\mathrm{1}}$ tends to increase linearly with increasing $E_{\mathrm{0}}$；however，$g_{\mathrm{p}\mathrm{2}}$ increases slowly in the high spatial-frequency domain under different $p/k$ conditions. In fact，under strong bias conditions the MI process should be treated with the local process and Eqs.（16） can be simplified to Eq.（21）.

Fig. 6 shows the dependence of the global gains $g_{\mathrm{g}\mathrm{l}}$ on $p/k$ with different $r$ in the case of $E_{\mathrm{0}}=\mathrm{1.8}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{V}/\mathrm{m}$. The maximum value of $g_{\mathrm{g}\mathrm{l}}$（i.e.，$g_{\mathrm{m}\mathrm{a}\mathrm{x}}$）appears in the low-spatial-frequency domain when $r$ is small；however，$g_{\mathrm{m}\mathrm{a}\mathrm{x}}$ arises in the high-spatial-frequency regime when $r$ is relatively large.

In contrast to the results obtained previously，we find that both $g_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and $p_{\mathrm{m}\mathrm{a}\mathrm{x}}$ follow the quadratic polynomials of $E_{\mathrm{0}}$，and the maximum MI gains increase rapidly with an increase in the external bias field according to Eqs.（24）and（25），respectively. In addition，the MI gain can be adjusted by altering the polarity of the external bias field，even for the same bias field strength. Moreover，further analysis of Eqs.（24）and（25）shows that $g_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and $p_{\mathrm{m}\mathrm{a}\mathrm{x}}$ reach a maximum when $r=\mathrm{0.88}$.

We investigated the one-dimensional modulational instability of a broad optical beam propagating in biased PR crystals with both linear and quadratic electro-optic effects under steady-state conditions. Both the one-dimensional global and local MI growth rates were obtained by treating the space-charge field equation globally and locally. It was shown that the global and local modulational instability gains were governed simultaneously by the strength and polarity of the external bias field and the ratio of the intensity of the broad beam to that of the dark irradiance. This means that the spatial period of spontaneously generated filaments can be controlled by adjusting the factors mentioned above. Moreover，under a strong bias field，the results obtained from these two methods were in good agreement in the low spatial frequency regime. The instability growth rates increase with the bias field，and the maximum instability growth occurs at $r=\mathrm{0.88}$.