Computational analysis of cooling dynamics in photonic-crystal-based thermal switches

The paper presents optical and thermal simulations of two kinds of thermally controlled silicon Photonic Crystal devices: an air hole in silicon slab that switches between two refractive behaviors, and a T-shaped circuit made by silicon rods in air that has an on-off behavior. Both effects have been obtained by increasing the device temperature, then exploiting the thermo-optic effect. Theoretical models of thermal dynamics, in particular for natural convective heat exchange on micrometric and sub-micrometric scale, have been validated by means of numerical simulations.


Background
The control of light propagation in a Photonic Crystal (PhC) slab can occur through the modulation of several parameters [1].Indeed, both structural parameters (as refractive index of constituent materials, array geometry, etc.) and light parameters (as incident angle, polarization status and wavelength) influence the light propagation inside the PhCs [1][2][3].Fixed the light parameters, the refractive index varies with the temperature, according to the thermo-optic effect [4,5].Therefore, it is possible to control light propagation in a PhC by means of temperature variations.This paper presents theoretical results about two different thermally modulated PhC devices: a PhC slab based on air holes in silicon that exploits the negative refraction behavior [6,7] and a T-shaped PhC circuit based on silicon rods in air that exploits the band gap properties [8].Both devices have been simulated thermally and optically using commercial software, COMSOL multiphysics® and RSoft™, respectively.In particular, a rigorous numerical analysis of the cooling phase of these structures has been performed in order to determine the dynamic and the thermal time constant at micrometric and sub-micrometric characteristic length.In [9], the trend of the natural convective heat exchange coefficient h c has been obtained analytically from the relations between Nusselt (Nu), Reynolds (Re) and Prandtl (Pr) numbers at characteristic length in the range of (0.01÷1000) mm.In this range, the h c coefficient is inversely proportional to the characteristic lengths below 10 mm and almost constant for higher dimensions.In [10], an experimental validation at the same characteristic length has been reported.However, the interesting range of PhC-based devices was almost unexplored.In this paper, the cooling dynamics of sub-micrometric devices has been investigated using a thermal-fluid-dynamic approach.Results emphasize that the dominant contribution is due to the conductive mechanism, in contrast to the statements of J. Peirs et al. in [9].Moreover, using thermooptical effect and optical Finite Difference Time Domain (FDTD) simulations, the optical efficiency of two analyzed devices and the value of the switching bandwidth of both geometries have been calculated.

Methods
Changing temperature is a simple way to obtain a refractive index change and then a photonic modulation of silicon based devices.However, for macroscopic devices, the relatively slow heat diffusion in solids strongly limits their performances.In this frame, PhC devices are promising to overcome the problem, because they are constituted by an array of elements of micrometric or sub-micrometric dimensions.Before to simulate light propagation and calculate efficiency of two specific devices, two preliminary analyses have to be done.The first is needed to establish if and which thermal exchange mechanism is predominant in a generic pillar-based device.The second has been done to simulate both thermal and optical behavior for specific PhC structures.

Thermal analysis
Preliminary Finite Element Analysis (FEA) simulations have been conducted on the elementary cell sketched in Fig. 1, in order to determine the predominant thermal exchange phenomena.Four different mechanisms of heat transfer have been separately considered, as follow: The first approach solves the total energy balance and the flow equations of the outside cooling air.In this case, a Computational Fluid Dynamics (CFD) analysis has been performed, producing detailed results for the flow field around the pillar, the temperature distribution and the cooling power.However, this approach is more complex and requires more computational resources than the others, but it can be considered as a rigorous model for the actual heat dissipation mechanism.The second approach assumes that heat flux from solid to air is not enough large to induce fluid density changes and consequently laminar convection cells initiation; in this frame, air is considered as a still medium where heat is dissipated only by conduction.The third approach describes the outside heat flux by means of a heat transfer coefficient function, that is automatically calculated by the software as function of surface area, geometry and orientation.This assumption is generally valid for macroscopic objects and has been numerically verified in this work for micrometric and sub-micrometric structures.The fourth approach is carried out with the aim to evaluate the relative influence of the radiative heat exchange, compared to the other mechanisms.
In all approaches, the basic cell geometry is parameterized in terms of pillar height hp, as specified in Table 1.Moreover, in cases A and B a cylindrical air volume surrounding the pillar is considered, while in cases C and D only the silicon pillar is considered, without the surrounding air volume.
The FEA analysis has been performed by means of Comsol® software.The thermo-physical characteristics of the investigated structure are reported in Table 2.  equations combined with realizable k-ε turbulence model [11].

Mathematical model and boundary conditions
where ρ is the fluid density, v → is the fluid velocity vector, p is the pressure, μ is the dynamic viscosity, I is the identity tensor, g → is the gravitational acceleration, T is the temperature, c p is the specific heat and k is the thermal conductivity.The above mentioned equations have to be particularized for the considered simulations.
In particular, for B, C and D simulations the terms that contain velocity vector v → are null.As initial conditions, the temperature T and the pressure p of the air (modeled as ideal gas) have been assumed equal to 288.15 K and 101325 Pa, respectively; moreover, the initial temperature of the silicon pillar T p has been assumed equal to 400 K. Considering a cylindrical coordinate system (r, φ, z), the boundary conditions, at the pillar/air interface for the considered cases, are summarized in Table 3.

FEA analysis results
The lack of reliable theoretical and numerical models for natural convective heat exchange on micrometric and sub-micrometric scale was overcome by means of parametric Finite Element Analysis (FEA) of the elementary PhC cell (sketched in Fig. 1).As previously explained, in this scaling analysis all dimensions are related to the pillar height h p (see Table 1).Considering a natural convective heat exchange by the vertical wall of the pillar, the dimensional scaling of the h p parameter down to 1 μm allows to extrapolate an equivalent convective coefficient h c from the resulting cooling dynamics.As shown in Fig. 2, FEA calculations confirm the literature results [9,10] in the size range from 10 μm to 100 μm.
If pillars height is smaller than 1 μm, the cooling dynamic changes dramatically.Figure 3 shows the resulting cooling dynamics for a singular sub-micrometric pillar structure with h p = 200 nm for all four different heat transfer mechanisms above described.The represented temperature, considered at the center of pillar, ranges from 400 K to ambient temperature.It is worth noting that the case A, the conjugate heat transfer model, can be assumed as the most rigorous approach which more closely reproduces the real dynamic of the structure cooling.At the same time, a comparison among the cooling dynamics obtained in the B, C and D cases indicates the predominant mechanism of heat exchange and dissipation from the pillar toward the surrounding environment.
First of all, Fig. 3 shows that the time required for pillar cooling, when only radiative heat exchange mechanism is considered (case D), is about one hundred times longer than other three cases.Hence, for the considered geometry and temperature, the radiative contribution to the cooling dynamics can be neglected.Moreover, the A and B curves overlap, at least when the pillar temperature is below 380 K; this implies that modeling air as a stationary medium and heat dissipation only by means of conduction is a  This delay confirms the hypothesis that a laminar convective heat exchange doesn't occur for the considered geometry and temperature.The well-known theory of convection, which well describes the heat exchange from hot macroscopic objects in air, underestimates the real cooling rate in the micron and sub-micron dimensional range.In fact, for sub-micron-scaled devices and relatively small thermal gradient, recent studies report that conduction is the dominant mechanism of heat transfer [12], even with surrounding air considered as static fluid.The presented fully coupled thermo-fluiddynamic FEA of a single silicon rod immersed in air (case A) confirms the preliminary results presented in [12]: for sub-micron range, air can be considered as a static medium where the heat is dissipated principally by conduction.Moreover, the simple scaling of classical convective coefficient produces incorrect results, underestimating the efficiency of the heat exchange.From these considerations it can be concluded that the model B is an enough accurate approximation of the real heat exchange mechanism which is rigorously described by the fully coupled thermo-fluidic model A. Therefore, in order to reduce the computational time without significantly compromising the accuracy of the simulation, for the thermal analysis performed on such sub-micrometric structures, we will adopt the model B, neglecting natural convective and radiative heat exchange mechanisms.

Band structure analysis
Band structure and equi-frequency surface analyses are fundamental research tools in the design of PhC devices.In this preliminary work, two kinds of PhC structures have been investigated through the band structure and the equi-frequency surfaces analyses, carried out with Bandsolve® software (https://optics.synopsys.com/rsoft/rsoft-passive-device-bandsolve.html).The first device follows the structure in [6] and [7], where a PhC with air holes in silicon switches between a normal and a negative refractive state, by changing the material refractive index.Exploiting the Pendellösung phenomenon in PhC slab, it is possible to choose suitably structural parameters in order to select an effective refractive status [3,6].Consider a PhC slab of air holes in silicon in square arrangement, with hole radius of 123 nm and period of 640 nm, fabricated in a slab, which is 1.5 μm thick and 12 periods wide, sketched in Fig. 4a.By increasing the silicon refractive index of Δn = 0.06, band structure of the PhC changes, as shown in Fig. 4b.Moreover, in the same Fig.4b the equi-frequency surfaces at wavelength of 1550 nm at lower and higher refractive index are reported.The sketch shows the different refraction behavior at the same frequency and incidence angle.
The second structure is an array of silicon pillars in square arrangement with radius of 123 nm and period 640 nm, as sketched in Fig. 5a.The band structure shows a photonic band gap useful to guide the light inside a linear defect: for frequencies inside the band gap range light cannot propagate inside the PhC, except in the linear defects.However, the band gap width and position are a function of the refractive index contrast between the two constituent materials (silicon and air).Figure 5b shows the changes in the band gap if the refractive index is increased of Δn = 0.08.Considering the wavelength on the top edge of the band gap, the light is not guided by the defect at the corner and it is possible to obtain a switching behavior.

Results and Discussions
As established in the preliminary analysis, the only conductive exchange has to be taken in account to calculate the cooling time and the bandwidth of two optical devices.In particular, performances and efficiency of the two PhC devices previously investigated have been   4a), by increasing the Si refractive index of Δn = 0.06, FDTD simulations show that the refractive status of the PhC structure changes, as shown in Fig. 6a and b.Exploiting the thermo-optical effect [5], and considering that at wavelength of 1.55 μm the silicon thermo-optic coefficient is δn/δT = 2.4 × 10 -4 K -1 [4], the required increase of refractive index is obtained with a temperature increasing of ΔT = 250 K. Assuming that the heating of a submicrometric device can be induced in a very efficient way, for instance using an integrated resistor conveniently polarized or a laser pulse of appropriate wavelength, the bottleneck of such component is the cooling phase, especially if not assisted by active systems (forced convection, Peltier cooling modules…).
Cooling dynamics simulations, according to the model described before, have been performed in order to calculate the typical thermal transient time and, as its consequence, the maximum switching frequency (Fig. 6c).In this case the whole PhC structure has been simulated, by using the mesh shown in Fig. 6d, where only four holes are depicted for sake of simplicity.
The second device with a Si pillars in air structure exploits the well knows band gap proprieties of photonic crystals.In a T-shaped circuit, as depicted in Fig. 5a, monochromatic light is strongly confined into the path if his wavelength is included into the band gap range [8].As previously explained, the band gap range is a function of the refractive index contrast between the two constituent materials (silicon and air).FDTD simulations show the influence on the optical path of a silicon refractive index increasing of Δn = 0.08 (see Fig. 7a and b).In this case, the increase of refractive index can be obtained with a temperature increase ΔT = 300 K.
In this case, using the condition of predominant conduction in the heat transfer, simulations show (Fig. 7c) that the second structure has smaller thermal transient time than the first (about 2 μs versus 10 μs), in despite of a greater required temperature increase.

Conclusions
PhC devices can be controlled by changing their temperature: it is not particularly difficult to control it in an extremely small region (few tens of μm 2 ).The optical efficiency of two analyzed devices is up to 80 % and the value of the bandwidth of both devices is in the range of hundreds of kHz.

A
Conjugate heat transfer (conduction in solids, modeling the convective flow of air surrounding the pillar to describe the thermal dissipation); B Conduction in solids and in fluid (neglecting convection mechanism, considering air as a stationary fluid); C Conduction in solids and convection in fluid (using heat transfer coefficients h c at surface to describe the thermal dissipation in fluid); D Conduction in solids and radiation (neglecting convection in fluid).

5 hFig. 1
Fig. 1 3D geometry of the elementary cell of a pillar-based PhC

Fig. 2
Fig. 2 Comparison between literature values of h c parameter and FEA results for micrometric structures

Fig. 3
Fig. 3 Comparison among cooling dynamics for different heat transfer mechanisms

Fig. 5 aFig. 6 Fig. 7
Fig. 5 a: Sketch and dimensions of a PhC made of silicon pillars in air.b: Band structures and detail of the band gap tuning

Table 1
dimensions of PhC elementary cell