Femtosecond laser direct-write waveplates based on stress-induced birefringence

The use of femtosecond lasers to introduce controlled stress states has recently been demonstrated in silica glass. We use this technique, in combination with chemical etching, to generate and control stress-induced birefringence over a well-defined region of interest, demonstrating direct-write wave plates with precisely tailored retardance levels. This tailoring enables the fabrication of laser-written polarization optics that can be tuned to any wavelength for which silica is transparent and with a clear aperture free of any laser modifications. Using this approach, we achieve sufficient retardance to act as a quarter-wave plate. The stress distribution within the clear aperture is analyzed and modeled, providing a generic template that can be used as a set of design rules for laser-machined polarization devices.

Femtosecond laser exposure in the non-ablative regime induces a confined change of volume [24], creating oriented tensile or compressive stress fields capable of reaching locally high levels of stress (up to ~ 2 GPa in compression) [25]. Interestingly, these laser-induced stress fields can be tailored and controlled by changing the orientation of the laser polarization with respect to the writing direction [25,26], and by tuning the pulse duration of the laser [27]. This not only allows for control of the stress orientation, but also the magnitude and state, i.e. switching from tensile to compressive stress. This ability to fine-tune stress fields finds its origins in the variety of nanostructures induced in laser-modified materials, as highlighted in [25][26][27]. This intriguing effect has already been applied to tensile testing procedures for investigating glass micromechanics [28], as well as for future concepts of optical component packaging [29].
Here, we explore the use of controlled stress-states for applications in polarization optics. Specifically, we present a technique for manufacturing polarization devices using the indirect action of laser-induced stress. Device geometry is defined using laser-machined silica substrates, limiting the action of the induced stress (and resulting birefringence) to a particular aperture. Unlike other techniques, which rely on direct light interaction with formbirefringent nanogratings [30], this clear aperture is free from laser-modifications, potentially affording higher power handling capabilities, higher overall device transmission, and a cleaner undistorted beam profile. Additionally, the use of fused silica substrates provides a broad transmission window, potentially enabling a wide tuning range for the operating wavelength. With fast fabrication times and full control over the device retardance, broadband-like devices consisting of many waveplates on a single substrate are possible. Furthermore, by using the same femtosecond laser for both machining and subsequent loading, the device manufacturing process can be greatly simplified.

Waveplate concept
The principle of a laser-machined waveplate is shown in Fig. 1. A rectangular clear aperture (labeled as zone c) is defined by two cuts made through a fused silica substrate using a wetetch-assisted laser machining process [31].
After etching, a set of lines (which we call 'stressors') are machined in the bulk of the material at opposing ends of the rectangle (b), using writing parameters sufficient to generate nanogratings within the laser-affected zone [5]. These nanogratings are oriented parallel to the laser writing direction, such that the principal component of the stress tensor [25,26] is directed perpendicular to the line orientation.
This arrangement of opposing stressors induces a quasi-uniaxial loading of the material in the center (as indicated by the red arrows in Fig. 1), generating a sizeable optical retardance within the clear aperture. With this simple arrangement, the cuts serve as a stress-free interface, confining the action of the stressors to the area of interest, in which the level of optical retardance can be tailored through careful control of the number, density, and exposure parameters of the stressors. It should be noted that, while the fabrication process presented here is done in multiple steps, the action of machining and device loading could be combined prior to etching, simplifying the production process further. Fig. 1. Microscope image of a laser-machined waveplate. In this arrangement, a clear aperture (c) is defined by a pair of horizontal cuts (a), fabricated using a wet-etch-assisted laser machining process [31]. Stress-induced birefringence is then generated in this aperture by writing a set of vertical lines (labeled 'stressor bar', (b)) within the bulk of the substrate at each end of the device, using exposure conditions sufficient to generate nanogratings. Here, the writing polarization, and hence the orientation of the nanogratings, is chosen such that the generated stress behaves like a quasi-uniaxial loading, as indicated by the red arrows.

Setup and measurement principle
Waveplates were fabricated in several steps, beginning with the machining of 20 blank devices (no stressors, each device 1x2 mm, staggered on a 4 mm grid) in a 25x25x0.5 mm fused silica substrate, using an amplified femtosecond laser system (Yb-fiber, Amplitude Systèmes), delivering 280 fs pulses at 1030 nm. Here, the laser was focused using a 20x objective (NA 0.4) optimized for 1046 nm. Following machining, the sample was etched in a 2.5% solution of hydrofluoric acid for ~ 8 hours. To maximize etch rate, the variable repetition rate of the laser was set to 760 kHz during machining [31].
Using the same writing configuration, each device was loaded by progressively forming stressors perpendicular to the cuts at both ends of the waveplate. The stressors were then arranged in blocks at each end of the waveplate at a fixed spacing of 5 or 10 µm, with the laser polarization set perpendicular to the writing direction in order to generate maximum stress along the long-axis of the waveplate (y-direction in Fig. 1) [25,26]. Here, we define a single 'stressor' as a series of stacked discrete laser-modified regions with a z-spacing of 15 µm, forming a modified 'sheet' of material in the xz plane, for a total of 29 lines per sheet.
Each group of stressors was embedded in the bulk of the substrate and spaced 20 µm from all exterior surfaces to limit stress concentration and consequent crack formation. Additionally, each block was inset 100 µm from the ends of the device. To ensure even loading, the repetition rate of the laser was reduced to ~ 80 kHz in order to maintain the required value of energy deposition for maximum stress, while keeping the writing velocity low, thereby limiting the influence of acceleration dynamics of the motion stages used during fabrication.
The resulting stress-induced birefringence was measured using an optical microscope (Olympus BX51) equipped with a liquid-crystal universal compensator for measuring birefringence (VariLC). The LC unit was controlled using freely available software (OpenPolScope [32][33]), which provided all the necessary functions for calibration of the LC compensator and retardance measurement acquisition. Further post-processing in MATLAB was then used to apply a false color-map to each measurement, as well as extract retardance data. For the purposes of this experiment, we only focus on the measured absolute retardance.

Experimental results
We first examine the evolution of stress within the clear aperture of the waveplates as a function of deposited energy [34]. A pulse energy of 250 nJ was selected, corresponding to a laser-modification which lies solidly within the nanograting regime for optimal stress generation [25,26]. The speed was varied from 1700 to 170 µm/s, resulting in a deposited energy range of 10 to 100 J/mm 2 . The results of this experiment are shown in Fig. 2 below.
For this experiment, the line spacing of individual stressors was fixed at 10 µm, with a total of 16 stressors per side (32 per device). From Fig. 2, the developed stress peaks around 20 J/mm 2 , giving a maximum retardance of just over 15 nm. Above this point, the retardance begins to decay, which we interpret as a consequence of stress-relaxation due to crack formation within the nanogratings. This measurement not only serves as a calibration to develop maximum stress during the writing process, but also confirms previous findings for stress evolution in the nanograting regime [24]. For reference, the phase shift at 546 nm (the wavelength used by the VariLC measurement system) is given on the right hand side of the graph. While the peak retardance developed in Fig. 2 is only on the order of ~ 15 nm, this plot serves as a calibration for process control to determine the best writing conditions for generating stress in a single line. For further optimization of the developed stress, we now turn our attention to the lateral spacing of the stressors.
Based on the results shown in Fig. 2, a deposited energy target of 20 J/mm 2 was chosen to fabricate three additional devices: one with a stressor spacing of 10 µm and 64 total stressors, and two with a reduced spacing of 5 µm, each with 64 and 128 total stressors respectively. The results are shown in Fig. 3 below. The device in Fig. 3(a) is the same displayed in Fig.  2(a), and is used as a reference (20 J/mm 2 , ~ 15 nm max retardance).
Here we measure the developed retardance at the center of the clear aperture as a function of the spacing and number of lines in each stressor. Note that the retardance scale has a maximum retardance of 100 nm. For the device shown in Fig. 3(d), a maximum retardance of ~ 55 nm is developed. This retardance map will be used as a comparison and calibration for the following section on device modeling.

Modeling and discussion
We now set out to describe the behavior of the developed stress in the waveplates in terms of a set of design parameters, such as the device geometry and the number and density of the stressors. These models may then be used as a set of design rules, giving full control over the stress field under different geometrical and loading conditions. As a first approximation, we consider a one-dimensional lumped model (illustrated in Fig. 4 below), following an approach reported in [28]. Here the modified and unmodified regions of the material are represented by springs, for which an equivalent stiffness is derived.

One dimensional analytical model
In this model, the geometry of the waveplate is broken down into several distinct regions, as shown in Fig. 4(a). These regions are labeled as the stressors (composite laser modified volume), the un-affected sidebars, and the clear aperture of the waveplate. Each of these regions may then be represented by an equivalent stiffness, as shown in Fig. 4(c).
The model is built by first considering the composite region containing each stressor block, which is treated as a parallel and series combination of springs and assigned an equivalent stiffness, K D , as illustrated in Fig. 4(c). This assembly is then placed in series with the clear aperture of the waveplate, K WP .
After formulation of the total stiffness of the system, and taking into account stressinduced birefringence, we arrive at the following expression for calculating the retardance in terms of the number of stressors, n: where the factor λ is given by: The parameters α, β and γ are dimensionless ratios, defined as α = L LAZ / L INT , β = t /t s and γ = w / w s , respectively, while κ is the ratio of the Young's modulus of pristine SiO 2 to that of laser modified SiO 2 . Note that in previous work [28], we estimated the Young's modulus for the laser-affected zones (nanograting modification regime) to be approximately one third of the modulus of pristine, un-modified silica. The quantities w s , t s , w, and t are the width and length of a stressor and surrounding material, respectively, while the quantities L INT and L LAZ correspond to the spacing between laser-modified regions (vertical sheets)and the thickness of these regions (along the y-direction), respectively. C is the stress-optic coefficient (see Αppendix), L WP is the initial length of the clear-aperture, E SiO2 is the elastic modulus of fused silica, and ε LAZ is the net volumetric expansion induced during laser exposure [25].
From a design point of view, Eq. (1) gives an overview of the tuning parameters available to the user. The choice of Young's modulus is fixed by the choice of fused silica (material component), while the induced strain (laser-induced component) is also relatively fixed for optimized stress generation, as demonstrated by the waveplates shown in Fig. 3. This leaves device geometry as the main parameter, and we note that, for large numbers of stressors, the parameter λ converges to zero, leaving only the ratio of the stressor spacing, material thickness, and the number of stressors.
It should be noted that this model is a simplification, as the interfacial energy between the various domains (laser-affected, pristine, etc.) is not considered. However, so long as we assume that the behavior of these regions is governed by linear elasticity, this approximation remains valid. Despite these subtleties, the mechanical behavior of the system is reasonably well predicted using this approach.
In addition to the simplifications mentioned, we make the following assumptions: 1/ Fine structures within the laser-modified composite are ignored, and instead treated as a homogeneous material with a lower density than the surrounding, unmodified bulk. An average elastic modulus of 34 GPa was chosen for this region, which was recently reported in [28].
2/ For the particular arrangement of embedded stressors used here, a certain portion of the device loading is transferred to the bulk as shear stress; in this case along the regions label 'unaffected sidebars'. We consider this shear low and negligible, as calculated based the shear-lag model used in [28].  Fig. 3), while the curves are the predicted values based on the model. For each measured data point, the error was estimated from the measurement noise, which on average was not greater than ± 1 nm. The reader should note that the measured data points shown in the above graph compare different sizes of clear aperture. The devices written with 32 stressors 10 µm spacing and 64 stressors 5 µm spacing have the same clear aperture size, while conversely the second two data points are also comparable. Figure 5 plots the developed retardance as a function of the number of stressors, comparing the predicted retardance based on the 1D model with the measured values from the center of each waveplate shown in Fig. 3.
Overall, the model is in good agreement for low and high numbers of stressors, however there is some deviation for higher numbers of stressors. While this model provides a reasonable estimation, the accuracy it provides is limited, as we do not consider the 2D interaction of stress present in the waveplate. For an improved model, we now examine the full 2D structure, which we cover in the next section.

Two-dimensional analytical model
While the simplified 1D model discussed in the previous section provides an adequate approach to predicting stress evolution from the number and density of machined stressors, it does little to describe the full contour of the induced retardance field. For analytical derivation of this full contour, we use a more refined model based on Airy stress functions [36][37][38][39]. In this model, the waveplate is treated as a two-dimensional body (finite rectangular domain) with a continuously distributed load along the width (provided by the action of the stressors), as depicted in Fig. 6.  While a full derivation on the use of Airy stress functions is beyond the scope of this work, we will briefly outline the process here, referring to the notation used in Fig. 6. Here, we assume a right-handed coordinate system, the origin of which is taken to be in the center of the rectangle defining the clear aperture. Following the work in [39], a general solution to the stress developed in the waveplate is found by defining a stress function φ(x,y), which satisfies the biharmonic equation: The equation φ(x,y) is of the form: sin ax The stress function then becomes, after integration of Eq. (5): where the coefficients C 1 -C 4 are computed using the boundary conditions of the substrate at the upper and lower edges (y=±c). In order to accommodate generalized representations of the load applied to one edge of the substrate, a Fourier series is used, as follows: with an additional loading term computed for the lower edge. For our case of symmetrical loading, the terms containing sin(αx) vanish from Eq. (7), and the coefficients A o and A' m may be computed by: Here, A o represents the uniform load applied to the right-hand edge. Finally, the 2D stress tensor is related to the stress function by: σ xx = ∂ 2 ϕ ∂y 2 , σ yy = ∂ 2 ϕ ∂x 2 , and τ xy = ∂ 2 ϕ ∂x ∂y (9) where σ xx , σ yy , and τ xy are the x, y, and shear components of stress, respectively. Computing the derivatives of Eq. (6), we arrive at the solutions that satisfy the biharmonic equation, giving a full tensorial description of the stress field in the clear aperture of the waveplate: Finally, using the stress components shown in Eq. (10), the birefringence field was computed using the photoelastic relationships outlined in the Appendix. Note that, in our calculations we ignore the contribution from shear stress, as this component plays a minimal role in the final device retardance. Finally, for the equations shown above, we compute the stress components for the first fifty orders of the Fourier series, which we find adequately describes the stress-field developed in the clear aperture. The coefficient D is a correction factor, which will be discussed in Section 4.4.

Finite element modeling
Finite element modeling was used as a means to verify the experimental results, as well as to provide an additional metric for comparison against the 2D analytical model discussed in the previous section. Of the many FEM software packages available, freeFEM++ [40] was chosen due to its low-level programmatic approach, giving a high degree of control over the modeling process. The simulation domain was modeled after the peak device shown in Fig. 3(d) (maximum center retardance of ~ 55 nm) using the displacement formulation [40]. As with the 1 and 2D analytical models, the stressors were treated as a homogeneous volume, and in this model, were represented as a rectangular region occupying the same area as the stressors of the original device. Furthermore, Dirichlet boundary conditions were used, specifying that all displacements must go to zero at the outer boundaries of the substrate. For this reason, the outer boundaries were grounded while all other boundaries internal to the domain were free to move.
In order to simulate the loading action of the stressors, an initial displacement was specified, directed along the long dimension of the device (see Fig. 1), for both lateral boundaries of each rectangle representing the stressor region. This value was computed by multiplying the size of the stressor region along the y-direction (here approximately 315 µm) by a factor of 0.03%, the average unit strain for nanograting-induced stress (experimentally measured in [24]).
Finally, stress-induced birefringence was calculated based on the photoelastic equation (outlined in the Annex), and the resulting field distribution was exported for processing in MATLAB. The results of this simulation, showing only the region of the clear aperture, are shown in Fig. 7. It should be noted that the value 0.03 % is an average strain for nanograting-induced stress, and as we have done here, was computed by treating the composite region of the lasermodified material as homogeneous. As such, the resulting calculated retardance in the FEM model (when using this value) is under-estimated. To achieve a retardance value closely matching that of the original device, an iterative approach was used, gradually increasing the 'virtual' strain in the model until the desired retardance at the center of the clear aperture was achieved, resulting in a final value of ~ 0.16 %. Figure 8 illustrates a comparison between experimental data, FEM simulation, and the analytical model for line profiles taken along the x and y directions through the center of the clear aperture. For reference, the retardance map for the peak device in Fig. 3(d) is shown in (a).

Comparison with experimental results
We first examine the prediction of the shape of the stress field, as illustrated by the contours in Fig. 8(e-g). Qualitatively, the 2D model (f) is in good agreement with the contour of the measured retardance (e), while the FEM model (g) has similar features, but deviates near the boundaries. While the 2D model performs well for the contour, the overall magnitude variation is not properly estimated, as shown by the line profile plots for the y and x directions in (c) and (d) respectively. The upper insets show the analytical model prediction, and while the shape of the profile closely follows the measurement, the overall magnitude is incorrect. We suspect that this may be caused by a slight out-of-plane deformation, suggesting that the loading induced by the stressors requires further optimization. Indeed, no compensation for spherical aberration was used during the writing process. This can account for a decrease of localized intensity with increased writing depth, resulting in a slightly unequal distribution of stress across the substrate thickness. This stress gradient generates a net moment, and in turn favors out-of-plane bending. Such effects have been reported in previous work [28] and can be easily corrected, either by adjusting the power or the stressor spacing with respect to depth within the sample. To adjust the analytical model, we have included a scaling factor (labeled as the variable D in Eqs. (10)), which is shown here for a value of 25, indicating that the overall prediction of the 2D model is in agreement with the measurement, provided that this adjustment is made. We call the reader's attention to the fact that this is an artifact for curve fitting. The exact physical meaning of the parameter D in terms of stress remains unclear at this stage.
The FEM model also deviates considerably from the measurement, which may be attributed to non-uniformity of the distributed load on the as-measured device over the thickness close to the boundaries of the applied load.

Conclusion
This work demonstrates the use of direct-write femtosecond laser processing to induce predictable, controlled stress states within a chosen region of a substrate through the indirect action of stress-induced birefringence. Here, we utilize this technique to fabricate a polarization device through the indirect action of stress-induced birefringence, resulting in a simple waveplate. In this example, a basic geometry is used, consisting of a series of vertical, laser-written planes arranged in a perpendicular fashion, creating a rectangular region under uniaxial stress.
Though the results presented here are not fully optimized, we have provided a generalized framework, with which further development can be made. The flexibility of the direct-write process allows for virtually infinite combinations of cuts and stressors, opening up new design opportunities for more complex polarization devices with clear apertures and potentially zero transmission losses. For example, one could consider devices based on twodimensional stress states (as demonstrated in [41] at the macro-scale), leading to the creation of optical vortices. The technique we have introduced here allows for such a concept to be implemented at the micro-scale, and may be tailored to a variety of complex stress states.