1.

INTRODUCTION

Optical microring resonators (ORRs) are circular waveguide structures with sub-wavelength thicknesses and diameters in the µm-range. High index resonator structures are embedded into low index dielectric claddings. Those devices enable a variety of biological [1,2], chemical [3,4] and physical [5,6] sensing applications with unique spatial resolution. A prominent physical application is the precise and robust measurement of temperature [7,8]. A typical configuration is the all-pass-filter [9], where light is evanescently coupled from a straight bus waveguide into the resonator ring [10]. Resonances in the ring lead to anti-resonances in the outcoupled light. As the spectral position of those resonances is shifted by temperature changes via the thermo-optical effect, probing this position gives access to the temperature of the ORR.

Refractive index pertubations, caused by optical absorption processes, lead to a superimposed resonance distortion. Created charge carriers couple to the dielectric function of the ring material by electro-optical effects such as bandgap shrinkage [11] and band filling [12]. Further, excess energies during the absorption and non-radiative recombinations insert heat into the resonator ring. Both these mechanisms lead to a superimposed resonance shift. Additionally, optical bistability effects [13], induced by the two-photon absorption (TPA) process [14], cause a strong asymmetry in the resonance shape [15] for very high intensities.

In this work, we derive analytical formulas for the self-heating effect in air-cladded ORRs. For that purpose, the heat equation with its corresponding boundary conditions and an absorption induced heat source is integrated on the ORR domain. Although not all occurring integrals can be expressed in terms of elementary functions, the explicit formulas for the spatially and time resolved temperature fields can be evaluated with substantially reduced computational costs compared to a rigorous numerical solution of the boundary value problem. Furthermore, physical dependencies manifest in a direct way here. For a comparison and verification of the analytical results, rigorous finite element simulations are performed.

2.

THE ANALYTICAL SELF-HEATING MODEL

Geometry of optical ring resonators

We investigate ORRs in the all-pass-filter configuration [9], as depicted in figure 1. A high-index resonator ring of width b, height d and diameter D is evanescently coupled to a straight bus waveguide through a gap g. The waveguides are placed on a low refractive index bottom cladding, typically made of silica glass, and are covered by a top cladding. In the scope of this work, the top cladding is air. For the formulation of the boundary value problem, see section 2.2, a Cartesian coordinate system (x₁, x₂, x₃) is utilized as shown in figure 1.

Figure 1.

All-pass-filter configuration of an optical ring resonator device. The ring and bus waveguides are placed on a low refractive index bottom cladding and are covered by a top cladding.

The boundary value problem

Although photon energies far below the band gap of the resonator material are typically used, guided light in optical ring resonators gets partially absorbed. Linear absorption processes that include defect states and TPA lead to the creation of free charge carriers [15]. The resulting change in the carrier concentration directly influences the resonance conditions by electro-optical effects [11,12]. While the excess energy of the absorption process is directly transferred into heat, recombinations occur after respective lifetimes τ. In indirect semiconductor materials like crystalline silicon, nearly all recombinations are non-radiative [16]. Thus, the annihilation energy is also transferred into heat. Charge carrier lifetimes in silicon microring resonators are in the order of τ = 1 ns [17], so that the time delay between absorption and heat creation is neglected here.

The following self-heating model is an extension to a previous calculation for thermal equilibrium, presented in Weituschat et al. [18], to the full heating dynamics. In order to formulate the boundary value problem that defines the temperature field, the following approximations are pre-determined:

• (I) The bottom cladding region Ω₁ is very large and set to be equal to the lower Euclidean half-space, i.e. Ω₁ = ℝ³_.

• (II) While the bottom cladding is typically made of silica glass [19], the resonator ring consists of crystalline materials such as silicon. Thus, the thermal diffusivity a in the former is much smaller than in the latter, i.e. a₁ ≪ a₂ (see table 1).

• (III) A critical self-heating is expected exclusively on resonance due to the high resonator finesse. Here, the intensity distribution I(x, t) in the ring is overall much higher than in the bus waveguide, which is consequently neglected. Further, I is assumed to be spatially homogeneous.

• (IV) We investigate air-cladded devices, making the heat flow into the top cladding negligible.

Table 1.

Input parameters for the heating calculations.

Quantity	Symbol	Unit	Value	Source
diameter	D	nm	19130	design
width	b	nm	500	design
height	d	nm	220	design
thermal conductivity cladding (SiO2)	λQ1		1.3	[28]
thermal conductivity ring	λQ2		130	[29]
linear attenuation coefficient	α		6.9	measured
TPA coefficient	β		8 × 10−12	[30]
thermal contact conductance	ℎ12		109	[31]
mass density cladding (SiO2)	ρ1		2200	[32]
heat capacity cladding (SiO2)	cm1		1000	[33]
loss factor laser → chip	C1	1	0.25	measured
coupling efficiency	C2	1	0.01	simulated
Q-factor	Q	1	2.5 × 105	measured
free spectral range	FSR	nm	9.575	measured
full width at half maximum	FWHM	pm	7	measured
temperature sensitivity	dλ/dT	pm/K	68	[18]
operation wavelength	λ	nm	1546.6	measured

In the whole device Ω, the temperature field T(x, t) satisfies the heat equation

where f is the local volumetric heat source [20]. According to approximation (IV), the whole inserted heat into the ring is transferred to the bottom cladding via the interface in between. Further, the heating time is predominantly determined by the cladding due its much lower thermal diffusivity, see table 1. Due to the very thin waveguide width compared to the ring diameter (b ≪ D, see table 1), the heat source f₁ utilized for the heating calculation in the bottom cladding is approximated as a delta distributed ring apart from the interface:

with the constant light intensity in the ring I, the linear attenuation coefficient α, the TPA coefficient β, the cladding mass density ρ₁, the cladding heat capacity c_m1 and cylindrical coordinates (ϱ, ϕ, x₃). Both linear absorption processes from defect states and TPA are considered here, leading to αI + βI² for the absorbed power density [21]. A very large bottom cladding ℝ³_, surrounded by air with negligible heat release, fulfills insulating boundary conditions (BCs) [22]. That leads to the Neumann problem:

determining the temperature field T₁(x, t) in the cladding region apart from the interface.

Due to the high heating rate in the ring region and the high thermal contact conductance (see table 1), a quasi-equilibrium temperature distribution is reached within few nanoseconds, as shown by the following estimation. The temporal heating behavior in the resonator ring is determined by the law of conservation of energy:

Here, ρ₂ is the ring mass density, C_m2 is the heat capacity and h₁₂ is the thermal contact conductance. Further, πDbd is the ring volume and T₁(x, t)|_Interface is the cladding temperature at the interface. Thus, the left side of equation (4) is the total absorbed power. That power is equal to the sum of the inserted heat in the ring per time and the heat flow through the interface per time, representing the right side. Considering a much slower heating rate in the large cladding, the solution of equation (4) is a saturating exponential increase of the ring temperature:

For typical silicon ORR systems (see table 1), the corresponding heating time t_rise = dρ₂C_m2/h₁₂ is equal to 0.4 ns. Thus, the fast heating period in the ring region is neglected and a quasi-equilibrium temperature field T₂(x, t) is calculated, following the much slower heating process (µs-range, see figure 4) in the cladding. For that purpose, the cylindrical coordinate ϱ is shifted by −D/2 + b/2, resulting in a new radial coordinate . Due to the cylindrical symmetry, the problem is effectively two-dimensional with the waveguide cross-section as its rectangular domain. On the lower boundary x₃ = 0, a Dirichlet condition T₂₀(x, t) as well as a Neumann condition are enforced by the cladding temperature distribution T₁(x, t)|_Interface and the thermal contact conductance at the interface [22]. On all other boundaries, exclusively insulating Neumann BCs are applied due to the air top cladding. The resulting mixed BC Poisson equation reads as

with the constant heat source

Solution for the temperature field

First, the bottom cladding temperature field T₁(x, t) according to equation (3) is calculated. The solution of that boundary value problem is a convolution of the heat source f₁(x, t) and a fundamental solution G(x, x′, t) for the respective BCs [20]:

Here, the fundamental solution is the Greens function for the lower Euclidean half-space and insulating Neumann BCs:

where is the reflection vector with respect to the -plane. Due to f₁(x) ∝ δ[x₃], see equation (2), exclusively values of G at contribute to the convolution (8). Consequently, the distance in cylindrical coordinates fulfills over the whole integration range. Using that distance in equation (8), the ϱ′- and x₃′-integration lead to

While the φ′-integral represents a hyperbolic Bessel function of first kind I₀ [23], the remaining time integral cannot be expressed by elementary functions:

Thus, the temperature increase in the bottom cladding ∆T₁ = T₁ − T₀ is proportional to the total absorbed power (αI + βI²)πDbd and anti-proportional to the thermal conductivity λ_Q1 = ρ₁c_m1a₁. It asymptotically reaches zero for large distances, i.e. . In thermal equilibrium, i.e. t → ∞, the cladding temperature field can be further simplified. For that purpose, the mathematical identity [24]

is utilized, where ℜ denotes the real part of the complete elliptic integral of first kind K. Thus, the thermal equilibrium temperature distribution T₁(x, ∞) reads as

Directly at the interface, an approximately homogeneous temperature distribution emerges:

In order to obtain the temperature field T₂(x t) in the resonator ring from equation (6), the lower boundary temperature T₂₀(x t) must be determined. For that purpose, an energy balance for the heat flow through the interface between ring and bottom cladding is formulated, considering a quasi-equilibrium distribution in the former one:

Here, the left side is the total absorbed power. In thermal equilibrium, that inserted heat completely flows through the interface to the bottom cladding. The corresponding heat transfer per time is represented by the right side of equation (15). Hereby, T₂₀(x, t) is determined by the bottom cladding temperature field T₁(x, t)|_x3=0 according to equation (14). Applying the superposition principle, the general solution of (6) is the sum of a particular solution of the inhomogeneous Poisson equation and the general solution of the homogeneous Laplace equation [25]. A particular solution that fulfills all Neumann BCs is

As holds at the interface here, the BCs for the Laplace problem of the homogeneous solution remain unchanged, except , compare equation (6). A product separation ansatz leads to eigenvalue problems for the radial and x₃-dependencies. The general homogeneous solution is a superposition of all eigenstates, expressed by the Fourier series [26]

The constant Dirichlet condition at the interface enforces c₀(t) = T₂₀(t), while all higher Fourier coefficients equal zero. Thus, the resulting temperature field within the ring resonator reads as

Here, the first summand describes the parabolic temperature distribution in x₃-direction, the second summand is the temperature step at the interface and the third summand is the cladding temperature.

3.

THE FINITE ELEMENT MODEL

For the numerical computation, the software Comsol Multiphysics [27] is utilized. Only thermal equilibrium temperature distributions, i.e. t → ∞, are determined here. Since the ring is the primary heat source when driven at the resonance wavelength, we can simplify the computation by using axis symmetry. Hereby, the geometric model consists of the silicon ring waveguide with width b, height d and diameter D (see table 1). The ring is surrounded by two half-spheres, as illustrated in figure 2: one for the upper air cladding and the one for the SiO_x bottom cladding beneath the waveguide. These two half-spheres are limiting the optical computation, using scattering boundary conditions. Another spherical domain is added, that is used to simulate an infinite SiO_x substrate as heat sink by imposing an infinite-element-domain condition with a fixed temperature boundary of T₀. The width of this domain is equal to 10 % of the radius R_clad of the entire system. The radius of both cladding domains is determined by sweeping the radius until the temperature of the ring waveguide converges to a constant value. This radius was found to be R_clad ≈ 3 · D. Thermal radiation and convection at the air-solid-interfaces are neglected. A triangular mesh is used with a maximum element size of λ/30 within and λ/12 outside the waveguide core with a growth rate of 1.1. As a first simulation step, an eigenfrequency study is conducted to determine the exact optical resonance wavelength of the system, which is found to be at λ = 1546.6 nm. Afterwards, a user-defined port combined with a frequency-domain study is used to calculate the electromagnetic field distribution for several input powers. The intensity I within the ring is calculated from the electric field E as

Figure 2.

Sketch of the axis-symmetric model for the FEM computation. The inner sphere with radius R_clad marks the optical computation area. The outer sphere with radius R_inf utilizes the infinite-element-domain property to simulate an infinite SiO_x substrate with a fixed temperature of T₀.

where n₂ is the refractive index of the ring, ε₀ is the vacuum permittivity and c is the vacuum speed of light. Cylindrical coordinates (ϱ, φ, x₃) are utilized. In thermal equilibrium, the time derivative in the heat equation vanishes, leading to

The heat source f(x) is defined just like in equation (7), but with a spatially dependent intensity distribution I(x) according to the electric field E(x).

4.

SELF-HEATING IN SILICON RESONATORS - RESULTS

Equations (11), (14) and (18) describe the spatially and time resolved temperature field in the whole ORR device analytically. The heating behavior is exemplarily calculated for a crystalline silicon resonator on silica glass with air as top cladding. All input parameters are summarized in table 1. This ORR was manufactured as part of a multi-project wafer run. Therefore, the geometric dimensions are the target design parameters without manufacturing design tolerances,

For a handy comparison between the measurable spectral shift ∆λ of the resonance and the calculated self-heating, both are determined as functions of the experimental accessible laser input power P₀ and wavelength sweep rate dλ/dt. The intensity I in the resonator is calculated as

where C₁C₂ represents the portion of light from the laser that is coupled into the resonator ring, and Q × FSR/λ is the resonator finesse. Further, the excitation time t is determined by the FWHM of the resonance and the wavelength sweep rate as

Figure 3 shows the analytically calculated radial temperature distribution in the bottom cladding in thermal equilibrium according to equation (13), i.e. ∆T₁(ϱ, 0, ∞). An input laser power P₀ = 0.4 mW is utilized here, leading to a maximal heating amplitude of ∆T₁(D/2,0, ∞) ≈ 19 mK at the interface. For comparison, the FEM computed result is additionally depicted in figure 3 (b) (dashed line). The associated heating amplitude is with about 22.5 mK slightly higher than the analytically calculated one.

Figure 3.

(a) Analytically calculated temperature field ∆T₁(x₁, x₂, 0, ∞) on the upper edge of the bottom cladding in thermal equilibrium. The input power P₀ is equal to 0.4 mW here. (b) Analytically calculated (solid line) and FEM computed (dashed line) radial temperature profile ∆T₁(ϱ, 0, ∞).

Figure 4 shows the analytically calculated resonator temperature increase as a function of time according to equation (18), i.e. ∆T₂(D/2,0, t). Additionally, the corresponding wavelength sweep rate and the resulting spectral resonance shift are depicted on the upper and on the right axis, respectively. Again, P₀ = 0.4 mW is utilized. It turns out that thermal response times of the ORR device under continuous wave excitation are in the µs-range: the excitation time for reaching half of the equilibrium temperature increase is t_0.5 ≈ 13 µs.

Figure 4.

Analytically calculated temperature increase on the lower edge of the resonator ring ∆T₂(D/2,0, t). The input power P₀ is equal to 0.4 mW here.

Figure 5 shows the analytically calculated (solid lines) power dependence of the resonator temperature increase ∆T₂(D/2,0, t) for three different wavelength sweep rates according to equation (18). Again, the resulting spectral resonance shift is depicted on the right axis. Slower wavelength sweeping leads to a stronger heating effect, as the heating period gets longer. For relatively low excitation powers, i.e. P₀ < P_th, linear absorption processes are the dominant heat source. Here, the temperature increase is proportional to the input power. For sufficiently high powers, i.e. P₀ > P_th, TPA gets dominant and hence ∆T₂ rises quadratically with P₀. The threshold between those two regimes is calculated from the heat source (2) and intensity (21) to be P_th,ana ≈ 22 mW.

Figure 5.

FEM computed (dashed line) and analytically calculated (solid lines) temperature increase on the lower edge of the resonator ring ∆T₂(D/2,0, t) as a function of input power P₀ for three different wavelength sweep rates dλ/dt.

For comparison, the respective FEM computed power dependence is depicted in figure 5 as well (dashed line). That simulation result is referring to thermal equilibrium, i.e. dλ/dt → 0. As already shown in figure 3 (b), the simulated heating amplitude is slightly higher than the analytically calculated one in the linear regime. That deviation increases in the TPA dominated regime. Thus, the FEM computed threshold is with about P_th,FEM ≈ 14 mW significantly lower than P_th,ana. The overestimation of P_th in our analytical model is a result of approximation (III) (see section 2.2): a homogeneous intensity distribution within the ring waveguide leads to lower peak intensities, and thus to a lower threshold. In general, P_th is determined by an equality of the total absorbed power due to linear and quadratic absorption:

The FEM computed electric field distribution E(ϱ, x₃) ∝ I^1/2(ϱ, x₃) is depicted in figure 6. Utilizing this simulation result, a homogeneous distribution for the analytical model and equation (23) leads to the ratio P_th,ana/P_th,FEM:

The threshold between the linear absorption and the TPA dominated regimes marks an important transition in ORR systems, as optical bistability effects are caused by TPA [13-15].

Figure 6.

FEM computed electric field distribution E(ϱ, x₃) on resonance within the resonator ring cross-section for an input power P₀ = 0.4 mW.

5.

CONCLUSION AND OUTLOOK

In this work, we investigate the dynamical self-heating effect in optically excited optical microring resonators (ORRs). For that purpose, the heat equation boundary value problem under continuous excitation is solved in air cladded ORRs. We derive completely analytical expressions for the spatially and time dependent temperature fields. The self-heating effect is explicitly calculated for a widely applied ORR technology platform [9,34,35]: a silicon ring resonator on a silica glass bottom cladding. More explicitly, this ORR was produced using a multi-project wafer prototyping process. Experimentally determined parameters (see table 1) were included in the calculations. It is shown that the heating amplitude, which is about 50 mK for an input power of 1 mW, is mainly determined by the bottom cladding due to its small thermal conductivity. Excitation times for reaching half of the equilibrium temperature increase are in the range of a few µs, and therefore about 10⁴ times larger than typical charge carrier lifetimes in the resonator ring [17]. As expected, the heating amplitude rises linearly with the input power in the linear absorption dominated and quadratically in the TPA dominated regime. In the latter power range (here above a threshold of about 14 mW), optical bistability effects cause significant asymmetries in the resonance shape [13-15].

The analytical results are verified by rigorous finite element simulations (FEM) in thermal equilibrium. For this purpose, the software Comsol Multiphysics [27] is utilized. Both analytical and numerical results are in good agreement down to a deviation of about 20 % in the heating amplitude and of about 55 % in the TPA threshold. Those deviations are due to a simplification in the analytical model, namely a homogeneous intensity distribution in the ring. The derived equations enable the determination of the whole temperature field with drastically reduced computational effort compared to FEM simulations, which enables fast scans over the whole parameter space. In future, the derived analytical model shall be extended by charge carrier dynamics and optical bistability effects. That extension would enable the modelling of absorption induced resonance distortions at low computational costs.