Analytical modeling of the static and dynamic response of thermally actuated optical waveguide circuits

Thermo-optic phase shifters allow one to dynamically tune and control the operation of integrated-optics interferometers. They have been demonstrated nowadays in different waveguide platforms, and their reliable functioning has enabled the realization of reconﬁgurable circuits of notable complexity. The design approach to such devices is often based on ﬁnite-element numerical simulations, which provide accurate descriptions of the underlying thermal phenomena, at the price of long computational times. Here, on the contrary, we devise an analytical model for the heat diffusion in a simpliﬁed geometrical conﬁguration. The model describes both static and dynamic regimes, and can be conveniently applied both to three-dimensional waveguide devices inscribed by femtosecond laser pulses and to planar lithographic circuits. The accuracy of the predictions of the model is validated with experimental measurements on Mach-Zehnder interferometers with different geometries, realized in both kinds of platforms.

Thermo-optic phase shifters are widely adopted to achieve an active control on integrated-optics interferometers. Such devices exploit the refractive index variations induced by local temperature changes, which in turn are typically produced by resistive heaters deposited on the surface of the chip.
These devices have been demonstrated in multiple platforms, including silica-on-silicon [11,14], silicon nitride [15,16], silicon-on-insulator [3,17], and femtosecond laser written circuits [18,19]. Such active components are often designed on the basis of an empirical approach, or by exploiting numerical finite-element simulations of the heat diffusion [20,21]. As a matter of fact, analytical models, developed in conveniently approximated settings, can be beneficial for a deeper understanding of the involved physical dynamics, and for laying out simple guidelines for the engineering process.
Analytical models of the heat diffusion have been indeed reported in a few papers studying thermo-optic phase shifters. However, to the best of our knowledge, these are limited to the analysis of the steady state [22,23] or to giving general and approximate indications on the dynamic behavior [24]. In addition, these models may not be suitable to describe femtosecond laser written devices, where crosstalk among waveguides of the same interferometer can be relevant due to the large substrate thickness.
In this paper, we delineate an analytical model for the heat diffusion that works in the geometrical setting of typical optical chips and is valid in both static and dynamic regimes. This model can be applied to three-dimensional waveguide devices inscribed by femtosecond laser pulses and, as well, to planar circuits fabricated by lithographic technologies. Experimental measurements, performed on both kinds of platforms, validate the accuracy of the model predictions.

II. PREMISE
The phase delay accumulated by coherent light, propagating in an arbitrary segment (possibly curved) of length L of a single-mode waveguide, is given by where λ is the wavelength and n( r) the effective refractive index of the optical mode in a specific point of the waveguide. If temperature variations are not too big, the effective refractive index can be considered a linear function of the local temperature, where n T is the thermo-optic coefficient. In this paper, we study essentially the temperature distributions, both in static and dynamic conditions, and we implicitly assume that the refractive index, and thus the induced phase shift, follow in- The method of images allows us to define a symmetrized system, which yields (in the region 0 z h) the same temperature distribution as the original one. Note that, for convenience of representation, in (b) the axes xz are rotated by 90 • with respect to (a). stantaneously and linearly the temperature trend according to the above equations.
In particular, we refer to the geometry depicted in Fig. 1. A planar dielectric layer of thickness h is placed in contact with a thermostat (or heat sink) at a fixed temperature T 0 . Above the dielectric, we consider the presence of air, assumed as an insulating medium for heat conduction. We neglect convection and radiative phenomena. The whole structure is infinite and uniform along the y direction (orthogonal to the picture). A wirelike heater is placed at the origin of the Cartesian system and runs parallel to the y axis.
Dealing with interferometric devices, it is useful to consider the difference in the phase acquired by the light propagating in two different waveguides. To this purpose, we consider two waveguides contained in the dielectric layer, which also run parallel to the y axis, and are characterized by their transversal coordinates (x 1 , z 1 ) and (x 2 , z 2 ) [ Fig. 1(a)]. If we consider waveguide segments of length L, exploiting the previous relations, we can write the relative phase as a function of the temperature difference between the two points, where T = T 1 − T 2 = T (x 1 , z 1 ) − T (x 2 , z 2 ). In particular, in the case of an ideal Mach-Zehnder interferometer (MZI), if coherent light is injected into one input port, the optical transmission on the opposite output port is written as where φ 0 is a phase term that is built in the interferometer.
We note that the geometry illustrated in Fig. 1 can conveniently model femtosecond laser written circuits: In this case waveguides are typically inscribed at arbitrary depths in the bulk of a thick (∼1-mm) glass substrate, which may be placed on a metallic support that works as a heat sink [18,25,26].
In addition, the same geometry can be employed to model planar lithographic waveguides, if the waveguides consist in doped regions or high-index insertions within a silica (dielectric) layer that is grown on top of the silicon substrate [14,15,27]. The silicon substrate has thermal properties that are similar to those of metals and can work as the heat sink [28].

III. STATIC TEMPERATURE DISTRIBUTION
For a general, arbitrary configuration of heat sources, the heat conduction within a uniform material region is described by where T ( r, t ) gives the temperature distribution, q( r, t ) is the thermal power generated per unit volume, κ is the thermal conductivity of the material, c is the specific heat capacity, and ρ is the density. To simplify the notation, in the following we will omit to write the dependencies on r and t. In the static case, the time derivative in Eq. (5) is zero: As a first step to develop our analytical model, we consider a single wirelike heater, infinitely long, immersed in a uniform medium. This geometry is characterized by cylindrical symmetry and all physical quantities only depend on the distance from the heater. For convenience, we can adopt a cylindrical coordinate system (r, θ, y) coaxial to the wirelike heater, where all quantities depend only on r.
In detail, in the uniform space without heat sources, Eq. (5) becomes and yields the general solution where C, r * , and T * are proper constants; in particular, T * is identified as the temperature at the coordinate r * . If the wirelike heater, placed on the y axis, dissipates statically a thermal power P per unit length, such thermal power will also have to flow through each cylindrical surface of unit length and of arbitrary radius. Choosing a surface with radius r 0 , the Fourier's law for heat conduction is written as Therefore, Eq. (8) becomes finally Note that, in the case of a source consisting in a resistive heater, it is indeed more appropriate to fix as a boundary condition the dissipated power, rather than the heater temperature. In fact, this is the quantity that is fixed by the Joule's law, when a certain voltage or current is imposed in the resistor.
As a second step to develop our model, we introduce the presence of the glass-air interface, with the wirelike heater lying exactly on that. We observe that the thermal conductivity of substrates such as silicon or fused silica (λ Si ∼ 130 W m −1 K −1 , λ SiO 2 ∼ 1 W m −1 K −1 ) is quite larger than that of air (λ air ∼ 0.02 W m −1 K −1 ). It can be reasonable to approximate air as a perfect heat insulator, with respect to the substrate. Thus, at the substrate-air interface the temperature gradient ∇T (which is indicative of the heat flux, because of the Fourier's law) can be assumed as parallel to the interface itself.
Such an interface corresponds geometrically to a symmetry plane of the cylindrical temperature distribution, along which the temperature gradient is parallel to the plane itself. Therefore, the upper half space can be replaced by a perfect heat insulator without affecting the solution of the equation in the lower half space. A solution of the kind of Eq. (10) is thus valid in the lower half space both in the uniform case and in the presence of the air-glass interface, under the reasonable approximation that the heat flux across such an interface is negligible. However, we need to drop a factor of 2 at the denominator, to take into account the fact that the power P is dissipated only on the lower semicylinder: As a third and final step, we add the thermostat. The presence of the thermostat implies the existence of a plane horizontal surface, at depth h, having a uniform temperature T 0 . To retrieve the temperature distribution in this configuration, we employ the so-called method of images. Namely, we look for an equivalent system, consisting in a uniform and infinite medium and a distribution of heat sources, which provides the same conditions at the boundaries of the original system.
We can consider in particular the configuration shown in Fig. 1(b), where the original system is symmetrized with respect to its top surface. Such a configuration consists of a slab of material with thickness 2h, delimited both from above and from below by surfaces with fixed temperature T 0 , with one wirelike heat source placed in the middle. Furthermore, we add an infinite series of sources at coordinates z = 2mh (with m integer, positive or negative), alternately dissipating a power −P or +P. For the linearity of the heat equation, the solution will be written as the sum of the contributions of these infinite sources.
Because of symmetry constraints, the gradient ∇T , along the whole plane z = 0, lies parallel to the plane itself. This provides the condition of vanishing heat flow across the top surface of the original system. Again by symmetry considerations, one observes that the series of sources is antisymmetric with respect to both planes at z = ±h. In other words, for each positive source placed at a certain distance from a given point of the plane z = h (or z = −h) there is another negative source placed at the same distance but in the opposite half space. This means that the summed contributions of the sources vanish across these whole planes: The temperature in these planes is uniform and can be fixed to an arbitrary value T 0 .
Writing the solution mathematically, one obtains with r (m) indicating the distance between the point of interest at coordinates (x, z) and the mth source: By collecting the contributions in couples of one positive and one negative source, in detail the (m + 1)th with the (−m)th, and by exploiting the mathematical properties of logarithms, Eq. (12) can be rewritten as It is easy to show that such an infinite sum is converging, using the Leibniz criterion. It is also easy to observe that Eq. (14) provides T = T 0 for every value of x if z = h, because each term of the series is vanishing. The phase shift produced in an interferometer composed by two waveguides, placed at coordinates (x 1 , z 1 ) and (x 2 , z 2 ), is proportional to the temperature difference T = T 1 − T 2 between those points, as per Eq. (3). From Eq. (14) we see that T , and then φ, is proportional to P. Modeling a realistic device, where both the waveguides and the heater are taken of length L, one can consider the full dissipated power on the thermal phase shifter as P ts = P · L and write φ = αP ts .
The proportionality coefficient α, which indicates somehow the efficiency of the device, is calculated from Eq. (14) to be Approximate expressions for α can be derived in several cases of significant experimental relevance.
In the case where the thickness of the substrate is much larger than the distances between the heater and the waveguides, i.e., for h x 1,2 and h z 1,2 , the terms of the sum in Eq. (16) become negligible already for m 1, and even in the term m = 0 we can approximate 1. Hence we produce a very simple expression, where R = r 2 r 1 = which is the ratio of the distances of the heater from the second waveguide and from the first waveguide. We may note that Eq. (17) coincides also with the result of Ref. [18].
In fact, such an approximation is appropriate in the case of femtosecond laser written waveguides inscribed at a depth of tens of microns inside a uniform glass substrate of millimetric thickness. Figure 2 plots the experimental measurements of α, conducted on MZIs realized by femtosecond laser waveguide writing, together with the theoretical predictions of Eq. (17). For each interferometer the waveguide depth is the same for both branches (z 1 = z 2 = z) and the first waveguide lies precisely below the heater (x 1 = 0). Different geometrical parameters are explored for z and x 2 , with z ranging from 15 to 53 μm and three possible values for x 2 = 50, 254, and 1000 μm. Experimental data points are acquired by coupling coherent light at a 1550 nm wavelength in one input port of each MZI, and by measuring the optical power simultaneously at both output ports, using distinct power-meter head sensors. A MATLAB script controls the power supply driving the resistive microheaters, and records the optical measurements for each value of the actuation current. The coefficient α is retrieved by a nonlinear fit, taking into account Eqs. (4) and (15). Further experimental details on the device fabrication and characterization are reported in Appendix A. Theoretical curves are calculated considering n T = 7 × 10 −6 K −1 and κ = 1.09 W m −1 K −1 [19,29]. Very good agreement between theory and experiment is observed for all devices, thus showing the potential of this simple model. (A dynamical characterization will be performed on a subset of these devices, i.e., devices A, B, and C in Table I, and will be discussed in Secs. IV and V.) On the other hand, in the case of lithographic devices it is common to have a very thin dielectric thickness h, a heater placed just above the first waveguide (x 1 = 0) and the second waveguide at a distance x 2 h. Waveguides are also usually placed at the same depth z = z 1 = z 2 in the substrate. In this case, for all values of m, which implies that the temperature of the second waveguide is approximately equal to T 0 and thus T T 1 . In this case As a further approximation, which allows a rapid estimation of the order of magnitude of the α coefficient, one could even truncate the series at its first term, i.e., m = 0. This is equivalent to employing Eq. (17) with  Table I. The experimental value α expt is TABLE I. List of the integrated-optics MZIs used to validate experimentally the dynamic model. Relevant geometrical parameters are listed for each device. In all cases, both waveguides are at the same depth z = z 1 = z 2 and the first waveguide is placed exactly below the heater (x 1 = 0). Devices A, B, and C are treated in the approximation of thick substrate and R is calculated according to Eq. (18). On the contrary, for device D the dielectric layer is thin and R follows Eq. (20).   Table I). The black line marks the value that is experimentally measured; measurement uncertainty is indicated by the line thickness. Dots represent successive approximations achieved truncating the series in Eq. measured with the same procedure used for the femtosecond laser written MZIs. In this lithographic device, waveguides are buried at z = 8.8 μm below the surface in a silica layer of thickness h = 19.8 μm. The silica layer is deposited on top of a silicon substrate, which can work as a heat sink. The first waveguide of the interferometer lies exactly below the heater (x 1 = 0) while the other one is placed at x 2 = 250 μm. Assuming as material properties n T = 11 × 10 −6 K −1 and κ = 1.4 W m −1 K −1 (standard values for pure silica), the estimate of Eq. (19) shows good matching with the experimental value, when the series is evaluated with a sufficient number of terms. However, even Eq. (17) catches the order of magnitude of this coefficient. Thus, it can provide useful indications, from a circuit design perspective, on how α scales with the geometrical parameters of the circuit and the thermal properties of the material.

IV. DYNAMIC RESPONSE: FREQUENCY ANALYSIS
To study analytically the dynamic response, i.e., the response of the system for a time-varying q, we can perform a frequency analysis and take the Fourier transform of Eq. (5): We can choose to solve the equation in cylindrical coordinates, assuming conditions of perfect cylindrical symmetry and a wirelike heat source. In particular, in the transformed space, we refer to a heater that dissipates an oscillatory power P =P (ω) per unit length, at a given frequency ω. The solution procedure, which is quite lengthy, is reported in Appendix B and results in the following expression forT (ω), where K 0 indicates the modified Bessel function of the second kind, of order 0.
As we did in the case of the static response, we can adapt that solution to the case of a semi-infinite medium interfaced with a perfectly insulating environment. Namely, we multiply it by a factor of 2, which takes into account the fact that P is fully dissipated in the lower half space. Then, to obtain the complete solution for the geometry considered in Fig. 1, we add up the contributions of an infinite series of sources with an alternate sign, where r (m) is defined as in Eq. (13). The Bessel function K 0 (ζ ) (where ζ is a complex number) yields the following two approximations: Using Eq. (24) it is easy to trace back Eq. (23) to its static limit Eq. (14) for ω → 0. The other limit, Eq. (25), shows instead that when ω is sufficiently high, the Bessel functions in Eq. (23) decrease faster than exponentially with |ζ |. On the one hand, this means thatT (ω) → 0 for ω → +∞. On the other hand, this also means that for increasing frequencies only the terms with the lower |r (m) | become significant. In particular, for the highest frequencies the m = 0 term (which corresponds to the lowest |r (m) |, i.e., r (0) = r = √ x 2 + z 2 ) tends to dominate over the other ones.
Considerations regarding the limits given by Eqs. (24) and (25) can be applied also to the expression ofα(ω). In particular, if r 1 < r 2 we can observe that r 1 = r (0) 1 is the lowest among all coefficients r (m) 1,2 . Application of the limit given by Eq. (25) implies that an interferometric device based on a thermo-optic phase shifter responds as a low-pass system, whose behavior at high frequencies is ultimately determined by the distance between the heater and the first waveguide.
Again in analogy with the static case, we can derive simplified expressions ofα(ω), which are noteworthy in practical occurrences. To this purpose, it is useful to get further insight into Eq. (26). The infinite sum in the latter equation is com-  31), which delimits the bandwidth with good approximation. (b) shows the same functions of (a), normalized to their limit value at low frequencies (which is 1/ ln R); the bandwidth reduction for increasing R is made more evident. posed of terms of the kind In detail,α with the parameters (m) and R (m) defined as The modulus off R , for different values of R and as a function of the normalized frequency ω/ , is plotted in Fig. 4. We observe that this function presents three different regions: a first flat region for low frequencies, a further flat region (with vanishing modulus) for high frequencies, and an intermediate region with |f R | monotonically decreasing with ω, in between the two. The first flat region, where the modulus is maximum, corresponds to the region where the limit given by Eq. (24) is valid for both K 0 functions, giving We could say that this region is delimited by the condition that the arguments of the Bessel functions are both smaller than 1. Specifically, if R > 1, this occurs for R ω < 1, Figure 4 shows that ω B can be indeed a good approximation for the bandwidth. The high-frequency region, conversely, corresponds to the other limit given by Eq. (25).
In the case of a very thick substrate (h x 1,2 and h z 1,2 ), which is relevant for femtosecond laser written devices, one observes that R (m) 1 for m = 0 and thus the functions f R (m) for m = 0 have a negligible modulus. It is then licit to approximate the sum in Eq. (26) to its zeroth term, Note that R here has the same definition as in static problem, Eq. (18). Using the latter definitions for and R the plots in Fig. 4 describeα approx (ω) apart from a proportionality constant. The operation bandwidth of the thermo-optic phase shifter may be identified with the width of the first flat region off R , i.e., the frequency ω B defined above. In this case, in particular, one reads which depends only on the distance r 2 of the farther waveguide from the heater. On the other hand, the distance r 1 of the first waveguide, which in turn governs the behavior at very high frequencies, is here of negligible influence. Figures 5(a), 5(c) and 5(e) report the experimentally measured frequency response |α(ω)| for three MZIs fabricated by femtosecond laser writing in a thick glass substrate. Geometrical details of the devices are given in Table I. To perform the measurements, small sinusoidal modulations at different frequencies, provided by an electronic function generator, were superimposed to a base current level dissipated on the heater, so that the interferometer phase is modulated around φ 0 = π/2. The modulation amplitude of the optical signal at the output of the device was recorded as a function of the frequency, using an IR-sensitive photodiode connected to a digital oscilloscope. In this small-signal regime, such an optical-modulation amplitude can be considered proportional to the phase modulation and thus proportional toα(ω). Further information on the experimental characterization procedure is provided in Appendix A.  Table I. Black crosses indicate experimental measurements; the red lines represent the frequency responses calculated according to the approximate expression α approx (ω), as in Eq. (32). In the case of device D, the response calculated using the more accurate formula given in Eq. (36) is also shown in green. All curves and experimental data are normalized to the static limit. Error bars on the experimental points are smaller or comparable to the marker size.  Table I. The black lines are experimentally acquired oscilloscope tracks; the red lines represent the step response calculated according to the approximate expression Eq. (43). In the case of device D, the response calculated on the basis of the antitransform of Eq. (36) is also shown as the green curve. All step responses, both experimental and theoretical ones, are normalized to their asymptotic value.

023094-7
As shown in the figures, the experimental frequency responses yield good overlap with the curves calculated with the approximate expression Eq. (32), where and R are given by (33). In particular, to calculate the theoretical curves, the following values of the parameters were employed: κ = 1.09 W m −1 K −1 , c = 768 J kg −1 K −1 , and ρ = 2380 kg m −3 [29].
Another case of interest, which has been already considered in Sec. III and is relevant for lithographic waveguide devices, is the case of a substrate that is much thinner than the interwaveguide distance, with the heater positioned just above the first waveguide (in detail, x 1 = 0, z 1 = z 2 = z, and t x 2 ). As for the static case, also here we neglect the contribution ofT 2 (ω), which means In order to perform further approximations, it may be useful to rewrite the sum in (35) collecting the (m + 1)th and the (−m)th terms together, thus involving again thef R functions, provided that the proper parameters are defined: We note that, since in any case z < h, the parameter R (m) is significantly different from 1 only for the lowest values of m and is strictly decreasing with m. Reminding that sup{| f R (m) |} = ln R (m) , it is reasonable to produce a rough estimation ofα(ω) by truncating the series in Eq. (36) just to the m = 0 term. Namely, also in the case of a thin substrate we can employ the approximate expression forα approx given by Eq. (32), provided that we define in analogy with the parameter R employed for the static case and defined by Eq. (20). Figure 5(g) reports the experimentally measured frequency response for one MZI realized in a SiON platform, which is device D in Table I. Measurements were performed with the same method used for the femtosecond laser written devices. In this case, the theoretical curve calculated according to Eq. (36) explains with remarkable accuracy the experimental data. The approximated model of Eq. (32), here employed with the parameters and R given by (38), is less accurate. However, it still follows the experimental points with a relative error better than 22%. The latter appears to be a surprisingly good figure, looking at the extreme simplicity and degree of approximation of our model. Here, to calculate the theoretical curves, employed parameters were κ = 1.4 W m −1 K −1 , c = 772 J kg −1 K −1 , and ρ = 2200 kg m −3 (standard values for pure silica).

V. DYNAMIC RESPONSE: STEP RESPONSE
After analyzing the frequency response we can pass to the time domain, studying the evolution of the phase delay φ(t ) under the action of a time-varying thermal power P ts (t ) dissipated on the wirelike heater. Given the relation φ (ω) = P ts (ω) ·α(ω), valid in the Fourier space, the function φ(t ) is calculated by means of the convolution, where α(t ) the Fourier antitransform ofα(ω): In particular, the response to a power step P ts (t )

is the Heaviside function] is calculated as
We will devote a specific analysis to the case in which the frequency response is described byα approx (ω), as given by Eq. (32). Indeed, we have discussed in Sec. IV how this simplified model can be very accurate in the case of femtosecond laser written devices and how it can describe with reasonable approximation also thermo-optic phase shifters fabricated by lithographic techniques. In that case, by defining the adimensional step response, one gets Figure 6(a) shows that, for increasing R and having fixed , the adimensional step response increases its asymptotic value. The latter actually coincides with the amplitude of the static responsef R (0) = ln R. The different curves, however, are practically overlapped in the first instants. We can interpret such a behavior in the light of the analysis performed on the frequency response: The early part of the step response is mainly determined by the high-frequency region off R ( ω ), which indeed has negligible dependence on R. On the other hand, the late part of the step response is strongly influenced by the low-frequency components, which are enhanced by increasing R. If we normalize the step response to the asymptotic value, as in Fig. 6(b), we see that the higher is R, the longer is the time required to reach any given fraction of such asymptotic value. In fact, we may associate to ω B a correspondent time constant: A numerical evaluation shows that F R (t · ) reaches in t = τ B about 90%-95% of its asymptotic value, with small variations depending on R.
In the case of the thick substrate, where R = r 2 /r 1 , the time constant reads which depends only on r 2 and the characteristics of the material. Note that, while r 2 may be seen as mainly responsible for a sort of settling time, namely of the late part of the step response, the distance r 1 of the first waveguide affects directly the timescale of the graphs in Fig. 6, because of the definition of . In particular, decreasing r 1 while leaving unvaried R = r 2 /r 1 causes a quadratic enhancement of the rapidity of the response. Figures 5(b), 5(d), 5(f), and 5(h) plot the experimental step response for the devices listed in Table I, whose frequency response was already discussed in Sec. IV. The step response was measured with the same experimental apparatus used for the frequency-response characterization. Here, the microheater is driven with a square wave (having a period of several seconds) instead of a sinusoid. As for the frequency analysis, measurements were carried out in a small-signal regime, i.e., in a condition in which the measured variation of the optical signal at the output can be considered proportional to φ(t ).
The red curves are theoretical predictions calculated using Eq. (43). Consistently with the observations made about the frequency response, the predictions of our simplest model overlap very well the experimental data for devices A, B, and C, i.e., the ones fabricated in a thick substrate by femtosecond laser waveguide writing. Such a model still captures the timescale of the response of the SiON device (device D), even if it is not accurate in predicting its details. However, a more refined calculation of φ(t ), performed using Eq. (41) with α(ω) given by Eq. (36), looks to be in excellent agreement with the observed response (green curve).
We note that our study has focused on the rise transients. Fall transients, in response to P (t ) = P 0 [1 − H (t )] should indeed follow the same dynamics, due to the linearity of heat equation. In our experimental measurements we have always measured both rise and fall transients in response to a squarewave signal (even if the fall transients are not shown in this work) and we never noticed significant differences between the two.
Differences in the rise and fall transients, which have been reported in several works on thermo-optic phase shifters [30,31], should be attributed to the presence of nonlinearities in the heat equation (e.g., because of the occurrence of heat irradiation), to insufficient accuracy in setting the operation point of the interferometers, or to other nonidealities that have not been taken into account in our thermal model.

VI. OPTICAL RESPONSE OF A MACH-ZEHNDER SWITCH
Thermal phase shifters may be employed as the dynamical element of integrated optical switches based on MZIs. Namely, the thermo-optic device acts on the optical phase between the two branches of the interferometer, to switch the normalized transmission on either output from 0 to 1 (or viceversa), with the optical response following Eq. (4) in the case of a perfectly balanced device.
It is worth analyzing briefly the dynamics of the largesignal optical response that originates from the thermal dynamics studied in the previous sections. In particular, we shall consider a further ideal situation in which φ 0 = 0, and a step of thermal power is applied to the heater in order to produce optimum contrast, i.e., The optical response of the MZI, given by Eq. (4), acts as a nonlinear transformation on the step responses studied in the previous section. Interestingly, in the condition given by Eq. (46) the nonlinearity may have the effect to delay the very first instants of the response, and to compress the late part of it. As a consequence, the effective optical settling time may be shorter than the thermal one, as seen in Fig. 7. . The continuous lines are experimentally acquired oscilloscope tracks. In particular, the red continuous line is the small-signal response (measured around φ 0 = π/2) normalized to 1, which is proportional, with good approximation, to the phase variation φ(t ). The blue continuous line describes the large-signal optical response of the devices, normalized to 1, measured by reaching experimentally the condition given by Eq. (46). The dashed lines are best fits of the two curves with exponential functions. Note that the timescale of the two panels is very different, with device D showing a response that is hundreds of times faster than device C.
We note that, in several works in the literature [17,32], the response of thermo-optic devices has been analyzed in terms of exponential time constants. However, as understood from our model and confirmed by the experimental data, these devices are far from being described by a single-pole dynamics, featuring a simple and classical exponential response of the kind F τ (t ) = 1 − exp(−t/τ ).
Deviations from the exponential curve can be diverse. For instance, in the case of a femtosecond laser written device realized in a thick substrate, with a large distance between the two waveguides, the settlement time of φ(t ) may be mostly determined by the need for the heat wave to reach the second waveguide. After a relatively steep rise in the first instants [see the red curve in Fig. 7(a)], the response curve slows down, yielding a poor adaptation to a curve of the kind F τ . On the other hand, the optical response, which is a nonlinear function of φ(t ), may yield a better similarity with an exponential best fit [see the blue curve in Fig. 7(a)].
Somehow different observations can be made regarding the response of the lithographic SiON device (device D in Table I). As discussed previously, in this device the dynamic response is mainly determined by the heat diffusion time in the region around the heater, while the farther waveguide has a negligible influence. This produces a fast rise and quick settlement of φ(t ) to its asymptotic value, following a curve that in this case may overall resemble an exponential [see the red curve in Fig. 7(b)].
However, in such a short timescale [ Fig. 7(b)] it becomes noticeable that the first derivative of φ(t ) in t = 0 is zero, which is indeed a nonexponential feature. In fact, a (short) time is needed also for the heat to diffuse across the few microns which separate the heater from the first waveguide. This feature is actually present also in the response of device C, but it may be less apparent in the long timescale of the response of that device.
We further note that, in the case of device D, the nonlinear transformation of Eq. (4) results in a worse fitting of the exponential curve in the case of the large-signal optical response (blue curve in Fig. 7), with respect to the case of the small-signal one.

VII. DISCUSSION AND CONCLUSIONS
We have developed an analytical model to predict the response of thermo-optic phase shifters realized in integrated optics, both in the static and dynamic regimes. Experimental validation of the developed theory is performed by measuring the optical response of MZIs, fabricated by femtosecond laser waveguide writing in glass, or by lithographic techniques in a SiON platform.
The main results of such model are summarized in Table II. Notably, we have worked out simple formulas that allow an approximate description of the dynamic response in relevant geometrical limits. In particular, we have shown that two parameters and R, defined on the thermal properties of the material and on geometrical features of the interferometer, are mostly determining the dynamic behavior. Therefore, the evaluation of these sole two figures permits us to predict in a simple and synthetic way how the device response changes along with the design variables.
Differently from an approach based on finite-element simulations, our analytical approach indeed enables the quick evaluation and comparison of many different configurations, with minimal computational effort. In addition, it allows us to point out and study with ease the peculiar features of the optical response, which are sometimes neglected in the literature, such as the nonexponential characteristics of the transients.
We note that our model explains well the experimental data even in the presence of rough approximations, such as considering the heater as infinitely long and wirelike, whereas in actual devices this may be of non-negligible lateral size. In fact, a large heater may be seen as the composition of many juxtaposed wires, which sum up their effects. From simple geometrical considerations it is easy to see that, on the one hand, the distance between these many wires and the farther waveguide does not change significantly. On the other hand, the closer waveguide sees an averaged effect of the many wires where, however, the closest ones play a major role. This may provide a justification for substituting them with a single one, placed in the closest position to the first waveguide. A deeper discussion on these aspects is reported in Appendix C.
The analytical model presented here cannot replace finiteelement simulations that take into account the precise geometrical and material features of the specific devices. This is especially true in cases of complex three-dimensional and microstructured devices [33,34]. However, we believe that this work provides useful insights in the physical operation of integrated thermo-optic devices, which can guide engineers in the design process. We further envisage that an analytical approach analogous to the one adopted here could be developed, in future work, to describe crosstalk dynamics in large optical meshes. Such a toolbox of approximated formulas that describe the thermal dynamics could be beneficial for a quick benchmark of complex circuit designs, before embarking in a time-consuming numerical optimization process.

Fabrication of femtosecond laser written circuits
Femtosecond laser micromachined devices are fabricated in a commercial aluminoborosilicate glass substrate (Corning EAGLE XG) by means of a two-step process. First, the optical circuit is inscribed in the substrate via femtosecond laser waveguide writing; then, the microheaters and electrodes are deposited and patterned on the substrate surface by maskless photolithographic lift-off.
In detail, waveguides are written using a Yb:KYW cavity-dumped femtosecond laser producing 300 fs pulses at 1030 nm wavelength and 1 MHz repetition rate. Pulses with 580 nJ energy are focused in the substrate by a 20× water immersion objective [numerical aperture (NA)=0.5], while the substrate is translated at the constant speed of 20 mm s −1 with respect to the laser focus.
The resulting waveguides allow for a single-mode operation at 1550 nm wavelength. Different MZIs are realized, composed of two balanced directional couplers and 3-mmlong straight waveguide segments in the central parts. Such devices are inscribed at depths ranging between 15 and 53 μm below the substrate surface, and the central waveguide segments are distanced 50, 254, or 1000 μm from each other. In the balanced directional couplers, waveguides are brought as close as 7 μm for an interaction length of 1 mm, to allow evanescent-field coupling. Curved waveguide sections follow circular arcs with a 45-mm curvature radius.
In order to fabricate the metallic heaters, a 1.5-μm-thick layer of photoresist (Microchemicals AZ 5214E) is deposited on the substrate as a sacrificial layer, using a Karl Süss RC8 spin coater. The photoresist is then exposed to UV light with a Heidelberg MLA100 maskless aligner, which allows for imprinting custom layouts on the surface of our sample without the need for a physical mask. After development of the photoresist the pattern is defined, and a 300-nm-thick chromium layer is deposited on the substrate surface using an e-beam evaporator (Evatec BAK 640). Finally, the photoresist is stripped with acetone and the desired metal pattern remains on the surface. Chromium was chosen because of its high thermal stability. Resistors are fabricated with a width of 15 μm and a length of 3 mm, resulting in a resistance of about 1 k (resistivity of chromium is 124.5 n m).
The terminal electrodes of the resistors are fabricated similarly, but evaporating copper instead of chromium, as it has a lower resistivity (17.2 n m) and it is compatible with wire bonding. In order to minimize the series resistance introduced by these electrodes, they are fabricated with a width of 100 μm and thickness of 1.5 μm. For this second step of lithography a thicker 4.5-μm layer of photoresist (Microchemicals AZ ECI 3027) is used.
After completing the two lift-off steps the devices are thermally annealed in vacuum at 420 • C for 60 min in the Evatec BAK 640 chamber. Annealing is required as it increases the breakdown power of the microheaters and significantly improves their stability in time.

Characterization of the device responses
Experimental characterization of the thermally actuated MZIs, fabricated in both platforms, is performed using coherent light at 1550 nm wavelength from a laser diode (Thorlabs L1550P5DFB), which is coupled into one input of the device by means of a single-mode fiber.
In detail, a Thorlabs SMF-28-J9 fiber is used to couple femtosecond laser written waveguides. In the case of the SiON device, in order to guarantee minimal losses, the SMF-28-J9 fiber is spliced to a Nufern UHNA4.
To characterize the static device response, a Keithley 2231A dc power supply is used to drive the resistive microheater, while output light from both ports of the MZI is collected with a 0.55-NA aspheric lens (Thorlabs C230TMD-C). The two optical signals are simultaneously monitored by two power meters (Ophir Nova II, equipped with PD300-IR photodiode heads). The electrical power supply and the optical power meters are computer controlled, and a MATLAB script manages the automated acquisition of the optical measurements for different values of the actuation current.
For the dynamic characterization (frequency response and step response) the driving current for the microheater is provided by a function generator (Tektronix AFG3011C), capable of supplying both sinusoidal and square-wave signals. In addition, a single output port of the MZI is monitored, using an IR photodiode (Thorlabs PDA20CS2). Waveforms of the input voltage and output optical signal are recorded with a digital oscilloscope (Tektronix DPO2024B).

APPENDIX B: SOLUTION OF THE DYNAMIC DIFFERENTIAL EQUATION
We consider a system with perfect cylindrical symmetry and we take the usual cylindrical coordinates system (r, θ, y), so that necessarilyT =T (r). In addition, if the heat source is wirelike and located in r = 0, outside of this pointq = 0. The heat equation Eq. (21) becomes homogeneous, where we have collected all constants as With the change of variables s = Br we can further rewrite Eq. (B2) as This is actually a modified Bessel equation of the form with α = 0. Such an equation yields the general solutioñ where C 1 and C 2 are proper constants to match the boundary conditions, and I 0 and K 0 are the modified Bessel functions of the first and second kind, respectively. Passing back to the coordinate r, we have We note that I 0 ( 1+ j √ 2 ωcρ κ r) diverges as r goes to infinite, while K 0 ( 1+ j √ 2 ωcρ κ r) tends to 0 in the same limit. Within a uniform and infinite medium, we are not putting any boundary at r → +∞. Therefore, the only function of the two that may retain some physical meaning is K 0 , and we can write the solution of the heat equation Eq. (B7) as whereT 0 is an arbitrary complex multiplicative constant. In order to find the value ofT 0 , we consider a cylindrical heat source of external radius r 0 that dissipates a certain power P (t ) per unit length through its external surface. The thermal flux across an elementary surface dS (oriented orthogonally to the vector u n ) is given by the Fourier law, −κ ∇T · u n dS = dP. (B9) By integrating over a cylindrical surface with radius r 0 and unit height, and by taking the Fourier transform, this gives Practically, by fixing the dissipated power through this cylindrical surface with radius r 0 we are fixing the derivative ofT with respect to r at r = r 0 . The derivative of Eq. (B8) is By equating Eqs. (B10) and (B11), the latter evaluated in r = r 0 , we obtaiñ If we consider a source such as the wirelike heater, which is infinitely thin, we can take r 0 → 0. In such a limit, K 1 (z) ∼ 1 z , thusT 0 simply becomesT and we finally obtain Eq. (22).

APPENDIX C: ROLE OF LARGER HEATERS
Our analytical model considers a heater with vanishing transverse dimension, i.e., a wirelike element. Experimentally, this is not always an obvious assumption, because resistors, which are patterned on thin metallic films covering the substrate, may extend laterally for several microns. In principle, our model could be expanded to the case of a non-wire-like heater by considering the latter as composed of infinitely many wires (of infinitesimal width), whose contributions are properly integrated.
The wirelike heaters, ideally composing the real one of width w, are placed at coordinates (x, y) = (ξ, 0) with −w/2 ξ w/2. Each of these wires alone is of width dξ and dissipates a power per unit length dP = P w w dξ (where P w is the power dissipated on the whole heater per unit length). Each wire then determines a contribution to the temperature distribution written adapting Eq. (23), where We note that we limit ourselves to discuss here the dynamic case, because the static condition coincides with its limit at ω → 0 [possibly up to an additive constant T 0 , for that regards the static distribution T (x, z)].
The contributions dT ξ (ω) are finally integrated over the width w to produce the resulting temperature distribution. The formula for determiningα(ω) can be adapted accordingly and results in the average of the differentα ξ (ω) coefficients, calculated considering the individual wires. Namely, one reads withα In the latter equation, r (m) ξ 1 and r (m) ξ 2 take the same definition as in Eq. (C2), adapted for (x, z) = (x 1,2 , z 1,2 ).
It is worth discussing how relevant is this refinement of the theoretical model, when it comes to practical purposes. Already looking at Eq. (C2), it is clear that the relative contribution of ξ rapidly becomes negligible with increasing m, and may be fundamentally irrelevant even for m = 0, if x or z are sufficiently large. In fact, Ref. [19] studied the influence of the heater width in the static regime, under the approximation of a very thick substrate, finding that an increasing width causes only a modest decrease in the amplitude of the α coefficient, unless unreasonably large values for w are considered.
As a notable example, here we report in Fig. 8 the modulus of the frequency response |α(ω)|, calculated for a device with the same geometrical and material parameters of device D in Table I, but assuming different widths w of the heater. In particular, the black line is calculated using Eq. (23) (wirelike heater, w → 0), while the other curves are calculated using Eq. (C3) and different values for w. One observes that increasing the heater width to several tens of microns, i.e., to a size that becomes comparable to or larger than the relevant dimensions of the device such as x 1,2 , h, z, causes an overall decrease of |α(ω)|, with further attenuation of the highest frequency. Namely, a very large heater makes the device less efficient and slower.
However, one should consider that the actual device, which was used in our experiments, featured a heater with w = 9 μm, corresponding to the red dashed curves in the figure. Even though w is not of negligible size, compared to z = 8.8 μm or h = 20 μm, the correspondent normalized frequency response in Fig. 8(b) is hardly discernible from the one calculated adopting the wirelike assumption. The amplitude of the response in static conditions [which can be evaluated from Fig. 8(a)] is about 96% of the amplitude in the case of a wirelike heater.  Table I, and analogous material properties. The black line corresponds to the assumption of a wirelike heater [response calculated using Eq. (23)]. The other lines correspond to different heater widths w, as indicated in the labels in (a) [Eq. (C3) was employed here]. In (a) the modulus of all responses is normalized to the static limit of the wirelike case, to facilitate a comparison of the amplitudes. In (b) each response is normalized to its own static limit.

023094-13
In conclusion, one could observe, first, that too large heaters are typically not desirable, because they slow down the response of thermo-optic MZIs and they reduce their efficiency. Indeed, in real devices the heater width is kept limited, compatibly with constraints on fabrication tolerances and on the maximum thermal power that can be dissipated on the unit area. In any case, the good agreement of the model predictions with the experimental data (shown in the main text), together with the considerations reported here, suggest that the assumption of a wirelike heater does not constitute a limitation of our analytical model in many cases of practical interest.