Three-Dimensional Scalar Time Dependent Photorefractive Beam Propagation Model

1. Introduction

In the latter part of the 20^th 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). 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 21^st century it was beginning to exceed the capabilities of optical systems. For example, optical compact disks and DVDs have now become practically obsolete, replaced by high-capacity solids state drives and cloud storage. There are, though, still applications where optics remains the best option. One example 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[1]. Ancora et al. have shown that pattern classification can be accomplished using only a multimode fiber coupled with a single logistic regression layer[2]

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[3]. Tegin et al. [4] have shown that nonlinear effects in multimode fibers can be used to design optical neural networks whose power requirements can be lower than that of conventional 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[5], 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 stimulated photorefractive scattering, known as fanning[6]. The results presented here show that fanning can be well modelled 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[7] or with resonant materials such as sodium vapor[8]. These included distortion compensating phase conjugate mirrors with reflectivity exceeding unity[9], photorefractive oscillators[10], and self-pumped phase conjugate mirrors[10,11].

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[6,12,13,14,15]. This amplified scattering is commonly known as fanning because of the 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 treat

• Image amplification with comparison to experimental data published by Fainman et al.[16]
• Beam fanning with comparison to prior theory by Zozulya et al[13,17] and new experiments.
• Time dependence with comparison to high gain beam amplification experiments by Fischer et al[18].
• Theoretical appearance of high order diffraction at high gain described by Brown and Valley[14]

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 angular spectrum of plane waves[19]. Linear absorption is neglected in the present realization but could be added easily. The photorefractive nonlinearity accumulated during the linear propagation step is then applied to the optical field as a thin transparency. 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 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 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, in that it results in the appearance of rings centered on the direction of the incident beams 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 5mm 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 unnecessarily introduce additional complexity. We have chosen to continue to use FFT methods with the understanding that wraparound effects need to be acknowledged 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 try to keep the nonlinear propagation model as simple as possible while maintaining the principle physical behaviors we have chosen in the following to formulate the model using scalar diffraction. Extension to the full vector case including use of the dielectric tensor for beam propagation and electro-optic tensor for nonlinear grating formation would lead to a significant increase in complexity and we leave that generalization to future research.

• The vector angular spectrum of plane waves method requires separate handling of the ordinary and extraordinary vector modes which are both dependent on the transverse spatial frequencies[20].
• 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 component 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.
• Nevertheless, the scalar quasi-3D method described below does match many experimental results.

2.1. Beam Propagation

The model describes the optical field in Cartesian coordinates with x and y transverse z longitudinal. The input field is given in terms of its scalar amplitude A(x,y) at the input plane z = 0. Propagation of the field amplitude A(x,y) by distance dz is accomplished using the angular spectrum of plane waves method[19]. 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)$$

(1)

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

(2)

where f₀ = n/λ, λ 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.

An inverse Fourier transformed is then performed to return the field to real space. The Fresnel approximation can be recovered by expanding the square root in the exponent to its first order Taylor expansion, 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.

2.2. Photorefractive Nonlinearity Model

The hopping model for the photorefractive effect, as originally published by Feinberg et al., 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[21]. It can be rewritten in terms of continuous variables if it is assumed that only nearest neighbor hopping is important, that trapping sites can be treated as being uniformly spaced, and 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 a time dependent differential equation for the space charge 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. Primes represent spatial derivatives in the transverse direction. Full details may be found in Appendix B.

$${\displaystyle \frac{\partial \mathit{E}}{\partial \mathit{t}}}=-\left(\mathit{E}\mathit{I}-{\mathit{I}}^{\mathit{\text{'}}}\right)\left(1+{\mathit{E}}^{\mathit{\text{'}}}\right)+{\mathit{E}}^{\mathit{\text{'}}\mathit{\text{'}}}\mathit{I}$$

(3)

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 [22]:

$$\begin{array}{c}{\displaystyle \frac{\mathit{d}{\mathit{A}}_1}{\mathit{d}\mathit{z}}}=-{\mathit{\gamma }}_{\mathit{p}}{\displaystyle \frac{{\mathit{A}}_1{\mathit{A}}_2^{\mathit{*}}}{{\mathit{I}}_0}}{\mathit{A}}_2\\ {\displaystyle \frac{\mathit{d}{\mathit{A}}_2^{\mathit{*}}}{\mathit{d}\mathit{z}}}=+{\mathit{\gamma }}_{\mathit{p}}{\displaystyle \frac{{\mathit{A}}_1{\mathit{A}}_2^{\mathit{*}}}{{\mathit{I}}_0}}{\mathit{A}}_1^{\mathit{*}}\end{array}$$

(4)

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:

$${\mathit{\gamma }}_{\mathit{p}}={\displaystyle \frac{\mathit{\pi }{{\mathit{n}}^3\mathit{r}}_{\mathit{e}\mathit{f}\mathit{f}}{\mathit{E}}_0}{2\mathit{\lambda }\mathbf{c}\mathbf{o}\mathbf{s}\mathit{\theta }}}{\displaystyle \frac{{\mathit{k}}_{\mathit{g}}}{1+{\mathit{k}}_{\mathit{g}}^2}}={\displaystyle \frac{\mathit{\gamma }}{\mathbf{c}\mathbf{o}\mathbf{s}\mathit{\theta }}}{\displaystyle \frac{2{\mathit{k}}_{\mathit{g}}}{1+{\mathit{k}}_{\mathit{g}}^2}}$$

(5)

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, both described in Appendix B. r_eff is the effective electrooptic coefficient[23]. γ is the part of the coupling constant that is independent of propagation angles. The analytical solutions of these equations[22,24] show that the output beam ratio r_out in terms of the input beam ratio r_in is given by:

$${\mathit{r}}_{\mathit{o}\mathit{u}\mathit{t}}={\mathit{r}}_{\mathit{i}\mathit{n}}\mathit{e}\mathit{x}\mathit{p}\left(2{\mathit{\gamma }}_{\mathit{p}}\mathcal{l}\right)$$

(6)

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

Equivalently, the coupling constant length product (gain) can be written as

$${\mathit{\gamma }}_{\mathit{p}}\mathcal{l}={\displaystyle \frac{1}{2}}\mathbf{l}\mathbf{n}\left({\displaystyle \frac{{\mathit{r}}_{\mathit{o}\mathit{u}\mathit{t}}}{{\mathit{r}}_{\mathit{i}\mathit{n}}}}\right)$$

(7)

Figure 1 shows the gain calculated using equation (7) from the output of the program when run with input gain $\gamma \mathcal{l}$ = 3 using both the generalized refractive index grating (equation(3)) and the plane wave sinusoidal refractive index grating (equation (B11)). 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 in our code 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 Tukey window is not used. Allowable symmetric beam input angles are given by

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

(8)

where the integer n is subject to $1\le \mathit{n}\le {\log }_2({\mathit{N}}_{\mathit{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[14]. 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 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 experiments in photorefractive barium titanate performed by Fainman et al. and treated by them using standard coupled wave theory[16]. They presented their results in terms of absolute signal amplification G₀:

$${\mathit{G}}_0={\displaystyle \frac{\left(1+\mathit{r}\right)\mathbf{e}\mathbf{x}\mathbf{p}\left(2\mathsf{\gamma }\mathcal{l}\right)}{1+\mathit{r}\mathbf{e}\mathbf{x}\mathbf{p}\left(2\mathsf{\gamma }\mathcal{l}\right)}}$$

(9)

where r = r_in is the ratio of the intensity of the signal beam to that of the pump beam. Figure 2 shows theoretical plots of normalized gain ${\mathit{G}}_0/{\mathit{G}}_{0\mathit{s}\mathit{a}\mathit{t}}$ where the saturated, small signal gain is ${\mathit{G}}_{0\mathit{s}\mathit{a}\mathit{t}}=\mathbf{exp}\left(2\mathit{\gamma }\mathcal{l}\right)$.

Figure 3. Plots of normalized gain in dB versus input beam intensity ratio in the format used in reference [16]. Coupled wave equation results are shown in the following lines: Dashed line; G_0sat=100, solid line; G_0sat=1000, dotted line; G_0sat=10,000. The black points are as calculated by the 3D program, and the red points are data from the experiment described in the reference. The results are so close that the data points sometimes overlap.

Figure 3. Plots of normalized gain in dB versus input beam intensity ratio in the format used in reference [16]. Coupled wave equation results are shown in the following lines: Dashed line; G_0sat=100, solid line; G_0sat=1000, dotted line; G_0sat=10,000. The black points are as calculated by the 3D program, and the red points are data from the experiment described in the reference. The results are so close that the data points sometimes overlap.

Reference[16] treats image amplification optimization as its main result and we have reproduced some of the fidelity tests of that paper. While the group labels on the resolution targets are different, the computer input image has been scaled to the same size as the section of the chart shown in the amplified image from ref[16]. 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 Tukey apodization at the grid edges.

Reference[16] 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 by our model as shown in Figure 6 and Figure 7.

3.3. Photorefractive Amplified Scattering

The model allows for the three-dimensional simulation of photorefractive amplified scattering (fanning) seeded by optical scatterers[25]. This fanning serves as the seed for various photorefractive oscillators such as unidirectional ring oscillators and self-pumped phase conjugate mirrors[23] 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[13,17]. This model used Crank-Nicholson integration of the coupled wave equations and a relaxation method for 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[17] 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^-1(σ/λ). 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[17] and Figure 9 shows the corresponding results from our program. The relevant parameters are given in Table 1

To provide a comparison of the transverse output of the code with experiment, we show in Figure 10 an image of fanning of a 10mW 488nm beam from a Toptica iBeam smart diode laser after passage through a 6.5mm long crystal of barium titanate in our laboratory. The corresponding output of the model using beam waist 3mm 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{gathered} \cos {\theta }_x=f_x\lambda \\ \cos {\theta }_y=f_y\lambda \end{gathered}$$

(10)

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 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 modelled in the undepleted pump approximation by Montemezzani et al.[26]. The development of these lobes can be seen in 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 modelling 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 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 numerical Euler integration of Equation Error! Reference source not found.

$$\mathit{E}\left({\mathit{t}}_{\mathit{n}}\right)=\mathit{E}\left({\mathit{t}}_{\mathit{n}-1}\right)+\mathsf{\Delta }\mathit{t}{\displaystyle \frac{\partial \mathit{E}}{\partial \mathit{t}}}\left({\mathit{t}}_{\mathit{n}-1}\right)$$

(11)

In all calculations that follow we use normalized time and normalized distance: see Appendix B).

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.[18], 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[27]:

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

(12)

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.

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

In the context of modelling 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 (B11)). The stability criterion of Euler integration for simple exponential decay is |Δt/τ-1| < 1 [28] so we can attempt to ensure 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 (B11) we choose the step size conservatively to be

$${4\mathsf{\Delta }\mathit{t}={\mathit{\tau }}_{\mathit{m}\mathit{i}\mathit{n}}={\left(1+{\left({\displaystyle \frac{{\mathit{k}}_{\mathit{g}\mathit{m}\mathit{a}\mathit{x}}}{{\mathit{k}}_0^2}}\right)}^2\right)}^{-1}=\left(1+{\left({\displaystyle \frac{\mathit{\pi }{\mathit{N}}_{\mathit{x}}}{\mathit{L}{\mathit{k}}_0^2}}\right)}^2\right)}^{-1}$$

(13)

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

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 has been 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 supplemental materials (Move 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.

4. Discussion

The split step beam propagation model has already been successfully applied to many photorefractive phenomena starting in the early 1990’s. 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 had been at least four gaps in knowledge about the method:

• 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 have given some guidelines for the mitigation of these artifacts, mainly by the judicious choice of propagation step sizes.
• Extension of the scalar model from two to three dimensions. This has been enabled by improvements in available computational power since the 1990’s, in terms of data storage, random access memory, the availability of multicore processors and GPUs for parallel processing.
• Extension of the model from scalar to vector fields and the ability to model propagation in birefringent crystals considering the full electrooptic tensor.
• Inclusion of thermal effects. These become important when the intensity of the beams becomes so large that temperature dependent refractive index changes important. This effect can be significant since many photorefractive crystals are ferroelectric crystals near their Curie temperatures.

This paper has addressed the first two of these gaps and has given 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, that 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 two beam coupling of 100μm waist gaussian beams.

The original 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[1,2]. I hope, however, that the extension of the scalar model to three-dimensions, 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[12]. 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 [29]. We have presented here a realization of the scalar split step method in three dimensions and have used it to revisit several prior results on fanning, image amplification, high gain effects, and the time dependence of the processes.