Three-Dimensional Scalar Time-Dependent Photorefractive Beam Propagation Model

Abstract

This paper presents an open-source time-dependent three-dimensional scalar photorefractive beam propagation model (PRProp3D) based on the well-known split-step method. The angular spectrum method is used for the diffractive steps, and the nonlinearities accumulated at the end of each diffractive step are applied using spatially varying phase screens. Comparisons with previously published experimental results are given for image amplification, photorefractive amplified scattering (fanning) and photorefractive screening solitons. Artifacts can be mitigated by use of step sizes less than 5~10 micrometers and by careful choice of the transverse computation grid size to ensure adequate sampling. Wraparound effects associated with the use of discrete Fourier transforms are mitigated by apodization and beam centering.

1. Introduction

In the latter part of the 20th century, optical computing became a subject of extensive research fueled in part by its inherent parallelism and the rapid development of laser technology. The new availability of bright coherent optical sources simplified phase-coherent optical image processing. Simple lenses could be used to provide Fourier transforms of optical images in parallel at the speed of light (notwithstanding the fact that wired electrical signals also travel at the speed of light) [1]. Pattern recognition using real-time holography became feasible, and it was not long before content-addressable memories and designs for neural networks began appearing. Electronic computational power continued to improve dramatically in terms of both processing speed and memory capacity. By the beginning of the 21st century, the capabilities of electronic neural networks exceeded even the potential of optical systems. With respect to data storage, optical compact disks and DVDs have now become practically obsolete, replaced by high-capacity solid-state drives and cloud storage. Nevertheless, optics remains vital to the infrastructure supporting modern computing and may still, with careful engineering, contribute to its capabilities. An interesting discussion of these issues may be found in a recent perspective on the physics of optical computing [2]. One of several instances where optics retains its hold is data transmission over optical fibers. Particularly relevant here is the continuing interest in the design of neural networks using multimode fibers. As early as 1990, Aisawa et al. described image classification using multimode fibers followed by a multilevel network [3]. Ancora et al. have shown that pattern classification can be accomplished using only a multimode fiber coupled with a single logistic regression layer [4].

Propagation distances in optical fibers can result in the appearance of nonlinear effects. While these effects are often deleterious and require compensation, they can be used to good effect, for example, in supercontinuum generation [5]. Tegin et al. [6] have shown that nonlinear effects in multimode fibers can be used to design scalable optical neural networks whose power requirements in the computation stage have the potential to be lower than that of electronic neural networks. Elements of a data set are represented optically using a spatial light modulator and transmitted through a multimode optical fiber. The authors showed that nonlinear interactions between the modes of the fiber can reformat the data into a “nearly linearly separable output data space” at the output of the fiber.

The original motivation behind the research presented in this paper was to explore the possibility of using another low-power nonlinear optical effect, the photorefractive effect [7], to perform similar transformations. This effect is a form of real-time holography using optically induced refractive index variations in a medium to enable spatial mode mixing for use in optical transformations and computation. One possibility to use photorefractive mixing in a way that might be applied for model training is amplified photorefractive scattering, known as fanning [8]. The results presented here show that fanning can be well modeled with step sizes of the order of 2~5 μm. The code can also be used to investigate other forms of photorefractive nonlinear coupling and generalized to other nonlinearities.

It was during the earlier period of intensive research into optical computing that most progress was made in the understanding and application of the photorefractive effect. A breakthrough took place with the discovery that ferroelectric crystals such as barium titanate, potassium niobate, and strontium barium niobate could support strong nonlinear optical interactions between low-power, milliwatt-level, continuous-wave laser beams. This enabled the realization of many nonlinear optical devices that could otherwise only be achieved with high-power pulsed lasers [9] or with resonant materials such as sodium vapor [10]. These included distortion-compensating phase conjugate mirrors with reflectivity exceeding unity [11], photorefractive oscillators [12], self-pumped phase conjugate mirrors [12,13] and the realization of photorefractive spatial solitons [14].

In the early 90s, we and others developed two-dimensional beam propagation codes to model photorefractive effects such as image amplification by two-beam coupling and amplified photorefractive scattering [8,15,16,17,18]. This amplified scattering is commonly known as fanning because it gives rise to a characteristic fan-shaped pattern observed in the far field. Over the past 30 years, as computing power has improved, it has become feasible to extend this model to three dimensions including time dependence. In this paper, I describe the development of code implementing such a model (PRProp3D) for scalar fields. It is available as open-source software running in Google COLAB using its GPUs at https://github.com/mcroning/PRProp3D. Limitations on the longitudinal step size are investigated in the context of tradeoffs between computation time and storage and accuracy. The model has been checked against experimental data on image amplification, beam fanning, and time dependence. Sections in the paper discuss the following:

• Image amplification with a comparison to the experimental data published by Fainman et al. [19].
• Beam fanning with a comparison to prior theory by Zozulya et al. [16,20] and new experiments.
• Time dependence with a comparison to high-gain beam amplification experiments by Fischer et al. [21].
• The theoretical appearance of high-order diffraction at high gain described by Brown and Valley [17].
• Photorefractive screening solitons.

The core of the program is a split-step fast Fourier transform beam propagation method combined with a dimensionless model of the photorefractive effect. The Fourier method was chosen over finite difference methods because of its simplicity and comparable accuracy and computing time requirements. Propagation is divided into equal longitudinal steps of length dz through the interaction region. Each computation step is divided into two parts. The first involves linear diffraction using the method of the angular spectrum of plane waves [1]. Linear absorption is neglected in the present realization but could be added easily. The second part is to apply the photorefractive nonlinearity accumulated during the linear propagation step as thin transparency. Its proper usage requires an understanding of the origin and management of several computational features, some of which are artifacts due to the periodic nature of the fast Fourier transform and the discrete step sizes involved in computational beam propagation. These artifacts are covered in detail in the appendix and include the following:

• The appearance of excess high-order diffraction over interaction lengths nominally in the Bragg regime: this effect decreases with decreasing longitudinal step size but does not tend toward to zero.
• Coupling to the longitudinal grating formed by the equally spaced nonlinear transparencies: This feature is especially apparent in models of the fanning effect resulting in the appearance of rings centered on the direction of the incident beam in the far field where the fanning pattern is usually observed. The diameter of these rings is approximately proportional to the inverse of the square root of the step size. Step sizes less than 5 μm will usually push the diameter of the innermost ring beyond the region of interest in the far field.
• Wraparound effects due to the intrinsic periodic nature of the discrete Fourier transform: This last effect could be controlled by using finite-difference beam propagation methods instead of FFT methods, but this would introduce additional complexity. We choose to continue to use FFT methods with the understanding that wraparound effects need to be recognized and minimized by choosing transverse apertures large enough to keep the main parts of the interacting beams away from the boundaries. We will see below that the periodic nature of the FFT is advantageous when seeking to verify the code against standard plane wave theory.

2. Methods

In keeping with the overall philosophy to keep the nonlinear propagation model as simple as possible while maintaining the principle physical behaviors, we have chosen to formulate the model using scalar diffraction. Extension to the full vector case including the use of the dielectric tensor for beam propagation and the electro-optic tensor for nonlinear grating formation would lead to a significant increase in complexity.

• The vector angular spectrum of the plane wave method requires the separate handling of the ordinary and extraordinary vector modes, which are both dependent on the transverse spatial frequencies [22].
• The use of the full electro-optic tensor would require knowledge of the longitudinal components of the photorefractive space charge field, which depend on the three-dimensional gradient of the optical intensity $\nabla I\left(z\right)$. The methods used in the past and in this paper ignore the longitudinal components of the intensity gradient and space charge electric field: the grating is calculated assuming that it is independent of the optical fields at neighboring planes.

We leave these modifications to future research. Nevertheless, the scalar quasi-3D method described below does match many experimental results.

2.1. Beam Propagation

The calculation of the propagation of an optical field in a photorefractive nonlinear medium requires the solution of Maxwell’s equations in a medium whose refractive index depends slowly (a millisecond-to-second time scale) on the optical intensity. At these time scales, the propagation of light through the medium (picosecond time scale) is practically instantaneous, so we can use the single-frequency version of the wave equation:

$${\nabla }^{2}E+{\left({\displaystyle \frac{\omega n\left(t\right)}{c}}\right)}^2=0$$

(1)

where ω is the radian frequency and c is the speed of light. The aim is to monitor the amplitude of the optical field as it responds slowly to the optically induced development of an inhomogeneous refractive index distribution.

The model describes the optical field in Cartesian coordinates with x and y transverse and z longitudinal. The input field is given in terms of its scalar amplitude A(x,y) at the input plane z = 0. The propagation of the field amplitude A(x,y) by distance dz is accomplished using the angular spectrum of the plane wave method [1]. The Fourier transform of the field is multiplied by the propagator h for propagation over distance dz:

$$h\left(f_x,f_y,dz\right)=\exp \left(2\pi i{f}_{z}dz\right)$$

(2)

$$f_z=\sqrt{f_0^2-f_y^2-f_x^2}$$

(3)

where f₀ = n/l, l is the wavelength of the light, n is refractive index, and f_x and f_y are the transverse spatial Fourier frequencies of the optical amplitude.

Inverse Fourier transformation returns the field to real space. The Fresnel approximation can be recovered by keeping the first term in a Taylor expansion of the square root in the exponent, but this leads to significant errors for the high angles of propagation that are commonly used. Discrete Fourier transform wraparound effects are minimized by using a Tukey apodization window at each step [23].

2.2. Photorefractive Nonlinearity Model

The most widely used theory of the photorefractive effect (Kukhtarev et al. [24]) is based on a standard drift and diffusion semiconductor model. Its spatial behavior can be shown to be equivalent to the hopping model (Feinberg et al. [25]). There are minor differences in the temporal behavior due to the presence of the charge mobility in the Kukhtarev model. The hopping model can be readily reduced to the single partial differential equation shown in Equation (4). I have chosen to use the hopping model because of its relative simplicity. The model was developed using the concept of charge carriers hopping between sites in a crystal, with the hopping probability dependent on the local optical intensity and electric field. The key is that it can be rewritten in terms of continuous variables if it is assumed (1) that only nearest neighbor hopping is important, (2) that trapping sites can be treated as being uniformly spaced, and (3) that there are many trapping sites in a single grating period so their distribution can be treated as a continuous variable. These approximations lead to Equation (4), which describes development of the internal electric field E in terms of the spatially varying optical intensity I, which includes the equivalent dark intensity I_d to account for thermal relaxation, as well as any additional uniform background intensity (in the treatment of spatial solitons, the addition of a uniform optical background and an external DC electric field E_app are essential. See Section 3.6). Primes represent spatial derivatives in the transverse x direction. All variables are normalized with respect to their characteristic quantities shown in

Table 1. Variable normalizations.

Variable	Normalization	Value
Distance	1/k0	$k_0=q\sqrt{{\displaystyle \frac{p_d}{\epsilon k_{B}T}}}$
Electric Field	E0	$E_0=\sqrt{{\displaystyle \frac{k_{B}T{p}_d}{\epsilon }}}$
Time	t0	$t_0={\left(D{s}^2{I}_0{k}_0^2\right)}^{-1}$
Intensity	I0	Sum of peak intensities of the input beams

$${\displaystyle \frac{\partial E}{\partial t}}=E_{app}{I}_d-\left(EI-I^{\prime }\right)\left(1+E^{\prime }\right)+E^{″}I$$

(4)

One is often only interested in finding solutions at a steady state to avoid the overhead of a time-dependent solution. In these cases, we can solve for the space charge field directly by setting the time derivative to zero in the above equation. This requires the solution of a nonlinear differential equation, leading to additional computational overhead. To avoid this, we can use a version linearized with respect to the electric field at the expense of moderately reduced accuracy:

$$E+E_{app}{E}^{\prime }-E^{″}={\displaystyle \frac{E_{app}{I}_d+I^{\prime }}{I}}$$

(5)

We have found that the errors resulting from this linearization are negligible in the absence of an applied DC field and small (few percent in intensity) with an applied field when the optical intensity at the computation boundary is zero. Except for Section 3.6 on solitons, all the results described below are without an applied DC electric field.

3. Results

3.1. Two-Beam Coupling

As a check on the validity of the program, we compared its output with predictions from the standard steady-state coupled wave equations for two-plane waves intersecting at half angle θ in the slowly varying envelope approximation [24,26]:

$$\begin{array}{c}{\displaystyle \frac{d{A}_1}{dz}}=-{\gamma }_p{\displaystyle \frac{A_1{A}_2^{\ast }}{I_0}}{A}_2\\ {\displaystyle \frac{d{A}_2^{\ast }}{dz}}=+{\gamma }_p{\displaystyle \frac{A_1{A}_2^{\ast }}{I_0}}{A}_1^{*}\end{array}$$

(6)

where A₁ is the electric field amplitude of the pump beam, A₂ is the amplitude of the signal beam and γ_p is the plane wave coupling constant in the absence of an applied DC electric field, first derived in ref. [24]:

$${\gamma }_p={\displaystyle \frac{\pi {n^{3}r}_{eff}{E}_0}{2\lambda \mathrm{c}\mathrm{o}\mathrm{s}\theta }}{\displaystyle \frac{k_g}{1+k_g^2}}={\displaystyle \frac{\gamma }{\mathrm{c}\mathrm{o}\mathrm{s}\theta }}{\displaystyle \frac{2k_g}{1+k_g^2}}$$

(7)

where n is the refractive index of the propagation medium, k_g the grating wavenumber normalized to the characteristic wavenumber k₀ and E₀ is the characteristic space charge field. r_eff is the effective electrooptic coefficient [27]. γ is the part of the coupling constant that is independent of propagation angles. The analytical solutions of these equations [24,26] show that the output beam ratio r_out in terms of the input beam ratio r_in is given by

$$r_{out}=r_{in}exp\left(2{\gamma }_p\mathcal{l}\right)$$

(8)

where the beam ratio is defined as the ratio of the intensity of beam 2 to that of beam 1.

This equation be used to find the coupling constant length product (gain) as follows:

$${\gamma }_p\mathcal{l}={\displaystyle \frac{1}{2}}\mathrm{l}\mathrm{n}\left({\displaystyle \frac{r_{out}}{r_{in}}}\right)$$

(9)

Figure 1 shows the gain calculated with no applied DC field using Equation (9) from the output of the program when run with input gain $\gamma \mathcal{l}$ = 3 using both the generalized refractive index grating (Equation (4)) and the plane wave sinusoidal refractive index grating (Equation (A11)). As can be expected, agreement with the plane wave coupled wave theory was better when using its specialized grating model. This infinite plane wave calculation is made possible by taking advantage the fact that the periodic nature of the discrete Fourier transform implies that the beams have infinite transverse profiles. To prevent wraparound effects, the beam crossing angles are chosen so that the grating is continuous across the boundaries and the peaks in Fourier space overlap from one Brillouin zone to the next. The apodization window is not used. Allowable symmetric beam input angles are given by

$${\theta }_1=-{\theta }_2={\sin }^{-1}\left(\lambda /\left(2^n\Delta x\right)\right)$$

(10)

where the integer n is subject to $1\le n\le {\log }_2(N_x-1)$. Δx is the transverse grid spacing and N_x is the number of grid samples in the transverse x direction.

In 1993, Brown and Valley showed that two-dimensional beam propagation calculations for two-beam coupling at large coupling constants predicted the appearance and amplification of beams (kinky beams) diffracted by grating harmonics [17]. The amplification of light scattered by these grating components would normally be considered forbidden by the Bragg condition and is not observed in practice. The authors hypothesized that the lack of the experimental confirmation of these higher-order beams is a result of their being overwhelmed by fanning. Our three-dimensional program reproduces these findings, as shown in Figure 2.

3.2. Image Amplification

One of the applications of photorefractive two-beam coupling is coherent image amplification. Several authors have studied the fidelity of the process under various conditions. To see how well our code can model existing image amplification data, we modeled the experiments in photorefractive barium titanate performed by Fainman et al. and treated by them using standard coupled wave theory [19]. They presented their results in terms of absolute signal amplification G₀:

$$G_0={\displaystyle \frac{\left(1+r\right)\exp \left(2\mathsf{\gamma }\mathcal{l}\right)}{1+r\exp \left(2\mathsf{\gamma }\mathcal{l}\right)}}$$

(11)

where r = r_in is the ratio of the intensity of the signal beam to that of the pump beam. Figure 3 shows theoretical plots of normalized gain $G_0/G_{0sat}$ where the saturated, small signal gain is $G_{0sat}=\exp \left(2\gamma \mathcal{l}\right)$.

Reference [19] treats image amplification optimization as its main result, and we reproduce some of the fidelity tests of that paper. While the group labels on the resolution targets are different, the computer input image is scaled to the same size as the section of the chart shown in the amplified image from ref. [19]. Figure 4 and Figure 5 show the comparisons for the small signal case treated in the Fainman paper. The dark stripe on the right of the output image in Figure 4 is caused by the apodization at the grid edges.

Reference [19] shows that when the ratio of the input signal intensity to the pump intensity increases, the amplified image experiences edge enhancement because the high-frequency parts of the image experience larger photorefractive gain than the low-frequency components. This effect is reproduced qualitatively by our model as shown in Figure 6 and Figure 7. The fact that the specifics of which edges are enhanced do not always correspond in theory and experiment may be due to differences in theoretical and experimental amplitude distributions of the beams illuminating the Air Force chart.

3.3. Photorefractive Amplified Scattering

The model allows for the three-dimensional simulation of photorefractive amplified scattering (fanning) seeded by optical scatterers [28]. This fanning serves as the seed for various photorefractive oscillators such as unidirectional ring oscillators and self-pumped phase conjugate mirrors [27] and contributes to noise in image amplification.

In 1994–5, Zozulya et al. published a two-dimensional computer model of the fanning effect and of several photorefractive phase conjugate mirrors [16,20]. This model used the Crank–Nicholson integration of the coupled wave equations and a relaxation method for the calculation of the photorefractive nonlinearity. We used their results as the basis for comparison for the output of our three-dimensional model. Following the method of Zozulya et al. [20] we modeled scattering noise in the crystal as a series of phase screens applied at each propagation step. The phase screens were represented as $S\left(x,y\right)=exp\left(i\phi \left(x,y\right)\right)$ where φ is a two-dimensional array of normally distributed random numbers filtered spatially with a Gaussian kernel with standard deviation σ. This standard deviation represents the scattering bandwidth in the transverse spatial spectrum. The corresponding angular scattering bandwidth is sin⁻¹(σ/λ). The mean square deviation of φ is given by $〈{\phi }^2〉=ϵ/N_s$ where $\mathit{ϵ}$ is a parameter representing the strength of the scattering and N_s is the number of scattering screens.

Figure 8 shows a fanning profile taken from Zozulya’s paper [20] and Figure 9 shows the corresponding results from our program. The relevant parameters are given in

Table 2. Parameters for static and dynamic fanning calculations.

Parameter	Static Calculation	Dynamic Calculation
Coupling $\gamma \mathcal{l}$	10.0	10.0
x aperture mm	3.0	1.5
y aperture mm	2.0	0.75
Number of x samples	16,384	4096
Number of y samples	4096	2048
Crystal length mm	5.0	5.0
Scattering corr. length μm	0.4	0.4
Longitudinal step size μm	2	8
Wavelength μm	0.488	0.488
Beam waist μm	0.3	0.3
Dark intensity	0.01	0.01
Angle of incidence radians	0.3	0.3
Azimuth of incidence rad	0.0	0.0
Time steps	NA	160
End time normalized	NA	10
Tukey window parameter	0.2	0.2
GPU batches	1	5
A100 calc. time hh:mm	0:05	1:57
A100 CPU RAM GB	7.2	71.9
A100 GPU RAM GB	31	38

.

To provide a comparison of the transverse output of the code with the experiment, we show in Figure 10 an image of the fanning of a 10 mW 488 nm beam from a Toptica iBeam smart diode laser (TOPTICA Photonics AG, Germany) after passage through a 6.5 mm long crystal of barium titanate (BAE Systems Electronics, formerly Sanders Associates, Nashua NH, USA) in our laboratory. The corresponding output of the model using beam waist 3 mm and longitudinal step size 2 μm and $\mathit{\gamma }\mathcal{l}$=10 is shown in Figure 11. The axes of these images are in terms of free-space propagation direction cosines:

$$\begin{array}{c}cos{\theta }_x=f_x\lambda \\ cos{\theta }_y=f_y\lambda \end{array}$$

(12)

These are the components of the individual wavevectors in the beam along the x and y directions, respectively. The direction cosines for propagation within the crystal may be found by dividing by the refractive index. In the small angle paraxial approximation, these direction cosines can be interpreted as propagation angles in radians.

A disadvantage of using a scalar model rather than a vector model with the full electrooptic tensor is we do not reproduce the characteristic three-lobed pattern of fanning in BaTiO₃ as modeled in the undepleted pump approximation by Montemezzani et al. [29]. The development of these lobes can be seen in the experiment video (Movie S1) corresponding to the static image of Figure 10. The central scattering angle calculated is less than that observed experimentally, possibly because of a mismatch in the characteristic grating wavenumbers. Notice that in the video the underlying speckle pattern remains fixed as opposed to in observations characteristic of high intensities where thermal effects result in a time-varying speckle pattern.

A question that arises about these fanning calculations is why large coupling constants of the order of $\gamma \mathcal{l}=10$ are required to reproduce experimental observations, while for the modeling two-beam coupling of gaussian beams or image amplification, $\gamma \mathcal{l}=3$ suffices. One possible explanation is that the effective interaction length for fanning from narrow beams is limited. Figure 12 shows how the calculated fanning efficiency increases with the beam waist. Because of computer limitations, it is often necessary to perform fanning calculations for beam waists that are smaller than those encountered in the lab.

3.4. Time Dependence

The model includes the option to follow the time dependence of the buildup of the interactions by using the numerical Euler integration of Equation (4)

$$E\left(t_n\right)=E\left(t_{n-1}\right)+\Delta t{\displaystyle \frac{\partial E}{\partial t}}\left(t_{n-1}\right)$$

(13)

In all calculations that follow, we use normalized time and normalized distance (Table 1).

One of the main features of the photorefractive effect is that when the intensity of the interacting beams is above the dark intensity I_d, the steady-state coupling strength is approximately independent of intensity. The tradeoff is that the response becomes slower as the intensity drops: when the intensity is larger than the equivalent dark intensity the response time is approximately inversely proportional to intensity. Additionally, when the coupling constant increases, the gain response time also increases. Figure 13 shows a comparison between the time behavior of two-beam coupling gain measured experimentally by Horowitz et al. [21] and that predicted by the present software. These results match those from the basic time-dependent model used in that paper. The dots represent the experimental data, and the solid lines represent our computed response. The response times are in accordance with the analytic estimate derived by Vachss [30]:

$$\begin{gathered} \tau \approx \gamma \mathcal{l}\left(1+\left({\left(1-e^{-1}\right)}^{{\displaystyle \frac{1}{2}}}-0.5\right){\left({\displaystyle \frac{4\pi }{\gamma \mathcal{l}}}\right)}^{{\displaystyle \frac{1}{2}}}\right) : \gamma \mathcal{l}>1.1 \\ \tau \approx 1+0.81\gamma \mathcal{l}+0.148{\left(\gamma \mathcal{l}\right)}^2 : \gamma \mathcal{l}\le 1.1 \end{gathered}$$

(14)

3.5. Time-Dependent Amplified Scattering

Here, we present a time-dependent model of photorefractive fanning. The parameters are the same as for the static study of Section 3.3, except the longitudinal step size had to be increased from 2 μm to 8 μm because of memory limitations. The calculation in the time-dependent case requires the retention of the full three-dimensional grating distribution at the previous time step. This increased demand on memory can be mitigated by batching the grating distribution into the CPU (see Appendix A).

The increased step size resulted in the appearance of ring artifacts whose origin is treated in Appendix C below.

In the context of modeling the time dependence of fanning, we found that the integration could become unstable. We believe this is due to the high spatial frequencies in the fanning grating: the time constant decreases as the grating spatial frequency increases (see Equation (A11)). The stability criterion of Euler integration for simple exponential decay is |Δt/τ − 1| < 1 [31], so we can attempt to ensure the stability of the integration by making the size of the time steps substantially less than the shortest time constant t_min of the system. The grating wavenumber cannot be larger than 2π times the maximum frequency of the model in Fourier space, so using Equation (A11), we choose the step size conservatively to be

$$\begin{array}{cc}4\Delta t={\tau }_{min}& ={\left(1+{\left({\displaystyle \frac{k_{g max}}{k_0^2}}\right)}^2\right)}^{-1}\\ & ={\left(1+{\left({\displaystyle \frac{\pi N_x}{L{k}_0^2}}\right)}^2\right)}^{-1}\end{array}$$

(15)

where N_x is the number of samples in the transverse direction, and L is the transverse aperture of the crystal.

These instabilities would be mitigated by using an implicit integration method such as the backward Euler method. However, the computational expense of the requisite root finding in this nonlinear scenario would be impractical.

Because we have made this empirical choice for time steps and because we are solving a space–time partial differential equation, not an ordinary differential equation as is usually the case for Euler integration, it is important to provide a monitor for instability in the program’s calculation. Such an instability manifests itself in the generation of exponentially growing values in the grating space charge field. We abort the calculation loop upon the occurrence of a NaN in the space charge field. For our purposes, the method used for setting time step sizes is found to be sufficient. The development of adaptive step sizes and better convergence criteria and integration methods could be addressed in future research. Figure 14 shows (a) the fanning far field computed by the time-dependent algorithm at the end of 10 characteristic time constants t₀ and (b) by the static algorithm. The experiment and time-dependent calculation videos can be found in the Supplementary Materials (Movies S1 and S2). The far field generated by the time-dependent model is not quite fully developed, as is evidenced by the more pronounced appearance of the additional artifact in the static model output. The artifacts can be a significant calculational sink of true fanning power and efficiency. The time dependence of the fanning power, shown in Figure 15 is consistent with the coupling constant-related lengthening of two-beam coupling response times described above.

3.6. Solitons

Photorefractive solitons have provided a rich area of research since their discovery in 1991 [32]. In this section, we show how the program can be used to model photorefractive screening solitons. To produce these solitons, a DC electric field E_app is applied across the crystal, usually of the order of 5 kV/cm for barium titanate, and a uniform background illumination is provided throughout the crystal. A tightly focused laser beam is introduced into the crystal and propagates along the longitudinal z direction. A typical beam waist is 10 μm. The beam excites charge carriers that screen the electric field from the propagating beam. Because the crystal is electrooptic, the refractive index in the screened region can be either higher or lower than the refractive index outside the beam, depending on the direction of application of the DC field. If the refractive index is higher, the region occupied by the beam can form a waveguide counteracting the tendency of the beam to diverge via diffraction. This optically induced waveguide supports the propagation of the laser beam as a soliton. Figure 16 shows the results of the model run in time calculation mode after 20 intrinsic time constants. The parameters chosen here match those found in Figure 1 of ref. [33], which studies the dynamics of screening solitons. The formation of the soliton compresses the transverse extent of the beam and gives it an ellipsoidal shape whose major axis alternates between the x and y axis as it propagates. The beam moves transversely by 32 μm in traveling 4 mm in the longitudinal z direction. These results match those reported in ref. [33] well. Supplementary Movie S3 shows a progression through 4 mm of z planes from the input plane to the output plane, showing the lateral shift of the beam and the changes in its transverse shape. Supplementary Movie S4 shows the temporal development of the soliton at the output face of the crystal over a period of 20 characteristic time constants, showing the approach to steady state at about 15 time constants. The program was modified so that the input beam focus was at the entrance to the crystal to make a more realistic model.

4. Discussion

The split-step beam propagation model has already been successfully applied to many photorefractive phenomena, starting in the early 1990s. These include image amplification, photorefractive oscillators, self-pumped phase conjugate mirrors, photorefractive surface waves and optical solitons. To the best of my knowledge, however, there are at least four gaps in knowledge about the method:

• The validity of the method and its computational limitations. Many of the prior results concentrated on analyzing the fields within the crystal, not the far-field output. It is in the far field that deficiencies of the model are most apparent through the appearance of scattering ring computational artifacts. In this paper, we gave some guidelines for the mitigation of these artifacts, mainly by the judicious choice of propagation step sizes.
• An extension of the scalar model from two to three dimensions. This has been enabled by improvements in available computational power since the 1990s in terms of data storage, random access memory, the availability of multicore processors and GPUs for parallel processing.
• An extension of the model from scalar to vector fields and the ability to model propagation in birefringent crystals considering the full electrooptic tensor. This is a topic for future research.
• The inclusion of thermal effects. These become important when the intensity of the beams becomes so large that temperature-dependent refractive index change is important. This effect can be significant since many photorefractive crystals are ferroelectric crystals near their Curie temperatures. This is a topic for future research.

This paper addressed the first two of these gaps and gave guidelines for proper usage. We found that propagation steps shorter than a few (~4) wavelengths are usually sufficient to push the artifacts’ direction cosines out of propagation range, which is greater than unity. On the other hand, longer steps can be used to reduce computation time for the interaction of beams with slow transverse variation. An example is the two-beam coupling of 100 μm waist gaussian beams.

The code presented here is unidirectional. A topic for future research would be to extend it to handle bidirectional beam propagation for modeling phase conjugate mirrors for example [16]. This extension would be straightforward if restricted to transmission gratings. The computation time and storage requirements would be significantly increased: one would have to either (a) store the full 3D grating and use a relaxation method on alternate forward and backward propagation to reach a solution or (b) store the full 3D optical intensity and use a relaxation method to converge to a solution for the full 3D space charge field. Reflection gratings are in principle treatable but would require a bidirectional propagation kernel.

The application that originally motivated this study, the use of photorefractive fanning mode mixing for deep learning, was not successful. While a photorefractive fanning layer followed by a single fully connected layer could indeed be used to classify audio and regular MNIST data sets, the results were no better than those previously obtained using linear mode mixing in multimode fibers [3,4]. I hope, however, that the extension of the scalar model to three dimensions, the discussion of its capabilities and limitations, and the associated open-source Python software will be of use to researchers in the future.

5. Conclusions

This paper is an elaboration of a paper published by my group at Tufts University of a “whole beam” method for photorefractive nonlinear optical beam propagation [15]. That paper was motivated by a need for a photorefractive theory that went beyond the standard slowly varying envelope coupled wave theory to address spatial effects on image bearing and multiwavelength beams and was subsequently used in a two-dimensional study of fanning in my group [34]. In this paper, we presented a realization of the scalar split-step method in three dimensions and used it to revisit several prior results on fanning, image amplification, high-gain effects, solitons and the time dependence of the processes.