Thermal inverse design for resistive micro-heaters

Recent years have witnessed a surge of interest in inverse design, which generally refers to algorithms that automate the design of physical systems based on pre-defined performance targets. This contrasts with conventional forward design approaches, where multiple candidates are evaluated to down-select an optimized design, either by exhausting the entire design space or through human-guided design iterations. The inverse design paradigm has transformed many fields in engineering, leading to photonic, mechanical, acoustic, and magnetic devices with non-intuitive configurations and/or unprecedented functions such as nanophotonic devices,^1–3 metasurface optics, mechanical metamaterials,⁴ phononic crystals,⁵ and micro-magnet arrays.⁶ The inverse design has also been applied to engineering the transport properties in thermal metamaterials.^7,8

In this work, we aim to develop an inverse design strategy for an important class of devices: resistive micro-heaters, whose applications are pervasive in electronics, photonics, microelectromechanical systems (MEMS), and microfluidics. The designs of these heaters have been implemented via numerical optimization in a trial-and-error process, where the heater geometry is adjusted iteratively based on finite element method (FEM) simulation results of the temperature map across the heater.^9–14 Alternatively, the general Cauchy-type thermal inverse design problem, which seeks to optimize a heating device (comprising a spatial distribution of heat sources) to produce a pre-defined temperature profile, has also been approached using variational methods and genetic algorithms.^15–17 However, all these methods suffer from one intrinsic limitation: they necessitate modeling the entire heater device in each iteration. Such “global” simulations are computationally expensive, which limits either the domain size or the number of iterations that can be executed. As a result, the above-mentioned methods will only work for very small heaters or simple heat source configurations that can be parameterized with a small set of variables. To circumvent this limitation, our approach uses only parameterized thermal properties of unit cells to deterministically infer the heater design, thereby curtailing the need for “global” simulation iterations. Using doped silicon resistive heaters as an example, we show that our objective-driven approach allows automated generation of micro-heater designs with on-demand steady-state temperature profiles.

We present our inverse design scheme in doped silicon micro-heaters fabricated on a silicon-on-insulator (SOI) platform, although the methodology can be readily generalized to other types of resistive heaters. Doped SOI heaters are fully CMOS-compatible and can be seamlessly integrated with frontend electronic and photonic devices for applications spanning thermo-optic switching, material processing, gas sensing,^18–20 and beyond. Electrical current paths in doped Si heaters can be defined with the doping profile without altering the physical structure, a feature that is particularly useful for photonic devices since thermal and optical design considerations can be decoupled. These features ensure that our work not only informs a general inverse design approach but also has significant practical relevance.

The inverse design problem is formulated as follows. A steady-state temperature distribution T(x, y), where x, y are the in-plane coordinates on the heater surface, is defined as the design objective. Figure 1 depicts the heater layout under consideration, consisting of a SOI slab connecting two constant-biased electrical contacts. The target temperature profile is specified in a domain encircled by the red solid lines. The design flow starts with defining a mesh to divide the domain into a set of unit cells with the condition that the curvature of the temperature profile across a single unit cell is sufficiently small so that the average temperature can approximately represent the entire unit cell's temperature. The doping regions have a uniform dopant concentration, and they form an array of filaments with varying widths W(x, y), each traversing a continuous row of unit cells connecting the two contacts. This binary doping pattern can be fabricated by lithography followed by ion implantation. As we shall see later, given a heat generation rate profile $g( x , y)$ (defined as the heat generated in each unit cell per unit time in Watts), one can deterministically solve a unique solution W(x, y) for each and every row of unit cells.

Similarly, C₂ corresponds to the product of the equivalent in-plane thermal conductance and the unit cell side area, and $∇ x / yT$ denotes the temperature gradient in the in-plane x or y direction. Heat dissipation from each unit cell is a combination of convective heat transfer to air and conductive heat transfer to the surrounding solids. The heat generation $g( x , y)$ can be solved by summing the heat transfer rates over all directions in the steady state.

We then use FEM simulations shown in Fig. 2 to validate Eqs. (1) and (2) and derive the constitutive parameters C₁ and C₂. It is noteworthy that the simulations are not performed over the whole heater, but only across a small domain that contains the target unit cell (marked by red dotted lines) and its neighbors. In our example, the FEM simulation domain contains 5 × 5 square unit cells, each having a size of 1 × 1 μm². The small simulation domain size used in the step to extract the thermal constitutive parameters implies much lower computational overhead compared to full-scale modeling of the entire heater device. In the more general case, for an arbitrary heater geometry, the meshing process leads to unit cells of parallelogram shapes. In this case, a pair of C₁ and C₂ parameters need to be computed for each unique unit cell shape to create a lookup table that can then be repeatedly used for different target temperature profiles. The heater stack comprises (from top to bottom) air, a 10 nm thick silicon dioxide layer, 220 nm SOI layer, 3 μm buried oxide (BOX), and the Si substrate that acts as the heat sink. Stripes of doped Si (labeled with blue-gray color) with a fixed doping concentration of 10²⁰ cm⁻³ traverse the unit cell columns. The constitutive parameters only weakly depend on the doping region width, and, therefore, we only perform the simulation for two widths and use linear interpolation to estimate the parameters for other width values.

To obtain the constitutive parameter C₁, we set all 25 unit cells to have the same temperature, i.e., with zero in-plane temperature gradient. The simulation allows us to validate Eq. (1), i.e., that the out-of-plane heat rates linearly scale with the temperature difference T – T_R between the unit cell's temperature and room temperature. Figures 2(b) and 2(c) show the dependence of the heat transfer rates on the temperature difference T – T_R, where the sum of the slopes of the two lines yields C₁. Because of the layer structure, heat dissipation primarily occurs through the substrate. Next, a linear temperature gradient across the simulation domain is defined by setting the temperatures of all 25 unit cells. The heat transfer rates along x- and y-directions are then computed [Figs. 2(d) and 2(e)] to infer C₂. We also verify that the out-of-plane fluxes barely change with varying in-plane temperature gradients, which corroborates our hypothesis that the heat flow into a unit cell can be decomposed into independent components in in-plane and out-of-plane directions, respectively. These results support our hypothesis that the thermal behavior of the unit cell can be characterized using two constitutive parameters C₁ and C₂ with adequate accuracy.

The conclusion above holds provided that the two assumptions we made hold: (1) the temperature curvature within a unit cell is sufficiently low; and (2) the in-plane temperature gradient is much smaller than the out-of-plane temperature gradient. The former assumption breaks down as the unit cell size increases. The latter assumption is in general valid in doped Si heaters, since the in-plane thermal conductance (given the exceptionally high thermally conductivity of Si) is far larger than that in the out-of-plane direction.

The procedures above allow the inverse design of heaters when the target temperature and in-plane temperature gradient (necessary to solve the in-plane heat transfer rate from Eq. (2)] are known for each and every unit cell. However, most heater design problems do not specify the temperature gradient at the boundaries. In other words, there is no immediate way to ascertain if an a priori assumed temperature gradient at the boundary is realistic. Therefore, we must identify a method to solve a temperature distribution outside the boundary such that the transition from the “hot” heater to the surrounding “cold” area is physical. To resolve this issue, we propose an additional optimization process as summarized in Fig. 3. We define three zones in a heater: a target zone, a transition zone, and a peripheral zone. The target zone is the actual domain within which the inverse design is performed; the peripheral zone is the area where the temperature drops considerably such that the second-order derivative of the temperature profile at its outer boundary becomes negligible, and the transition zone has a smoothly varying temperature profile that bridges the two other zones. Initially, a temperature map T_target is defined. The T_target in the target zone is the design target, and the T_target in the other two zones is some smooth function as an initial guess. Next, we perform an inverse design step across all three zones to obtain an initial heater design. An FEM simulation is then executed across the entire domain to generate a temperature map T_response. The T_response in the peripheral zone is set as the new T_target for the next iteration. The T_target in the target zone remains the same in the subsequent iteration and an exponential function is used as the T_target in the transition zone to produce a smooth transition. The process reiterates until the T_response within the target zone meets the design specifications (more on this later). We found that the algorithm is robust and converges independent of the initial T_target heuristic in the peripheral zone, as we will demonstrate in the next section. The process also converges to a design with satisfactory accuracy quite quickly: in the examples we give in the following section, the maximum number of iterations used is only 5. Therefore, the computational load is well manageable. Even though not implemented in this paper, it is possible to further reduce the computation load of the optimization process by simulating only the transition zone and the peripheral zone. The temperature distribution within the target zone is used to define the boundary condition for such “partial” simulations.

Here, we present a set of heater design examples derived from our inverse design method. The first example is a square doped Si heater targeting uniform temperature distributions. Figure 4 compares the temperature profile from a reference heater with a uniform doping concentration [Figs. 4(a) and 4(c)] and that of an inverse-designed heater [Figs. 4(b) and 4(d)]. In both cases, the heater has an area of 38 × 38 μm². The target temperature profile is set to be 700 K across the center 30 × 30 μm² area. Figure 4(a) plots the temperature map of a reference heater with a uniform doping concentration, showing a significant temperature change of 102 K from the center to the edge within the 30 × 30 μm² target zone. In comparison, the inverse-designed heater exhibits a uniform temperature averaged at 699.7 K across the target zone with a root mean square (RMS) variation of only 0.18 K. The width map in nm is provided in the supplementary material (Table S1).

In this example, our design figure-of-merit (FOM) adopted during the boundary optimization process illustrated in Fig. 3 is given as 1/[(T – T_target)T_std], where T – T_target is averaged over the target zone and T_std is the standard deviation of temperature across the target zone. This FOM prioritizes temperature uniformity over the absolute temperature. As a result, the optimization produces a heater design with uniform a 705.5 K temperature. This deviation from the initial 700 K target can be easily corrected by slightly decreasing the applied voltage to obtain the results shown in Figs. 4(b) and 4(d). The optimization procedure is also largely agnostic to the choice of initial T_target. As an example, two overly simplistic T_target options [gray dashed lines in Fig. 4(d)] were tested: one is a uniform 700 K temperature profile throughout all three zones and one assumes a linear temperature gradient in the transition and peripheral zones that drops the temperature to ambient at the boundary. Both yield identical outcomes after the optimization procedure.

We further test uniform temperature heater designs with different sizes of target zones, 10 × 10 and 40 × 40 μm² [Figs. 5(a) and 5(b)]. It is worth pointing out that instead of executing the same optimization procedures in Fig. 3, we directly adopted the final T_target for the transition and peripheral zones from the prior heater design [the blue dot curve in Fig. 4(d)]. This important feature provides a facile means to scale heater designs while completely circumventing the need for full-heater-scale FEM simulations.

To demonstrate the versatility of our method, we apply the inverse design scheme to implement two non-uniform temperature profiles. The first example seeks to achieve a linear gradient temperature along the diagonal direction of a square heater as shown in Figs. 6(a) and 6(b). The target zone is 20 × 20 μm² while the whole heater size is 40 × 40 μm². Since the linear gradient temperature profile imposes a challenge that the local temperature gradient has to switch its sign from positive to negative as it moves from the target zone to the transition and peripheral zone, we have to utilize a larger transition zone so there is sufficient space for the temperature to decrease gradually.

The second example is a T_target map that features areas with depressed/raised temperatures assuming the shape of an MIT logo within a 60 × 40 μm² target zone. Figures 6(c) and 6(d) display the FEM-simulated heater temperature distributions that closely match our design target. These examples indicate that our method can be broadly applied to produce heater designs with on-demand temperature distributions¹² provided that the design target is physically viable.

The method is not limited to planar heater structures and can be applied to heaters with uneven surface morphologies as well. As an example, we consider a doped Si heater incorporating a ridge waveguide across its center as shown in Fig. 7. Similar heater configurations are commonly adopted in photonics, for instance, to switch the structure of phase change materials (PCMs),¹⁵ where the ability to precisely engineer temperature profiles is essential for controlling the materials' phase transition kinetics and avoiding damage due to hot spots.²¹ In this example, we assumed a ridge waveguide structure fabricated in a 220 nm thick silicon-on-insulator layer with an etch depth of 110 nm and a ridge width of 450 nm. The waveguide is covered by a 100-nm thick silicon dioxide cladding [Fig. 7(c)]. Two sets of constitutive parameters were computed corresponding to two types of unit cells: unit cells on the waveguide and unit cells outside the waveguide. With a design target of uniform temperature distribution at 700 K, the inverse design scheme leads to an average temperature of 700.3 K with a standard deviation of 0.16 K, as shown in Fig. 7(a), where the waveguide is highlighted by the white solid line and the target zone is within the red dashed line.

The examples we cited here all assume heaters with a rectangular shape. Our approach, however, can be readily extended to heaters with irregular geometries. In the latter case, the unit cells will no longer be rectangular but rather parallelogram in shape and might not be uniform in size. Thermal constitutive properties for each size and shape of unit cells need to be calibrated following protocols similar to that in Fig. 2 and compiled into a lookup table. They can then be used in conjunction with Eqs. (1) and (2) to determine the heat generation map and the subsequent design process flow is similar.

The accuracy of the inverse design approach is presently limited by the mesh size of the FEM simulations. Since the FEM model has a minimum mesh size which all structures must snap to, this inevitably results in the deviation of the local filament width from the analytically derived values from Eqs. (3) and (4). For the design illustrated in Figs. 4–7, the minimum mesh size (constrained by the computer memory capacity) is approximately 5 nm. This is comparable to, and in some cases considerably larger than the typical filament width variations between neighboring unit cells. This limitation partially accounts for the observed deviations from design targets in our examples.

In conclusion, we proposed and validated a generic inverse design framework to engineer doped Si micro-heaters' temperature profiles. The approach can be fully automated with minimal-to-none human intervention. Moreover, it circumvents the need for extensive full heater simulations and instead uses the thermal constitutive properties of small unit cells to construct the heater design. While we demonstrated the inverse design approach in doped Si heaters, the concept can be adapted to other resistive heater systems as well. Given the ubiquity of resistive heaters in microsystems, the ability to generate on-demand temperature distributions foresees many exciting potential applications in photonics, electronics, and micro-mechanics.

See the supplementary material for the details on the filament width map of the inversed designed doped silicon heaters showcased in Figs. 4–7.

The funding support was provided by the National Science Foundation (NSF) under Award Nos. 2329088 and 2225968. We cordially acknowledge Dr. Kiumars Aryana and Dr. Yifei Zhang for insightful discussions and technical assistance.

The authors have no conflicts to disclose.

Khoi Phuong Dao: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Juejun Hu: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.