On the role of thermo-electro-ionic dynamics in hysteresis and transient performance of perovskite solar cells

Abstract

Understanding the thermal origins of performance instabilities and hysteresis in perovskite solar cells (PSCs) is essential for advancing their long-term stability and reliable operation. In this perspective, we develop a novel coupled multiphysics mathematical framework that integrates layer-resolved optical absorption, non-isothermal electronic–ionic transport, and a layer-resolved, self-consistent energy balance with explicit bulk and interfacial heat-generation pathways. These pathways include hot-carrier thermalization, Joule heating, Peltier effects, radiative/non-radiative recombination, and parasitic optical absorption. This mathematical framework extends PSC characterization beyond conventional J–V analysis by introducing temperature–voltage (T–V) and heat–voltage (P–V) characteristics curves, enabling quantitative tracking of transient self-heating and its interactions with electronic–ionic dynamics. It is shown that PSCs develop internal thermal inertia that evolves on timescales comparable to ionic relaxation under bias sweeps, leading to strongly scan-rate-dependent heating. At intermediate scan rates, this thermo-electro-ionic coupling produces non-monotonic temperature evolution with dual-peak profiles during forward sweeps and pronounced T–V hysteresis, coinciding with S-shaped J–V distortions that are shifted to lower scan rates relative to isothermal predictions. The scan-rate-dependent heating can be resolved into interfacial- and bulk-dominated regimes: interfacial heating governs the temperature evolution at low-to-intermediate scan rates, bulk heating controls the profile at intermediate-to-high rates, while rapid sweeps leave insufficient time for heat to accumulate, ultimately driving the response toward an ion-frozen, quasi-isothermal limit. These distinct thermal regimes reshape carrier extraction asymmetry, internal field screening, and mobile-ion redistribution, thereby aggravating or mitigating hysteresis relative to isothermal electronic–ionic transport predictions. Neglecting thermo-electro-ionic effects underestimates transient temperature rises by >10 K and misidentify the scan-rate window associated with S-shaped J–V distortions. By integrating these multiphysics effects, the framework provides a diagnostic tool for next-generation PSC characterization and identifies design strategies such as interface engineering, nanostructural thermal management, and scan-protocol optimization to enhance device performance and stability under real-world operating conditions.

---

Broader context

Perovskite solar cells (PSCs) have achieved record efficiencies and are considered one of the most promising candidates for low-cost, next-generation photovoltaics. However, their real-world deployment is hindered by instability and inconsistent performance under fluctuating outdoor conditions. A central challenge is hysteresis: device output depends strongly on measurement history and scan rate, making it difficult to assess efficiency and stability with confidence. Conventional models typically neglect the role of heat generation and thermal feedback in hysteresis, even though sunlight absorption inevitably raises device temperature and alters ionic and electronic processes. In this work, we present a multiphysics mathematical framework that explicitly couples optical absorption, charge transport, ion migration, and heat generation across all layers. This approach reveals how heating at interfaces and within the bulk reshapes internal fields and carrier dynamics, producing asymmetries and transient losses not captured by conventional isothermal analysis. By introducing temperature–voltage and heat–voltage diagnostics, we provide practical tools to identify stability bottlenecks and clarify the regimes where interfacial engineering or bulk thermal management are essential. Together, these insights advance the predictive design of thermally robust PSCs and help bridge the gap between laboratory testing and reliable operation in real environments.

Introduction

Perovskite solar cells (PSCs) have rapidly emerged as high-efficiency thin-film photovoltaics, but ensuring their operational stability remains a critical hurdle. In particular, thermal instabilities in PSCs—exacerbated by device self-heating and spatial temperature gradients—are increasingly recognized as key degradation pathways.¹ During operation, PSCs generate heat through multiple loss mechanisms,² reaching temperatures that accelerate material decomposition, nanostructural evolution, and interfacial degradation. Such self-heating, driven by intrinsic loss channels and modulated by environmental exposure, can disrupt charge transport,³ accelerate ion migration,^4,5 and aggravate J–V hysteresis,⁶ thereby impairing photovoltaic performance.^7–9 Yet, compared with well-studied stressors such as moisture or UV light, the origin and dynamics of in situ thermal effects remain insufficiently resolved. Importantly, it is not only the magnitude of heating and temperature-gradient formation that matters but also their temporal evolution across the stack. These spatiotemporal dynamics reshape J–V characteristics, charge transport and extraction, field screening, ionic redistribution, and hysteresis, which become particularly important under real-world operating conditions.

Under real-world operation, PSCs experience dynamic environmental conditions that impose fluctuating and often asymmetric optical and thermal boundary conditions, thereby reshaping the spatiotemporal evolution of self-heating and the coupled electronic–ionic processes. For example, a device exposed to solar illumination on the front surface, while simultaneously subjected to convective and radiative heat fluxes on both sides, develops transient and spatially non-uniform heating across the stack. These boundary-driven optical-thermal loads act in concert with internal self-heating mechanisms, including hot-carrier thermalization, Peltier effects at interfaces and contacts, Joule heating in bulk layers, and radiative/non-radiative recombination losses. These heating mechanisms can evolve on timescales comparable to both electronic and ionic dynamics, depending on the bias sweep rate. They create temperature gradients across the stack that modulate charge transport, electric field distribution, and ion migration through the thermoelectric (Seebeck) effect, where temperature gradients induce a Seebeck field that superimposes on the electrostatic field and enables carriers in hotter regions to drift or diffuse toward cooler regions. This thermally driven redistribution alters the internal electric field landscape and can either enhance or hinder carrier extraction depending on the gradient direction and the locally dominant carrier type. Furthermore, because the predominant recombination pathway (bulk versus interfacial) critically determines where self-heating localizes, the intrinsic asymmetry of mixed electronic–ionic processes further modulates self-heating mechanisms and temperature-gradient formation and reinforces feedback with electronic–ionic transport. The overall effects establish thermo-electro-ionic feedback loops that span disparate timescales during bias sweeps—from ultrafast electronic responses to slower ionic redistribution and the finite thermal relaxation time of the stack—a phenomenon that cannot be captured under the conventional assumption of uniform device temperature. Ultimately, the interplay between device optoelectronic properties, mixed electronic–ionic transport, recombination pathways, self-heating mechanisms, and the thermal properties of individual layers—together with asymmetric boundary conditions—governs both PSC transient performance and J–V hysteresis.

Despite their practical importance, the fundamental mechanistic understanding of thermo-electro-ionic dynamics in PSCs remains in its infancy. Experimental evidence on thermal gradients across PSCs under operating conditions is still scarce. For instance, Wang et al.¹ showed that even a modest temperature differential across a device can trigger significant chemical instability in perovskites via increased ion migration, even when absolute temperature is well below the decomposition threshold. This finding highlights that spatial temperature gradients, independent of absolute temperature, can accelerate degradation, yet systematic investigations remain limited. From the modeling perspective, prior studies^2,4,10–13 have employed classical multiphysics toolboxes (e.g., finite-element approaches in COMSOL Multiphysics) to probe nanoscopic thermodynamics in PSCs, but these frameworks typically neglect critical features governing thermo-electro-ionic interactions—including non-linear field screening, mobile-ion redistribution, Debye layers at cross-layer interfaces, and spatially resolved heat generation. These limitations arise because the electric potential profile in PSCs differs fundamentally from that in classical p–n or p–i–n junctions due to the presence of mobile ions,¹⁴ producing nearly linear potential variations across bulk layers and sharp changes across narrow Debye regions. Standard semiconductor drift-diffusion solvers fail to reproduce this behavior, limiting their applicability to perovskite devices. Therefore, there is a need for self-consistent multiphysics models that explicitly couple thermal, electrical, and ionic domains with layer- and interface-resolved, time-dependent heat-charge-ion generation/transport, without relying on idealized isothermal assumptions. Although PSC-specific solvers such as IonMonger,^10,11 Driftfusion,¹² and SIMsalabim¹³ incorporate coupled electrical-ionic physics, they inherently assume fixed device temperature and thus cannot capture the transient, multiphysics essence of thermo-electro-ionic dynamics under real-world operating conditions.

In the present study, we develop a coupled optical-electrical-thermal (OET) multiphysics framework for planar PSCs—self-consistently solving the drift-diffusion-Poisson-energy equations across all layers—to understand thermo-electro-ionic interactions that reshape transient performance and hysteresis. The OET framework is structured into three interconnected models—optical, electrical, and thermal—that are solved self-consistently. (i) The optical model consists of a spectrally- and polarization-resolved partial-coherence formulation that simulates light absorption and photogeneration across the device stack. (ii) The electrical model extends a standard electronic–-ionic transport model, coupled with Poisson's equation, to a non-isothermal formulation by incorporating thermal-diffusion terms into the transport equations and including temperature-dependent material properties such as bandgap and ionic diffusivity. (iii) The thermal model treats layer-resolved self-heating sources, with separate temperature equations for each layer, distinguishing interfacial contributions (Peltier heating and interfacial recombination heat) from bulk processes (hot-carrier thermalization, Joule heating, bulk radiative/non-radiative recombination, and parasitic absorption in non-photoactive layers). This fully coupled approach enables genuinely thermo-electro-ionic transport simulations, thereby reshaping J–V characteristics, charge transport/extraction asymmetry, internal field screening, and ionic redistribution in ways absent from isothermal models. The self-consistent coupling of the energy balance with optical and electrical domains ensures rigorous energy conservation and allows spatiotemporal tracking of individual heat-generation pathways alongside electrical output during voltage sweeps.

The developed OET framework is applied to fulfill two objectives: (i) to quantify how transient self-heating reshapes mixed electronic–ionic response across scan rates, and (ii) to introduce T–V (temperature-voltage) and P–V (heat-voltage) diagnostic tools. By treating temperature as a self-consistent dynamic state variable rather than a prescribed parameter, we unveil how scan rate (reflecting electronic and ionic timescales) and thermal relaxation time jointly govern key photovoltaic behaviors, including photovoltaic performance metrics, ion migration and hysteresis, internal field screening, and charge transport and extraction between isothermal and non-isothermal treatments. We then go beyond standard J–V analysis by resolving thermal hysteresis and the operating-state cell temperature evolution during bias sweeps and demonstrate that T–V and P–V diagnostics provide quantitative, bias-resolved metrics to interpret the observed phenomena. Importantly, the strong scan-rate dependence of device temperature predicted by our model fundamentally alters interpretations of thermal degradation. For example, the downshift of the scan-rate window for S-shaped J–V curves relative to isothermal predictions indicates that thermally driven instabilities can distort transient device response in ways invisible to isothermal models. As a result, performance or stability conclusions drawn from scan-rate-independent temperature measurements risk overlooking coupled thermo-electro-ionic effects. In addition, our model further elucidates the origin of these behaviors by resolving scan-rate-dependent interfacial and bulk heating contributions within the perovskite absorber. We identify how transient performance and hysteresis are governed by the interplay between non-isothermal charge transport and extraction, thermally induced field screening, thermally accelerated iodide vacancy redistribution, and heat-generation pathways. By explicitly including thermo-electro-ionic effects, the framework provides both a diagnostic tool for next-generation PSC characterization and practical design guidelines for optical, electrical, and thermal engineering strategies aimed at enhancing PSC performance and stability under real-world conditions.

Materials and methods

The simulated device is an n–i–p planar perovskite solar cell with the Glass/SnO₂/SiO₂/FTO/TiO₂/MAPbI₃/Spiro-OMeTAD/Au architecture (Fig. S1a), with the reduced one-dimensional simulation domain used for coupled electrical-thermal modeling presented in Fig. S1b. SnO₂ and SiO₂ are introduced as front-side optical interlayers to suppress reflection and parasitic absorption, as suggested by experimental optimization,¹⁵ while the FTO/TiO₂ interface texture is represented using modified ellipsometry-derived effective optical constants for FTO. In this configuration, glass serves as the substrate, SnO₂ and SiO₂ as optical spacers, FTO as the transparent cathode, compact TiO₂ as electron transport layer (ETL), MAPbI₃ as the absorber, Spiro-OMeTAD as the hole transport layer (HTL), and Au as the anode. Nominal thicknesses and optical constants are taken from Ball et al.¹⁵ and are detailed in Section S1.4 in the SI.

A fully coupled multiphysics mathematical framework is developed that integrates spectrally-resolved optical absorption, non-isothermal electronic–ionic transport, and a self-consistent, layer-resolved energy balance with explicit bulk and interfacial heat-generation pathways. The framework consists of three models for the optical, electrical, and thermal domains, and resolves the spatiotemporal dependence of electrostatic potential, charge-carrier populations, iodide vacancy density, and lattice temperature across the stack. The optical model employs the transfer-matrix method with partial coherence to account for thin-film interference, incoherent propagation in glass, and multiple internal reflections.¹⁵ The electrical model extends the isothermal electronic–ionic transport formulation of IonMonger^10,11 to a non-isothermal framework by incorporating thermal diffusion terms into the drift-diffusion-Poisson equations and explicitly including temperature-dependent bandgap and ionic diffusivity, with temperature updated self-consistently through direct coupling to the energy balance equation. Carrier statistics in the transport layers follow Fermi–-Dirac distributions to ensure accuracy for organic and heavily doped transport materials. The thermal model solves transient, layer-resolved energy balance equations, incorporating all major heating mechanisms: hot-carrier thermalization, Joule heating, Peltier exchange, radiative and non-radiative recombination, and parasitic absorption in non-photoactive layers. The overall energy flow captured by our multiphysics framework is summarized schematically in Fig. 1, which shows how incident solar power is partitioned into optical losses, useful electrical output, and bulk and interfacial heat-generation pathways, all coupled to charge generation and transport, ion migration, and polarization, built-in, and thermally induced fields. By explicitly integrating these energy flows with the governing electrostatic potential, coupled electronic–ionic transport, and temperature field, the framework enables term-by-term spatiotemporal energy tracking during transient and directional voltage sweeps. This explicit mapping of loss channels and field interactions provides a physically rigorous basis for interpreting how self-heating modulates electronic–ionic transport and thereby reshapes device transient performance and hysteresis. Full details of model assumptions and model validation and verification are provided in Section S1 of the SI, heat-generation pathways and the complete list of symbols and units (Table S1) in Section S2, and a term-by-term proof of energy conservation together with the physical input parameters and material properties (Table S2) in Section S3.

Fig. 1 Schematic of the optical-electrical-thermal (OET) energy flows and field profiles in a perovskite solar cell, illustrating the partitioning of incident optical power into reflection/transmission losses, parasitic absorption, useful electrical output, and resolved bulk/interfacial heat-generation pathways, together with the self-consistent temperature field and electrostatic band diagram including polarization, built-in, and thermally-induced (Seebeck) fields. P^Eg: below-bandgap optical loss, P^el: electrical output, Pⁱⁿ: input solar irradiance, P^J: Joule heating, P^par: parasitic optical absorption, P^Pl: Peltier heating/cooling, P^rec: recombination heat, P^ref: reflectance power, P^th: thermalization loss, P^tra: transmittance power, E_pol: polarization field, E_bi: built-in electric field, E_S: Seebeck field, E_c(v): conduction(valence) band edge, E_fn(fp): quasi-Fermi levels, T: temperature field. The illustrated band-bending and temperature profiles correspond to a representative operating steady state near the maximum power point (MPP) and are shown to clarify principles of self-heating.

Numerically, the coupled drift-diffusion-Poisson-energy equations are formulated as a differential-algebraic system and solved in MATLAB with the ode15s solver for temporal evolution. The solver employs the finite-element/method-of-lines scheme of Courtier et al.¹⁶ with hyperbolic grid refinement and adaptive node clustering near Debye layers and interfaces, ensuring uniform error control without overrefinement. For transient J–V sweeps, the device is initialized at the built-in potential and preconditioned at 1 V for 5 s, after which a reverse scan from 1 V to 0 V is applied, followed by a forward scan from 0 V back to 1 V. The scan rate is systematically varied across several orders of magnitude to probe the crossover between electronic, ionic, and thermal relaxation timescales. This computational protocol enables consistent resolution of scan-rate-dependent J–V, T–V, and P–V curves while ensuring strict instantaneous energy conservation and numerical closure across optical, electrical, and thermal domains within the solver tolerances.

The local absorbed power density per unit volume and wavelength in each layer, A_j(x, λ), is computed from the Poynting vector formalism, expressed as the product of the incident spectral irradiance, the absorption coefficient, and the normalized electric-field intensity:

A_j(x, λ) = Pⁱⁿ(λ)α_j(λ)η_j(λ)(|E^x_j|² + |E^z_j|² + |E^y_j|²) (1)

where Pⁱⁿ(λ) is the spectral solar irradiance, η_j(λ) is the real part of the refractive index of layer j and represents the phase velocity of light, and α_j(λ) = 4πk_j/λ is the absorption coefficient of layer j. The Euclidean electric field components (E^x_j, E^y_j, and E^z_j) in each layer are computed for both forward- and backward-propagating waves as a function of the ambient incident angle θ₀ and the layer-specific wavevector ξ_j = 2πñ_j cos θ_j/λ, where ñ_j is the complex refractive index. This formalism ensures a consistent treatment of oblique incidence and angular dispersion effects in multilayer stacks. Parasitic absorption in non-photoactive layers is obtained by integrating A_j(x, λ) over wavelength and is applied as a volumetric heat source in the corresponding energy-balance equations.

The depth-resolved photogeneration rate in the photoactive perovskite layer is obtained as:¹¹

where A_A is the local power density absorbed in the absorber layer, λ_g = hc/E_g, E_g is the bandgap, and h and c are Planck's constant and speed of light, respectively. This photogeneration profile serves as the explicit source term in the perovskite charge-transport equations. In the following, the governing non-isothermal heat-charge-ion transport equations are presented for each layer of the planar perovskite solar cell. For clarity, all model variables and their definitions are summarized in Table S1.

Electron transport layer (−l_E < x < 0)

Accounting for only electrons as the majority charge carriers, the governing equations for electron density (n_E), electrostatic potential (ϕ_E), and temperature field (T_E) are:Here, k_B is the Boltzmann constant, q is the elementary charge, jⁿ is charge current density, μ is charge mobility, ε is permittivity, d is the effective doping density, ρ is density, c_p is specific heat capacity, and k_c is thermal conductivity. The inverse statistical function follows a Fermi–Dirac distribution, as detailed by Clarke et al.¹⁶ The Joule heating P^J_E and parasitic optical absorption P^par_E terms are defined in eqn (S1) and (S2) of the SI, respectively.

Perovskite absorber layer (0 < x < l)

The perovskite absorber is modeled with electrons, holes, and mobile-ion vacancies as active species. The governing drift-diffusion-Poisson-energy equations are given by:where N̂₀ is the equilibrium ion vacancy density, P is the halide ion vacancy density, and F^P is its flux. The heat sources in eqn (10) are hot-carrier thermalization heating P^th, Joule heating P^J, and bulk recombination heating P^rec (eqn (S3)–(S5)). In eqn (8), D_I is the ion vacancy diffusivity, given by:¹⁷where D₀ is the pre-exponential factor and E_ac is the ion activation energy.

Hole transport layer (l < x < l_H)

In the HTL, where holes dominate transport, the governing equations for hole density (p_H), electrostatic potential (ϕ_H), and temperature field (T_H) are:The Joule heating term P^J_H and parasitic optical absorption heating P^par_H are given in eqn (S6) and (S7).

Electrodes (−l_ca < x < −l_E & l_H < x < l_an)

The FTO cathode and Au anode are treated only in the thermal domain, with heating arising exclusively from parasitic optical absorption:The parasitic optical absorption heating in electrodes (P^par_ca and P^par_an) are given by eqn (S8) and (S9).

Results

The results are organized to first quantify how finite thermal inertia (finite time for heat accumulation or dissipation) and directional bias sweeps generate distinct scan-rate-dependent cell temperature trajectories, and to explain these trends by decomposing self-heating into interfacial- and bulk-dominated regimes. We then use the resulting thermal response to interpret the associated T–V and P–V characteristics, and investigate the effects of these signatures on mixed electronic–ionic device dynamics by analyzing the scan-rate dependence of photovoltaic figures of merit, ion migration and hysteresis, internal potential redistribution, and carrier populations, relative to an isothermal reference held at the steady-state temperature (310 K). Overall, the results show that these responses are governed by thermo-electro-ionic coupling—wherein bias-dependent self-heating and Seebeck contributions reshape ionic polarization and the internal potential landscape—rather than by ionic screening or fixed-temperature effects in isolation.

Photovoltaic performance characteristics

Fig. 2a illustrates the simulated cell temperature as a function of scan rate during forward and reverse sweeps. Because the stack is thin, through-thickness temperature variations are negligible, and we therefore report the temperature at the perovskite mid-plane. At very low scan rates (quasi-steady state), the cell has ample time to equilibrate at each bias, thereby reaching the quasi-steady-state temperature. At very high scan rates (ion-frozen regime), the cell cools toward ambient because the cumulative heating deposited per bias interval is greatly reduced and thermal relaxation outpaces heat accumulation. Between these limits, the response differs strongly between forward and reverse scans: during the forward sweep the temperature is non-monotonic, with two local maxima separated by a pronounced mid-rate minimum, whereas the reverse-scan temperature decreases quasi-monotonically. These forward-reverse differences reveal strong scan-rate-dependent thermal hysteresis across a broad mid-rate window, with the forward sweep remaining hotter over much of this range.

Fig. 2 Forward/reverse scan-rate dependence of (a) cell temperature and perovskite heating partitioned into (b) interfacial and (c) bulk components. The dual forward-scan temperature peaks reflect a crossover from interfacial-dominated heating at low-to-intermediate rates to bulk-dominated heating at intermediate-to-high rates. At rapid scan rates, device cools toward ambient. Acronyms “Interf” and “PVK” denote interfacial and perovskite, respectively.

Fig. 2b and c decompose perovskite heating into interfacial (Peltier and interfacial recombination) and bulk (thermalization, bulk radiative/non-radiative recombination, and Joule) components, providing insight into the thermal response observed in Fig. 2a. In the forward sweep, the interfacial term (Fig. 2b) displays a pronounced peak (with a maximum hysteretic heating reaching ∼40 W m⁻²) over a narrower low-to-intermediate scan-rate window (centered around 0.01–10 V s⁻¹) than the bulk term, whereas the reverse sweep remains comparatively flat. This behavior is consistent with scan-direction-dependent ionic redistribution and interfacial charge-carrier accumulation, which increase Peltier and interfacial recombination heating during low-to-intermediate forward scans. Bulk heating (Fig. 2c) also shows hysteresis over a broader rate window but with a smaller maximum hysteretic heating (∼25 W m⁻²) than the interfacial term, and it varies more slowly with scan rate.

The dual-peak temperature profile (Fig. 2a) follows from the competition between interfacial and bulk heating sources under finite thermal inertia. The first peak aligns with the low-rate regime where interfacial heating dominates. The second peak occurs at intermediate-to-high rates where bulk heating variation within the perovskite becomes predominant. The intervening minimum corresponds to the crossover point, where interfacial heating has already subsided while bulk heating has not yet fully risen, resulting in reduced total heat generation. Beyond the second peak, although instantaneous heating increases modestly, the dwell time per bias step becomes shorter than the thermal relaxation time, and boundary heat removal outpaces further heat accumulation; as a result, the cell temperature declines. Overall, the scan-rate-dependent thermal hysteresis stems from the interplay between interfacial and bulk heat sources and the cell's finite thermal inertia: interfacial processes dominate hysteretic heating at low-to-intermediate rates, bulk heating at intermediate-to-high rates, and rapid sweeps leave insufficient time for heat accumulation.

Fig. 3 compares the scan-rate dependence of J_sc, V_oc, FF, and PCE predicted by a thermally coupled OET model (self-consistent T) and an isothermal reference. The results show that the photovoltaic figures-of-merit change profoundly when the lattice temperature is resolved self-consistently and coupled to mixed electronic–ionic transport. At the two extremes—quasi-static and ion-frozen—J_sc and FF are nearly identical in both treatments, indicating that neither transient self-heating nor ionic redistribution limits photocurrent extraction in these two regimes. Small V_oc offsets persist even at these limits because V_oc is intrinsically temperature-sensitive through the bandgap and recombination statistics, and PCE inherits these differences. Across intermediate scan rates, the differences between treatments in V_oc, FF, and PCE become decisive, while J_sc is only weakly affected: relative to the isothermal reference, the thermally coupled model exhibits a deeper forward-scan V_oc dip and shifts both the V_oc and FF minima to lower scan rates (Fig. 3b and c). Consistent with this, the forward J–V develops an S-shaped distortion under thermal coupling already at ∼0.1 V s⁻¹ (Fig. S10), with a concave kink near mid-bias that sharply reduces FF; the isothermal case shows a comparable distortion only at higher scan rates (Fig. S11). This indicates that thermo-electro-ionic feedback lowers the critical rate for extraction-barrier formation: self-heating accelerates ion kinetics (Arrhenius), increases ion-induced screening of the built-in field, and elevates non-radiative recombination, thereby jointly depressing V_oc and FF at rates where the isothermal treatment appears innocuous. Accordingly, the V_oc (and hence PCE) comparison between treatments is strongly scan-rate dependent: at very slow (≲0.03 V s⁻¹) and very fast (≳0.3 V s⁻¹) sweeps, the isothermal reference tends to underestimate V_oc because a single setpoint with fixed temperature cannot reproduce the actual quasi-Fermi level splitting trajectory as the device self-heats or self-cools during the sweep; within the intermediate window (≈0.03–0.3 V s⁻¹) the trend reverses in forward scan—the isothermal model over-predicts V_oc while the thermally coupled calculation incurs a stronger temperature-voltage dependence (Fig. 3b). In essence, self-heating combined with mobile-ion redistribution produces a non-monotonic voltage response in PSCs that propagates to V_oc, FF, and PCE. This behavior can be viewed as a quasi-isothermal response at an effective temperature T_eff = T_ss + ΔT (Fig. S3), where ΔT depends on the device optical/thermal coupling to the environment. Consequently, fixed-temperature optoelectronic analyses misplace and misestimate the V_oc/FF/PCE losses over practical scan rates, whereas a fully coupled OET model is required for faithful scan-rate-dependent prediction.

Fig. 3 Scan-rate dependence of (a) J_sc, (b) V_oc, (c) FF, and (d) PCE for forward (solid) and reverse (dashed) sweeps under non-isothermal (black) and isothermal (orange) models. The isothermal temperature is fixed to the non-isothermal steady state at 0.001 V s⁻¹. Thermo-electro-ionic coupling shifts the hysteresis window in perovskite solar cells.

The companion T–V curves, presented in Fig. 4, provide deeper insight into how the cell temperature evolves with bias at different scan rates, offering a physical basis for the differences observed in the non-isothermal J–V curves. Fig. 5 complements this analysis by decomposing perovskite heating into interfacial (Peltier and interfacial recombination) and bulk (thermalization, Joule heating, and bulk recombination) contributions, thereby attributing the T–V response to explicit, bias-resolved heat dissipation terms. Consistent with Fig. 2a, and 4 demonstrates that the T–V trajectory reflects the interplay between bias-dependent heat generation and thermal inertia: the former dominates at low-to-intermediate scan rates, whereas the latter controls higher scan rates. At very low scan rates (i.e. ∼0.001 V s⁻¹), the T–V traces converge to the steady-state profile and decrease monotonically from J_sc to V_oc because Joule heating collapses as current decreases, while interfacial heating remains localized and is further suppressed by two-sided cooling. The quasi-steady-state T–V curve exhibits two distinct bias-dependent regimes: from J_sc to MPP and from MPP to V_oc. In the second regime, beyond MPP, the cell temperature decreases more sharply because Joule heating diminishes rapidly as the current drops steeply with increasing voltage. As the scan rate increases, however, thermal inertia (lag) becomes apparent, and the cell temperature deviates markedly from the quasi-steady-state trajectory. At intermediate scan rates (∼0.1–10 V s⁻¹), ionic polarization repartitions heat generation across interfacial and bulk pathways. The most pronounced re-partitioning occurs at 0.1 V s⁻¹ (Fig. 5c), coinciding with the S-shaped J–V distortion (Fig. S10c). At this scan rate, the reverse-scan T–V profile remains close to the quasi-steady-state trajectory, whereas the forward-scan temperature drops significantly with increasing bias as ion-induced field screening suppresses current and electrical power extraction, redistributing dissipation toward resistive loss channels. Consequently, the forward-scan MPP shifts to lower voltages, and its power contribution decreases from ∼15.3% to ∼8%, while Joule heating increases nearly fourfold. At scan rates >0.1 V s⁻¹, the device starts cooler in both scan directions but can exceed the quasi-steady-state temperature during the forward scan near V_oc, because thermal inertia dominates and heat generated earlier in the sweep accumulates at high bias. The subsequent cooling at the end of the sweep is then governed by boundary heat dissipation, consistent with the post-overshoot decline in Fig. 2a. At much faster scan rates (10–100 V s⁻¹), the temperature response collapses toward a quasi-isothermal trajectory and becomes weakly bias-dependent. In this limit, J–V hysteresis is expected to arise primarily from ionic polarization rather than thermal memory, highlighting the tight coupling between electronic and ionic timescales that governs the dynamic response of perovskite solar cells. Under this quasi-isothermal condition, both Joule heating (∼12.7 W m⁻²) and interfacial recombination heating (∼9.8 W m⁻² forward and ∼9.5 W m⁻² reverse) are reduced relative to the quasi-steady state, while the electrical power density increases by ∼5.8 W m⁻² (forward) and ∼5.2 W m⁻² (reverse). This corroborates the higher ion-free PCE in Fig. 3d under both treatments and explains why it typically exceeds the quasi-steady-state value—as discussed in Section S4.1 of the SI and previously observed experimentally.^18,19 In summary, direction- and bias-dependent heat-generation pathways—together with finite thermal inertia and ionic polarization—govern the forward-reverse divergence in key photovoltaic parameters (J_sc, V_oc, FF, and PCE) of perovskite solar cells. These results highlight that only a thermally coupled treatment can faithfully capture the intertwined thermo-electro-ionic dynamics.

Fig. 4 Forward and reverse T–V curves illustrating the transition from quasi-steady-state cooling at (a) low scan rates to thermally lagged and direction-dependent heating at (b–e) intermediate scan rate, and finally to (f) a quasi-isothermal limit at high scan rates governed by finite thermal inertia and boundary cooling.

Fig. 5 (a–f) Voltage-dependent power densities associated with electrical output and different heat-generation pathways (thermalization, Joule heating, Peltier exchange, and bulk/interfacial recombination) across scan rates (SRs). Insets show the MPP energy budget (in W m⁻²) for the forward sweep (reverse in brackets): optical losses—reflection/transmission (P^opt), parasitic absorption (P^par), and below-bandgap losses (P^Eg) evaluated under the AM 1.5G spectrum—together with electrical output and thermal losses. Joule heating and interfacial recombination heat re-partition the losses at intermediate rates, biasing the forward MPP toward lower voltages and increasing forward-reverse disparity. At high scan rates, short dwell times suppress cumulative heating despite large instantaneous interfacial heat-generation peaks.

By tracking the device thermal response during a bias sweep, the T–V and P–V curves provide a complementary diagnostic tool alongside established optoelectronic techniques such as photoluminescence (PL) and electroluminescence (EL) spectroscopy. Whereas PL/EL primarily probe radiative emission and quasi-Fermi level splitting and thus diagnose recombination- and voltage-loss pathways, the T–V curve directly reveals thermal memory and lag (thermal hysteresis)—distinguishing thermally quasi-steady, rate-limited, and quasi-isothermal regimes. The P–V curve, in turn, provides a voltage-resolved dissipation budget by decomposing self-heating into its constituent mechanisms, thereby identifying dominant irreversible loss channels and the bias- and direction-dependent localization of heat dissipation under ionic polarization—thereby identifying operating conditions prone to hotspot formation. This information is not directly accessible from PL/EL alone, because optical emission metrics do not uniquely constrain non-radiative dissipation or resistive losses, nor do they encode the thermal boundary conditions and heat-flow time constants that govern self-heating, thermal gradients, and thermal hysteresis under dynamic operation. Consequently, the combined T–V/P–V diagnostics provide a practical route to distinguish purely electrical hysteresis from thermo-electrical hysteresis and to identify when explicit thermal coupling is required for faithful interpretation of the J–V response.

Ion migration and hysteresis

Fig. 6 shows the spatiotemporal evolution of the iodide-vacancy distribution across the perovskite Debye layers during forward and reverse sweeps at a scan rate of 0.1 V s⁻¹. The corresponding scan-rate dependence is presented in Fig. S12, while Fig. S5 shows the fixed-temperature isothermal trends. In both scan directions, the bulk vacancy density remains close to its mean value, whereas pronounced, bias-dependent deviations are confined to the narrow Debye layers. The largest divergence between non-isothermal and isothermal vacancy profiles occurs at intermediate scan rates (∼0.1–1 V s⁻¹), where thermal, electronic, and ionic timescales are comparable (see Fig. 6 and Fig. S12). At 0.1 V s⁻¹ and near J_sc, both sweep directions exhibit weaker ionic polarization under non-isothermal conditions, manifested as shallower depletion at the ETL/perovskite interface and reduced accumulation at the perovskite/HTL interface, relative to the isothermal case. During the forward sweep (Fig. 6a), as the scan progresses toward the MPP, the vacancy profiles under both treatments converge, indicating that the competing influences of the Seebeck field and ionic polarization become comparable. Toward V_oc, however, the non-isothermal case diverges again, showing stronger ionic polarization relative to the isothermal reference. In contrast, during the reverse sweep (Fig. 6b) the non-isothermal profile remains less polarized than the isothermal reference throughout the bias sweep. This bias-dependent crossover during the forward sweep demonstrates that thermal coupling does not monotonically increase or suppress ion-vacancy migration; rather, thermo-electro-ionic dynamics redistribute the net field and thereby modulate ionic polarization.

Fig. 6 Spatiotemporal evolution of iodide vacancy density across the perovskite Debye layers during (a) forward and (b) reverse scans at a scan rate of 0.1 V s⁻¹ under isothermal (dashed) and non-isothermal (solid) treatments. Thermal coupling modifies ionic polarization depending on scan rate and sweep direction: near J_sc, polarization is suppressed in both sweep directions; during the forward sweep the profiles nearly overlap around the MPP and then diverge toward V_oc, where the non-isothermal case exhibits stronger polarization; during the reverse sweep, the non-isothermal profile remains less polarized than the isothermal reference across the bias sweep. Color gradient (gray → black) denotes the bias sweep from J_sc to V_oc (forward) and from V_oc to J_sc (reverse).

At very slow scan rates (Fig. S12a and d), ion vacancies have sufficient time to approach quasi-equilibrium Debye-layer profiles, leading to convergence between the two treatments. At very fast scans, ions become effectively frozen on the sweep timescale and the Debye layers contract, causing the non-isothermal and isothermal results to converge again toward the ion-frozen, quasi-isothermal limit. Overall, Fig. 6 and Fig. S12 show that quantitative prediction of ionic polarization in PSCs requires treating thermal, electrical, and ionic dynamics jointly, particularly in the intermediate-rate regime where their feedback is strongest.

Fig. 7 quantifies the scan-rate dependence of the hysteresis factor (HF) under isothermal and non-isothermal treatments. Both treatments exhibit the characteristic bell-shaped HF curve, with negligible hysteresis at very low and high scan rates, and a distinct maximum at intermediate rates. While in the isothermal case, the HF peak shifts to higher scan rates with increasing temperature due to enhanced ionic mobility (Fig. S6), under non-isothermal conditions the HF maximum is systematically shifted to lower scan rates (compared to the isothermal reference), while maintaining a comparable peak amplitude. This shift indicates that the scan rate at which forward-reverse mismatch is maximized is governed by thermo-electro-ionic coupling rather than solely by ion migration at a prescribed uniform cell temperature. Overall, Fig. 7 shows that accurate modeling of hysteresis in perovskite devices under realistic dynamic operation requires an explicit non-isothermal treatment, as fixed-temperature electronic–ionic descriptions cannot reproduce the thermally induced scan-rate shift in hysteresis.

Fig. 7 Hysteresis factor as a function of scan rate under isothermal and non-isothermal treatments. The non-isothermal treatment shifts the hysteresis maximum to lower scan rates, demonstrating that the scan-rate dependence of hysteresis is governed by thermo-electro-ionic coupling rather than solely by ion migration at a prescribed uniform cell temperature.

Thermally induced field-screening

Fig. 8 demonstrates the evolution of the internal electric potential (band bending) across the device during forward and reverse scans at 0.1 V s⁻¹, with scan-rate dependence summarized in Fig. S14, and Fig. S7 showing the corresponding fixed-temperature isothermal trends. Compared with the isothermal reference, thermo-electro-ionic dynamics reshape the potential landscape through the coupled superposition of (i) ionic polarization, which redistributes the electrostatic potential drops between the perovskite bulk and the Debye layers, and (ii) Seebeck (thermoelectric) field arising from local temperature gradients (E_S ∝ −∇T). The Seebeck field may either oppose or reinforce the electrostatic field depending on the sign of ∇T and carrier polarity, explaining the spatiotemporal divergence between isothermal and non-isothermal potential profiles.

Fig. 8 Spatiotemporal evolution of the electric potential across the device during (a) forward and (b) reverse scans at a scan rate of 0.1 V s⁻¹ under non-isothermal (solid) and isothermal (dashed) treatments. Local temperature gradients generate a Seebeck field that, together with ionic polarization, dynamically reshapes the potential landscape and re-partitions the voltage drop between the perovskite bulk and the interfacial Debye layers. The largest treatment-dependent deviations occur near the sweep endpoints—most prominently near J_sc and V_oc in the forward scan and near J_sc in the reverse scan. Color gradient (gray → black) denotes the bias sweep from J_sc to V_oc (forward) and from V_oc to J_sc (reverse).

During the forward sweep (Fig. 8a), at J_sc, Joule heating dominates and non-isothermal ionic polarization is weaker than that in the isothermal reference. In this regime, the ETL → HTL Seebeck field is also co-aligned with the built-in field, reinforcing the net internal field relative to the isothermal reference. Consequently, less of the non-isothermal potential drop is confined to the Debye layers and a larger fraction occurs across the absorber bulk. As the bias progresses toward MPP, the isothermal and non-isothermal profiles partially converge, indicating a reduced perturbation of the potential by thermo-electro-ionic dynamics at this operating point. Approaching V_oc, the non-isothermal case retains stronger ionic polarization than the isothermal reference, while the Seebeck field partially offsets this increased ionic screening, yielding a higher net field relative to the isothermal case.

When the scan is reversed (Fig. 8b), the device begins near V_oc in a pre-polarized state, so the field is intrinsically small and the two profiles nearly coincide. As the bias is reduced toward MPP, ionic redistribution progressively screens the electric field for both treatments; however, the isothermal case exhibits a stronger ionic screening in this region, whereas in the non-isothermal case the Seebeck contribution partially compensates for this screening, producing a steeper net bulk slope. Toward J_sc, this compensation becomes more pronounced, yielding the largest discrepancy between the two treatments.

At slow and fast scans (Fig. S14), the potential distribution differs from those observed at 0.1 V s⁻¹. At slow scan rates, ions have sufficient time to equilibrate and heat can accumulate, so the applied bias is preferentially redistributed into the Debye-layer regions and the bulk field is reduced for both isothermal and non-isothermal treatments. At fast scan rates, ionic motion and thermal buildup are limited, leading to more linear potential drops across the perovskite bulk and convergence of both treatments toward a quasi-isothermal potential profile. Overall, Fig. 8 and Fig. S14 show, compared to the isothermal reference, the non-isothermal net internal field is set by superposition of ionic polarization and the Seebeck field during bias sweep, with the largest deviations between the two treatments occurring at intermediate rates where neither ions nor self-heating fully equilibrate with charge carriers.

Thermo-electro-ionic dynamics of charge transport

Fig. 9 shows the spatiotemporal evolution of electron and hole densities across the device at an intermediate scan rate of 0.1 V s⁻¹, with Fig. S15 presenting the scan-rate dependence and Fig. S8 illustrating the corresponding fixed-temperature isothermal trends. During the forward sweep at 0.1 V s⁻¹ (Fig. 9a), non-isothermal electron and hole densities at and near J_sc are noticeably lower than their isothermal counterparts, consistent with the weaker ionic polarization and reinforced net field discussed for the same operating region in Fig. 6a and 8a. As a result, carriers are extracted more efficiently early in the forward sweep, suppressing charge accumulation in the perovskite absorber despite strong photocurrent and Joule-heating dissipation. As the bias increases toward MPP, the isothermal and non-isothermal profiles partially converge, reflecting minimal differences in both the internal field and ionic polarization near this operating point. Toward V_oc, the divergence increases again, with the non-isothermal treatment exhibiting a more pronounced charge imbalance—characterized by increased electron accumulation toward the HTL side and stronger hole depletion toward the ETL side—as the Seebeck field reinforces the net field relative to the isothermal case (consistent with Fig. 8a).

Fig. 9 Spatiotemporal evolution of electron and hole density profiles at a scan rate of 0.1 V s⁻¹, showing divergence between the isothermal reference and the non-isothermal treatment during (a) forward and (b) reverse sweeps. During the forward sweep, the non-isothermal treatment suppresses carrier accumulation near J_sc relative to the isothermal reference and amplifies interfacial charge imbalance toward V_oc. During the reverse sweep, carriers are more effectively extracted under non-isothermal conditions. Overall, the deviations between the two treatments largely originate from the interfacial regions. Color gradients (gray → black for electrons, light orange → dark orange for holes) denote the bias progression from J_sc to V_oc (forward) and from V_oc to J_sc (reverse).

During the reverse sweep (Fig. 9b), the device begins near V_oc in a pre-polarized state, resulting in nearly identical bulk carrier distributions for both treatments and deviations confined primarily to the interfaces. As the bias approaches MPP and then J_sc, stronger ionic screening in the isothermal case suppresses the bulk electric field, which limits drift-assisted carrier redistribution, whereas in the non-isothermal case the Seebeck contribution partially offsets this screening and accelerates recovery of the bulk field. Consequently, carriers are more effectively swept across the absorber under non-isothermal conditions.

At very slow and very fast scans (Fig. S15), isothermal and non-isothermal carrier profiles largely overlap during forward and reverse sweeps, with only minor differences in majority-carrier populations accumulated at interfaces. These differences become most visible near V_oc, where the non-isothermal treatment yields systematically lower interfacial majority-carrier densities. This observation, consistent with Fig. 9, indicates that thermo-electro-ionic coupling mostly affects interfacial charge-carrier extraction and recombination. Overall, Fig. 9 and Fig. S15 show that thermo-electro-ionic effects on carrier populations peak in the same intermediate scan-rate window where the electric field and ionic polarization landscapes diverge most strongly between the two treatments, whereas outside this window the device approaches either quasi-steady operation (slow scans) or the ion-frozen, and quasi-isothermal limit (fast scans), with comparatively minor differences at interfaces.

Discussion

This section explores how variations in device architecture, thermal boundary conditions, and material properties can mediate the interplay between self-heating dynamics, ion migration, internal field, and charge transport, thereby affecting the conclusions drawn from the previous section.

Impacts of transport layers’ electrical and thermal properties

Transport layers (TLs) not only control carrier selectivity but also strongly modulate the internal electrostatic distribution, recombination pathways, and interfacial heating. Courtier et al.⁵ identified TL permittivity and doping as key factors governing field screening, ionic redistribution, and recombination losses in perovskite solar cells: high-permittivity or heavily doped TLs reshape the internal field by shifting the major potential drops from the transport layers into the perovskite Debye layers. Consequently, replacing an inorganic ETL such as TiO₂ with an organic counterpart like PCBM not only redistributes interfacial band bending and the depth of electronic–ionic accumulation and depletion but also alters the localization (bulk vs. interfaces) and magnitude of self-heating. This can intensify thermo-electro-ionic dynamic asymmetry relative to the isothermal reference and may worsen hysteresis artifacts in J–V behavior. Beyond their electrical role, the thermal properties of transport layers determine whether heat is efficiently dissipated or instead confined near interfaces. Thermally resistive HTLs such as Spiro-OMeTAD are particularly prone to generating localized hot spots, which in turn exacerbates ion migration and accelerates interfacial degradation. Future work should therefore focus on how key electrical (e.g., permittivity, doping, mobility) and thermal (e.g., conductivity, diffusivity) properties of transport layers jointly govern photovoltaic metrics, charge transport and extraction, internal field distribution, and ionic redistribution—processes that ultimately determine whether hysteresis is aggravated or mitigated.

Cooling flux and boundary effects

The external cooling environment is a decisive yet often overlooked factor in shaping thermo-electro-ionic dynamics of PSCs under real-world conditions. In this study, we assume symmetric cooling fluxes on the front and rear sides as a generic boundary condition to demonstrate the feasibility of our model and capture intrinsic thermo-electro-ionic interactions without the added complexity of external asymmetries. However, real devices rarely operate under such idealized environmental conditions. Asymmetric convective or radiative boundary conditions—arising from differences in airflow, encapsulation, device orientation, or surface optical and thermal emissivity—can significantly reshape temperature gradients and heat-generation pathways across the device stack, thereby altering the direction and magnitude of internal field screening, ionic redistribution, charge transport and extraction asymmetry, and recombination. These alterations, in turn, not only shift the magnitude of peak operating temperatures but also redefine where self-heating localizes within the stack—whether at the interfaces or in the bulk—and ultimately change hysteresis and key photovoltaic performance metrics. Taken together, these considerations indicate that conclusions drawn under symmetric cooling flux conditions, as assumed in this work, should not be uncritically generalized to devices operating under asymmetric environments. Future studies should therefore systematically investigate cooling boundary conditions to establish how external fluxes reshape T–V and J–V characteristics, scan-rate-dependent self-heating profiles, internal temperature gradients, ionic redistribution, field screening, and ultimately the severity or mitigation of hysteresis and long-term stability.

Temperature-dependent electronic properties

Although our model accounts for temperature-dependent bandgap and ion diffusivity, several other semiconductor properties—most notably charge-carrier mobilities, recombination kinetics, trap densities, and density of states—are treated as temperature-invariant to preserve generality and demonstrate the core thermo-electro-ionic coupling. This simplification enables us to isolate the role of thermal gradients on field screening, ion migration, and hysteresis without the confounding effects of material-specific dependencies. However, in real perovskites these parameters are strongly temperature-dependent and could reshape device behavior. Experiments^20,21 and prior simulations³ have indicated that higher temperatures typically enhance carrier mobility but simultaneously accelerate non-radiative and trap-assisted recombination, both of which directly affect open-circuit voltage, fill factor, and stability. Local heating can also activate additional defects and shift recombination pathways, further complicating the balance between mobility gains and recombination losses. If such effects were explicitly included, they could either mitigate (through faster carrier extraction) or exacerbate (through stronger recombination and trap activation) the thermo-electro-ionic dynamic effects revealed by our model. Future efforts should therefore incorporate empirically measured temperature dependencies of mobility, recombination coefficients, trap densities, and density of states for specific PSC architectures. Such studies would not only refine the quantitative conclusions drawn here under constant-property assumptions but also extend the applicability of thermo-electro-ionic dynamic analysis to a broader class of perovskite devices.

Dominant recombination pathways

In this study we assume that the ETL/perovskite interfacial recombination is the predominant loss mechanism, as is often the case in well-optimized planar PSCs where bulk recombination in the perovskite absorber is suppressed by defect passivation and long carrier diffusion lengths. This leads to pronounced local heating near the ETL, producing a directional thermal gradient from front to back and asymmetric ionic redistribution. For inverted p–i–n architecture, for example, the device exhibits different recombination pathways compared to the standard structure,^22,23 which can invert or redistribute the heat-generation pathways relative to those observed here, thereby producing distinct thermo-electro-ionic coupling effects. Further, many perovskite solar cells—particularly those with high trap densities or short carrier lifetimes—are bulk-recombination-limited, in which case heat generation is distributed more uniformly across the absorber. In such cases, a more uniform heating profile within the perovskite would likely drive a correspondingly more uniform ionic redistribution, rather than the strong ion accumulation at one interface observed under interfacial recombination dominance. This bulk-dominated heating scenario would also diminish the strong front-back thermal gradient and instead produce more center-hot profiles, thereby altering T–V and J–V characteristics, internal field, charge transport and extraction, ionic redistribution, and hysteresis. Prior work has already shown that device J–V hysteresis depends strongly on whether recombination occurs at interfaces or in the bulk,⁵ suggesting that our thermo-electro-ionic dynamic conclusions are not universal but conditional on this assumption. Future studies should therefore incorporate bulk-recombination scenarios—by introducing trap densities or shortened carrier lifetimes into the model—to determine whether the effects reported here persist or are fundamentally reshaped under bulk-dominated conditions.

Applicability to inverted p–i–n architectures

Although the present study is parameterized for a regular n–i–p device in which ETL/perovskite interfacial recombination dominates, the coupled OET framework is architecture-independent at the governing-equation level. Extension to inverted p–i–n architectures primarily requires reassigning electrode-contact boundary conditions and re-parameterizing interface and contact parameters (e.g., interfacial recombination velocities and band offsets, and transport-layer thermal and electrical parameters) to reflect the inverted layer sequence and its dominant loss channels. Inverted devices can experience different thermo-electro-ionic dynamics than those observed here through several features. These include: (i) a shift in the dominant recombination pathway and its spatial localization (interface versus bulk), (ii) reversal of the built-in electric field polarity and interfacial band bending, which can affect field screening, extraction asymmetry, and hysteresis; (iii) differences in transport-layer chemistry (e.g., organic HTLs and fullerene-based ETLs) that reshape ionic accumulation and migration kinetics; and (iv) a different exposure to optical (i.e., illumination) and heat fluxes that reshapes photogeneration profiles, and therefore the spatial distribution of thermalization and recombination heating. Accordingly, while well-optimized p–i–n devices often exhibit reduced electrical hysteresis compared to n–i–p architectures,^24,25 the OET framework remains fully applicable and predicts that distinct thermo-electro-ionic signatures may still emerge under realistic operating conditions. These considerations clarify how architectural inversion alters, rather than invalidates, the coupled thermo-electro-ionic dynamics captured by the proposed model.

Conclusions

This work presents a comprehensive investigation of thermo-electro-ionic dynamics in perovskite solar cells, revealing how self-heating, ion migration, and electronic transport interact dynamically during bias sweeps. Using a fully coupled optical-electrical-thermal multiphysics framework, we show that the cell temperature evolves nonlinearly with scan rate, applied bias, and sweep direction—behavior that cannot be captured by conventional isothermal electronic–ionic transport models. It is found that internal heat generation dynamically reshapes the internal electric field, ion migration, and charge accumulation and extraction, particularly at intermediate scan rates where ionic, thermal, and electronic timescales overlap. This dynamic self-heating produces a dual-peak temperature profile in forward sweeps, governed by a crossover from interfacial to bulk heating under finite thermal inertia, and introduces a scan-rate-dependent thermal hysteresis. The resulting temperature gradients drive a Seebeck field that actively reshapes internal potentials, influencing both ionic and electronic transports. These effects distort the J–V response, particularly under intermediate scan rates, where S-shaped distortions, fill factor losses, and a systematic shift of the hysteresis maximum to lower scan rates emerge. By modulating the internal field, ion migration, and carrier distribution, self-heating acts as an active feedback pathway that couples thermal memory to electrical hysteresis, offering a previously unexplored explanation for anomalous thermal hysteresis.

Overall, our findings demonstrate that J–V hysteresis in PSCs is not solely governed by electro-ionic effects at a prescribed temperature, but can be significantly amplified and reshaped by internal thermal gradients evolving during operation. Neglecting thermo-electro-ionic coupling risks misinterpreting device transient behavior, especially when hysteresis appears or disappears depending on scan rate or measurement direction. Beyond conventional J–V analysis, the introduced T–V and P–V characteristics provide bias-resolved diagnostics that track temperature memory and attribute self-heating to specific loss channels, helping distinguish purely electrical hysteresis from thermo-electrical hysteresis. Therefore, self-heating dynamics are not peripheral but central to dynamic PSC response. Accurate modeling and mitigation of self-heating and ion migration—through materials with tailored thermal and electronic properties, nanostructural interface engineering to control heat-generation pathways, and optimized device architectures—are essential for achieving robust, hysteresis-free, and thermally stable PSCs.

Author contributions

Hadi Rostamzadeh: conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing – original draft. Hamid Montazeri: conceptualization, funding acquisition, project administration, resources, supervision, writing – review & editing.

Conflicts of interest

The authors declare no competing interests.

Data availability

A part of data supporting this study includes the version of the IonMonger code that is available online (https://github.com/PerovskiteSCModelling/IonMonger).

Supplementary information (SI): Document SI. Fig. S1–S16 and Tables S1, S2, and Supplemental references. See DOI: https://doi.org/10.1039/d5ee05840f.

References

1. X. Wang, H. Liu, F. Zhou, J. Dahan, X. Wang, Z. Li and W. Shen, Temperature Gradient-Induced Instability of Perovskite via Ion Transport, ACS Appl. Mater. Interfaces, 2018, 10(1), 835–844.
2. A. Shang and X. Li, Photovoltaic Devices: Opto-Electro-Thermal Physics and Modeling, Adv. Mater., 2017, 29(8), 1603492.
3. R. L. Milot, G. E. Eperon, H. J. Snaith, M. B. Johnston and L. M. Herz, Temperature-Dependent Charge-Carrier Dynamics in CH₃NH₃PbI₃ Perovskite Thin Films, Adv. Funct. Mater., 2015, 25(39), 6218–6227.
4. Z. Ai, T. Ma, Y. Zhang, Y. Bao, L. Shi, Z. Yang, Y. Zhan, L. Qin, G. Cao and X. Li, Unveiling Energy Conversion Mechanisms and Regulation Strategies in Perovskite Solar Cells, Small, 2024, 20(49), 2404012.
5. N. E. Courtier, J. M. Cave, J. M. Foster, A. B. Walker and G. Richardson, How Transport Layer Properties Affect Perovskite Solar Cell Performance: Insights from a Coupled Charge Transport/Ion Migration Model, Energy Environ. Sci., 2019, 12(1), 396–409.
6. L. K. Ono, S. R. Raga, S. Wang, Y. Kato and Y. Qi, Temperature-Dependent Hysteresis Effects in Perovskite-Based Solar Cells, J. Mater. Chem. A, 2015, 3(17), 9074–9080.
7. R. A. Belisle, W. H. Nguyen, A. R. Bowring, P. Calado, X. Li, S. J. C. Irvine, M. D. McGehee, P. R. F. Barnes and B. C. O'Regan, Interpretation of Inverted Photocurrent Transients in Organic Lead Halide Perovskite Solar Cells: Proof of the Field Screening by Mobile Ions and Determination of the Space Charge Layer Widths, Energy Environ. Sci., 2017, 10(1), 192–204.
8. W. Tress, N. Marinova, T. Moehl, S. M. Zakeeruddin, M. K. Nazeeruddin and M. Grätzel, Understanding the Rate-Dependent J–V Hysteresis, Slow Time Component, and Aging in CH₃NH₃PbI₃ Perovskite Solar Cells: The Role of a Compensated Electric Field, Energy Environ. Sci., 2015, 8(3), 995–1004.
9. M. Diethelm, T. Lukas, J. Smith, A. Dasgupta, P. Caprioglio, M. Futscher, R. Hany and H. J. Snaith, Probing Ionic Conductivity and Electric Field Screening in Perovskite Solar Cells: A Novel Exploration through Ion Drift Currents, Energy Environ. Sci., 2025, 18(3), 1385–1397.
10. Y. Zhang, Y. Bao, Y. Zhao, T. Ma, L. Shi, C. Zhang, Z. Ai, Y. Zhan, L. Qin, C. Wang, G. Cao, X. Li and Z. Yang, Tracing the Energy Losses in All-Perovskite Tandem Solar Cells From Opto-Electro-Thermal Perspectives, Adv. Funct. Mater., 2025, 2503408.
11. A. Kamppinen, H. Palonen and K. Miettunen, Self-Heating of Planar Perovskite Solar Cells Depending on Active Material Properties, ACS Appl. Energy Mater., 2024, 7(10), 4324–4334.
12. Y. An, C. Sheng and X. Li, Radiative Cooling of Solar Cells: Opto-Electro-Thermal Physics and Modeling, Nanoscale, 2019, 11(36), 17073–17083.
13. Y. An, C. Wang, G. Cao and X. Li, Heterojunction Perovskite Solar Cells: Opto-Electro-Thermal Physics, Modeling, and Experiment, ACS Nano, 2020, 14(4), 5017–5026.
14. N. E. Courtier, Interpreting Ideality Factors for Planar Perovskite Solar Cells: Ectypal Diode Theory for Steady-State Operation, Phys. Rev. Appl., 2020, 14(2), 024031.
15. J. M. Ball, S. D. Stranks, M. T. Hörantner, S. Hüttner, W. Zhang, E. J. W. Crossland, I. Ramirez, M. Riede, M. B. Johnston, R. H. Friend and H. J. Snaith, Optical Properties and Limiting Photocurrent of Thin-Film Perovskite Solar Cells, Energy Environ. Sci., 2015, 8(2), 602–609.
16. W. Clarke, L. J. Bennett, Y. Grudeva, J. M. Foster, G. Richardson and N. E. Courtier, IonMonger 2.0: Software for Free, Fast and Versatile Simulation of Current, Voltage and Impedance Response of Planar Perovskite Solar Cells, J. Comput. Electron., 2023, 22(1), 364–382.
17. S. Tammireddy, S. Reichert, Q. An, A. D. Taylor, R. Ji, F. Paulus, Y. Vaynzof and C. Deibel, Temperature-Dependent Ionic Conductivity and Properties of Iodine-Related Defects in Metal Halide Perovskites, ACS Energy Lett., 2022, 7(1), 310–319.
18. V. M. Le Corre, J. Diekmann, F. Peña-Camargo, J. Thiesbrummel, N. Tokmoldin, E. Gutierrez-Partida, K. P. Peters, L. Perdigón-Toro, M. H. Futscher, F. Lang, J. Warby, H. J. Snaith, D. Neher and M. Stolterfoht, Quantification of Efficiency Losses Due to Mobile Ions in Perovskite Solar Cells via Fast Hysteresis Measurements, Sol. RRL, 2022, 6(4), 2100772.
19. J. Thiesbrummel, S. Shah, E. Gutierrez-Partida, F. Zu, F. Peña-Camargo, S. Zeiske, J. Diekmann, F. Ye, K. P. Peters, K. O. Brinkmann, P. Caprioglio, A. Dasgupta, S. Seo, F. A. Adeleye, J. Warby, Q. Jeangros, F. Lang, S. Zhang, S. Albrecht, T. Riedl, A. Armin, D. Neher, N. Koch, Y. Wu, V. M. Le Corre, H. Snaith and M. Stolterfoht, Ion-Induced Field Screening as a Dominant Factor in Perovskite Solar Cell Operational Stability, Nat. Energy, 2024, 9(6), 664–676.
20. M. C. Gélvez-Rueda, N. Renaud and F. C. Grozema, Temperature Dependent Charge Carrier Dynamics in Formamidinium Lead Iodide Perovskite, J. Phys. Chem. C, 2017, 121(42), 23392–23397.
21. S. H. Lee, K. Moon, M. Shoaib, C. N. B. Pedorella, K. O’Brien, M.-J. Sher, S. Kim and T. L. Cocker, Temperature-Dependent Recombination Dynamics of Photocarriers in CsPbBr₃ Microcrystals Revealed by Ultrafast Terahertz Spectroscopy, Adv. Opt. Mater., 2024, 12(30), 2401162.
22. F. Ye, S. Zhang, J. Warby, J. Wu, E. Gutierrez-Partida, F. Lang, S. Shah, E. Saglamkaya, B. Sun, F. Zu, S. Shoaee, H. Wang, B. Stiller, D. Neher, W.-H. Zhu, M. Stolterfoht and Y. Wu, Overcoming C₆₀-Induced Interfacial Recombination in Inverted Perovskite Solar Cells by Electron-Transporting Carborane, Nat. Commun., 2022, 13(1), 7454.
23. J. Warby, F. Zu, S. Zeiske, E. Gutierrez-Partida, L. Frohloff, S. Kahmann, K. Frohna, E. Mosconi, E. Radicchi, F. Lang, S. Shah, F. Peña-Camargo, H. Hempel, T. Unold, N. Koch, A. Armin, F. De Angelis, S. D. Stranks, D. Neher and M. Stolterfoht, Understanding Performance Limiting Interfacial Recombination in Pin Perovskite Solar Cells, Adv. Energy Mater., 2022, 12(12), 2103567.
24. M. Saliba, J.-P. Correa-Baena, C. M. Wolff, M. Stolterfoht, N. Phung, S. Albrecht, D. Neher and A. Abate, How to Make over 20% Efficient Perovskite Solar Cells in Regular (n–i–p) and Inverted (p–i–n) Architectures, Chem. Mater., 2018, 30(13), 4193–4201.
25. I. Levine, P. K. Nayak, J. T.-W. Wang, N. Sakai, S. Van Reenen, T. M. Brenner, S. Mukhopadhyay, H. J. Snaith, G. Hodes and D. Cahen, Interface-Dependent Ion Migration/Accumulation Controls Hysteresis in MAPbI₃ Solar Cells, J. Phys. Chem. C, 2016, 120(30), 16399–16411.

---