Designable transverse thermoelectricity through interlayer thermal coupling for conformal heat harvesting

Materials and reagents

Polydimethylsiloxane (PDMS, Sylgard 184) was purchased from Dow Corning Co., Ltd. Polyimide (PI) films and polytetrafluoroethylene (PTFE) films were commercially obtained from RAYITEK Hi-Tech Film Co., Ltd. Silica aerogel ultrathin films were purchased from Nanuo Technology Co., Ltd. Thermal rubbers with different thermal conductivities were purchased from Kafuter New Materials Technology Co., Ltd. Bi₂Te₃-based thermoelectric plates (15 mm × 15 mm × 1 mm) were purchased from Quanzhou Qijin New Material Co., Ltd. and characterized before assembly [Supplementary Figure 1]. Flexible printed circuit (FPC) boards were custom-fabricated by JLCPCB (Shenzhen Jialichuang Technology Co., Ltd.) and adopted as the circuit substrates. Anhydrous ethanol was purchased from Sinopharm Chemical Reagent Co., Ltd. Low-temperature curing conductive silver paste was purchased from Kaitesen Co., Ltd. The detailed parameters of the films and the corresponding coupling parameters are listed in Supplementary Table 1.

Fabrication of flexible thermal metamaterial devices

The Al₂O₃ skeleton was fabricated by stereolithography (SLA) 3D printing using a photosensitive resin mixed with high-purity alumina powder, followed by sintering at 1,400 °C for 4 h to obtain a tilted multilayered framework. The skeleton was regionally infused with silicone rubbers of different thermal conductivities according to the designed lamination configuration. The composite was vacuum-degassed for 30 min to remove residual air bubbles and fully cured at 80 °C for 24 h. Detailed structural specifications and parameters of all thermal conductive layers are summarized in Supplementary Table 2. To prepare the Bi₂Te₃ thermoelectric plate for electrode fabrication, the sample surface was sequentially polished with fine sandpaper and briefly etched in 5% diluted hydrochloric acid to remove native oxide layers, then rinsed with anhydrous ethanol and dried naturally. Four nickel pads (1 mm × 1 mm) were evaporated onto the corners of the processed Bi₂Te₃ plate as diffusion barriers; copper leads were subsequently soldered and reinforced with conductive silver paste, followed by thermal curing at 80 °C for 2 h. The bilayer thermoelectric chip was assembled by sequentially stacking the tilted thermal layer, an interchangeable coupling film, and the processed Bi₂Te₃ thermoelectric plate. The peripheral gaps were filled with uncured PDMS and cured at 80 °C for 2 h to form an integrated device package. For flexible device applications, the assembled thermoelectric chip was pre-soldered onto an FPC board. Three-dimensional schematic illustrations of the device architecture were rendered using KeyShot Studio 2025.

The thermal measurement platform comprised a positive temperature coefficient (PTC) heating plate and a Peltier cooler positioned on the left and right sides, respectively, to establish a one-dimensional heat flow along the x-direction^[31]. Two wedge-shaped copper blocks were machined and inserted between the heat source/cold sink and the sample to concentrate heat transfer and ensure uniform thermal contact. The hot and cold terminal temperatures were monitored by K-type thermocouples and stabilized by a PID temperature controller to maintain a constant temperature difference across the device.

Finite element simulation

The finite-element simulations were performed using COMSOL Multiphysics 6.2. The lower thermal layer was modeled as an anisotropic medium with a uniform off-diagonal thermal conductivity tensor determined by equivalent medium theory^[32]. Heat convection, radiative heat loss, and interlayer thermal contact resistance were neglected. The temperatures of the hot side and cold side were set as T_h = 353 K and T_c = 293 K, respectively, with the remaining boundaries adiabatic. Detailed material properties were provided in the main text. Theoretical calculations and post-processing were performed using MATLAB R2022a.

Mechanical characterization

Mechanical durability of the flexible bilayer devices was evaluated through cyclic bending tests at a fixed radius of 50 mm and static bending tests at various radii down to 20 mm. Cross-sectional micrographs before and after 150 bending cycles were acquired using a stereomicroscope (Olympus SZX16) to assess interfacial integrity.

Thermal and electrical characterization methods

The interlayer coupling strength was varied by replacing the coupling layer material, while the temperature difference was precisely regulated by a temperature control driver board. Thermal measurements were performed using an infrared camera (Fotric 287+) to capture the surface temperature distribution, and the acquired data were imported into OriginPro 2024 for subsequent processing, including the plotting of isotherms and temperature contour maps with color bars. Electrical measurements were conducted by a two-pair shorting method: corresponding electrodes were shorted to measure transverse and longitudinal thermoelectric voltages. Reference electrodes were connected through a preamplifier (100×), and the voltage and current were recorded by a 5½-digit digital multimeter (Keysight U3606A). Each data point was measured five times (n = 5) and averaged. Error bars in the figures denote the standard deviation (SD) of the five repeated measurements. Photographs and circuit wiring diagrams of the voltage measurement system are provided in Supplementary Figure 2. The optical image of the experimental setup for characterizing the transverse thermoelectric performance of the flexible device is shown in Supplementary Figure 3.

Thermoelectric conversion efficiency evaluation

The conversion efficiency η was evaluated by the direct heat-balance method. The maximum output power P_max was extracted from the measured current–voltage characteristics under matched-load conditions. The input thermal power Q_in was evaluated via Fourier’s law from the longitudinal temperature gradient |∂T/∂x| measured along the edge of the upper thermal conductive filling, where the heat flux is predominantly horizontal. The conversion efficiency follows as η = P_max/Q_in.

Theoretical framework of transverse thermoelectricity via interlayer thermal coupling

Consider a bilayer system stacked along the z-direction with a temperature difference between high temperature T_h and low temperature T_c applied along the x-direction, and the planar dimension of the system is L × a [Figure 1B]. The lower layer (Layer 2) of thickness d₂ is an anisotropic thermal conductor with an off-diagonal thermal conductivity tensor κ₂, while the upper layer (Layer 1) of thickness d₁ is an isotropic thermoelectric layer with scalar thermal conductivity κ₁. The two layers exchange heat through a thin coupling layer of thickness d_c and equivalent thermal conductivity κ_c. The geometric configuration of this bilayer thermal coupling system is illustrated in Supplementary Figure 4. Under steady-state conditions and neglecting the through-thickness temperature gradient within each layer [Supplementary Figure 5]. The governing heat balance equations read^[33]:

where T₁ and T₂ denote the temperature fields of the upper thermoelectric layer and the lower conductive layer, respectively. The temperature distribution within the lower layer rotates clockwise by an angle θ about the z-axis to form an inclined thermal field [Figure 1C]. The substantial thermal conductivity mismatch between the two layers drastically suppresses thermal feedback from the thermoelectric layer to the lower thermal layer. Accordingly, the thermal field of the lower layer is simplified as a predefined oblique field: T₂(x,y) = T₀ - $$ \frac{\Delta T}{L} $$(xcosθ + ysinθ), where T₀ = (T_h + T_c)/2 is the central temperature, and ΔT = T_h - T_c is the applied total temperature difference. Introducing the normalized coordinates ξ = x/L, η = y/a and the dimensionless temperature τ = (T₁ - T₀)/ΔT, the upper-layer Equation (1) reduces to an inhomogeneous Helmholtz-type equation^[34,35]:

where the dimensionless coupling parameter Γ = κ_cL²/κ₁d₁d_c represents the coupling strength and β = a/L refers to the width-to-length ratio. Physically, Γ describes the degree of interlayer thermal coupling, defined as the ratio of the characteristic in-plane thermal resistance of the upper layer L²/κ₁d₁ to the coupling layer thermal resistance d_c/κ_c that governs heat exchange between the two layers.

Under adiabatic lateral boundary conditions (∂τ/∂η = 0 at η = ±1/2) and longitudinal Dirichlet boundary conditions (τ = ±1/2 at ξ = ±1/2), the analytical solution is obtained via separation of variables [Supplementary Note 1]:

where λ_n = $$ \sqrt{\Gamma +(n\pi /\beta )^2} $$. This solution explicitly reveals the dependence of the upper-layer thermal field on the coupling parameter Γ and the deflection angle θ of the lower layer. The effective deflection angle ϕ of the upper-layer thermal field is defined by the local temperature-gradient direction as tanϕ = (∂T/∂y)/(∂T/∂x) = β^-1(∂τ/∂η)/(∂τ/∂ξ). As Γ → 0, the upper-layer temperature field degenerates into a linear distribution along the x-axis and ϕ → 0. When Γ → ∞, the thermal deflection angle of the upper layer infinitely approaches that of the lower layer, namely ϕ → θ [Figure 1D]. Therefore, for any fixed deflection angle θ of the lower layer, continuous adjustment of the interlayer coupling parameter Γ allows the thermal deflection angle ϕ of the upper layer to be tuned smoothly from 0 to θ. The parametric evolution of ϕ with respect to Γ and θ verifies that the full range of deflection angles can be realized by regulating the interlayer thermal coupling [Figure 1E].

Under this oblique thermal excitation, the originally isotropic upper-layer thermoelectric material exhibits an effective anisotropic transport behavior. Within the linear Onsager framework^[36], denoting the longitudinal and transverse temperature differences across the upper layer as ΔT_x and ΔT_y, respectively, which can be evaluated from the analytical temperature field as: ΔT_x = $$ \int_{-L/2}^{L/2}\frac{\partial T_1}{\partial x}dx=\Delta T \int_{-1/2}^{1/2}\frac{\partial \tau }{\partial \xi }d\xi $$ and ΔT_y = $$ \int_{-a/2}^{a/2}\frac{\partial T_1}{\partial y}dy=\Delta T \int_{-1/2}^{1/2}\frac{\partial \tau }{\partial \eta}d\eta $$, its effective Seebeck tensor can be written as:

where S is the intrinsic Seebeck coefficient. Notably, the off-diagonal component S_yx^* = S·ΔT_y/ΔT_x is not an intrinsic material property, but an extrinsic response purely induced by the interlayer geometric thermal coupling^[8,25]. Consequently, the effective transverse power factor and figure of merit obey [Supplementary Note 2]:

Here PF₀ = S²σ and ZT₀ = S²σT/κ₁ are the intrinsic power factor and figure of merit, respectively, and σ is the electrical conductivity. Notably, ZT₀ reflects only the intrinsic properties of the upper thermoelectric layer; the influence of the lower and coupling layers is captured by the transverse-to-longitudinal temperature ratio ΔT_y/ΔT_x rather than by modifying ZT₀ itself. The above results demonstrate a hierarchical amplification mechanism: the off-diagonal Seebeck coefficient S_yx^* is activated linearly by the transverse-to-longitudinal temperature ratio ΔT_y/ΔT_x, while the transverse power factor and figure of merit are amplified quadratically by the same ratio. As plotted in Figure 1F, the blue straight line confirms linear transverse voltage output over ΔT_y/ΔT_x ∈ [-1,1], and the symmetric red parabola reflects quadratic power-factor variation. By tuning both the interlayer coupling parameter Γ and the thermal deflection angle θ of the lower thermal layer, the gradient ratio can be set to arbitrary values within the range [-1,1], enabling flexible manipulation of transverse thermoelectric output performance. This mechanism creates an innovative design route for reconfigurable flexible thermoelectric devices.

Numerical simulation of interlayer thermal coupling

To validate the theoretical framework, we performed finite-element simulations in COMSOL Multiphysics with preset parameters matching the analytical model: L = 15 mm, a = 15 mm, d₁ = d_c = d₂ = 1 mm, θ = 45°, κ₁ = 2.25 W/(m·K), κ₂ =30 W/(m·K), T_h = 353 K and T_c = 293 K. The bilayer geometry is illustrated in Figure 2A. As a reference, the uncoupled upper layer (Γ → 0) exhibits a purely longitudinal temperature distribution T_ref with isotherms perpendicular to the x-axis. The dimensionless coupling parameter Γ = κ_cL²/κ₁d₁d_c can be decomposed into two independent factors: L²/κ₁d₁, an intrinsic structural factor determined by geometry and upper-layer thermal conductivity, and κ_c/d_c, an extrinsic experimental knob governed by the coupling layer’s material and thickness. Physically, Γ represents the ratio of interlayer thermal conductance to the characteristic in-plane conductance of the upper layer, thereby quantifying the degree of thermal communication across the stack. The distribution of Γ over the parameter plane spanned by κ_c/d_c and L²/κ₁d₁ is plotted, with white dashed isolines at Γ = 10 and Γ = 80 partitioning the domain into weak (Γ < 10), moderate (10 ≤ Γ ≤ 80), and strong (Γ > 80) coupling regimes for subsequent analysis [Figure 2B].

Figure 2. Thermal manipulation via interlayer coupling and corresponding simulation results. (A) Schematic illustration of the influence of bottom thermal field rotation by θ on the thermal reference value T_ref of the upper layer; (B) The dimensionless coupling parameter Γ as a function of κ_c/d_c and L²/κ₁d₁. The white dashed lines denote the contours at Γ = 10 and Γ = 80, delineating the boundaries between weak (Γ < 10), moderate (10 ≤ Γ ≤ 80), and strong (Γ > 80) coupling regimes; (C-K) Temperature distributions for three typical coupling strengths, all simulated with identical parameters: L = 15 mm, a = 15 mm, d₁ = d_c = d₂ = 1 mm, θ = 45°, κ₁ = 2.25 W/(m·K), κ₂ = 30 W/(m·K), T_h = 353 K, T_c = 293 K. Only the coupling-layer thermal conductivity κ_c is varied: κ_c = 0.02 W/(m·K) for weak coupling, κ_c = 0.2 W/(m·K) for moderate coupling, and κ_c = 1 W/(m·K) for strong coupling. (C-E) Temperature contour of the upper layer (C), transverse and longitudinal temperature profiles comparing simulated and expected results (D), and temperature difference relative to the reference value (E) under weak coupling. Corresponding analyses for moderate and strong coupling cases are presented in (F-H) and (I-K), respectively; (L) Dependence of the thermal deflection at the center on the aspect ratio L/a under different coupling strengths. The red, black, and blue curves correspond to Γ= 40, 10, and 2, respectively; (M) Continuous evolution of the thermal deflection transfer efficiency with increasing Γ. The blue, black, and red curves correspond to θ = 30°, 45°, and 60°, respectively; (N) Y-direction temperature difference (black line) and magnitude of heat flux (red line for q_x, blue line for q_y) as functions of Γ, where the heat flux is normalized to that in the x-direction.

Figure 2C-E illustrates the simulation results under the weak-coupling condition [κ_c = 0.02 W/(m·K), Γ = 2]. The upper-layer temperature contour [Figure 2C] shows nearly vertical isotherms, indicating that the thermal field is only weakly perturbed by the 45° lower-layer field. In good agreement with simulation and theoretical curves, the transverse and longitudinal temperature profiles along the center axes reveal a longitudinal temperature span close to the full 60 K applied bias, whereas the transverse difference ΔT_y remains below 5% of ΔT_x [Figure 2D]. The deviation field T - T_ref exhibits a faint twist with peak deviations less than 2 K [Figure 2E]. In this regime, the coupling layer acts effectively as a thermal insulator, and heat transfer is dominated by one-dimensional longitudinal diffusion. Under moderate coupling [κ_c = 0.2 W/(m·K), Γ = 20], the isotherms visibly tilt toward the 45° orientation, and the temperature field develops a pronounced curvature near the lateral boundaries [Figure 2F]. The center-axis profiles show that ΔT_y rises to approximately 33% of ΔT_x, reflecting a competition between in-plane longitudinal diffusion and transverse heat injection from the lower layer [Figure 2G]. The deviation field quantifies the progressive redistribution of thermal energy from the x-direction to the y-direction, with peak deviations exceeding 8 K [Figure 2H]. In the strong-coupling regime [κ_c = 1 W/(m·K), Γ = 100], the upper-layer isotherms effectively lock onto the lower-layer oblique angle [Figure 2I]. The center-axis profiles show that ΔT_y exceeds 52% of ΔT_x [Figure 2J]. The deviation field displays a fully developed twist with peak deviations surpassing 18 K [Figure 2K]. The consistency between simulation and analytical theory is further corroborated by off-center line scans at x = ±L/4 and y = ±a/4, where the maximum relative deviation remains below 3% across all cut-lines [Supplementary Figure 6].

The central thermal deflection angle is maximized at an aspect ratio of unity (L/a = 1) and gradually declines as the aspect ratio increases further [Figure 2L], indicating that a square active area provides the most efficient geometric channel for transverse thermal transfer. Figure 2M shows that the thermal deflection transfer efficiency tanϕ/tanθ increases monotonically with Γ, converging toward unity at strong coupling. Saturation at large Γ arises because, once interlayer resistance becomes negligible, further increases in κ_c/d_c yield diminishing returns in deflection angle; the upper-layer field cannot exceed the lower-layer driving field. Figure 2N demonstrates that the x-direction heat flux q_x remains virtually invariant (varying by less than 1%) across all Γ, confirming that interlayer coupling acts primarily as a directional reallocating mechanism rather than a dissipative loss channel. Meanwhile, the y-direction temperature difference and the associated transverse heat flux q_y rise monotonically with Γ. This amplified transverse thermal bias directly corresponds to an enhanced effective off-diagonal thermal conductivity κ_yx^* in the upper layer. As derived in Supplementary Note 3, κ_yx^*/κ₁ = ΔT_y/ΔT_x, which scales linearly with the transverse-to-longitudinal temperature ratio. Supplementary Figure 7 further confirms that this linear scaling holds across the full angular range from -90° to +90°, including negative deflection angles, validating the all-azimuthal capability of the interlayer coupling mechanism. For Γ = 100, κ_yx^* reaches approximately 60% of the intrinsic κ₁, thereby establishing the thermodynamic link between interlayer thermal engineering and the emergent anisotropic transport that underpins the subsequent transverse thermoelectric response.

Experimental demonstration of interlayer thermal coupling

To experimentally validate the interlayer thermal coupling mechanism, anisotropic thermal conductive layers with engineered oblique thermal fields were fabricated. Notably, we measured the material properties of the experimental thermoelectric plate initially [Supplementary Figure 1], and corrected the simulated and designed model based on the actual characteristics. The key thermal conductivity adopts κ₁ = 2.0 W/(m·K). Then, by engineering the off-diagonal thermal conductivity to satisfy tanθ = κ_xy/κ_yy, a uniform oblique thermal field is obtained; however, in a finite rectangular domain, this anisotropic conduction inevitably suffers from edge distortion^[37]. To eliminate this boundary effect, we designed gradient-index boundary layers using coordinate transformation theory^[38], constructing a virtual parallelogram space and extracting a central rectangular active area while mapping the lateral regions into impedance-matching transition zones [Supplementary Figure 8]. Crucially, the virtual spaces for all deflection angles share identical projected lengths along the x-direction, ensuring consistent heat flux densities regardless of the deflection angle [Supplementary Note 4]. The spatially varying thermal conductivity tensors in the transition regions are continuous and smooth [Supplementary Figure 9], and simulations confirm that the integrated boundary architecture successfully eliminates edge distortion throughout the central domain [Supplementary Figure 10].

To reconcile the continuous theoretical tensor with practical fabrication constraints, the boundary transition zones were discretized into piecewise-homogeneous elements via the finite-element midpoint method^[35,39] [Supplementary Figure 11], where the Jacobian matrix is evaluated at the edge midpoints of each quadrilateral element to yield a uniform effective conductivity [Supplementary Note 5]. Supplementary Figures 12 and 13 present the calculated tensor component distributions under progressively refined mesh densities for both 30° and 60° configurations, and Supplementary Figure 14 demonstrates that the central active region consistently exhibits the target thermal gradient deflection regardless of mesh refinement, validating that the boundary layers function primarily as impedance-matching transition zones where moderate discretization is tolerated. Each discrete element was physically realized by laminated composites with tailored orientation angles and volume fractions [Supplementary Note 6]. We sequentially laminate the coupling and thermoelectric layers onto the thermal layer to fabricate the thermoelectric chip, then load predefined thermal boundaries to carry out interlayer thermal coupling tests [Figure 3A]. On the test bench, cured thermal conductive silicone gel fills the side areas next to the thermoelectric layer to complete boundary construction and guarantee intimate thermal contact. The corresponding effective thermal conductivity tensors are computed via transformation theory, enabling consistent thermal impedance matching throughout the whole device stack.

Figure 3. Experimental characterization of interlayer thermal coupling. (A) Schematic of the experimental setup. The inset presents an exploded-view illustration of the thermally coupled transverse TE chip; (B) Optical images of the thermal conductive layers with engineered deflection angles of 0°, 30°, and 60° (without thermal grease filling; scale bar, 10 mm) and the corresponding infrared thermal images showing the temperature fields; (C) Simulated (top) and experimental (bottom) temperature distributions of the thermoelectric layer at Γ = 135 for the three deflection angles (θ = 0°, 30°, and 60°); (D) Transverse temperature profiles along the y-axis of the upper layer for three different lower-layer deflection angles θ, comparing experimental data with simulation results. Red, blue, and black circles denote experimental data for θ = 60°, 30°, and 0°, respectively; dashed lines represent the corresponding simulated results; (E) The thermal deflection ratio |$$ \nabla $$_yT/$$ \nabla $$_xT| vs. coupling parameter. Red and blue lines denote simulated results for θ = 60° and 30°, respectively; red and blue diamonds denote the corresponding experimental data. Insets show temperature distributions at representative observation points; (F) Bottom panel: longitudinal (Δ_xT) and transverse (ΔT_y) temperature differences and transverse thermal power (Q_y) of the thermoelectric layer under different temperature gradients; top panel: schematic of the calculated thermal deflection ratios. Red and blue colors correspond to θ = 60° and 30°, respectively. TE: Thermoelectric; IR: infrared.

Optical photographs of the fabricated thermal layers and the corresponding infrared thermal images are presented [Figure 3B]. Within the central 15 mm × 15 mm square area, the thermal layers with deflection angles successfully establish a uniform oblique thermal field consistent with the designed angle, featuring parallel isotherms of constant slope, whereas the 0° control sample maintains straight longitudinal isotherms. Near the edges, the isotherms smoothly transition from the oblique orientation in the bulk to a nearly vertical alignment at the outer boundaries, evidencing that the transformation-based boundary layers effectively suppress edge distortion and maintain the intended heat flux topology throughout the active region. We subsequently characterized interlayer coupling by placing the upper thermoelectric layer atop the substrates separated by a 0.25 mm PI coupling film [κ_c = 0.3 W/(m·K)], yielding Γ = 135 based on κ₁ = 2.0 W/(m·K). Simulated and measured upper-layer temperature fields for thermal layer deflection angles of 0°, 30°, and 60° are compared in Figure 3C; the experimental infrared contours closely match the predicted oblique patterns. Side-view infrared thermography further confirms that the temperature remains nearly uniform across the thickness of both the thermoelectric and thermal layers, with the dominant thermal resistance confined to the thin coupling interlayer [Supplementary Figure 15].

Figure 3D presents quantitative line scans along the central y-axis of the upper thermoelectric layer. The experimental transverse temperature profiles exhibit excellent agreement with simulation results across all three deflection angles (0°, 30°, and 60°), with a maximum deviation of less than 5% over the entire scanning range, thereby validating the predictive capability of the theoretical model. To systematically tune the dimensionless coupling parameter Γ across the full dynamic range, we selected a series of coupling media spanning a broad thermal conductance range (see Supplementary Table 2 for detailed parameters). Figure 3E plots the measured thermal deflection ratio |$$ \nabla $$_yT/$$ \nabla $$_xT| as a function of the coupling parameter Γ. For the 30° thermal layer, the ratio increases steadily with Γ: it reaches 0.08 at Γ = 5, 0.24 at Γ = 25, 0.28 at Γ = 61, 0.29 at Γ = 90, and ultimately approaches saturation near 0.3 at Γ = 135. The 60° thermal layer exhibits a similar monotonic rise, tending toward its own saturation plateau under strong coupling. Supplementary measurements at varied Γ further confirm that the upper thermal field gradually shifts from longitudinal to tilted alignment with growing coupling strength [Supplementary Figure 16]. To assess non-ideal effects, complementary finite-element simulations were performed incorporating interfacial contact thermal resistance (R_s), natural convection [h = 5 W/(m²·K)], and surface radiation (ε = 0.4) [Supplementary Note 7 and Supplementary Figure 17]. The results demonstrate that convective and radiative losses produce temperature deviations below 0.3 K under laboratory conditions and are therefore negligible. The dominant factor is the interfacial contact thermal resistance, which acts in series with the coupling layer and reduces the effective coupling parameter to Γ_eff = Γ/[1 + (ΓR_sκ₁d₁/L²)]. This correction indicates that contact thermal resistance is the primary source of the systematic deviation between experimental data and ideal theoretical predictions in Figure 3E. In addition, the impedance-matching boundary transition zones surrounding the active area account for the saturation values of the thermal deflection ratio lying below the ideal strong-coupling limit tanθ. Figure 3F presents the thermal response of the thermoelectric layer under various applied temperature gradients. All measurements were performed on the same device and with the same 0.25 mm PI coupling film (nominal Γ = 135) as in Figure 3C-E, with only the applied temperature difference varied. The bottom panel summarizes the longitudinal and transverse temperature differences together with the transverse thermal power of the thermoelectric layer, which is defined by the product of the heat flux density and the cross-sectional area. For the 60° thermal layer configuration, the longitudinal temperature difference ΔT_x increases from 6 to 11.5 K, while the transverse temperature difference ΔT_y rises from 3.1 to 5.3 K as the overall thermal bias is enlarged. The top panel shows the theoretically calculated heat flux deflection, which remains nearly constant across the measured range. Notably, the thermal deflection ratio |$$ \nabla $$_yT/$$ \nabla $$_xT| stays approximately 0.5 (0.52 at ΔT_x = 6 K and 0.46 at ΔT_x = 11.5 K), indicating that the thermal field deflection ratio is preserved despite the increasing temperature span. This constancy confirms that the interlayer thermal coupling mechanism is robust against fluctuating thermal loads and that the oblique thermal field is reliably transmitted to the upper layer for subsequent thermoelectric conversion. Collectively, the systematic variation of deflection angle θ and coupling parameter Γ demonstrates reconfigurable multi-field manipulation of the upper thermoelectric layer, where the thermal field orientation, the induced transverse temperature gradient, and the resulting heat flux distribution are all continuously tunable without altering the intrinsic material properties.

Interlayer thermal coupling induced transverse thermoelectric response

To isolate the extrinsic transverse thermoelectric response from the intrinsic longitudinal Seebeck effect, thermoelectric chips with a four-electrode layout were fabricated^[40,41] [Figure 4A]. Numerical simulations of the electric potential distribution reveal that, as the lower-layer thermal deflection increases from 0° to 60°, the equipotential lines in the upper thermoelectric layer rotate from a strictly longitudinal orientation to progressively oblique configurations, while the transverse electric field component evolves from a negligible value to a pronounced single-peak profile with increasing amplitude [Figure 4B].

Figure 4. Transverse thermoelectric conversion driven by interlayer thermal coupling. (A) Schematic of the thermoelectric chip with four-probe electrode layout for oblique thermal coupling characterization; (B) Simulated electric potential contours (upper panels) and corresponding vertical-center-axis electric field profiles (lower panels). The red and blue lines denote the electric potential V and the transverse electric field E_y, respectively; (C) Optical photograph of the fabricated four-probe device mounted on the measurement stage; (D) Effective off-diagonal Seebeck coefficient induced by interlayer thermal coupling plotted against the coupling parameter Γ. The colored curves represent different lower-layer deflection angles: red: θ = 60°, black: θ = 45°, blue: θ = 30°, pink: θ = 15° and green: θ = 0°; (E) Calculated distribution of the transverse effective figure of merit ZT_yx^*/ZT. White dashed lines mark the representative deflection angles θ = 30° and θ = 60° adopted for experimental characterization; (F) Experimental open-circuit voltages and numerical predictions for longitudinal (V_x, red) and transverse (V_y, blue) configurations at θ = 30° under an applied temperature difference of 40 K. Error bars denote the SD of experimental replicates; (G) Corresponding experimental and theoretical voltages for θ = 60°; (H) Experimentally measured transverse thermoelectric conversion efficiency η (bar charts) and maximum output power P_max (red line) as functions of the applied temperature difference. SD: Standard deviation.

The fabricated device mounted on the measurement stage is shown in Figure 4C. To eliminate the influence of contact and lead resistances, a two-pair shorting method was adopted: for the longitudinal measurement, electrodes V₁ and V₂ were shorted together to serve as one terminal while electrodes V₃ and V₄ were shorted to serve as the other, thereby probing the voltage drop along the x-direction; for the transverse measurement, electrodes V₁ and V₃ were shorted as one terminal and electrodes V₂ and V₄ as the other, capturing the voltage drop across the y-direction. This four-probe configuration ensures that the measured voltages reflect purely the thermoelectric potential generated by the oblique thermal field, free from spurious ohmic drops at the interfaces. The effective off-diagonal Seebeck coefficient S_yx^* induced by interlayer thermal coupling increases monotonically with the coupling parameter Γ and eventually saturates at strong coupling, consistent with the theoretical prediction in Equation (4) [Figure 4D]. The physical origin of this saturation lies in the asymptotic behavior of the transverse-to-longitudinal temperature ratio ΔT_y/ΔT_x. Furthermore, varying the coupling-layer thickness of silica aerogel and PDMS films at θ = 60° yields a monotonic decrease in S_yx^*/S with increasing thickness [Supplementary Figure 18], consistent with Γ $$ \propto $$ 1/d_c. The corresponding thermoelectric figure of merit ZT_yx^*/ZT indicates that the transverse thermoelectric performance is maximized under the joint action of strong coupling and large deflection angles [Figure 4E].

Experimental open-circuit voltages for longitudinal and transverse configurations at lower-layer deflection angles of 30° and 60° under an applied temperature difference of 40 K are presented [Figure 4F and G]; complementary data for 15° and 45° are provided in Supplementary Figure 19. The longitudinal voltage remains nearly constant across all deflection angles, as it is governed primarily by the applied longitudinal thermal bias ΔT_x. In contrast, the transverse voltage rises substantially from approximately 0.2 to 1.0 mV at 30°, and from approximately 0.5 to 1.5 mV at 60°. The monotonic growth of the transverse voltage with θ directly validates that the transverse thermoelectric response is extrinsically engineered by the interlayer thermal coupling rather than by any intrinsic material anisotropy. The experimentally measured transverse thermoelectric conversion efficiency η and the maximum output power P_max under different applied temperature gradients are summarized in Figure 4H. Bar-chart data for additional deflection angles are provided in Supplementary Figure 20. Under the 60°, Γ = 135 configuration, the conversion efficiency increases monotonically from 0.19% at ΔT = 20 K to 0.53% at ΔT = 60 K, while the maximum transverse output power grows from 0.011 to 0.107 mW over the same span. The modest efficiencies stem from parasitic dissipation at interlayer interfaces and boundary transition zones, which attenuates the central thermal gradient and causes the experimentally determined conversion efficiency to fall below the standard thermoelectric efficiency relation under identical boundary conditions. Systematic analysis of the normalized off-diagonal Seebeck coefficient across the full parameter space reveals that interfacial contact thermal resistance constitutes the dominant baseline loss, while boundary-layer fabrication tolerances become increasingly pronounced at larger deflection angles [Supplementary Note 8 and Supplementary Figure 21]. Finite-element simulations incorporating interfacial contact thermal resistance R_s reveal that both the transverse power factor and output power density decrease monotonically with increasing R_s [Supplementary Figure 17E and F], confirming that interfacial thermal resistance is the primary performance bottleneck. To verify that this bottleneck is not specific to the choice of thermoelectric material, complementary simulations were performed using PbTe as the upper layer. The resulting normalized off-diagonal Seebeck coefficient S_yx^*/S falls precisely onto the theoretical curve derived from the Bi₂Te₃-based model [Supplementary Figure 22], demonstrating that the transverse response is governed solely by interlayer thermal coupling geometry rather than by intrinsic material properties. Consequently, performance gains can be achieved by optimizing interlayer impedance matching and minimizing parasitic leakage through refined coupling-layer engineering, without relying on exotic material selection.

Conformal thermoelectric device for non-planar surfaces

Rigid thermoelectric modules inherently struggle to establish conformal thermal contact on curved surfaces, leaving substantial heat-harvesting opportunities untapped on non-planar geometries such as human skin, wearable electronics, and cylindrical industrial piping. The present bilayer design overcomes this limitation by physically separating thermal field engineering and thermoelectric conversion into independently optimizable layers, allowing each layer to be thinned, patterned, and made flexible without compromising the interlayer coupling mechanism. To extend the interlayer-coupled architecture beyond planar geometries, the assembled thermoelectric chip was integrated onto an FPC board to form a mechanically compliant conformal device^[42,43]. The FPC substrate, pre-patterned with copper traces and contact pads, enabled seamless integration of the Bi₂Te₃ thermoelectric layer, the coupling film, and the tilted thermal conductive layer into a single stack. A preamplifier circuit was co-mounted on the same FPC board to minimize lead length and electromagnetic interference. The total thickness of the flexible bilayer device was approximately 2.5 mm, and the bending stiffness was dominated by the thin PI substrate rather than the rigid ceramic components, ensuring that the stack remained mechanically compliant [Figure 5A]. Conformal attachment was demonstrated by mounting the device onto a cylindrical surface (radius of curvature, 50 mm) and by wearing it directly on a human arm [Figure 5B]. The flexible stack readily conformed to the curved geometry without delamination or cracking of the thermoelectric layer, confirming that the mechanical compliance was preserved after integration.

Figure 5. Conformal attachment and energy harvesting on curved surfaces. (A) Schematic illustration of the flexible device integrated on an FPC substrate. Inset: optical photograph; (B) Optical photograph of the conformally attached device on a cylindrical surface and worn on a human arm (scale bar, 3 cm). The photograph of the device worn on the human arm was taken with the informed consent of the volunteer; (C) Time-dependent transverse open-circuit voltage (V_y) measured at a fixed temperature difference of ΔT = 40 K. The red and blue curves correspond to θ = 60° and θ = 30°, respectively; (D) Thermal cycling stability of the conformal devices. Transverse open-circuit voltage (V_y) of the θ = 30° (blue) and θ = 60° (red) flexible devices recorded over 60 min under periodic thermal cycling at ΔT = 40 K; (E) Current–voltage characteristics and corresponding power densities (normalized to the 15 mm × 15 mm active area) for the 60° and 30° flexible devices. Red and blue colors correspond to θ = 60° and θ = 30°, respectively. Solid lines and scattered points denote simulation and experimental results, respectively, while dashed lines correspond to fitted curves; (F) Proof-of-concept demonstration of an angle-resolved heat-flux indicator. (Left) Time-dependent transverse output voltage of the device under a fixed thermal bias (ΔT = 40 K) for θ = 60° (red), 30° (blue), and 15° (pink), driving a dual-LED threshold circuit. (Right) Corresponding optical states of the indicator: dual-LED ON (top), red-LED ON (middle), and dual-LED OFF (bottom). FPC: Flexible printed circuit.

Two sets of flexible devices with lower-layer deflection angles of 30° and 60° were characterized. A uniform temperature difference of ΔT = 40 K was applied using the same PTC heater and Peltier cooler configuration. The transverse open-circuit voltage V_y was recorded continuously over a 1,200-second interval. As shown in Figure 5C, the voltage stabilized within 360 s at approximately 0.937 and 1.346 mV for the 30° and 60° configurations, respectively, with SDs below 0.08 and 0.12 mV. Thermal cycling over 60 min yields reproducible periodic output without obvious degradation [Figure 5D], confirming robust interlayer coupling under cyclic thermal loads. To address mechanical durability, bending-cycle tests were conducted at a fixed radius of 50 mm and static bending tests were performed at various radii down to 20 mm [Supplementary Figure 23]. Cross-sectional micrographs confirm intact layer interfaces without delamination after 150 bending cycles [Supplementary Figure 24]. The output current–voltage characteristics and corresponding power densities normalized to the 15 mm × 15 mm active area of the thermoelectric chip are presented in Figure 5E; solid lines and scattered points denote simulation and experimental results, respectively, while dashed lines correspond to fitted curves. Although the measured power densities fall approximately 10% below theoretical predictions owing to additional contact thermal resistance at the curved interface and lateral heat leakage through the PI substrate, the 60° device consistently outperforms the 30° counterpart, confirming that the interlayer thermal coupling mechanism remains effective in the flexible, non-planar geometry.

A side-by-side comparison with representative transverse thermoelectric strategies [Supplementary Table 3] highlights that the present architecture offers full-range angular tunability and conformal adaptability without requiring cryogenic temperatures or external magnetic fields, distinguishing it from topology-driven and magneto-thermoelectric alternatives. Furthermore, by integrating an amplification stage and a dual-LED threshold indicator onto the FPC substrate, the device directly converts heat-flux deflection angle into discrete optical states [Figure 5F]. Under a fixed thermal bias of ΔT = 40 K, the three configurations (θ = 60°, 30°, and 15°) yield distinct outputs: dual-LED ON above ~1.0 mV, red-LED ON above ~0.65 mV, and dual-LED OFF below the threshold, respectively. This proof-of-concept validates that the geometrically induced transverse voltage can directly drive low-threshold microelectronics, establishing a feasible route toward conformal, angle-resolved thermal sensors.

Limitations and outlook

This study demonstrates that interlayer thermal coupling can decouple thermal field engineering from thermoelectric conversion, enabling a continuously tunable transverse response in intrinsically isotropic materials. While the predictive capability of the dimensionless coupling parameter Γ has been validated across the full dynamic range, three principal limitations should be acknowledged. First, interfacial contact thermal resistance and fabrication tolerances of the thermal conductive layer cause the experimentally measured transverse thermoelectric performance to fall below ideal theoretical predictions, particularly at large deflection angles. Second, the present devices were demonstrated at preset, discrete deflection angles and coupling strengths; in-situ tunability of the thermal field within a single device remains to be experimentally demonstrated. Third, the current work is restricted to a single-chip architecture; scaling to large-area, multi-module conformal arrays remains an open engineering challenge.

Looking forward, several directions merit exploration. The integration of actively tunable thermal interfaces - such as thermally responsive polymers or pressure-adjustable coupling layers - offers a pathway to liberate the architecture from static film-based configurations^[44-47], enabling dynamic modulation of the transverse output without physical replacement of the coupling medium. For the thermal metamaterial layer, the prescribed anisotropic conductivity can be realized in monolithic thin films via scalable processing routes such as blade coating or roll-to-roll assembly of aligned two-dimensional nanofillers in polymer matrices, liberating the design from discrete lamination^[48,49]. Beyond wearable electronics, the conformal nature of the device architecture suggests immediate applicability to curved industrial surfaces: conformally wrapping the bilayer stack around exhaust pipes or cylindrical heat exchangers opens the possibility of harvesting otherwise dissipated thermal energy from non-planar geometries where conventional flat modules fail to make adequate thermal contact^[50,51]. Extension to organic and thin-film thermoelectric systems would further reduce weight and bending stiffness^[52], broadening the design space toward epidermal patches and soft implantable energy harvesters. Ultimately, tiling the architecture into distributed multi-module arrays with integrated energy management circuits establishes a scalable route toward large-area, shape-adaptive thermoelectric harvesting on arbitrarily curved surfaces, translating the present proof-of-concept into deployable soft energy systems^[53,54].