Combined thermomechanical and optical simulations of planar-optical polymer waveguides

In this work, we describe a theoretical approach for combined thermal, mechanical and optical simulation and analysis of planar polymer waveguides. We consider a finite element approach for thermal and stress/deformation simulation. Also, a Crank-Nicholson finite difference beam propagation method (CN-BPM) is implemented to perform the optical simulation. The results of the finite element (thermo-mechanical) analysis are coupled with the CN-BPM results to carry out the optical simulation of poly(methyl methacrylate) (PMMA) waveguides as function of temperature. For thermal simulation, a model was designed where a polysilicon microheater was added to the upper cladding of the PMMA waveguides to vary the temperature between 20 °C and 200 °C. Thus, the impact of the induced temperature gradients on the refractive index modulation of the PMMA waveguides and the corresponding change in numerical aperture are obtained. In addition, the temperature gradients influence the beam intensity profiles and the movement of the primary eyes within the optical waveguides, thus, impacting the optical properties. Furthermore, the thermally induced mechanical stress and deformation were calculated for transverse and axial directions. In the next step, validation of the model by systematic experimental studies will be performed. In general, our approach provides a toolbox for more comprehensive multi-physics theoretical analysis of polymer-optical waveguides which, in future, can be extended to more complex and functional structures as required for flexible sensor networks, as example.


Introduction
Polymer waveguide structures lend themselves as promising components for broadband telecommunication, i.e. for optical network and computing systems, or sensor networks. All these systems rely on the advancement of optical waveguide manufacturing technology to produce low-loss waveguides and interconnect as well as methods for their comprehensive analysis taking into account all relevant physical effects. So far, many photonic devices such as thermal optical switches (TOSs) [1,2], variable optical attenuators (VOAs) [3,4], optical couplers/splitters [5][6][7], and arrayed waveguide gratings (AWGs) [8,9] were fabricated by employing polymer waveguides. The design and fabrication of a polymer transmission path for communication and sensing application is demonstrated in [10]. Also, the polymer based photonic components integrated to photonic circuits for next generation datacenter systems were reviewed [11]. A simple and reliable technique to fabricate optical interconnects also exploits the self-writing of the waveguides in photosensitive and photopolymer media [12][13][14]. Another simple and cost efficient fabrication technique for planar polymer optical waveguide through hot embossing and doctor blading techniques was investigated by Rezem et al [15,16]. A detailed fabrication of waveguide biosensors by the dip-floating technique was demonstrated in [17]. The poly(methyl methacrylate) (PMMA) grating structures are realized using the femtosecond laser writing technique in [18]. The investigation on photochemical mechanism and the nonlinear photo-absorptive effect on dry acrylamide/polyvinyl alcohol based photopolymer was analyzed by Li et al [19]. Polymers are, in general, attractive for integrated optical devices as they enable rapid processability, cost-effectiveness, high yield, and efficient optical performance, i.e. relatively low optical losses and smaller birefringence compared to their silica counterparts. Furthermore, polymers allow for power-efficient thermal actuation due their large negative thermo-optic coefficients ( dn dT = −(1 to 4 * 10 −4 • C −1 ) which can be 10-40 times higher than in more conventional optical materials such as glass [20]. The thermo-optic coefficients play a vital role in determining and assessing the performance of photonic polymer devices, as detailed by Zhang et al [21]. On the other hand, the operating conditions may easily drift and lead to a change in refractive index, induced by the change in the environmental conditions during the field service, causing thermal and functional instabilities in such devices. The change in refractive index with respect to the temperature mostly depends on two counteracting effects, i.e. the increase in polarizability with temperature and the change in density of the polymer caused by the positive thermal expansion coefficient. A small thermally induced change of refractive index can significantly alter the light intensity distribution within waveguides, especially in coupled systems [22]. By using microheaters, temperature gradients can practically be induced within the waveguides, which in turn is used to switch the light paths [22]. The service temperature of polymers is limited to the range between 60 • C and 250 • C, which is substantially less than for optical glasses (400 • C-700 • C) [23]. Also, the heat and stress generated in the polymer waveguides will eventually lead to mechanical deformations as well. Compressive and tensile stress are introduced in the polymer waveguide either parallel or perpendicular to the beam propagation directions. This is mainly caused due to the difference in thermal expansion coefficients between the core and cladding of the material. The change in the environmental conditions thus represents a major drawback for the usability of polymer materials in this field.
In this article, we focus on combining thermal, mechanical and optical simulations of polymer waveguide structures as basic elements for polymer based photonic devices to more comprehensively analyse their properties and performance and take into account the relevant physical effects. We theoretically investigate the variation of refractive index of optically transparent PMMA waveguides within the temperature range from 20 • C to 200 • C which is a reasonable [24,25]. Alongside, optical simulations of beam intensity profiles within the PMMA waveguides at different temperatures induced by the microheaters are presented. The intensity losses during the transmission and the normalized output power from the exit port are calculated to quantify the combined influence on waveguide geometry, refractive index and numerical aperture. Furthermore, the mechanical deformation generated by induced thermal loads in transverse and axial beam direction are determined. In general, we investigated the relatively simple case where a change in refractive index of the polymer waveguide was induced by the change in the temperature of the microheater. This, in turn, lead to the variation in beam intensities carried by waveguide. Also, the mechanical stress built up due to the heat producing elongations of the waveguide dimensions was studied. This approach to multi-physics simulations can thus be applied to test the influence of environmental instabilities on similar systems, i.e. thermo-optic switches or complex polymer-optical structures for sensor networks.

Finite element thermo-mechanical simulation
The finite element thermo-mechanical simulation and analysis of PMMA polymer waveguides are carried out by using the FEATool that combines a multi-physics and partial differential equation (PDE) simulation environment and allows for modeling of heat transfer, fluid dynamics, chemical engineering, structural mechanics and electromagnetic effects [26]. The modeling process is divided into sequential sub-processes such as definition of the geometry, the generation of the mesh/grid, and the specification of the physical effects and boundary conditions, among others [27]. These processes are often iterative in nature. Also, one may repeat some of the procedure steps to achieve an optimal performance for the designed model. In our approach, we simultaneously used the heat transfer module and plain stress analysis module of the tool under the multiphysics mode to carry out the finite element simulations of the waveguides. The model structures were designed in two dimensions so that they can be easily coupled to our 2-dimensional finite difference based beam propagation method (BPM) used for the optical simulation. Thus, the simulations do not require large memory space. We selected a fast and efficient triangular grid for mesh generation as it allows more flexibility in specifying the grid size for different subdomains of the geometry and boundary regions. In our case, the heat transfer by the conduction mechanism is considered. With respect to thermal analysis, the effect of radiation and convection in the polymer waveguide can be neglected [28]. For the thermal simulations, we used the following equation [29]: Here, ρ is the density, C p the specific heat and k the thermal conductivity of the PMMA material. T is the temperature distribution within the waveguide and T ′ its derivative with respect to time. The detailed values of the parameters used for the simulations are provided in table 1. Furthermore, to calculate the deformation (elongation) in the waveguides the plane stress equation was used: Here, E is the Young's modulus, ν is the Poisson's ratio and α is the thermal expansion coefficient. and v are the displacements along the x-axis and y-axis, i.e. material density, u perpendicular to the optical axis of the waveguide, respectively.

Beam propagation based optical simulation
To calculate the light distribution within PMMA waveguide core, a scalar two-dimensional beam propagation method was implemented, inspired by the work of Pedrola [30]. The slowly varying field envelope is assumed and then the paraxial approximation to the Helmholtz equation is applied to obtain the beam propagation equation as shown below: A Crank-Nicholson based finite difference scheme is applied to equation (4) to obtain the numerical solution, which is non-local in time. Here, n 0 is the initial average value for the refractive index of the homogenous media. The current value of refractive index is n, which is a dynamic value and is calculated from the thermal distribution within the PMMA sample, as shown below. The free space wavenumber is k 0 and calculated from the wavelength (λ) of the input Gaussian beam via k 0 = 2 * π/λ.The electromagnetic field amplitude is denoted by u (x, z). x and z are the transverse and axial beam directions, respectively. The waist of the Gaussian beam is assumed to illuminate the PMMA sample from the left input facet and the output beam is collected from the right exit facet. Along with the CN-BPM transparent boundary conditions (TBC) are employed to equation (4) to carry out the optical simulations [31]. Because of lower computational complexity and reduced computational time, the CN-BPM was considered to carry out the optical simulations of the waveguides. In general, our approach can also be combined with the more accurate finite difference time domain (FDTD) method for the optical part. Figure 1 shows the schematics of the cross-sectional view of the planar PMMA waveguide that was further considered for the thermo-mechanical analysis along with the optical simulations. The core part of the waveguide is assumed to be PMMA polymer material, whereas the upper/lower cladding of the waveguide is considered to be from monomer in liquid form. In general, polymers can be used for thermo-optical switches, due their large negative thermo-optic coefficients. To mimic this case, we added a polysilicon microheater to the model. As the position and the dimension of the microheater have a significant impact on the efficiency of the thermo-optic (TO) waveguide switch and the switching time, the microheater was placed at two different positions, on the top and inside the upper waveguide cladding, respectively. In the first case, the polysilicon microheater was placed on the top of the upper cladding; see figure 1(a). Here, the width of the core and cladding part of the waveguide along the transverse x-(and y-) axis were 40 µm and 20 µm, respectively. The length along the propagation z-axis of the PMMA waveguide was 1600 µm ( figure 1(a)). The polysilicon microheater in figure 1(a) had a width of 20 µm along the x-axis and a length of 1200 µm along the z-axis. In the second case, the microheater with the same dimensions was buried in the middle part of the upper cladding. In figure 1(b), the waveguide core was again considered to be PMMA material with a width of 40 µm along the transverse x-(and y-) axis and a length of 1600 µm along the z-axis. The upper/lower claddings were again monomer material with a width of 80 µm along the x-(and y-) axis. The microheater was placed at a distance of 30 µm from the core part in figure 1(b). The width of the waveguide core is the same for both models shown in figures 1(a) and (b). The cladding width in figure 1(b) is different compared to figure 1(a) in order to fit the top microheater within the upper cladding of the waveguide. The optical axis of the waveguides, the central line L1, passes through the middle section of the core along the z-direction and was considered further to obtain the induced elongation. Using these pre-defined geometries, the thermal simulations were performed by using the already described heat transfer module. The initial temperature distribution was T 0 = 20 • C for the whole waveguide. Then, the temperature of the polysilicon heater was varied from the initial temperature T 0 to the maximum temperature T max = 200 • C. Thereby, the output intensity from the exit port of the waveguide dropped below 50% of the initial input beam intensity. Appropriate boundary The temperature along the line L1 was increased from 20 • C to 110 • C for the top microheater on the upper-cladding model. A larger increase in temperature from 20 • C to 148 • C was achieved along the line L1 for the buried microheater within the upper cladding model. For both these cases, the microheater was at the temperature of 200 • C. The solver was static for all simulations. However, in one case it was selected time dependent to find the spreading time for the heat throughout the waveguide. In that case, the spreading time was 7 ms. The microheater placed on top of the waveguide cladding provided better and controlled increase in temperature while the embedded microheater could possibly be used to obtain quicker switching times. Also, the temperature distribution along the reference line L1 shows a slightly smoother increase for the top microheater compared to the buried microheater case; see figure 3. It is also more complex to implement a heater within the cladding of a waveguide in real world situations. Also, the modulation in temperature along the line L1 for the top microheater increased gradually when the microheater temperature was increased from 20 • C to 200 • C; figure 4. The temperature through the line L1 in the PMMA core part of the waveguide never exceeded the value 110 • C, when the top microheater was at a high temperature of 200 • C. The increase in temperature along the line L1 when the temperature of the buried microheater was varied between 20 • C and 200 • C is provided in figure 5. The increase in temperature along L1 was proportional to the temperature of the microheater. The buried microheater induced a larger increase in temperature along the line L1 than the heater on the top of the waveguide.

Results and discussions
Zhang et al evaluated many polymers including PMMA and demonstrated a linear relationship between the thermo-optic coefficients (dn/ dT) and the coefficients of thermal expansion (α) [21]. The change in refractive index as a function of temperature was measured for PMMA material, which is also considered in our case [21].
Here, the refractive indices for the core and cladding of the polymer waveguide are calculated by considering the temperature distributions obtained from the thermal simulations. Equation (5) was utilized for this [32].
Here, dn dT is the thermo-optic coefficient for PMMA and is taken as dn dT = −1.3 × 10 −4 • C −1 . T 0 = 20 • C is the initial room temperature and T 1 is the induced increased in temperature due to the heating by the microheater, as obtained  from the solution of the finite element heat analysis. n 0 is the initial value of the refractive index at room temperature and n temp is the new/updated value of refractive index calculated from the change in temperature within the waveguide due to the heating process and the thermo-optical coefficient.
The initial values of refractive indices of the PMMA core and monomer cladding were 1.48 and 1.41, respectively. The polymer material possesses a negative thermo-optical coefficient, so there is a decrease in the refractive index value with an increase in temperature, as also observed in the simulations. The observed decrease in refractive index of the core of the waveguide was from 1.4885 to 1.4600, while the temperature of microheater was varied between 20 • C and 200 • C. The corresponding refractive index of the cladding decreased from 1.412 to 1.392 during the heating process. Also, the associated change in numerical aperture was calculated for  the temperature range considered; see figure 6. The numerical aperture (NA) of the PMMA waveguide decreased from 0.471 to 0.457 during this process. Here, the NA for the waveguide was calculated by using √ n 2 core − n 2 clad . There was a refractive index gradient across the core of the waveguide. However, it was not significant (the size of the gradient across the core was 0.00015) due to the smaller core diameter. The refractive index gradient was considered for performing the optical simulations.
Subsequent to the finite element analysis for the thermal distribution, the computational waveguide models were transferred to the CN-BPM simulation (implemented in MATLAB) to carry out the optical simulations. The initial input beam with Gaussian beam profile was coupled into the PMMA waveguide core. The intensity of the input beam was normalized to the maximum peak value of one. The optical simulation was then carried out at different temperature distribution profiles of the PMMA waveguides imported from the FEATool. The obtained beam intensity distribution profiles are shown in figure 7, for selected temperatures of the microheater between 60 • C and 180 • C, respectively. The decrease in the output intensity of waveguide was not significant until the temperature of microheater reached 60 • C. Afterwards, the decrease became more prominent. The primary eyes which are the localized maximum intensity points that indicate the guiding properties of PMMA waveguides tend to move within the waveguide as the light beam propagates within them. Back and forth movements of the primary eyes within the PMMA waveguide cores were observed when the microheater was at selected temperatures of 60 • C, 80 • C, 100 • C and 150 • C, as shown in figure 7. This is an indication of changed beam propagation conditions and, in this case, decreased transmitted beam intensity at the output of the waveguide core. The primary eyes almost disappeared when the microheater was at 180 • C. For this temperature, the output beam intensity from the PMMA waveguide dropped to 50% of the initial beam intensity. The losses in the waveguide are due to the change in NA as well as due to the residual stresses build up during the heating process. However, we assume that the residual stress is less significant due to the small Young's modulus of the polymer material. As the current simulation is a multi-step process, the accumulated error could be larger than the actual physical process. The change in refractive index of the material depends on the thermo-optic coefficient. In addition, the change in refractive index depends on the size of the microheater as well as the size of the wave-guiding element. In the current simulations, the dn dT value of PMMA was considered for both the core and the cladding part. This could be different for different materials  and could vary in the actual physical processes. Experimental validation of the simulation is currently being performed to validate the approach. The output power at the exit port of the waveguide was then calculated for each case by considering the last slice of the PMMA waveguide and integrating the intensity values over the area. The obtained output power from the exit port of the waveguide is shown in figure 8. The normalized output power decreased as the heating by the microheater was increased. A sharper decrease in the output power was observed between 60 • C and 80 • C. At this point, the PMMA material approaches the glass transition temperature [33]. The output power continued to decrease further and at 180 • C the P output = 0.54 * P input . When the temperature of the heater is increased to even higher temperatures, the decrease in refractive indices of core and cladding will be larger and finally the beam would not be confined to the waveguide core anymore, but diffuse to the surrounding media.
The deformation or elongation produced in polymer based integrated photonic device is an important consideration to access their performance in varying environmental conditions. Thus, for the integrity of the performance, a link was established between the thermal distribution and stress analysis module in the FEATool. The corresponding elongations/deformations generated in the waveguide are quantified. The finite analysis solver of the Multiphysics FEATool was unitized for the stress and deformation analysis. In particular, the plane stress module with proper boundary conditions was employed for the deformation analysis of PMMA waveguides. The left and right extremities of the waveguide were fixed for the model in figure 1(a). The boundary conditions such as fixing one end or both ends in deformation analysis are crucial. Fixing both ends of the waveguide is a better practical approach as optical waveguides and interconnects are usually placed between fixed active components (LDs or waveguides) in integrated photonic circuits. The heat transfer module was coupled to the plane stress module of the FEATool to perform the thermo-mechanical multiphysics simulations.  Here, the initial temperature of the waveguide was again considered as 20 • C. The parameters used for the stress analysis are provided in the table 2.
In table 2, ρ is the material density, E is its Young's modulus and ν is the Poisson's ratio signifying the ratio between the transverse contraction strain to the longitudinal extension strain in the direction of the stretching force. The thermal expansion coefficient is denoted by α. As can be expected, the maximum deformation was obtained when the microheater was at the temperature of 200 • C. The initial waveguide without deformation (when the microheater was at 40 • C) is shown in figure 9(a), whereas the deformed waveguide (when the microheater was at 200 • C) is given in figure 9(b). The deformation analysis for the microheater at 200 • C is presented in figure 10. The elongations along the transverse x-axis and the beam propagation z-axis are shown in figures 10(a) and (b), respectively. We obtained a maximum elongation (∆L x ) of 831 nm along x-axis and a maximum elongation along z-axis ∆L z = 1397nm when the microheater was at 200 • C. The elongations produced along the x-and z-axis are small compared to the original waveguide dimensions. The residual stress build-up in the waveguide is assumed to be less severe because of the small Young's modulus of the polymer material. The effect of residual stress is thus not included in the current BPM simulation.

Conclusion
A combined simulation approach for the thermo-mechanical and optical study of polymer waveguide structures was carried out by incorporating a polysilicon microheater into a waveguide model. The microheater was located on top of or inside the upper cladding of the polymer PMMA waveguide. A finite element analysis was performed by using the FEATool box to obtain the heat distribution and mechanical deformation within the waveguides. The finite difference approach was adapted and implemented into a CN-BPM model in MAT-LAB. The simulation results obtained from the FEATool box were coupled back to the MATLAB module. The obtained modulation in temperature distribution was found to affect the light beam intensity distribution within the PMMA waveguide significantly, as the polymer material is very sensitive to environmental changes. The refractive indices of core and cladding were computed from the thermal distributions by considering their influence on the thermo-optic coefficients. Due to the heating, the NA of the polymer PMMA waveguide was observed to decrease from 0.471 to 0.457. The obtained beam intensity profiles within the PMMA waveguide at various microheater temperatures were determined. The movement of maximum intensity primary eyes signifying the beam guiding properties of the waveguide was no longer visible when the temperature of the microheater was increased beyond 180 • C. The normalized output power from the exit port of the waveguide decreased to 0.54 * P input when the microheater reached 180 • C. Consequently, the normalized output power from the exit port of the waveguide is also varied during the heating. No significant change in the output power occurred until 60 • C. Furthermore, for the model dimensions considered here, the deformation of the PMMA waveguide was found to be 831 nm along the transverse x-direction (∆L x ) and 1397 nm along the beam propagation z-direction (∆L z ) when keeping the waveguide end facets at fixed locations. The study performed helps us to understand the relevant effects and dynamics of PMMA waveguides as a function of relevant parameters and conditions. The modulation of the refractive index with respect to the varied temperature load is a critical aspect for thermo-optical switches or future polymer-based planar sensor networks. More complex PMMA thermo-optic switches and PMMA photonic components will be designed and simulated for thermomechanical analysis and optical simulations in future. Also, systematic experimental studies will be initiated to verify the theoretical findings and eventually extend and optimize the simulations.