Microresonator Effective Thermal Parameters Deﬁnition via Thermal Modes Decomposition

: High-Q optical microresonators are particularly efﬁcient practical tools of modern applied optics and photonics. Using them, one inevitably faces the problem of thermal effects. Accurate determination of effective thermal parameters of high-Q microresonators (effective thermal relaxation rate and optical absorption rate) is of particular importance for developing microresonator-based devices. Our investigation looks into diverse methodologies to estimate these effective parameters for such systems, ultimately revealing a divergence between the commonly employed simpliﬁed model, the direct numerical approach, and classical analytical formulas. We introduce a novel approach to calculate effective parameters based on the decomposition of the thermal ﬁeld into microresonator thermal modes, which inherently considers the intricate geometry and material anisotropy inherent in microresonators, as well as the inﬂuence of external conditions. The method for the accurate determination of the effective thermal parameters of the microresonator for corresponding thermal modes is developed. As a result of applying this method, we modiﬁed the classical approach for the simulation of thermal effects in optical microresonators for better agreement with the numerical simulations. By accounting for the complexities of microresonator shapes, material properties, and external factors, our proposed method contributes to a more accurate understanding of thermal dynamics and enhances the predictive capabilities of simulations for these systems. We demonstrated the application of this method on the example of integrated microring resonators, but it can be used to analyze thermal effects in other microresonator platforms.

Thus, for accurate simulation of linear and nonlinear processes in optical microresonators and the development of effective microresonator-based devices, it is essential to estimate the effective microresonator thermal parameters, such as effective thermal relaxation rate and effective optical absorption rate (and, consequently, effective thermal nonlinearity value), defining the dynamics of thermal processes in such structures.Especially, this issue is of paramount importance for the compact on-chip structures widely studied and developed nowadays.
Integrated photonic platforms find their foundation in a variety of materials, including LiNbO 3 [50][51][52][53], Si [54][55][56], AlGaAs [57][58][59][60], GaP [61], TiO 2 [62], and others.However, the noticeable progress in recent decades in the technology of manufacturing integrated silicon nitride (Si 3 N 4 )-based photonic structures has led to increasing utilization of this platform in addressing actual problems in optics and photonics.Such a variety of devices based on Si 3 N 4 integrated technology is beneficial for silicon nitride low losses in the telecommunications wavelength range, the absence of two-photon absorption, manufacturing reproducibility, and high nonlinearity [63][64][65].It is also shown in [66] that irradiation with high-energy protons does not have a long-term effect on the linear optical losses of silicon nitride structures, which opens up prospects for space applications.
One of the most promising types of integrated structures are high-Q optical microresonators, which are widely used to create compact, highly stable lasers [67][68][69] and compact generators of optical frequency combs [70][71][72][73].However, light absorption leads to inevitable heating and, consequently, to thermal effects [11][12][13]50], such as the thermo-optical and thermal expansion effects.Previously, it was shown that these effects strongly influence the dynamics of various processes, including nonlinear processes [24,74] and thermomechanical fluctuations in optical microresonators [14,18].Therefore, in order to accurately model these processes, it is necessary to describe the dynamics of thermal processes accurately.
In numerous studies, the rate equation is employed to represent microresonator thermal dynamics, consequently determining the thermal shifts of microresonator eigenfrequencies [21,22,24,50].Nevertheless, our prior investigation [75] highlighted that such an approach might need more accuracy under specific conditions.For a more precise description of the thermal dynamics, we propose a universal approach based on the decomposition of the thermal field into microresonator thermal modes [76,77].We demonstrated the application of this method on the example of integrated microring resonators, but it can be used to analyze thermal effects in other microresonator platforms as well.

Modeling of Thermal Effects in Optical Microresonators
Let us consider the description of thermal effects in optical microresonators using the example of an on-chip microring resonator (see Figure 1).The electromagnetic field in microresonators is confined to a small volume compared to the entire microchip.The heat conduction equation for the integrated microresonator can be written as [13] where D T = k C p ρ is the thermal diffusion coefficient, k is the thermal conductivity, C p is the heat capacity, ρ is the density of the microresonator material.Q( r) is the power of the heat sources per unit volume for the m-th microresonator mode, which is equal to the optical mode power density absorbed in the waveguide Q |e m ( r)| 2 dV , P is the part of the optical power coupled into the m-th mode that is absorbed in the microresonator material, leading to its heating (see color map in Figure 1), T 0 is the ambient temperature, and h air is the heat transfer coefficient ≈ 10 W/(m 2 •K) for convective air cooling.S t and S b are the top and bottom surfaces of the microchip, S s is the side surface, and e m ( r) is the microresonator m-th mode profile.The complexity of this system, represented by Equation ( 1), makes it challenging to effectively simulate the dynamics of an electromagnetic field while considering the influence of the thermal effects in microresonators.Thus, usually, simplified semi-empirical rate equations are used.In the following, we briefly consider their introduction and thermal parameters for determination.Color map shows heat source distribution Q(r) inside the waveguide.R 0 is the radius of the microresonator.We simulated microresonators with a radius from 24 µm to 840 µm.The transverse dimensions dr and dz of the waveguide were varied in the range from 20 nm to 2000 nm.
As an example, we performed finite element modeling (FEM) using COMSOL Multiphysics of temperature dynamics during constant heating for the configuration presented in Figure 1 with default material and geometrical parameters specified in Tables 1 and 2. Obtained results for P = 10 mW heating of 24 µm ring resonator with waveguide width dr = 1 µm and height dz = 0.8 µm are shown in Figure 2a.

Analytical Formulas
The temperature change δT = (T(t) − T 0 ) from the room temperature T 0 caused by the intracavity heating power P is usually calculated with the rate equation [24,27,50].If we average (1) by multiplying by the square modulus of the field of the electromagnetic mode of the microresonator |e m ( r)| 2 and integrating by the volume of the microresonator, using Green's formula for the second term of the equation, we can obtain the following expression (also known as the rate equation) [13]: where γ = 1 . Effective thermal relaxation rate δ θ and effective optical absorption rate γ depend on the microresonator configuration.Effective thermal relaxation rate determines the time of establishment of thermal equilibrium.The combination of δ θ and γ determines the magnitude of the temperature change and, consequently, thermal mode shift value and hence the magnitude of the effective thermal nonlinearity.
When heating occurs due to the absorption of the electromagnetic field of the m-th mode excited in the microresonator by external laser, P can be written as where A m is slowly varying amplitude of m th mode, Q abs is a part of Q-factor, determined by the heat-related absorption, ω m is microresonator circular eigenfrequency of the pumped mode, and n 0 is refractive index of the microresonator material and 0 is vacuum dielectric constant.In order to convert the temperature change δT into a pumped mode resonance frequency shift δω, it is possible to use the relation: Then, passing from the temperature change in (2) to the thermal frequency shift, we can write the following equation: This equation is similar to the one usually used for the analysis of nonlinear processes in Kerr medium.Thus, it allows to compose an effective thermal nonlinearity coefficient n 2T depending on δ θ and γ: where α is absorption coefficient in 1/m (Q abs = n 0 ω m αc , c is speed of light).It is convenient to compare this parameter with the value of the Kerr nonlinearity n 2 in order to understand which nonlinear mechanism will prevail [21,22,27,46,74].Note that additional spectral components may appear in Kerr microresonators due to the process of four-wave mixing [6,7,22] that should be taken into account in (3).
Using Gaussian approximation of optical fields in δ θ and γ from Equation (2), one can obtain the estimations for thermal parameters for the considered microresonator: where r r and r z are mode transverse profile half-width and half-height, respectively, V ≈ 2π 2 r r r z R 0 is mode volume, R 0 is the microresonator radius.In the case of a ring waveguide, the mode dimensions can be well-approximated with the guide cross-section half-widths r r ≈ dr/2 and r z ≈ dz/2.

Thermal Dynamics Approximation by the Sum of Exponents
We considered the model of microresonator thermal relaxation due to thermal conductivity.However, suppose the power changes or the pump frequency is tuned very slowly.In that case, besides thermal conductivity in Si 3 N 4 , microring other processes may also appear, in particular thermal relaxation of the entire microchip as a whole due to air convection and thermal conductivity of the silicon substrate [11].Thus, there may be several thermal relaxation times associated with processes in different thermal subsystems of microchip.This leads to the appearance of the additional thermal effective parameters, and, as a consequence, for a more accurate description of the thermal dynamics (see Figure 2), it is necessary to use several rate equations [24,27,50]: so that δT = ∑ δT (k) is the sum of partial modal contributions δT (k) .Then, δ θ and γ (k) are, respectively, the effective thermal relaxation rate and the effective optical absorption rate of the k th thermal mode.
Further, we compare two theoretical methods for the determination of the effective microresonator thermal parameters from the data obtained by the numerical solution of (1) for the setup shown in Figure 1.Within the first method, we obtain the effective parameters by approximating the temperature dependence obtained numerically (solid line in Figure 2a) by the sum of the formal solution of ( 9) with a step-like excitation: where initial temperature deviation δT 0 = 0.As result, we obtain the effective γ (k) and δ θ .The coefficient of determination r 2 is used to estimate the accuracy of the approximation method: T i (t i ) is the dataset of the temperature from numerical simulation, where each T i is associated with the time moment t i , and T is the mean value of dataset T i (t i ) used for the normalization.The fitted value of the temperature T(t i ) was obtained using a sum of the expressions with effective parameters (10).For the range of parameters where the approximation method is in good agreement with the numerical simulation data, r 2 tends to 1 (see Figure 2c), and, for the range with poor agreement, r 2 decreases.This effect is explained by the fact that the actual dependence of the temperature on time during heating is a more complex function than the solution of (2) that is an exponential dependence.

Thermal Mode Decomposition Method
Within the second method, we use the fact that the general solution of the system (1) can be represented as decomposition into the microresonator thermal modes [77]: where T (k) 0 and C (k) ( r) are amplitude and spatial shape of the thermal mode, δ θ is eigenvalues corresponding to the times of thermal relaxation of various "thermal subsystems".
We should note that formal mode eigenvalues δ (k) θ introduced here and thermal relaxation rates introduced in previous subsection are not the same.However, as we obtain the same form of equations for the same system, it looks like they should tend to each other with increasing the number of exponents in the fit.In the following, we study this assumption.
In order to obtain C (k) ( r) and δ (k) θ , using COMSOL Multiphysics, we numerically solved the eigenvalue problem in the thermal domain for the microring resonator shown in Figure 1: As a result, we obtained mode shapes C (k) ( r) (normalized to 1 Kelvin) and eigenvalues θ (see Figure 3).We substituted the solution (12) into Equation ( 1) and multiplied the resulting expression by the normalized mode shape and integrated over volume.We numerically checked that thermal modes form an orthogonal basis such that where δ km is the Kronecker symbol.Using the orthogonality of thermal modes, we obtained the rate equations for the thermal modes similar to (9): where It can be seen that the expression on the right side ( 14) coincides in dimension with γ (k) in Equation ( 9): We calculated the expression (15) for each thermal mode shape C (k) ( r).Examples of shapes C (k) ( r) in Figure 3a,b correspond to the thermal mode with the maximum overlap integral , and (c), (d) correspond to the slowest thermal mode of the cladding subsystem.To obtain the exponent amplitude, we normalized the expression (15) to δ (k) θ .The resulting spectrum of exponent amplitude γ k /δ k obtained from (15) is shown by the stem line in Figure 4. We studied the dependence of the effective parameters on the number N of exponents.In Figure 4, the number of exponents and the corresponding effective parameters are indicated by different symbols (on the horizontal axis-effective thermal relaxation rate δ θ , on the vertical axis-the amplitude γ/δ θ of the exponent from Equation (10), δT 0 = 0, P = 1 W).
It can be seen in Figure 4 that the effective thermal relaxation rate corresponding to the maximum of the exponent amplitude spectrum γ (k) /δ θ obtained by the thermal mode decomposition method (stem lines) is very close to the effective thermal relaxation rate for the case of one-exponent approximation (N = 1, lime square).However, in such a spectrum, several other peaks are observed that are comparable in amplitude with the maximum one (for example, at δ θ ≈ 10 6 ).Nonetheless, the sum of all amplitudes γ (k) /δ (k) θ for all thermal modes is equal to the amplitude for the one-exponent approximation.This fact demonstrates the fulfillment of the energy conservation law.Note that, when we decompose into thermal modes, the coefficient of thermal nonlinearity n 2T is also split into the set of n  θ obtained by the thermal mode decomposition method for microring resonator with default material and geometrical parameters specified in Tables 1 and 2; different symbols represent the amplitudes for exponential approximation (10) with the certain amount of exponents from 1 to 15, red line with red circle-theoretical effective thermal relaxation rate and exponent amplitude obtained from Equations ( 7) and ( 8), red dotted line with red rhombus-theoretical effective thermal relaxation rate and exponent amplitude obtained from Equations ( 16) and ( 8), red dashed and dash-dotted lines-theoretical lowest effective thermal relaxation rate of the cladding (Equation ( 17)) and silicon substrate (Equation ( 18)).
Increase in the number of exponents (from 1 to 15) leads to a gradual decomposition of the temperature dynamics into thermal modes, where the eigenvalues correspond to effective thermal relaxation rates in exponents, and the amplitudes of the spectral components (Figure 4, stem lines) correlate with the amplitude of the exponent.However, after N = 9 exponents (Figure 4, green stars), the values of the effective thermal relaxation rates δ (k) θ and exponent amplitudes stop changing.Thus, convergence in the number of exponents can be a mathematically rigorous criterion for the optimal number of exponents.
In order to test our thermal mode decomposition, we reconstructed the temperature dynamics from the exponent amplitude spectrum (see stem lines in Figure 4), substituting into Equation (10) thermal mode eigenvalues δ (k) θ and γ (k) (see Equation ( 15)) for these thermal modes.As a result, we obtained a good agreement between the reconstructed dynamics and the simulation data (see Figure 5).Thus, two independent approaches to thermal mode decomposition converge to the same result.

Thermal Modes Existence Range
We studied analytically and numerically thermal modes existence range.It was revealed that one can limit the range of the effective thermal relaxation rates δ (k) θ using analytical formulas.First, it was found numerically that possible values of δ (k) θ for microresonator thermal modes (stem lines in Figure 4) are limited from above by the value δ theor θ (from Equation ( 7)) describing thermal processes in the waveguide material (solid red line with red circle in Figure 4).Also, it was noticed that the effective thermal relaxation rates used for multi-exponent approximation do not agree well with the analytical effective thermal relaxation rate (7) and corresponding exponent amplitude.It can be seen in Figure 4 that thermal modes (stem lines in Figure 4) near the theoretical value (7) (solid red line with red circle in Figure 4) have low amplitudes (less than 10 0 K/W and significantly smaller than exponent amplitude defined by Equations ( 7) and (8) shown by the red circle) and, accordingly, low contributions to the temperature dynamics.This effect is associated with the presence of thermal conductivity of the cladding, which limits the effective thermal relaxation rate.Thus, if the thermal conductivity of the cladding is substituted into the analytical Formula (7) instead of the thermal conductivity of the waveguide, then we obtained a limitation on δ (k) θ when the amplitudes of the exponents are still comparable (red dotted line with red rhombus, Figure 4): Only for the case of four exponents (purple triangles, Figure 4), the fastest effective thermal relaxation rate coincides with the analytical Formula ( 16).
Another factor influencing the limitation of the effective parameters at high thermal conductivity is the thickness of the cladding and the presence of a silicon substrate, which efficiently removes heat due to the high thermal conductivity.We found that the boundary with silicon is approximately equivalent to the boundary with a fixed temperature.We also found that the slowest effective thermal relaxation rate of the cladding subsystem is related to the cladding thickness.To calculate the slowest thermal mode of the cladding subsystem, we temporarily removed the silicon substrate from the model and replaced silicon substrate with a fixed temperature boundary T 0 .We further generalize Formulas ( 7) and ( 16) to find analytical approximation of this cut-off, substituting corresponding subsystem parameters into the formula.In this way, for cladding, substituting the cladding thickness H clad instead of r 2 r and r 2 z , we obtain (Figure 4, dashed red line): where H clad is SiO 2 cladding thickness (see Figure 1).Doing the same for the substrate, we obtain the lowest relaxation rate of the whole system (Figure 4, dash-dotted red line): This also implies that putting the system on a bigger mount can bring more low-frequency modes.However, they will not make high impact due to low amplitudes, and the lowest significant frequency will be defined the cladding subsystem.
As a result, we obtained the region of existence of relaxations δ ; δ theor θ ] and the region where the amplitudes of the exponents are comparable (the same order of magnitude) δ

Comparison with Experiment
Our results are in good agreement with the experimentally measured values for a similar integrated microresonator [78], which indicates the correctness of the chosen model and method.In [78], effective thermal relaxation rates for different times during the cooling were measured.The data of such a measurement are fitted with the sum of two exponents in order to take two different thermal time scales into account.The faster of the two time scales is 760 ns, and the slower is 14 µs.At the same time, Figure 2c shows that further addition of the number of exponents quickly reaches 1 in approximation accuracy (r 2 ).Hence, for simple applications, where only the shape of the relaxation curve in a particular heating process is of interest, two exponents are sufficient, judging by r 2 .However, we can see that the 2-exponent case does not capture the thermal mode with the highest exponent amplitude, which is there for all other exponent numbers N. So, for more complex processes, where individual thermal eigenfrequencies can respond, this may not be enough, and exponent number convergence should be searched for.

Applicability Range and Accuracy Analysis of Microresonator Thermal Parameters Estimation Methods
As mentioned in the previous section, the effective parameters can be estimated using three methods-analytical formulas, approximation of numerical simulation, and thermal mode decomposition.We compared these methods in a wide range of parameters, including the thermal conductivity of the microresonator materials, the geometrical dimensions of the waveguide, and the position of the waveguide relative to the cladding.
In Figure 6, the dependencies of the effective thermal relaxation rate and effective optical absorption rate on the thermal conductivity of the waveguide calculated using three methods are presented.The red line indicates the analytical dependence using Formulas (7) and (8) in Figure 6a and Figure 6b, respectively.It can be seen that, with an increase in the thermal conductivity of the waveguide, the value of the effective thermal relaxation rate increases linearly.At low values of the waveguide thermal conductivity (k wg < 10 −1 W/(m•K)), all methods provide the same result that is well-described by analytics (7).However, the effective thermal relaxation rates obtained by both numerical methods, approximations by sum of exponents (Equation (10), case N = 1-squares, N = 2-crosses, and N = 7-circles) and maximal-amplitude relaxation from the thermal modes decomposition (blue stars), in contrast to the analytical ones, tend to saturate with increasing thermal conductivity and are close to the values obtained from one-exponent approximation.17), squares-one-exponent approximation (( 10), k = 0), crosses-two-exponent approximation (( 10), k = 0, 1), circles-seven-exponent approximation (( 10), N = 7).The coefficient of determination r 2 for the approximations is shown by color.Blue stars correspond to the maximum of exponent amplitude spectrum (stem lines in Figure 4); grey stars stand for the slowest thermal mode of the cladding subsystem.Vertical lines correspond to the thermal conductivity coefficients of Si 3 N 4 , SiO 2 , and Si (see Table 1).
The coefficient of determination r 2 indicated by the color code in Figure 6 was used to estimate the accuracy of the fitting method.However, with an increase in k wg , it can be seen that the effective thermal relaxation rate for the case of two (N = 2) exponents fit (crosses, Figure 6) is split at k wg ≈ 10 −1 W/(m•K) into two effective parameters, δ (1) θ and δ (2) θ .Moreover, the maximum of the exponent amplitude spectrum (blue stars, Figure 6) from stem lines in Figure 4 is between the fast and slow effective thermal relaxation rate branches and coincides with the empirical exponent.Figure 6 shows that adding the second exponent (N = 2, Equation ( 9)) increases r 2 for a real microresonator from 0.97 to 0.995.Note that numerical eigenvalues for the slowest thermal modes of the cladding subsystem-i.e., calculated without silicon substrate-are well-described by Equation ( 17) (grey stars and dashed line in Figure 6a) over a large range of the waveguide thermal conductivity k wg .Also, it is clearly seen that almost all effective thermal relaxation rate values are in the range between red dotted (Equation ( 16)) and red dashed (Equation ( 17)) lines in Figure 6a.
To confirm our proposed method and formulas, we further investigated the dependence of the effective parameters on the cladding thermal conductivity k clad .It can be seen from Figure 7 that the effective thermal relaxation rate depends almost linearly (dotted red line, Figure 7a) on the thermal conductivity k clad , which is in good agreement with the Formula ( 16) we proposed earlier.One can see that numerical eigenvalues for the slowest thermal modes of the cladding subsystem over a large range of the cladding thermal conductivity k clad are well-described by Equation ( 17) (see dashed line and grey stars in Figure 7a) and are close to the smallest relaxation rate values for the multi-exponent approximation (compare the positions of grey stars and lowest circles in Figure 7a).Also, it is clearly seen that almost all effective thermal relaxation rate values are in the range between red dotted (Equation ( 16)) and red dashed (Equation ( 17)) lines in Figure 7a.17), squares-one-exponent approximation ((10), k = 0), crosses -two-exponent approximation ((10), k = 0, 1), circles-seven-exponent approximation ((10), N = 7).The coefficient of determination r 2 for the approximations is shown by color.Blue stars correspond to the maximum of exponent amplitude spectrum (stem lines in Figure 4); grey stars stand for the slowest thermal mode of the cladding subsystem.Vertical lines correspond to the thermal conductivity coefficients of Si 3 N 4 , SiO 2 , and Si (see Table 1).
Next, to test our methods and formulas, we studied the temperature dynamics for different microchip geometries.Figure 8 shows the dependence of the effective parameters on the microresonator radius R 0 .It can be seen that the use of more exponents provides the same approximation accuracy over the entire range of radii R 0 , in contrast to the classical approach (N = 1).From Figure 8, we can conclude that the effective thermal relaxation rates do not depend on the radius of the integrated microresonator.
It should be noted that the numerical eigenvalues for the slowest thermal modes of the cladding subsystem (grey stars, Figure 8a) are well-described by Equation ( 17) over a large range of microresonator radius R 0 .Also, it is clearly seen that almost all effective thermal relaxation rate values are in the range between red dotted (Equation ( 16)) and red dashed (Equation ( 17)) lines in Figure 8a.
At the same time, for narrow waveguides (dr < dz), there is a strong discrepancy with the theoretical curve (red dotted line, Figure 9).However, when adding exponents, the discrepancy with the theoretical curves for narrow waveguides decreases.Moreover, with a decrease in the waveguide width dr, a larger number of exponents is needed for better agreement with the theory (7).The slowest effective thermal relaxation rate of the cladding subsystem is in good agreement with the first thermal mode of the cladding (grey stars, Figure 9a) and the Formula (17) (red dashed line, Figure 9a).Also, it is clearly seen that almost all effective thermal relaxation rate values are in the range between red dotted (Equation ( 16)) and red dashed (Equation ( 17)) lines in Figure 9.  17), squares-one-exponent approximation ((10), k = 0), crosses-two-exponent approximation ((10), k = 0, 1), circles-seven-exponent approximation ((10), N = 7).The coefficient of determination r 2 for the approximations is shown by color.Blue stars correspond to the maximum of exponent amplitude spectrum (stem lines in Figure 4); grey stars stand for the slowest thermal mode of the cladding subsystem.Red vertical line in (a) and (b) corresponds to the default microresonator radius (see Table 1) R 0 = 24 µm.17), squares-one-exponent approximation ((10), k = 0), crosses-two-exponent approximation ((10), k = 0, 1), circles-seven-exponent approximation ((10), N = 7).The coefficient of determination r 2 for the approximations is shown by color.Blue stars correspond to the maximum of exponent amplitude spectrum (stem lines in Figure 4); grey stars stand for the slowest thermal mode of the cladding subsystem.Red vertical line in (a) and (b) corresponds to square waveguide dr = 0.8 µm (see Table 1).(10) with N = 1 and 9 exponents, respectively, for waveguides with a width of 1000 nm, 500 nm, and 2000 nm (red, green, and gray symbols, respectively); stem line-exponent amplitude spectrum γ (k) /δ (k) θ obtained by the thermal mode decomposition method for microring resonator with default material and geometrical parameters specified in Tables 1 and 2 for different waveguide widths dr: red lines-1000 nm, green lines-500 nm, grey-2000 nm, red/green/grey solid line-theoretical effective thermal relaxation rates from Equation (7) with dr = 1000 nm, 500 nm, and 2000 nm, respectively, red/green/grey dotted line-theoretical effective thermal relaxation rate from Equation ( 16) with dr = 1000 nm, 500 nm, and 2000 nm, respectively, red dashed line-theoretical effective thermal relaxation rate from Equation (17), red/green/grey dash-dotted line-rough estimation via averaging of effective thermal relaxation rates Equation (16) and Equation (17) ) with dr = 1000 nm, 500 nm, and 2000 nm, respectively, red dashed line-Equation (17).
We also studied the exponent amplitude spectrum for different waveguide widths dr. Figure 10 shows that changing the width dr of the waveguide over a wide range does not lead to a significant change in the spectrum.We tried different ways of averaging analytical Formulas (7) (red/green/grey solid line in Figure 10), ( 16) (red/green/grey dotted line in Figure 10), and (17) (red dashed line in Figure 10) to obtain a rough approximation for one-exponent relaxation rate.The best results were obtained with the geometric mean (red/green/grey dash-dotted line in Figure 10).We combined the two proposed analytical formulas, Equations ( 16) and (17), and obtained values δ av θ = δ mod θ • δ slow θ with a similar order to the case of a one-exponent fit (see dash-dotted lines in Figure 10).This proves that the effective thermal relaxation rate from one-exponent fit is essentially a combination of the effective thermal relaxation rates of two thermal subsystems (mode and cladding).However, we note that the difference between the proposed average and the numerical simulation is quite large, so this can only serve as a rough estimate for a particular design.
It is also worth noting that, for a wide waveguide, there are modes (gray stars in Figure 10) with effective thermal relaxation rates greater than that provided by the theoretical Formula (7) (gray solid line in Figure 10).This is due to the fact that, for wide waveguides, estimating the mode transverse profile half-width as the waveguide rosssection half-widths may not be accurate enough.
It is worth paying attention to the fact that Formula (7) does not consider anything other than the resonator by derivation.Furthermore, the resonator itself is considered only in the sense of the mode form, not the resonator form.We found that the distance between the waveguide and the silicon substrate h has a significant effect on the effective parameters.Figure 11 shows the dependence of the effective parameters on the distance h.It can be seen that, in the range from 10 −6 m to 4 × 10 −5 m, the effective thermal relaxation rate changes by more than an order of magnitude.
One can see that the numerical eigenvalues for the slowest thermal modes of the cladding subsystem (grey stars, Figure 11a) are well-described by Equation ( 17) (red dashed line) over a large range of distance between the waveguide and the silicon substrate h.Also, it should be noted that almost all the effective thermal relaxation rate values are in the range between red dotted (Equation ( 16)) and red dashed (Equation ( 17)) lines in Figure 11a.
We also studied the influence of the microchip boundary with air at various heat transfer coefficients in a wide range and did not find a significant effect on the effective parameters.Thus, it is permissible to consider the boundary with air as thermally insulated.17), squares-one-exponent approximation ((10), k = 0), crosses-two-exponent approximation ((10), k = 0, 1), circles-seven-exponent approximation ((10), N = 7).The coefficient of determination r 2 for the approximations is shown by color.Blue stars correspond to the maximum of exponent amplitude spectrum (stem lines in Figure 4); grey stars stand for the slowest thermal mode of the cladding subsystem.Red vertical line in (a) and (b) corresponds to h = 4 µm.
To check the correctness of the boundary conditions and material parameters of the integrated microresonator model, the experimentally measured heater-induced frequency shifts were compared with the results of the numerical simulations [79,80].In the model, the source of thermal power was specified as P = U 2 /R in the heater region, where U is the voltage that is supplied to the contacts of the heater in the experiment, R is the resistance of the heater material, specified by the manufacturer R = 65 ohm.A typical metal (for example, aluminum) was chosen as the heater material because the real parameters of the heater are unknown and are a trade secret of the microchip manufacturer.A voltage U in the range from 1 to 5 volts was applied to the contacts of the heater integrated into the microchip.In this case, the stationary shift in resonances during laser scanning was measured.The measurement points (frequency shifts at the corresponding voltage on the heater) were approximated by a quadratic function: Figure 12 shows that the experimental data (red line) and numerical simulation (blue dots) are in good agreement, which verifies the developed integrated microresonator.

Experimental Limitations
The proposed theory appeals to subtle processes described in detail by thermal mode decomposition.We were interested in how accurately thermal parameters can be measured in experiments with respect to fundamental noise mechanisms and whether it is possible to see the predicted dynamics.The proposed experimental setup that might allow to obtain necessary accuracy is shown in Figure 13.A continuous wave laser heats the microresonator by exciting its mode.The pump laser frequency is bounded to the eigenfrequency of the same microresonator using a self-tuning system based on the Pound-Drever-Hall (PDH) technique [81].We propose to measure the error signal module proportional to the microresonator frequency shift; further, it can be approximated by the sum of the exponents according to the developed approach discussed in the previous section.Note that pump power should be small enough to avoid unwanted nonlinear effects (threshold power depends on the microresonator quality factor) and to minimize Kerr frequency shift.
There are fundamental limitations that add noise to the experimentally measured dependencies, which can limit the number of experimentally observed exponents.In particular, we consider such fundamental factors as thermal noise in the integrated microresonators, frequency and amplitude noise of the pump laser, photodetector noise, and digital noise.We made estimations of the contributions of such factors.In an optical resonator, there are thermodynamic fluctuations of temperature, δT noise , whose variance is where T is the temperature of the heat bath, k B is the Boltzmann constant, ρ is the density, C is the specific heat of the waveguide material, and V is the optical mode volume.Using the material parameters of Si 3 N 4 (see Tables 1 and 2) and the optical mode volume of a typical 24 µm microresonator, we can obtain a value of the standard deviation of the temperature as < δT 2 noise > ≈ 60 µK [18].It can be seen in Figure 2 that the fundamental temperature noise is much smaller (maximum temperature deviation from data (T i (t i ) − T(t i )) ≈ 3 K) than the heating of the microresonator due to the absorption of the pump power for the standard range of the microresonator parameters.Frequency noise for a laser with a linewidth of 10 kHz and a microresonator with a resonance linewidth of 100 MHz provides a relative error only of 5 × 10 −5 .Figure 2 also shows that the relative error reaches T 0 ≈ 14%.The noise of the PDH-stabilizing system, including the thermal noise of the detector and digital noise of the system, can be quite small [81].Thus, considering fundamental limiting mechanisms, it can be argued that such a detailed thermal response distribution can be measured.

Conclusions
In our work, we compared various methods for the determination of the effective thermal parameters for a wide range of material and geometric parameters of an on-chip microring resonator.In particular, we showed that, for the range of standard microresonator parameters, the accuracy of the conventional one-exponent approximation obtained from the single rate equation describing microresonator heating can be relatively low.However, as we showed, the addition of several exponents for the description of this process corresponding to different thermal modes of the microresonator can significantly improve the accuracy of the approximation.We investigated the dependence of the approximation accuracy on the number of exponents (corresponding number of the rate equations) for a wide range of realistic parameters of the considered structures and found the optimal number of the thermal modes that should be taken into account.For the considered microresonator, convergence in the number of exponents was achieved at N = 9 thermal modes.
As a result, we propose an original and effective method for the accurate calculation of the effective thermal parameters in high-Q optical microresonators based on the decomposition of the thermal field into microresonator thermal modes.
The proposed method agrees better with direct numerical simulation for the considered microring resonators and well-known theoretical formulas over the entire range of microresonator parameters, in contrast to the classical approach.The developed method is universal and can be used for a wide range of microresonator platforms.The obtained results also provide deep insight into the thermal dynamics of microresonator-based platforms.

Figure 1 .
Figure 1.The integrated microresonator model used in FEM simulation: (a) top view, (b) side cut.Color map shows heat source distribution Q(r) inside the waveguide.R 0 is the radius of the microresonator.We simulated microresonators with a radius from 24 µm to 840 µm.The transverse dimensions dr and dz of the waveguide were varied in the range from 20 nm to 2000 nm.

Figure 2 .
Figure 2. (a): Blue line is data obtained by the numerical solution of (1) for P = 10 mW heating of 24 µm ring resonator with waveguide width dr = 1 µm and height dz = 0.8 µm; red dashed line shows approximation (fitting) by one exponent; green dashed line depicts approximation (fitting) by two exponents; (b): discrepancy (T i (t i ) − T(t i )) of one-, two-, and three-exponent fit (red, green, and blue dashed lines, respectively) from numerical simulation data.(c): Coefficient of determination r 2 for different number of exponents N from 1 to 15.

Figure 3 .
Figure 3. (a): Examples of the spatial shapes C (k) ( r) normalized to 1 Kelvin (colorbars on the right side of each panel) of the thermal modes of the microring resonator with default material and geometrical parameters specified in Tables 1 and 2, shown in Figure 1: (a) side cut and (b) top view for the thermal mode with the maximum overlap integral C (k) ( r)•Q( r)dV C (k) ( r)•C (k) ( r)dV normalized to δ (k) θ with Equation (6)), and the sum n (k) 2T is equal to the initial one (∑ n (k) 2T = n 2T ).Each coefficient corresponds to the effective thermal relaxation rate δ (k) θ .

Figure 5 .
Figure 5. (a): Blue line-numerical simulation data, red dashed line-approximation by one exponent, green dots-reconstructed dynamics from the exponent amplitude spectrum (stem lines in Figure 4) and (b): discrepancy (T i (t i ) − T(t i )) of one (red dashed line) from numerical simulation data for P = 10 mW heating of 24 µm ring resonator with waveguide width dr = 1 µm and height dz = 0.8 µm.

Figure 6 .
Figure 6.Comparison of the effective (a): thermal relaxation rate and (b): optical absorption rate calculated using different methods for a wide range of the waveguide thermal conductivity k wg : red solid line corresponds to Equation (7) in panel (a) and to Equation (8) in panel (b), red dashed line-Equation (17), squares-one-exponent approximation ((10), k = 0), crosses-two-exponent approximation ((10), k = 0, 1), circles-seven-exponent approximation ((10), N = 7).The coefficient of determination r 2 for the approximations is shown by color.Blue stars correspond to the maximum of exponent amplitude spectrum (stem lines in Figure4); grey stars stand for the slowest thermal mode of the cladding subsystem.Vertical lines correspond to the thermal conductivity coefficients of Si 3 N 4 , SiO 2 , and Si (see Table1).

Figure 7 .
Figure 7.Comparison of the effective (a): thermal relaxation rate and (b): optical absorption rate calculated using different methods for a wide range of the cladding thermal conductivity k clad : red solid line corresponds to Equation (7) in panel (a) and to Equation (8) in panel (b), red dashed line-Equation (17), squares-one-exponent approximation ((10), k = 0), crosses -two-exponent approximation ((10), k = 0, 1), circles-seven-exponent approximation ((10), N = 7).The coefficient of determination r 2 for the approximations is shown by color.Blue stars correspond to the maximum of exponent amplitude spectrum (stem lines in Figure4); grey stars stand for the slowest thermal mode of the cladding subsystem.Vertical lines correspond to the thermal conductivity coefficients of Si 3 N 4 , SiO 2 , and Si (see Table1).

Figure 8 .
Figure 8.Comparison of the effective (a): thermal relaxation rate and (b): optical absorption rate calculated using different methods for a wide range of the microresonator radius R 0 : red solid line corresponds to Equation (7) in panel (a) and to Equation (8) in panel (b), red dashed line-Equation (17), squares-one-exponent approximation ((10), k = 0), crosses-two-exponent approximation ((10), k = 0, 1), circles-seven-exponent approximation ((10), N = 7).The coefficient of determination r 2 for the approximations is shown by color.Blue stars correspond to the maximum of exponent amplitude spectrum (stem lines in Figure4); grey stars stand for the slowest thermal mode of the cladding subsystem.Red vertical line in (a) and (b) corresponds to the default microresonator radius (see Table1) R 0 = 24 µm.

Figure 9 .
Figure 9.Comparison of the effective (a): thermal relaxation rate and (b): optical absorption rate calculated using different methods for a wide range of the waveguide width dr: red solid line corresponds to Equation (7) in panel (a) and to Equation (8) in panel (b), red dashed line-Equation (17), squares-one-exponent approximation ((10), k = 0), crosses-two-exponent approximation ((10), k = 0, 1), circles-seven-exponent approximation ((10), N = 7).The coefficient of determination r 2 for the approximations is shown by color.Blue stars correspond to the maximum of exponent amplitude spectrum (stem lines in Figure4); grey stars stand for the slowest thermal mode of the cladding subsystem.Red vertical line in (a) and (b) corresponds to square waveguide dr = 0.8 µm (see Table1).

Figure 10 .
Figure 10.Squares and stars represent the amplitudes for exponential approximation (10) with N = 1 and 9 exponents, respectively, for waveguides with a width of 1000 nm, 500 nm, and 2000 nm (red,

Figure 11 .
Figure 11.Comparison of the effective (a): thermal relaxation rate and (b): optical absorption rate calculated using different methods for a wide range of the distance between the waveguide and the silicon substrate h: red solid line corresponds to Equation (7) in panel (a) and to Equation (8) in panel (b), red dashed line-Equation (17), squares-one-exponent approximation ((10), k = 0), crosses-two-exponent approximation ((10), k = 0, 1), circles-seven-exponent approximation ((10), N = 7).The coefficient of determination r 2 for the approximations is shown by color.Blue stars correspond to the maximum of exponent amplitude spectrum (stem lines in Figure4); grey stars stand for the slowest thermal mode of the cladding subsystem.Red vertical line in (a) and (b) corresponds to h = 4 µm.

Figure 12 .
Figure 12.(a) Temperature distribution in the integrated chip for a voltage of U = 4.3 V. (b) The blue dots are the results of the numerical simulation of the frequency shift by the heater as a function of voltage, and the red line is the approximation of the experimentally measured frequency shift:

Figure 13 .
Figure 13.Proposed experimental setup for the validation of the method of the effective thermal parameters determination via thermal mode decomposition.P LS and ω LS -laser source optical power and generation frequency, ω 0 -microresonator resonant frequency, δT = ∑ δT (k) from Equation (10).

Table 1 .
Default chip parameters used for FEM simulations.

Table 2 .
Physical properties of waveguide material used for FEM simulations.