## Abstract

Diffraction of a 3D optical beam on a multilayer phase-shifted Bragg grating (PSBG) is considered. It is shown that the PSBG enables optical computation of the spatial Laplace operator of the electromagnetic field components of the incident beam. The computation of the Laplacian is performed in reflection at normal incidence. As a special case, the parameters of the PSBG transforming the incident Gaussian beam into a Laguerre–Gaussian mode of order (1,0) are obtained. Presented numerical results demonstrate high quality of the Laplace operator computation and confirm the possibility of the formation of Laguerre–Gaussian mode. We expect the proposed applications to be useful for all-optical data processing.

© 2014 Optical Society of America

## 1. Introduction

Optical devices that perform temporal and spatiotemporal transformations of optical signals are of great interest for a wide range of applications including ultrafast all-optical information processing and analog optical computations [1–12]. Among the most important operations of the analog processing of optical signals are the operations of temporal and spatial differentiation. The operation of temporal (spatial) differentiation assumes the differentiation of the envelope of an optical pulse (the spatial profile of an optical beam). Various types of Bragg gratings [2–6] and resonant diffraction gratings [7,8] were proposed for the temporal differentiation of optical pulses. Spatial differentiation of 2D optical beams was proposed in the authors’ recent work [9]. In [9], it was shown that a phase-shifted Bragg grating (PSBG) can be used for the optical computation of the first derivative of the 2D incident beam profile at oblique incidence. At normal incidence, the structure computes the second derivative of the beam profile. At the same time, spatial transformations of 3D optical beams were not considered in [9].

In the present work, the possibility of the optical computation of the Laplace operator in spatial coordinates using a PSBG is theoretically and numerically demonstrated for the first time. It should be noted that sophisticated techniques are currently used to implement this operation which are based on the subtraction of two holograms recorded at different distances or with different wavelengths [13, 14]. As an important practical application of PSBG, the transformation of a 3D Gaussian beam into a Laguerre–Gaussian beam is considered. The proposed applications of PSBG are useful for various all-optical image processing applications including edge detection, image sharpening, pattern recognition and motion estimation.

## 2. Diffraction of a 3D optical beam by a multilayer system

Consider a 3D beam normally incident upon a multilayer system (Fig. 1). Let us represent the incident beam through its angular spectrum as a superposition of TE- and TM-polarized plane waves with different propagation directions. Since the *z* axis is perpendicular to the considered multilayer system, it is contained in the planes of incidence of the plane waves. Consequently, the electric field vector of the TE-polarized waves is perpendicular to the *z* axis (*E _{z}* ≡ 0). Similarly, the TM-polarized waves have the magnetic field vector perpendicular to the

*z*axis (

*H*≡ 0). Assuming that the time dependence has the form exp(−i

_{z}*ωt*), the following representation of the TE- and TM-polarized plane wave field can be obtained from Maxwell’s equations [15]:

**r**= (

*x*,

*y*,

*z*),

**Φ**(

*r*) = [

**E**(

**r**)

**H**(

**r**)]

^{T}is a column vector of the electric (

**E**= [

*E*

_{x}*E*

_{y}*E*]) and the magnetic (

_{z}**H**= [

*H*

_{x}*H*

_{y}*H*]) field components,

_{z}**p**= (

*α*,

*β*,

*γ*) is the wave vector, where $\gamma =\pm \sqrt{{k}_{0}^{2}{n}_{\text{sup}}^{2}-({\alpha}^{2}+{\beta}^{2})}$,

*k*

_{0}= 2

*π/λ*is the wave number,

*λ*=

*ω/c*is the wavelength,

*n*

_{sup}is the refractive index of the medium,

**Φ**

_{TE}(

*α*,

*β*) and

**Φ**

_{TM}(

*α*,

*β*) are the vectors of the following form:

*γ*component of the wave vector define the waves propagating in the positive and negative directions of the

*z*axis, respectively. Equations (2) are not defined at

*α*=

*β*= 0. We assume that

**Φ**

_{TE,TM}(0, 0) = lim

_{α→0}

**Φ**

_{TE,TM}(

*α*, 0). With this definition, the vectors

**Φ**

_{TE}(0, 0) and

**Φ**

_{TM}(0, 0) correspond to the plane waves with the electric field vector directed along the

*y*axis and the

*x*axis, respectively.

In the general case, the incident light beam can be represented as a sum of TE- and TM-polarized beams, which in turn correspond to the superpositions of TE- and TM-polarized plane waves propagating in the negative direction of the *z* axis:

**A**

_{inc}(

**r**) = [

**E**

_{inc}(

**r**)

**H**

_{inc}(

**r**)]

^{T}is a column vector containing the field components of the incident beam,

*G*

_{TE}(

*α*,

*β*) and

*G*

_{TM}(

*α*,

*β*) are the spatial frequency spectra of the TE- and TM-polarized beams, and

**G**

_{inc}(

*α*,

*β*) is the spatial frequency spectrum of the incident field components defined by the following expression:

*G*

_{TE}(

*α*,

*β*) and

*G*

_{TM}(

*α*,

*β*) properly, two arbitrary components of the electromagnetic field of the incident beam in the plane

*z*= 0 can be specified. As an example, let us consider a beam polarized along the

*x*axis [16, 17]. In this case,

*G*(

_{x}*α*,

*β*) is the spatial frequency spectrum of the

*x*-component. Using the expressions for the

*E*and

_{x}*E*field components in Eqs. (2), (3), (5), one can obtain the following system of linear equations:

_{y}*G*

_{TE}and

*G*

_{TM}allows us to obtain the coefficients for the expansion of the linearly polarized beam in TE- and TM-polarized waves:

**Φ**

_{TE,TM}vectors at

*α*=

*β*= 0, we obtain

*G*

_{TE}(0, 0) = 0, ${G}_{\text{TM}}(0,0)=-{G}_{x}(0,0)/({k}_{0}\sqrt{\epsilon})$.

For a Gaussian beam, *E*_{x,inc}(*x*, *y*,0) = exp{−(*x*^{2} + *y*^{2})/(2*σ*^{2})}, and the function *G _{x}*(

*α*,

*β*) in Eqs. (5) and (7) takes the following form:

We will henceforth assume that the beam incident on a multilayer system is defined by Eqs. (3)–(7). In this case, the reflected beam has the form

*R*

_{TE}(

*α*,

*β*),

*R*

_{TM}(

*α*,

*β*) are the reflection coefficients of the TE- and TM-polarized plane waves. Let us note that for the studied multilayer structure the reflection coefficients depend only on the angle between the surface normal and the wave vector, so they can be considered as the functions of the in-plane wave vector component:

*R*

_{TE,TM}(

*ξ*) =

*R*

_{TE,TM}(

*ξ*, 0)]. Equation (9) implies that the TE- and TM-components

**A**

_{refl,TE}(

*x*,

*y*, 0) and

**A**

_{refl,TM}(

*x*,

*y*, 0) of the reflected field in the plane

*z*= 0 correspond to the results of transformation of the incident field components

**A**

_{inc,TE}(

*x*,

*y*, 0),

**A**

_{inc,TM}(

*x*,

*y*, 0) by linear time-invariant systems with transfer functions (TF) ${H}_{\text{TE}}(\alpha ,\beta )={R}_{\text{TE}}\left(\sqrt{{\alpha}^{2}+{\beta}^{2}}\right)$ and ${H}_{\text{TM}}(\alpha ,\beta )={R}_{\text{TM}}\left(\sqrt{{\alpha}^{2}+{\beta}^{2}}\right)$. Note that according to Eq. (10),

*R*(−

*α*) =

*R*(−

*α*, 0) =

*R*(

*α*, 0) =

*R*(

*α*). This means that the reflection coefficients

*R*

_{TE,TM}(

*ξ*) are even functions of the argument

*ξ*. As a consequence, the functions

*H*

_{TE}(

*α*,

*β*),

*H*

_{TM}(

*α*,

*β*) are quadratic in the first approximation.

## 3. Computation of the Laplace operator

Let us assume that the angular (spatial frequency) spectrum of the incident beam is nonzero at $\sqrt{{\alpha}^{2}+{\beta}^{2}}<\delta $. It follows from Eqs. (3), (9) that if

*R*

_{TE,TM}(0) = 0) is necessary for the computation of the Laplace operator. In this case, the following relation holds up to the second order:

*H*

_{TE, TM}(

*α*,

*β*) ∼

*R″*

_{TE,TM}(0)

*H*

_{Δ}(

*α*,

*β*). Thus, the magnitude of the second derivative

*R″*

_{TE, TM}(0) determines the energy efficiency of the computation of the Laplacian. To provide a large value of the second derivative

*R″*

_{TE,TM}(0), rapid variation of the reflection coefficient in the vicinity of normal incidence is necessary. Such a reflection spectrum is typical for PSBGSs [5, 6, 12] which are widely used as narrow-band spectral filters.

PSBG provides zero reflectance at the incidence angle *θ*_{0} (at *α* = *α*_{0} = *k*_{0}*n*_{sup} sin*θ*_{0}) simultaneously for TE- and TM-polarizations of the incident wave. PSBG consists of two symmetric Bragg gratings (BGs) separated by a phase-shift (or defect) layer (Fig. 1). In the simplest case, the BG layers have equal optical thickness:

*i*= 1, 2;

*n*,

_{i}*h*are the refractive indices and layer thicknesses of the BG,

_{i}*λ*

_{B}is the Bragg wavelength. If the defect layer has the optical thickness

*ñh*

_{def}=

*λ*

_{B}/2, where $\tilde{n}=\sqrt{{n}^{2}-{n}_{\text{sup}}^{2}{\text{sin}}^{2}{\theta}_{0}}$,

*n*being the refractive index of the defect layer, the BG reflection coefficient vanishes at wavelength

*λ*

_{B}and angle of incidence

*θ*

_{0}[6]. Let us note that this reflectance zero is located at the center of the first photonic band gap of the BG. The appearance of a zero in the reflectance spectrum has resonant nature and is associated with the excitation of the eigenmodes localized in the defect layer.

A resonant representation of the reflection coefficient as a function of the spatial frequency can be written [7, 10]. According to Eq. (10), the reflection coefficient of the studied multilayer structure can be considered as a single-argument function of the *α* component of the wave vector (*β* ≡ 0) at fixed polarization (TE or TM). The reflection (or transmission) coefficient of a multilayer structure having several resonances can be represented as a sum of rational functions [10]:

*r*is the non-resonant reflection coefficient,

*M*is the number of eigenmodes supported by the structure,

*b*is the coupling coefficient of the mode [18],

_{j}*p*is the complex propagation constants of the eigenmodes of the multilayer structure. The value of Im

_{j}*p*characterizes the mode extinction coefficient. For the PSBG this quantity depends essentially on the number of BG periods: a greater number of periods leads to lower scattering of the mode, which corresponds to the smaller value of Im

_{j}*p*. The quantities

_{j}*r*,

*b*and

_{j}*p*depend on the polarization of the incident wave.

_{j}Due to the fact that *R*(*α*) is an even function, the *R* function can be approximately represented as a rational function with respect to *α*^{2}. The simplest case assumes the presence of two eigen-modes with the propagation constants ±*p*. Moreover, given the condition of zero reflectance at *α* = *α*_{0} = *k*_{0}*n*_{sup} sin*θ*_{0}, the reflection coefficient can be written as:

*α*

_{0}= 0 allows us to write the TFs

*H*

_{TE,TM}(

*α*,

*β*) in the form

*H*

_{Δ}(

*α*,

*β*) providing the computation of the Laplace operator and a smoothing filter

*H*

_{err}(

*α*,

*β*). Thus, if the function

*H*

_{err}(

*α*,

*β*) is slowly varying in the neighborhood of the point

*α*=

*β*= 0, the PSBG will perform the optical computation of the Laplacian of TE- and TM-components of the incident beam.

## 4. Simulation results and discussion

As an example, let us consider a PSBG consisting of two BGs having *N* periods separated by a defect layer (total number of the layers is 4*N* + 1). We use the following values of the refractive indices of the BG layers, defect layer, superstrate and substrate: *n*_{1} = 1.5, *n*_{2} = 2.25, *n* = 1.5, *n*_{sub} = *n*_{sup} = 1. Layer thicknesses were calculated from Eq. (13) at *λ*_{B} = 1500 nm and *α*_{0} = 0. To simulate the beam diffraction we used the following rigorous approach. First, we solved the plane wave diffraction problem using the scattering matrix method [19]. Then, we calculated the reflected field distribution by evaluating the integrals (9) numerically with the use of the Fast Fourier Transform.

Figure 2 shows the TFs of the PSBGs consisting of 13, 17, and 21 layers. According to Fig. 2, the TFs are approximately quadratic in the vicinity of *α* = 0 (dashed lines show quadratic TF providing exact computation of the Laplace operator). Let us note that the increase in the layer quantity (number of the BG periods *N*) leads to the decrease in Im*p*_{TE,TM}, which results in the decrease in the dip width of the reflection coefficient and, consequently, to the decrease in the *α/k*_{0} range where the TF can be considered quadratic. It is easy to show that the decrease in the dip width corresponds to the decrease in
$\left|{p}_{\text{TE},\text{TM}}^{2}\right|$ and, consequently, the increase in the error function *H*_{err}(*α*, *β*) values in Eq. (16). At the same time, the decrease in |*p*_{TE,TM}| results in the increase in the
${R}_{\text{TE},\text{TM}}^{\u2033}(0)~1/{p}_{\text{TE},\text{TM}}^{2}$ value, which defines the amplitude of the reflected beam (or energy efficiency). Thus, the number of layers governs the trade-off between energy efficiency and differentiation quality.

According to Eq. (16), the Laplace operator of the TE- and TM-components of the incident beam is in general case computed with different weights
${w}_{\text{TE},\text{TM}}={r}_{\text{TE},\text{TM}}/{p}_{\text{TE},\text{TM}}^{2}$. At the same time, Fig. 2 shows that for the considered example the reflection coefficients for TE- and TM-polarizations are close in magnitude. Provided that *R*_{TE}(*α*, *β*) ≈ *R*_{TM}(*α*, *β*), the reflected beam [Eq. (9)] takes the form

**G**

_{inc}(

*α*,

*β*) is the spatial frequency spectrum of the incident field defined by Eq. (4). According to Eq. (17), in this case the PSBG performs the optical computation of the Laplace operator of all components of the incident field.

Let us study the performance of the 13-layer PSBG for a Gaussian incident beam defined by Eqs. (5)–(8). Figure 3 shows the transverse distributions of the *x*-component magnitudes of the electric field of the incident beam *E*_{x,inc}(*x*, *y*,0) = exp {−(*x*^{2} + *y*^{2})/(2*σ*^{2})} (the half-width at the 1/e^{2} amplitude level amounts to 2*σ* = 10*μ*m) [Fig. 3(a)] and of the reflected beam [Fig. 3(b)]. Figure 3 also shows the field profiles of the incident [*E*_{x,inc}(*x*, 0, 0)] and the reflected [|*E*_{x,refl}(*x*, 0, 0)|] beams. The field profile of the reflected beam is shown together with the Laplacian of the incident beam |*χ*Δ*E*_{x,inc}(*x*, *y*, 0)|_{y=0}, where *χ* = Δ*E*_{x,inc}(0, 0, 0)/*E*_{x,refl}(0, 0, 0) is the scaling factor. Comparison of the profiles indicates that the reflected beam corresponds to the Laplacian of the incident beam with high accuracy. Indeed, the normalized half-width of the incident beam spectrum amounts to *δ/k*_{0} = 2/(*σk*_{0}) = 0.095. According to Fig. 2, at |*α/k*_{0}| ≤ 0.095 the amplitudes of the TFs *H*_{TE,TM}(*α*, 0) almost coincide (and are nearly quadratic) for the considered 13-layer PSBG, and the TF phase differences do not exceed *π*/40.

It is worth noting that the transverse field distribution of the beam corresponding to the Laplacian of the Gaussian function is not conserved during the propagation. A similar beam with conserved profile is the Laguerre–Gaussian mode of order (1,0). The transverse field distribution of this mode at the beam waist has the following form [20]:

*L*

_{1}(

*x*) = 1 −

*x*is the Laguerre polynomial. This mode can be represented as a superposition of the incident Gaussian beam and its Laplacian:

*H*(

*α*,

*β*) =

*σ*

^{−2}+

*H*

_{Δ}(

*α*,

*β*). It follows from Eqs. (15) and (19) that the PSBG designed for

*α*

_{0}=

*σ*

^{−1}has a similar TF:

*H*(

*α*,

*β*) ∼ (

*σ*

^{−2}+

*H*

_{Δ})

*H*

_{err}(

*α*,

*β*). Hence, the PSBG allows to perform the aforesaid mode transformation. Note that in this case, the reflection coefficient of the PSBG vanishes at the angle of incidence

*θ*

_{0}= asin[1/(

*k*

_{0}

*n*

_{sup}

*σ*)]. Figure 4 shows the reflected beam generated by the PSBG calculated for ${\alpha}_{0}^{2}={\sigma}^{-2}=0.04\mu {\text{m}}^{-2}$. The incident beam is the same as in the previously considered example [Fig. 3(a)]. The reflected beam profile matches the profile of the Laguerre–Gaussian mode

*E*

_{x,LG}(

*x*, 0) with high accuracy. Figure 4(b) shows the profile of the reflected beam at distances

*z*= 0,

*z*= 50

*μ*m and

*z*

_{refl}= 100

*μ*m. It is evident that the beam profile is conserved up to scale during propagation.

## 5. Conclusion

In the present work, we have shown that the PSBG enables optical computation of the Laplace operator in spatial coordinates of the electromagnetic field components of the incident beam. The computation of the Laplacian is performed in reflection at normal incidence. The parameters of the PSBG transforming the incident Gaussian beam into a Laguerre–Gaussian mode of order (1, 0) were also derived analytically. Presented numerical results demonstrate high quality of the Laplace operator computation and confirm the possibility of the formation of Laguerre–Gaussian mode. We believe that the proposed applications of PSBG could be useful in real-time all-optical image sharpening, edge detection and pattern recognition.

## Acknowledgments

The work was funded by the Russian Science Foundation grant 14-19-00796.

## References and links

**1. **A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science **343**(6167), 160–163 (2014). [CrossRef] [PubMed]

**2. **R. Slavík, Y. Park, M. Kulishov, and J. Azaña, “Terahertz-bandwidth high-order temporal differentiators based on phase-shifted long-period fiber gratings,” Opt. Lett. **34**(20), 3116–3118 (2009). [CrossRef] [PubMed]

**3. **L. M. Rivas, S. Boudreau, Y. Park, R. Slavík, S. LaRochelle, A. Carballar, and J. Azaña, “Experimental demonstration of ultrafast all-fiber high-order photonic temporal differentiators,” Opt. Lett. **34**(12), 1792–1794 (2009). [CrossRef] [PubMed]

**4. **M.A. Preciado, X. Shu, P. Harper, and K. Sugden, “Experimental demonstration of an optical differentiator based on a fiber Bragg grating in transmission,” Opt. Lett. **38**(6), 917–919 (2013). [CrossRef] [PubMed]

**5. **N. K. Berger, B. Levit, B. Fischer, M. Kulishov, D. V. Plant, and J. Azaña, “Temporal differentiation of optical signals using a phase-shifted fiber Bragg grating,” Opt. Lett. **15**(2), 371–381 (2007).

**6. **M. Kulishov and J. Azaña, “Design of high-order all-optical temporal differentiators based on multiple-phase-shifted fiber Bragg gratings,” Opt. Express **15**(10), 6152–6166 (2007). [CrossRef] [PubMed]

**7. **D. A. Bykov, L. L. Doskolovich, and V. A. Soifer, “Single-resonance diffraction gratings for time-domain pulse transformations: integration of optical signals,” J. Opt. Soc. Am. A **29**(8), 1734–1740 (2012). [CrossRef]

**8. **D. A. Bykov, L. L. Doskolovich, and V. A. Soifer, “Temporal differentiation of optical signals using resonant gratings,” Opt. Lett. **36**(17), 3509–3511 (2011). [CrossRef] [PubMed]

**9. **L. L. Doskolovich, D. A. Bykov, E. A. Bezus, and V. A. Soifer, “Spatial differentiation of optical beams using phase-shifted Bragg grating,” Opt. Lett. **39**(5), 1278–1281 (2014). [CrossRef] [PubMed]

**10. **V. Lomakin and E. Michielssen, “Beam transmission through periodic subwavelength hole structures,” IEEE Trans. Antennas Propag. **55**(6), 1564–1581 (2007). [CrossRef]

**11. **G. D. Landry and T. A. Maldonado, “Gaussian beam transmission and reflection from a general anisotropic multilayer structure,” Appl. Opt. **35**(30), 5870–5879 (1996). [CrossRef] [PubMed]

**12. **N. Q. Ngo, “Design of an optical temporal integrator based on a phase-shifted fiber Bragg grating in transmission,” Opt. Lett. **32**(20), 3020–3022 (2007). [CrossRef]

**13. **C.-S. Guo, Q.-Y. Yue, G.-X. Wei, L.-L. Lu, and S.-J. Yue, “Laplacian differential reconstruction of in-line holograms recorded at two different distances,” Opt. Lett. **33**(17), 1945–1947 (2008). [CrossRef] [PubMed]

**14. **J. P. Ryle, D. Li, and J. T. Sheridan, “Dual wavelength digital holographic Laplacian reconstruction,” Opt. Lett. **35**(18), 3018–3020 (2010). [CrossRef] [PubMed]

**15. **V. A. Soifer, *Diffractive Nanophotonics* (CRC Press, 2014). [CrossRef]

**16. **S. M. Sepke and D. P. Umstadter, “Exact analytical solution for the vector electromagnetic field of Gaussian, flattened Gaussian, and annular Gaussian laser modes,” Opt. Lett. **31**(10), 1447–1449 (2006). [CrossRef] [PubMed]

**17. **G. Zhou, “The analytical vectorial structure of a nonparaxial Gaussian beam close to the source,” Opt. Express **16**(6), 3504–3514 (2008). [CrossRef] [PubMed]

**18. **S. Fan, W. Suh, and J. D. Joannopoulos, “Temporal coupled-mode theory for the Fano resonance in optical resonators,” J. Opt. Soc. Am. A **20**(3), 569–572 (2003). [CrossRef]

**19. **L. Li, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. A **13**(5), 1024–1035 (1996). [CrossRef]

**20. **A. E. Siegman, *Lasers* (University Science Books, 1986).