## Abstract

Understanding and modeling the propagation of polarized light through thick, space variant birefringent media is important in both fundamental and applied optics. We present and experimentally evaluate two methods to model the off axis propagation of polarized light through a thick stress-engineered optic (SEO). First, we use a differential equation solving method, which utilizes the analytic expression for the Jones matrix of the SEO leading to a numerical solution for the output electric field. Then we present a geometric method to obtain similar results with much less computational complexity. Finally, a comparison is done between the data and the simulations.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Stress-related effects are ubiquitous in optical materials; while stress is often minimized through suitable optomechanics, there have been occasions when either static or dynamic stress is used to engineer useful optical properties that could be difficult to achieve otherwise. The early applications of this included both waveplates and modulators [1]. Stress engineering can be particularly useful for optical wavelengths for which no suitable (or suitably inexpensive) birefringent material is available. In recent years, papers have been published that report the use of spatially inhomogeneous stress for applications in pupil filtering to create unconventional polarization states such as cylindrical vector beams [2,3] and full Poincaré beams [4,5]. These symmetrically stressed elements (stress-engineered optics, or SEOs) imprint unique topological features (eg. phase vortices) on an incident laser beam. Some of these have also been shown to be an interesting option for single shot polarimetry [6–8].

Such elements are often optically thick, with complex spatial dependence to the birefringence albeit with a maximum birefringence that is generally orders of magnitude less than anisotropic materials such as Calcite. Polarization ray tracing has been shown to accurately model general cases, especially with large birefringence [9–13]. However, weakly birefringent media such as stressed materials are important but often don’t require rigorous ray tracing since there is a little deviation between the energy and the ray direction.

SEO models presented and analyzed in previous papers [4,5] have mostly considered a strict paraxial approximation; when valid, the paraxial approximation allows the entire window to be treated as a single space-variant Jones matrix that locally transforms the incident polarization as if it were a thin element. However, thick stressed elements show easily visible angular dependence to the polarization distribution in the transmitted light [2]. For continuously varying optically thick media, if the birefringence distribution is constant in the z direction but varying in a plane perpendicular to the face of the window, an off axis ray will experience a birefringence gradient as it passes through the medium.

In this paper, we present and experimentally evaluate two methods to model the off-axis propagation of polarized light through a thick stress-engineered optic and compare the results with an experimental test that makes use of an SEO.

## 2. Stress-engineered optics

The SEOs studied here are optical quality windows with symmetrically distributed peripheral forces. Fabrication of these SEOs are described elsewhere [2, 3] and can be summarized as follows. A BK7 optical window (diameter 12.7 *mm*, thickness 8 *mm*) is stressed using a thermal compression method where the glass window is inserted into a metal ring with inner diameter 25 *μm* smaller than the outer diameter of the window. Here the metal ring was heated to about 300°*C* and the higher expansion coefficient of metal allowed the window to be inserted into the ring. Material has been removed in the metal ring symmetrically so that it created three points of contact on the edge of the window once the ring is cooled. This process is illustrated in Fig. 1.

When used in an optical system, the SEO is generally apodized with a circular aperture. If used as a pupil filter, such as has been done in polarimetry, the aperture size is limited to less than half the physical size of the window [6].

## 3. Foundational theory

When an electric field, **E*** _{in}*, is incident on a retarding element, the output electric field,

**E**

*, is given by*

_{out}**E**and 𝕁 are referred to as the Jones vector of the electric field, and the Jones matrix of the retarder respectively. For a generic non-absorbing waveplate with

*δ*phase retardance with respect to a certain wavelength, and the waveplate orientation

*α*with respect to the horizontal axis, the Jones matrix can be written as

*α*) is the

*pseudo rotation matrix*, which can be written in linear and circular bases as follows (

*l*and

*c*subscripts correspond to linear and circular bases)

For SEOs with *m* number of rotationally symmetric force points with one of the force points aligned with the x axis the form of the Jones matrix near the center of the window takes the form

*m*being the number of stress points,

*ϕ*being the azimuthal coordinate, and

*ρ*being the normalized radial coordinate of the pupil. The unitless parameter

*c*(the stress parameter [2–4,6]) represents the magnitude of the retardance in radians at the edge of the pupil compared to the center. Although SEOs can be made with an arbitrary number of stress points (

*m*≥ 3), it is evident that when

*m*= 3 the retardance in the center region is linearly varying as a function of the radial coordinate. For this case the Jones matrix of the SEO in circular basis can be written as

When light with a certain polarization is incident perpendicularly on an SEO, it undergoes a transformation that gives rise to beams having all possible polarization states (Full Poincaré beams). These beams exhibit a variety of phase vortices in the scalar field components as well as polarization singularities [2, 4, 14, 15]. These vortices are similar to those used in Q-plates and related technology with the important difference that an SEO is a thick element whereas Q-plates are typically thin. Also, it has been shown that the point spread function arising from an imaging system in which the SEO is placed at a pupil plane is sensitive to the input polarization, enabling an efficient method of single-shot Stokes polarimetry [6,7]. So far all the modeling and experimentation have been focused on the perpendicular propagation of light through the SEO. It would be interesting to look at the off-axis effects when light incident on the SEO with an angle.

The geometry is as follows: When the light is incident on the SEO with angle *θ*, it refracts with refraction angle *θ′* given by Snell’s law (Fig. 2). A ray experiences a varying amount of birefringence as it travels obliquely inside the material, and exits with a transverse deviation Δ. In this work we use two methods to model the oblique propagation through an SEO: 1) Numerically solve two coupled differential equations for the two field components (continuous model); 2) Discretize the SEO into many thin space-variant waveplates (discrete waveplate model).

#### 3.1. Continuous model

The Jones matrix for a single layer of a general SEO of thickness *δz* and m-fold symmetry, as depicted in Fig. 2, is given by

*δ*=

*cρ*

^{m−2}/2 has been scaled by

*δz*/

*L*. Thus, the output electric field vector,

**E**(

*z*+

*δz*) can be written, as a product of

**E**(

*z*) and the corresponding Jones matrix of the thin slab, 𝕁

*, as where*

_{δ}_{z}*k*is the propagation constant. It should be noted that since 𝕁

*is space-variant,*

_{δ}_{z}**E**(

*z*) can be parameterized by the radial coordinate

*ρ*and the azimuthal angle

*ϕ*. In the following derivation, this parameterization is implied.

**E**(

*z*) can be written as

**e**(

*z*)

*exp*(

*ikz*), and by substituting this in equation 7, one gets

*δz*to be small, Eq. 10 can be written as a difference quotient:

Using the circular basis and taking the limit as *δz* → 0, one would obtain

*e*) would have an azimuthal phase term, implying a spin orbit coupling of light when passing through the SEO. In this paper, however, we will not explore this phenomenon.

_{R}For the light incident on the SEO with a non-zero angle, Eq. 13 was solved numerically. For the input electric field, two matrices were defined where each matrix being the corresponding element of the Jones vector of the initial polarization distribution. Each light ray would undergo a deviation determined by the incident angle and the SEO thickness. Hence, when numerically solving the two differential equations the birefringence that a particular light ray would “see” at a particular longitudinal distance z changes as a function of z. Thus, making *ρ* and *ϕ* functions of z. Assuming the shift is only in the x direction,

These expressions for *ϕ*(*z*) and *ρ*(*z*) were used in Eq. 13 with *m* = 3 to solve for the electric field **e**(*z* = *L*) using the ode45() function in Matlab, which utilizes an explicit Runge-Kutta(4,5) formula, with the input electric field as the initial condition. If *M* × *M* matrices are used, then it is required to solve *M*^{2} number of differential equation pairs (one per each pixel) for a particular incident angle, making this method computationally expensive for large *M*.

#### 3.2. Discrete retarder approximation

The SEO can be thought as a collection of thinner “slabs”,Fig. 3 where each slab has the birefringence of the SEO divided by the number of slabs, *N*. Then the Jones matrix 𝕁* _{n}* for such single thin slab with three stress points (i.e.

*m*= 3) becomes

Thus, for normal incidence, the output Jones vector of the electric field after passing through the *N* slabs is given by

To model the oblique propagation of light, the electric field is shifted after propagating through each slab using a linear phase applied to the Fourier transform. At each step, the amount of shift is equal to the total shift of the field divided by *N*.

For each method, the *x* and *y* components of the electric field are computed and used to obtain both a polarization map of the electric field and an irradiance profile following a suitably chosen analyzer. The images were compared on the basis of irradiance profiles, polarization maps, and contours of the relevant Stokes parameters.

## 4. Experimental procedure

In order to image the birefringence pattern of the SEO a setup as shown in Fig. 4 was used. Light from a doubled Nd:YAG laser (*λ* = 532 nm) was spatially filtered and directed to a rotating diffuser to reduce spatial coherence prior to illuminating the SEO. A 10 mm aperture following the diffuser limited the range of incident illumination angles to ±1.5°. The SEO was attached to a rotation stage with 1° precision. Using an afocal relay (4f) system the central region was imaged (0.5× magnification) to a CCD camera (Imaging Source DMK33G445) connected to a computer. A right-hand circular analyzer was placed before the camera to block one component of the Jones vector. The rotation stage was rotated from 0° to 24° in 2° steps and images were acquired using a Matlab script. This was repeated for different input polarizations H, V, P, M, RHC and LHC. A dark frame was taken at the end of the data acquisition by covering up the camera and subtracted from each of the images of the SEO.

Then a quarter wave plate was added following the SEO, and the circular analyzer was replaced with a vertical analyzer. Using input LHC polarization, for each tilt position of the SEO the quarter wave plate was rotated from 0° − 180° in 10° steps. These images were analyzed using another Matlab script to compute the Stokes parameters.

## 5. Results and discussion

Figure 5 compares the continuous model and discrete retarder model with experimental results, showing irradiance profiles in the SEO center region for three incident angles 0°, 10°, and 20°, assuming LHC input polarization. The top right image shows the experimental irradiance image for a maximum angle of 20° with the predicted contours of equal irradiance shown in red.

Figure 6 shows polarization ellipse maps and *S*_{3} (the third component of the Stokes vector) contours of the simulated outputs using the continuous model. At normal incidence, when going outward from the center (where there is zero birefringence and the polarization is unchanged), the ellipticity and the handedness changes periodically, resulting in a beam containing all fully polarized polarization states (the so-called Full Poincaré Beam). As the incident angle increases, the lines of circular polarization (sets of points in the pupil where *S*_{3}, becomes ±1, corresponding to the poles on the Poincaré sphere) transform from concentric circles to singular points. Contours of linear polarization morph from circles to the deformed (dashed) contours seen in figure 6.

As seen in Fig. 7, in the central region of the window, the experimental and theoretical contours are nearly indistinguishable, with some differences appearing near the edge of the region under test. The contours of equal S3 and of linear polarization shown side by side for three example angles of incidence (see Visualization 1 and Visualization 2).

We now address the convergence of the discrete retarder method. Simulated irradiance patterns from both methods were normalized by dividing each image by the sum of the pixel values of each image. Then the error, *E*, was computed as the root sum square of the difference *I _{CM}* −

*I*, where

_{DRM}*I*and

_{CM}*I*are the irradiance patterns from the continuous method and the discrete retarder method respectively.

_{DRM}Figure 8 depicts the behavior of this error as a function of the incident angle used in the simulations. It agrees with the general intuition that the error should increase with increasing angle; furthermore, for any particular angle the error decreases as the number of slabs increases. We found that the computation time taken to simulate an irradiance pattern using the differential equation method was approximately 30 times higher than that when the discrete retarder method was used with 16 slabs (without any parallelizing routine).

It should be emphasized that, even though this discussion focuses on the SEO, one can model any space-variant medium with the knowledge of the Jones matrix at each point using these methods. When the birefringence distribution or the propagation geometry becomes complicated, numerical methods can still be used to solve for the output electric field. When the material is rapidly varying in space, however, it is necessary to include backward traveling waves (i.e. as is the case with CLC structures).

## 6. Conclusion

We have shown that the numerical methods could be used in successfully modeling the off axis propagation of light through birefringent masks. The discrete retarder method was shown to be very efficient and computationally less complex compared to the differential equation method. Both methods generated very good results. It would be interesting to investigate the behavior of the points of circular polarization and lines of linear polarizations more closely, as this could potentially be another way of generating polarization Möbius strips [16,17].

## Funding

National Science Foundation (PHY-1507278).

## Acknowledgments

The authors gratefully acknowledge helpful conversations and advice from Prof. Miguel A. Alonso and Anthony Vella.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References and links

**1. **L. Mollenauer, D. Downie, H. Engstrom, and W. Grant, “Stress plate optical modulator for circular dichroism measurements,” Appl. optics **8**, 661–665 (1969). [CrossRef]

**2. **A. K. Spilman and T. G. Brown, “Stress birefringent, space-variant wave plates for vortex illumination,” Appl. Opt. **46**, 61–66 (2007). [CrossRef]

**3. **A. K. Spilman and T. G. Brown, “Stress-induced focal splitting,” Opt. Express **15**, 8411–8421 (2007). [CrossRef] [PubMed]

**4. **A. M. Beckley, T. G. Brown, and M. A. Alonso, “Full poincaré beams,” Opt. Express **18**, 10777–10785 (2010). [CrossRef] [PubMed]

**5. **T. G. Brown and A. M. Beckley, “Stress engineering and the applications of inhomogeneously polarized optical fields,” Front. Optoelectronics **6**, 89–96 (2013). [CrossRef]

**6. **R. D. Ramkhalawon, T. G. Brown, and M. A. Alonso, “Imaging the polarization of a light field,” Opt. Express **21**, 4106–4115 (2013). [CrossRef] [PubMed]

**7. **B. G. Zimmerman, R. Ramkhalawon, M. Alonso, and T. G. Brown, “Pinhole array implementation of star test polarimetry,” (2014).

**8. **S. Sivankutty, E. R. Andresen, G. Bouwmans, T. G. Brown, M. A. Alonso, and H. Rigneault, “Single-shot polarimetry imaging of multicore fiber,” Opt. Lett. **41**, 2105–2108 (2016). [CrossRef] [PubMed]

**9. **R. Chipman, “Polarization ray tracing,” in “Proceedings of SPIE - The International Society for Optical Engineering,” , vol. 766 C. Londono and R. Fischer, eds. (SPIE, 1987), vol. 766, pp. 61–68.

**10. **R. Chipman, “Mechanics of polarization ray tracing,” Opt. Eng. **34**, 1636–1645 (1995). [CrossRef]

**11. **G. Yun, K. Crabtree, and R. Chipman, “Three-dimensional polarization ray-tracing calculus i: Definition and diattenuation,” Appl. Opt. **50**, 2855–2865 (2011). [CrossRef] [PubMed]

**12. **G. Yun, S. McClain, and R. Chipman, “Three-dimensional polarization ray-tracing calculus ii: Retardance,” Appl. Opt. **50**, 2866–2874 (2011). [CrossRef] [PubMed]

**13. **R. Chipman, “The polarization ray tracing calculus,” in “Frontiers in Optics, FiO 2014,” (Optical Society of America (OSA), 2014).

**14. **M. R. Dennis, K. O’Holleran, and M. J. Padgett, *Chapter 5 Singular Optics: Optical Vortices and Polarization Singularities* (Elsevier, 2009), vol. 53 of Progress in Optics, pp. 293–363.

**15. **E. J. Galvez and S. Khadka, “Poincare modes of light,” in “Proc.SPIE,” , vol. 8274 (2012), vol. 8274.

**16. **E. J. Galvez and S. Zhang, “Multitwist mobius polarization in crossed complex light beams,” in “Proc.SPIE,” , vol. 10549 (2018), vol. 10549.

**17. **T. Bauer, P. Banzer, E. Karimi, S. Orlov, A. Rubano, L. Marrucci, E. Santamato, R. W. Boyd, and G. Leuchs, “Observation of optical polarization möbius strips,” Science. **347**, 964–966 (2015). [CrossRef] [PubMed]