Evolutionary Generator and Path-Integral Control Framework for Underground Heat Recovery: Unifying Local and Fractional Transport

1. Introduction

Medium-to-long-term underground heat recovery systems, such as ground source heat pumps (GSHPs) and underground thermal energy storage (UTES), play a foundational role in large-scale decarbonization and sustainable energy utilization [1]. However, under long-term cyclic operations spanning multiple seasons, these systems frequently encounter global thermal imbalance and progressive thermal field degradation [2]. This irreversible performance deterioration not only reduces thermodynamic efficiency but also poses a critical sustainability bottleneck in contemporary subsurface thermal engineering.

Conventional engineering designs and numerical simulations primarily rely on the classical parabolic heat conduction equation formulated by Fourier’s law. However, these macroscopic models frequently exhibit severe deviations from long-term monitoring data when predicting sustained, multi-cycle thermal evolutions over decades [3]. The physical root of this discrepancy lies in the fact that the classical heat equation is governed by a local Gaussian generator, which assumes a Markovian, memoryless mechanism. This assumption fundamentally fails to capture the non-local long-range correlations and intrinsic memory effects induced by complex porous geological media under prolonged cyclic loading [4]. Consequently, a central physical question arises: when traditional local operators fail to encompass multi-scale anomalous migration, what underlying physical mechanism governs this large-scale evolutionary process?

To mitigate the aforementioned long-term thermal imbalance, researchers have explored the integration of phase change materials (PCMs) to modulate thermal dissipation. Recent engineering advances focus on two paradigms: PCM ground backfilling and structural encapsulation. For instance, Ahmed et al. [5] investigated GSHP systems utilizing PCM-enhanced backfill, revealing that latent heat buffering significantly stabilizes the circulating fluid temperature during peak loads. Subsequently, Aljabr et al. [6] conducted a comprehensive numerical study of PCM-integrated boreholes, identifying that increasing PCM thickness enhances heat storage capacity but follows a diminishing-return law under continuous extraction. To optimize thermal response under complex loads, recent studies such as Mousa et al. [7] proposed multi-layer PCM configurations, demonstrating that graded melting temperatures can effectively bridge the gap between regional thermal demand and subsurface dissipation rates.

Building upon these local interventions, recent studies have approached system-level thermal management through advanced composite structures. Goderis et al. [8] provided a systematic analysis of gradient PCMs, highlighting that multiple gradient PCMs can maintain a more stable heat extraction rate compared to single-PCM systems. Similarly, recent experimental work by Žirgulis et al. [9] on PCM-geothermal heat exchangers confirmed that phase-front morphology directly influences global temperature distribution over multi-year horizons. Furthermore, theoretical analysis by Feng et al. [10] on the effective thermal conductivity of porous media with phase change indicated that the non-linear storage term alters effective thermal resistance. However, a critical theoretical gap persists: virtually all PCM mitigation strategies remain grounded in heuristic optimization or parameter sweeping. Heat transfer engineers utilize PCMs to passively buffer temperature, yet they lack a rigorous mathematical framework to describe how gradient PCMs fundamentally alter the underlying long-term evolutionary geometry of the global heat field across vast spatiotemporal scales.

In non-equilibrium statistical physics, by contrast, the foundational mechanisms of non-local transport and structural memory have been systematically elucidated through anomalous diffusion. Representative work led by Liang’s group [11] decoded anomalous sub-diffusion, showing that power-law relaxation and non-Gaussian states naturally emerge in heterogeneous environments with complex trapping potentials. Advancing the computational frontier, fractional Langevin equations [12] have enabled characterization of memory-dependent trajectories in fractal-like media. In the context of spatial disorder, Song et al. [13] demonstrated that fractional diffusion operators intrinsically embed spatial long-range correlations that classical local operators miss. Most critically, regarding control, Zhang [14] derived analytical results on the optimal control of diffusion processes through time-inhomogeneous generators, providing a mathematical pathway for entropy-production minimization.

This contrast exposes a methodological separation between engineering application and fundamental physics. Thermal engineering uses PCM arrays to pragmatically modulate cross-cycle behavior without a full fractional stochastic framework; meanwhile, statistical physics formulates fractional evolutionary equations but seldom maps these theories onto actionable thermodynamic-conjugate engineering structures, such as gradient PCMs.

To bridge this disciplinary gap, this paper introduces a unified framework based on the concept that the evolutionary generator intrinsically determines the dissipative geometry of the system. We posit that gradient PCMs act as physical realizations of structural modifiers for the heat-field generator. Based on this premise, we propose a hybrid evolutionary generator with intrinsic non-locality, incorporating both a local Laplacian and a non-local Lévy operator. By mapping this mechanism onto the Onsager variational principle and the Martin-Siggia-Rose (MSR) formalism, we reinterpret cyclic heat conduction as a gradient flow intertwined with a thermodynamic conjugate information field. We further demonstrate mathematically that gradient PCMs reshape the field topology to ensure intrinsic sustainability. Subsequently, based on this structurally modulated stable field, we deploy a path-integral optimal-control strategy for the thermal source, pushing the thermodynamic efficiency of medium-to-long-term underground heat recovery toward its theoretical limit.

The remainder of this paper is organized as follows: Section 2 reviews experimental phenomena, revealing the failure of the classical local Gaussian equation under macroscopic memory effects. Section 3 rigorously derives the hybrid fractional generator, validates the MSR field-theoretic framework, and establishes the optimal-control architecture. Section 4 provides concluding remarks on the physical and engineering implications.

2. Experimental Observations and Intrinsic Limitations of the Classical Heat Equation

Our previous experimental findings [15] show that the dynamical behavior of medium-to-long-term underground heat recovery fundamentally differs from classical thermal diffusion. Specifically, the differences manifest in three dimensions:

(i) Emergence of structured steady/cyclic states: As underground heat-recovery cycles continue, the global thermal-state function g_appr. in sustainable Scenario B exhibits harmonic-like oscillations driven by the periodic heat source. Conversely, in unsustainable Scenario A, it exhibits an irreversible exponential decay superimposed on periodic oscillations (as shown in Figure 1), indicating failure of the system to converge to a bounded periodic regime.

(ii) Cross-cycle structural memory: The thermal-field structural entropy S exhibits macro-cyclic evolutionary behavior exclusively in sustainable Scenario B (as shown in Figure 2), characterizing the system’s ability to maintain organizational stability after multiple heat-injection/extraction shocks.

(iii) Cross-scale correlations: The structural-entropy vectors across scales and the global thermal-state vectors remain strictly parallel (as shown in Figure 3), indicating that the subsurface thermal structure possesses high cross-scale coherence and intrinsic self-similarity.

Why does the classical heat equation inherently fail to accurately predict medium-to-long-term underground heat-recovery behavior? The following section starts from the evolutionary generator of classical thermal diffusion and dissects the underlying physical limitations of the classical heat equation from first principles.

In Scenario A, the global thermal baseline exhibits periodic oscillations superimposed on an exponential downward drift, indicating progressive degradation and a fundamental failure to converge toward a bounded periodic state under cyclic operation.

Structural entropy quantifies the distribution of thermal structures across spatial scales. Bounded, small-amplitude, and macrocycle entropy oscillations (Scenario B) indicate resilient and sustainable thermal-field organization governed by robust structural memory.

High directional similarity across scales confirms cross-scale coherence and existence–uniqueness properties of diffusion-dominated thermal fields, supporting the use of compact structural descriptors to represent long-term subsurface thermal evolution.

2.1. The evolutionary Generator of Classical Thermal Diffusion

The standard diffusion equation is parabolic (Eq. 1), which is the natural consequence of Fourier’s constitutive relation (Eq. 2) within the framework of the first law of thermodynamics. The fundamental solution of this equation is the Gaussian kernel, which describes the spatiotemporal evolution of macroscopic observables such as temperature in physical space (Eq. 3). Under a uniform initial temperature distribution and natural boundary conditions, the temperature field driven by heat source f can be strictly expressed as the convolution of f and the Gaussian kernel (Eq. 4).

$$\left\{\begin{array}{l}\left({\partial }_t-a^2\Delta \right)\varphi \left(x,t\right)=0\hfill \\ \varphi \left(x,0\right)=\mathsf{\delta }\left(x\right)\hfill \\ \underset{\left|x\right|\to \infty }{\lim }\varphi \left(x,t\right)=0\hfill \end{array}\right.$$

(1)

$$\mathrm{q}\propto -\nabla \varphi \left(x,t\right)$$

(2)

$$\varphi \left(x,t\right)=K\left(x,t\right)={\displaystyle \frac{1}{\sqrt{\mathsf{\pi }}}}{\displaystyle \frac{1}{2\sqrt{a^{2}t}}}{\mathrm{e}}^{-{\displaystyle \frac{x^2}{4a^{2}t}}}$$

(3)

$$\begin{array}{c}\varphi \left(x,t\right)={\displaystyle {\int }_0^t\mathrm{d}\tau {\displaystyle {\int }_{{ℝ}^n}f\left(\xi ,\tau \right)K\left(x-\xi ,t-\tau \right)\mathrm{d}\xi }}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{1}{\sqrt{\mathsf{\pi }}}}{\displaystyle {\int }_0^t\frac{1}{2\sqrt{a^2\left(t-\tau \right)}}\mathrm{d}\tau {\displaystyle {\int }_{{ℝ}^n}f\left(\xi ,\tau \right){\mathrm{e}}^{-\frac{{\left(x-\xi \right)}^2}{4a^2\left(t-\tau \right)}}\mathrm{d}\xi }}\hfill \end{array}$$

(4)

where ϕ represents the dissipative diffusion field function (e.g., temperature); δ is the Dirac delta function; q is the heat-flux density vector within the thermal field; K denotes the Gaussian kernel; f is the external heat-source function; ℝⁿ denotes n-dimensional real space; ∆ is the classical Laplacian operator; ∂ is the partial differential operator; d is the differential operator; t and τ represent time; x and ξ are continuous variables in n-dimensional real space; and a² is the apparent thermal diffusivity of the thermal-field medium.

At the microscopic level, classical diffusion is described by a stochastic differential equation (Eq. 5). Here, μ is the drift coefficient vector determining the macroscopic translational trend of evolution; σ is the diffusion coefficient matrix determining the amplitude and correlation of random perturbations; and W_t is the standard Wiener process. The infinitesimal generator A of this diffusion process is a second-order elliptic spatial partial differential operator acting on a smooth function g (Eq. 6). Through a linear combination of the first derivative (drift) and the second derivative (diffusion), the generator characterizes the expected instantaneous rate of change of any observation function g at the current state x.

Based on the current value of an arbitrary function g, superimposing the cumulative expectation of the generator over time yields the expected value of the function after a certain period (Eq. 7). If ϕ denotes this expected value, then ϕ naturally satisfies a Cauchy problem (Eq. 8), the solution of which is the convolution of the fundamental solution and the initial distribution g (Eq. 9). Comparing Eq. 9 with Eq. 4 reveals that: (i) The generator serves as the core topological bridge connecting the microscopic stochastic process with the macroscopic deterministic partial differential equation. In the evolution of observables, it acts as the spatial evolutionary generator, jointly formulating the spatiotemporal evolution rules of ϕ with the temporal metric ∂_t. (ii) This rule is mathematically represented as a convolution structure. The convolution kernel (Gaussian kernel) is entirely dominated by the second-derivative term of the generator, whereas the first-derivative term only accounts for spatial translation (e.g., the −ct term in the exponential function) and has no effect on the geometric topology of diffusion.

$$\mathrm{d}{X}_t=\mu \left(X_t\right)\mathrm{d}t+\sigma \left(X_t\right)\mathrm{d}{W}_t$$

(5)

$$Ag\left(x\right)={\displaystyle \sum _{i=1}^n{\mu }_i\frac{\partial g}{\partial x_i}\left(x\right)+\frac{1}{2}{\displaystyle \sum _{i,j=1}^n{\left(\sigma {\sigma }^{\mathsf{{\rm T}}}\right)}_{ij}\frac{{\partial }^{2}g}{\partial x_i\partial x_j}\left(x\right)}}$$

(6)

$${\mathrm{E}}^x\left[g\left(X_{\tau }\right)\right]=g\left(x\right)+{\mathrm{E}}^x\left[{\displaystyle {\int }_0^{\tau }Ag\left(X_s\right)\mathrm{d}s}\right]$$

(7)

$$\left\{\begin{array}{l}\left({\partial }_t-A\right)\varphi \left(x,t\right)=0\hfill \\ \varphi \left(x,0\right)=g\left(x\right)\hfill \end{array}\right.$$

(8)

$$\begin{array}{c}\varphi \left(x,t\right)=g\left(x\right)\ast K\left(x,t\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{1}{\sqrt{\mathsf{\pi }}}}{\displaystyle \frac{1}{2\sqrt{a^{2}t}}}{\displaystyle {\int }_{{ℝ}^n}g\left(\xi \right){\mathrm{e}}^{-\frac{{\left(x-\xi -ct\right)}^2}{4a^{2}t}}\mathrm{d}\xi }\hfill \end{array}$$

(9)

where E^x[·] denotes the mathematical expectation; and * represents the convolution operation.

In one-dimensional space, the stochastic differential equation simplifies to Eq. 10. The corresponding generator A_B (Eq. 11) reveals the essence of Brownian motion: the classical Laplacian term a²∆ dictates that the root-mean-square displacement of Brownian motion strictly obeys (a²t)^1/2, thereby fundamentally locking the system into a continuous, memoryless Gaussian process at the lowest level. Consequently, the classical heat equation with the Gaussian generator a²∆, under the uniformly flowing temporal metric ∂_t, is fundamentally deprived of the capability to produce non-local long-range phenomena such as cross-cycle memory and stable cyclic states.

$$\mathrm{d}{B}_t=-c\cdot \mathrm{d}t+\sqrt{2}a\cdot \mathrm{d}{W}_t$$

(10)

$$\begin{array}{c}{A}_{\mathrm{B}}g\left(x\right)=\left(-c\right)\cdot g^{\prime }\left(x\right)+{\displaystyle \frac{1}{2}}\cdot {\left(\sqrt{2}a\right)}^2\cdot g^{″}\left(x\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=a^2\cdot g^{″}\left(x\right)-c\cdot g^{\prime }\left(x\right)\hfill \end{array}$$

(11)

where −c represents the migration velocity along the positive direction of x.

2.2. Physical Implications Within the Gradient Flow Framework

According to the Onsager variational principle, the actual evolutionary path of a system must extremize the functional representing the sum of the dissipation potential and the rate of free-energy dissipation. Even an irreversible process such as pure diffusion strictly adheres to this minimization principle. We define the Rayleighian functional ℜ as shown in Eq. 12, where Ξ is the dissipation function characterizing the quadratic form of the evolutionary rate (Eq. 13), and ϒ is the free energy driving evolution (i.e., entropy potential or Dirichlet energy, which is also a quadratic form, Eq. 14). The actual evolutionary path minimizes this functional (Eq. 15), which is precisely equivalent to the classical heat equation.

$$\Re \left[{\partial }_t\varphi ;\varphi \right]=\Xi \left({\partial }_t\varphi ,{\partial }_t\varphi \right)+{\displaystyle \frac{\mathrm{d}\gamma }{\mathrm{d}t}}\left({\partial }_t\varphi ;\varphi \right)$$

(12)

$$\Xi \left({\partial }_t\varphi ,{\partial }_t\varphi \right)={\displaystyle \frac{1}{2}}{\displaystyle {\int }_{{ℝ}^n}{‖{\partial }_t\varphi ‖}^2\mathrm{d}x}$$

(13)

$$\gamma \left[\varphi \right]={\displaystyle \frac{1}{2}}〈\varphi ,-A\varphi 〉={\displaystyle \frac{1}{2}}{\displaystyle {\int }_{{ℝ}^n}\varphi \cdot \left(-A\varphi \right)\mathrm{d}x}$$

(14)

$$\begin{array}{c}{\displaystyle \frac{\mathsf{\delta }\Re }{\mathsf{\delta }\left({\partial }_t\varphi \right)}}\left[{\partial }_t\varphi ;\varphi \right]={\displaystyle \frac{\mathsf{\delta }}{\mathsf{\delta }\left({\partial }_t\varphi \right)}}\left[{\displaystyle \frac{1}{2}}{\displaystyle {\int }_{{ℝ}^n}{‖{\partial }_t\varphi ‖}^2\mathrm{d}x}\right]+{\displaystyle \frac{\mathsf{\delta }}{\mathsf{\delta }\left({\partial }_t\varphi \right)}}\left[{\displaystyle {\int }_{{ℝ}^n}\frac{\mathsf{\delta }\gamma }{\mathsf{\delta }\varphi }{\partial }_t\varphi \mathrm{d}x}\right]\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\partial }_t\varphi +{\displaystyle \frac{\mathsf{\delta }\gamma }{\mathsf{\delta }\varphi }}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\partial }_t\varphi +\left(-A\varphi \right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=0\hfill \end{array}$$

(15)

where δ denotes the variational operator.

Eq. 15 shows that the classical heat equation is essentially the Euler−Lagrange equation of the Onsager variational principle. Within this framework, heat conduction is not merely the empirical intuition of “heat flowing from high to low temperatures”; rather, driven by the free-energy functional, the macroscopic state of the system evolves toward equilibrium along the steepest-descent path of free energy within the L² metric defined by the Onsager operator. The macroscopic state under the action of the generator is the variational derivative of this free energy, and the Gaussian kernel is the analytical solution of this gradient flow in free space. Importantly, the external heat source f does not alter the intrinsic evolutionary properties of the system—dissipation is embedded in the underlying spatiotemporal geometry and is the thermodynamic consequence of descent along the free-energy gradient.

2.3. Physical Implications Within the Field Theory Framework

By leveraging the Martin-Siggia-Rose (MSR) formalism, the classical heat equation can be incorporated into a field-theoretic framework: a functional Dirac delta function is introduced to constrain the deterministic dynamical equation (Eq. 16), where the integration measure corresponding to the auxiliary real field ψ is the functional Lebesgue measure. The evaluation of any functional of the system on the solution manifold can be expressed as Eq. 17, and the action S (Eq. 18), together with its Lagrange density L, becomes the generating functional of the primary equation system. In this thermodynamic context from microscopic to macroscopic levels, the action S characterizes the “total informational energy” of the system.

$$\mathsf{\delta }\left[{\partial }_t\varphi -A\varphi -f\right]={\displaystyle \int D\psi \exp \left(\mathrm{i}{\displaystyle \int \mathrm{d}t\mathrm{d}x\psi \left[{\partial }_t\varphi -A\varphi -f\right]}\right)}$$

(16)

$$〈F\left[\varphi \right]〉={\displaystyle \int D\varphi F\left[\varphi \right]\mathsf{\delta }\left[{\partial }_t\varphi -A\varphi -f\right]}\propto {\displaystyle \int D\varphi D\psi F\left[\varphi \right]{\mathrm{e}}^{\mathrm{i}S\left[\varphi ,\psi \right]}}$$

(17)

$$S\left[\varphi ,\psi \right]={\displaystyle {\int }_0^{\infty }\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\psi \left({\partial }_t\varphi -A\varphi -f\right)\left(x,t\right)}}={\displaystyle {\int }_0^{\infty }\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}L\mathrm{d}x}}$$

(18)

$$\left\{\begin{array}{l}\left({\partial }_t\varphi -A\varphi -f\right)\left(x,t\right)=0\hfill \\ \left(-{\partial }_t\psi -A\psi \right)\left(x,t\right)=0\hfill \end{array}\right.$$

(19)

where ψ is the conjugate information field accompanying the physical field ϕ; D represents the functional integration measure; and i is the imaginary unit.

In mathematical terms, the MSR formalism is a constrained variational principle, wherein the auxiliary field ψ acts as a Lagrange multiplier in spacetime. Taking variation with respect to ϕ yields the adjoint equation that ψ must satisfy (Eq. 19; detailed derivations are provided in Appendices A1 and A2). In physical terms, the adjoint equation of the thermal field possesses time-reversal character. ψ has two roles: (i) as the thermodynamic conjugate field, it enforces consistency of the irreversible process with the Second Law of Thermodynamics; and (ii) as an order parameter for the organizational structure of the physical field ϕ, it governs the backward propagation mechanism of structural information in the thermal field.

Therefore, the classical Gaussian generator generates a dissipative structure conforming to the (a²t)^1/2 scaling law. Classical diffusion based on a²Δ thus leads to field-theoretic information indicating structural collapse. This elucidates why mainstream heat equations fundamentally fail to guide “high-resilience, long-term” deep-decarbonization engineering.

3. The Evolutionary Generator with Intrinsic Non-Locality

The core question addressed in this section is: What evolutionary generator is intrinsically embedded within sustainable medium-to-long-term underground heat recovery?

3.1. Discovery of the Hybrid Generator

Based on the experimental observations, we arrive at the following conclusions:

1. There exists no singular jump or bifurcation in macroscopic observables between the unsustainable and sustainable heat-recovery scenarios. This suggests a generalized, unified theoretical framework that naturally encompasses both diffusion modes: exhibiting long-range correlation features under critical sustainability conditions, while degenerating to Gaussian diffusion under conventional conditions.

2. The observed cross-cycle structural memory indicates that the system has evolved into a non-equilibrium stationary state (NESS). Maintaining such a resilient dissipative structure necessarily relies on a delicate microscopic balance within the system.

3. Local and non-local mechanisms coexist in the subsurface medium, dominating evolution in different spatiotemporal regimes; consequently, a critical crossover scale must exist where they compete.

4. Structural memory is an emergent phenomenon of macroscopic power-law long tails, arising from the dominance of non-local mechanisms. It is not an independent entity and cannot be artificially mimicked by introducing an ad hoc relaxation time, a procedure that violates first principles.

Guided by this physical logic, we construct a generalized hybrid generator capable of accommodating both unsustainable and sustainable behaviors: namely, the weighted superposition of the classical local Laplacian operator and the non-local fractional Laplacian operator (Eq. 20). Under the physically observed temporal metric ∂_t, these two operators independently generate two distinct spatial scaling laws, (a²t)^1/2 and (b²t)^1/β, respectively. Their counterbalance naturally gives rise to the critical length scale ℓ_cr and the prolonged polarization transition period t_cr (Eq. 21), while simultaneously constraining the permissible ranges of the apparent thermal diffusivities a² and b², as well as the fractional jump dimension β (Eq. 22). Experimental validation [15] confirms that all parameters a², b², and β can be engineered and modulated via gradient PCMs, thereby enabling sustainable medium-to-long-term underground heat recovery in practice.

$$\begin{array}{cc}A=-a^2\left(-\Delta \right)-b^2{\left(-\Delta \right)}^{{\displaystyle \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}},\hfill & 0<\beta <2\hfill \end{array}$$

(20)

$$\left\{\begin{array}{l}{\mathcal{l}}_{\mathrm{cr}}\sim {\displaystyle \frac{1}{\left|k_{\mathrm{cr}}\right|}}\sim {\left({\displaystyle \frac{a^2}{b^2}}\right)}^{{\displaystyle \frac{1}{2-\beta }}}\hfill \\ t_{\mathrm{cr}}\sim {\displaystyle \frac{{\left({\mathcal{l}}_{\mathrm{cr}}\right)}^{\beta }}{b^2}}\sim {\left({\displaystyle \frac{{\left(a^2\right)}^{\beta }}{{\left(b^2\right)}^2}}\right)}^{{\displaystyle \frac{1}{2-\beta }}}\sim {\left({\displaystyle \frac{{\left|a\right|}^{\beta }}{{\left|b\right|}^2}}\right)}^{{\displaystyle \frac{2}{2-\beta }}}\hfill \end{array}\right.$$

(21)

$$\left\{\begin{array}{l}{\displaystyle \frac{a^2}{b^2}}\sim 1\hfill \\ 0<\beta \ll 1\hfill \end{array}\right.$$

(22)

where (−∆)^β/2 denotes the fractional Laplacian operator; β is the fractional jump dimension governing the spatial non-linear damping of the medium; b² is the second apparent thermal diffusivity coefficient; and k_cr is the critical scale in Fourier space.

Re-substituting Eqs. 20−22 into the gradient-flow variational framework (Eqs. 12−15) reveals a key insight: although the fundamental dissipative nature of the system remains unchanged, its trajectory toward equilibrium in L² space is substantially reconfigured (as illustrated in Figure 1). The extremely small fractional order β induces strong long-range correlations, resulting in a high initial dissipation rate at system startup. Furthermore, during the protracted transition period t_cr, the dynamic equilibrium between the two mechanisms steers the system onto a prolonged power-law tail trajectory. Under periodic excitation f, the coefficients characterizing the competing mechanisms must be close in magnitude and significantly less than unity to prevent the system from diverging into instability.

Subsequently, we subject this hybrid generator to field-theoretic symmetry (conservation-law) tests. Eq. 23 confirms invariance under time translation, and Eq. 24 verifies spatial translation invariance. Eq. 25 demonstrates that the hybrid generator exhibits approximate scale invariance only when parameters a², b², and fractional order β are all exceedingly small. This condition explains the parallel evolution of structural entropy across scales.

In summary, the sustainability of medium-to-long-term heat recovery is a stringent dynamic balance between Gaussian local thermal dissipation and fractional long-range interactions. The action functional S must be driven to the vicinity of the scale-invariant critical boundary. From an engineering perspective, the role of gradient PCMs is reinterpreted: they are not merely sensible/latent-heat reservoirs, but physical intervention agents introduced into the subsurface structure to reshape this fractional evolutionary generator, thereby regulating the propagation rates and pathways of both thermal energy and its conjugate information field.

$$\begin{array}{cc}\because H\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\left(\frac{\partial L}{\partial \left({\partial }_t\varphi \right)}{\partial }_t\varphi +\frac{\partial L}{\partial \left({\partial }_t\psi \right)}{\partial }_t\psi -L\right)}\hfill \\ & ={\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\left(\psi {\partial }_t\varphi -L\right)}\hfill \\ & ={\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\left(\psi {\partial }_t\varphi -\psi {\partial }_t\varphi +\psi A\varphi \right)}\hfill \\ & ={\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\left(\psi A\varphi \right)}\hfill \\ \therefore {\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}}\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\left({\partial }_t\psi A\varphi +\psi A{\partial }_t\varphi \right)}\hfill \\ & ={\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\left(-A\psi A\varphi +A\psi A\varphi \right)}\hfill \\ & =0\hfill \end{array}$$

(23)

$$\begin{array}{ll}\because {\rm P}\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\left(\frac{\partial L}{\partial \left({\partial }_t\varphi \right)}{\partial }_i\varphi +\frac{\partial L}{\partial \left({\partial }_t\psi \right)}{\partial }_i\psi \right)}={\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\psi \nabla \varphi }\hfill \\ \therefore {\displaystyle \frac{\mathrm{d}{\rm P}}{\mathrm{d}t}}\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\psi \nabla \varphi }\hfill \\ & ={\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\left({\partial }_t\psi \nabla \varphi +\psi \nabla {\partial }_t\varphi \right)}\hfill \\ & ={\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\left(-A\psi \nabla \varphi +\psi \nabla A\varphi \right)}\hfill \\ & ={\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\left(-A\psi \nabla \varphi +A\psi \nabla \varphi \right)}\hfill \\ & =0\hfill \end{array}$$

(24)

$$\begin{array}{c}\because \left\{\begin{array}{l}{x}^{\prime }=\lambda x\hfill \\ t^{\prime }={\lambda }^{s}t\hfill \end{array}\right.\hfill \\ \therefore \left\{\begin{array}{l}{\varphi }^{\prime }\left(x^{\prime },t^{\prime }\right)={\lambda }^{-{\Delta }_{\varphi }}\varphi \left(x,t\right)\hfill \\ {\psi }^{\prime }\left(x^{\prime },t^{\prime }\right)={\lambda }^{-{\Delta }_{\psi }}\psi \left(x,t\right)\hfill \end{array}\right.\hfill \\ \therefore \mathrm{d}{t}^{\prime }\mathrm{d}{x}^{\prime }={\lambda }^{n+s}\mathrm{d}t\mathrm{d}x\hfill \\ \because \left\{\begin{array}{l}{\partial }_t\to {\lambda }^{-s}{\partial }_t\hfill \\ \nabla \to {\lambda }^{-1}\nabla \hfill \end{array}\right.\hfill \\ \therefore \left\{\begin{array}{l}\left(-\Delta \right)\to {\lambda }^{-2}\left(-\Delta \right)\hfill \\ {\left(-\Delta \right)}^{{\displaystyle \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}\to {\lambda }^{-\beta }{\left(-\Delta \right)}^{{\displaystyle \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}\hfill \end{array}\right.\hfill \\ \therefore L^{\prime }={\lambda }^{-{\Delta }_{\psi }}\psi \left({\lambda }^{-s}{\partial }_t\left({\lambda }^{-{\Delta }_{\varphi }}\varphi \right)-a^2{\lambda }^{-2}\left(-\Delta \right)\left({\lambda }^{-{\Delta }_{\varphi }}\varphi \right)-b^2{\lambda }^{-\beta }{\left(-\Delta \right)}^{{\displaystyle \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}\left({\lambda }^{-{\Delta }_{\varphi }}\varphi \right)\right)\hfill \\ \hspace{1em}\hspace{1em}={\lambda }^{-{\Delta }_{\psi }-{\Delta }_{\varphi }-s}\psi {\partial }_t\varphi -{\lambda }^{-{\Delta }_{\psi }-{\Delta }_{\varphi }-2}{a}^2\left(-\Delta \right)\varphi -{\lambda }^{-{\Delta }_{\psi }-{\Delta }_{\varphi }-\beta }{b}^2{\left(-\Delta \right)}^{{\displaystyle \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}\varphi \hfill \\ \therefore S^{\prime }={\displaystyle {\int }_0^{\infty }\mathrm{d}{t}^{\prime }{\displaystyle {\int }_{{ℝ}^n}\mathrm{d}{x}^{\prime }{L}^{\prime }}}={\lambda }^{n+s}{\displaystyle {\int }_0^{\infty }\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x{L}^{\prime }}}\hfill \\ \hspace{1em}\hspace{1em}={\lambda }^{n+s}{\displaystyle {\int }_0^{\infty }\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\left({\lambda }^{-{\Delta }_{\psi }-{\Delta }_{\varphi }-s}\psi {\partial }_t\varphi -{\lambda }^{-{\Delta }_{\psi }-{\Delta }_{\varphi }-2}{a}^2\left(-\Delta \right)\varphi -{\lambda }^{-{\Delta }_{\psi }-{\Delta }_{\varphi }-\beta }{b}^2{\left(-\Delta \right)}^{\raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\varphi \right)}}\hfill \\ \therefore \left\{\begin{array}{l}n+s-{\Delta }_{\psi }-{\Delta }_{\varphi }-s=n-\left({\Delta }_{\psi }+{\Delta }_{\varphi }\right)=0\hfill \\ n+s-{\Delta }_{\psi }-{\Delta }_{\varphi }-2=s-2=0\hfill \\ n+s-{\Delta }_{\psi }-{\Delta }_{\varphi }-\beta =s-\beta =0\hfill \end{array}\right.\hfill \\ \therefore \left\{\begin{array}{l}0<a^2\sim b^2\ll 1\hfill \\ 0<\beta \ll 1\hfill \end{array}\right.\hfill \end{array}$$

(25)

where L is the Lagrange density; H is the conserved quantity (Hamiltonian), reflecting the coupling strength between thermal field ϕ and information field ψ; Ρ is the generalized momentum; λ is the scale transformation coefficient; s is the dynamic exponent; and ∆_ϕ, ∆_ψ are the scaling dimensions of the thermal and information fields, respectively.

At the stochastic-dynamics level, employing the Itô−Lévy calculus framework, we construct a composite random process by concatenating the standard Brownian-motion term with a Lévy jump process governed by a Poisson intensity measure (Eqs. 26−29). A further application of Itô’s formula proves that this hybrid generator is the macroscopic mathematical projection of the underlying stochastic dynamics (Eqs. 30−35).

$$\begin{array}{c}\mathrm{d}{X}_t=\mathrm{d}{X}_t^{\left(\mathrm{B}\right)}+\mathrm{d}{X}_t^{\left(\mathrm{N}-\mathrm{L}\right)}\hfill \\ \hspace{1em}\hspace{1em}=\sqrt{2}a\mathrm{d}{W}_t+{\displaystyle {\int }_{\left|y\right|<1}y\tilde{N}\left(\mathrm{d}t,\mathrm{d}y\right)}+{\displaystyle {\int }_{\left|y\right|\ge 1}yN\left(\mathrm{d}t,\mathrm{d}y\right)}\hfill \end{array}$$

(26)

$$\begin{array}{c}\mathrm{E}\left[N\left(\mathrm{d}t,\mathrm{d}y\right)\right]=\mathrm{d}t\nu \left(\mathrm{d}y\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\mathrm{d}t\mathrm{d}y{\displaystyle \frac{1}{{\left|y\right|}^{1+\beta }}}{\displaystyle \frac{b^2}{{\displaystyle {\int }_{ℝ}\frac{1-\cos w}{{\left|w\right|}^{1+\beta }}\mathrm{d}w}}}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\mathrm{d}t\mathrm{d}y{\displaystyle \frac{1}{{\left|y\right|}^{1+\beta }}}C\left(\beta \right)b^2\hfill \end{array}$$

(27)

$$\mathrm{E}\left[{\mathrm{e}}^{\mathrm{i}k{X}_t^{\left(\mathrm{N}-\mathrm{L}\right)}}\right]={\mathrm{e}}^{-b^2{\left|k\right|}^{\beta }t}$$

(28)

$$\tilde{N}\left(\mathrm{d}t,\mathrm{d}y\right)=N\left(\mathrm{d}t,\mathrm{d}y\right)-\mathrm{d}t\nu \left(\mathrm{d}y\right)$$

(29)

$$\begin{array}{c}g\left(X_t\right)=g\left(X_0\right)+{\displaystyle {\int }_0^t{g}^{\prime }\left(X_{s-}\right)\mathrm{d}{X}_s}+{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^t{g}^{″}\left(X_{s-}\right)\mathrm{d}{\left[X^{\mathrm{c}}\right]}_s}\hfill \\ \hspace{1em}\hspace{1em}=g\left(X_0\right)+{\displaystyle \sum _{0<s\le t}\left(g\left(X_s\right)-g\left(X_{s-}\right)-g^{\prime }\left(X_{s-}\right)\Delta X_s\right)}\hfill \\ \Delta X_s=X_s-X_{s-}\hfill \end{array}$$

(30)

$$\begin{array}{l}\because \mathrm{d}{X}_t^{\mathrm{c}}=\sqrt{2}a\mathrm{d}{W}_t\hfill \\ \therefore \mathrm{d}{\left[X^{\mathrm{c}}\right]}_t=2a^2\mathrm{d}t\hfill \\ \therefore g\left(X_t\right)=g\left(X_0\right)+{\displaystyle {\int }_0^t\sqrt{2}a{g}^{\prime }\left(X_{s-}\right)\mathrm{d}{W}_s}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\int }_0^t{a}^2{g}^{″}\left(X_{s-}\right)\mathrm{d}s}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\int }_0^t{\displaystyle {\int }_{\left|y\right|\ge 1}\left[g\left(X_{s-}+y\right)-g\left(X_{s-}\right)\right]N\left(\mathrm{d}s,\mathrm{d}y\right)}}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\int }_0^t{\displaystyle {\int }_{\left|y\right|<1}\left[g\left(X_{s-}+y\right)-g\left(X_{s-}\right)\right]\tilde{N}\left(\mathrm{d}s,\mathrm{d}y\right)}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\int }_0^t{\displaystyle {\int }_{\left|y\right|<1}\left[g\left(X_{s-}+y\right)-g\left(X_{s-}\right)-y{g}^{\prime }\left(X_{s-}\right)\right]\nu \left(\mathrm{d}y\right)\mathrm{d}s}}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{{\rm M}}_t+{\displaystyle {\int }_0^{t}Ag\left(X_{s-}\right)}\mathrm{d}s\hfill \end{array}$$

(31)

$$Ag\left(x\right)=a^2{g}^{″}\left(x\right)+{\displaystyle {\int }_{ℝ}\left[g\left(x+y\right)-g\left(x\right)-y{1}_{\left|y\right|<1}{g}^{\prime }\left(x\right)\right]\nu \left(\mathrm{d}y\right)}$$

(32)

$$\begin{array}{l}\because -\Delta =-{\displaystyle \frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}}\hfill \\ \therefore -a^2\left(-\Delta \right)g\left(x\right)=a^2{g}^{″}\left(x\right)\hfill \end{array}$$

(33)

$$\begin{array}{l}\because {\rm I}g\left(x\right)={\displaystyle {\int }_{ℝ}\left[g\left(x+y\right)-g\left(x\right)-y{1}_{\left|y\right|<1}{g}^{\prime }\left(x\right)\right]\frac{\mathrm{d}y}{{\left|y\right|}^{1+\beta }}}\hfill \\ \therefore F\left\{{\rm I}g\right\}\left(k\right)=\widehat{g}\left(k\right){\displaystyle {\int }_{ℝ}\frac{{\mathrm{e}}^{\mathrm{i}ky}-1-\mathrm{i}ky{1}_{\left|y\right|<1}}{{\left|y\right|}^{1+\beta }}\mathrm{d}y}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\widehat{g}\left(k\right){\displaystyle {\int }_{ℝ}\frac{\cos \left(ky\right)-1}{{\left|y\right|}^{1+\beta }}\mathrm{d}y}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=-{\left|k\right|}^{\beta }\widehat{g}\left(k\right){\displaystyle {\int }_{ℝ}\frac{1-\cos w}{{\left|w\right|}^{1+\beta }}\mathrm{d}w}\hfill \\ \because F\left\{{\left(-\Delta \right)}^{{\displaystyle \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}g\right\}\left(k\right)={\left|k\right|}^{\beta }\widehat{g}\left(k\right)\hfill \\ \therefore {\rm I}g=\left(-{\displaystyle {\int }_{ℝ}\frac{1-\cos w}{{\left|w\right|}^{1+\beta }}\mathrm{d}w}\right){\left(-\Delta \right)}^{{\displaystyle \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}g\hfill \\ \therefore {\displaystyle {\int }_{ℝ}\left[g\left(x+y\right)-g\left(x\right)-y{1}_{\left|y\right|<1}{g}^{\prime }\left(x\right)\right]\nu \left(\mathrm{d}y\right)}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{b^2}{{\displaystyle {\int }_{ℝ}\frac{1-\cos w}{{\left|w\right|}^{1+\beta }}\mathrm{d}w}}}{\displaystyle {\int }_{ℝ}\left[g\left(x+y\right)-g\left(x\right)-y{1}_{\left|y\right|<1}{g}^{\prime }\left(x\right)\right]\frac{\mathrm{d}y}{{\left|y\right|}^{1+\beta }}}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{b^2}{{\displaystyle {\int }_{ℝ}\frac{1-\cos w}{{\left|w\right|}^{1+\beta }}\mathrm{d}w}}}{\rm I}g\left(x\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{b^2}{{\displaystyle {\int }_{ℝ}\frac{1-\cos w}{{\left|w\right|}^{1+\beta }}\mathrm{d}w}}}\left(-{\displaystyle {\int }_{ℝ}\frac{1-\cos w}{{\left|w\right|}^{1+\beta }}\mathrm{d}w}\right){\left(-\Delta \right)}^{{\displaystyle \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}g\left(x\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=-b^2{\left(-\Delta \right)}^{{\displaystyle \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}g\left(x\right)\hfill \end{array}$$

(34)

$$Ag\left(x\right)=-a^2\left(-\Delta \right)g\left(x\right)-b^2{\left(-\Delta \right)}^{{\displaystyle \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}g\left(x\right)$$

(35)

where B denotes Brownian motion; N−L indicates Lévy jumps; N is the Poisson random measure; E[·] represents the intensity measure or characteristic function; ⸟ denotes the compensated measure; X^c is the continuous part; Μ_t is a local martingale (the sum of the Brownian-motion and compensated Poisson integral); 1_{·} is the indicator function; y1 is the truncation function; Ι is the integral operator; and F{·} is the Fourier transform.

3.2. Asymptotic Behaviors in Different Spatiotemporal Regimes

The generator defined by Eqs. 20−22 corresponds to the evolution equation given in Eq. 36. It can also be derived naturally within the first-law framework from constitutive relation Eq. 37. In Eq. 37, Ι^1−β denotes the (1−β)-order integral representing long-range correlation, as defined in Eq. 38.

For simplicity, we consider the fundamental solution in one dimension, i.e., the Green’s function for the Cauchy (initial-value) problem (Eq. 39). In Fourier space, Eq. 39 transforms into Eq. 40, whose exact solution is given by Eq. 41. Taking the inverse Fourier transform of Eq. 41, we confirm that the Green’s function for Eq. 39 is the convolution of two independent propagators (Eq. 42): the Gaussian kernel K (the fundamental solution of standard Brownian motion) and the fractional kernel G_β (the fundamental solution of pure fractional diffusion). This verifies that the macroscopic evolution law converges to an intermediate scaling regime between (a²t)^1/2 and (b²t)^1/β under competition between the dual mechanisms.

$$\left\{\begin{array}{l}\left({\partial }_t-A\right)\varphi \left(x,t\right)=f\left(x,t\right)\hfill \\ \varphi \left(x,0\right)=0\hfill \\ \underset{\left|x\right|\to \infty }{\lim }\varphi \left(x,t\right)=0\hfill \end{array}\right.$$

(36)

$$\mathrm{q}\propto -\nabla \left(A+B{I}^{1-\beta }\right)\varphi \left(x,t\right)$$

(37)

$$\begin{array}{c}{I}^{1-\beta }\varphi \left(x\right)={\displaystyle {\int }_{{ℝ}^n}\varphi \left(y\right)k\left(x-y\right)\mathrm{d}y}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{1}{\mathsf{\Gamma }\left(1-\beta \right)}}{\displaystyle {\int }_{{ℝ}^n}\frac{\varphi \left(y\right)}{{\left(y-x\right)}^{\beta }}\mathrm{d}y}\hfill \end{array}$$

(38)

$$\left\{\begin{array}{l}{\partial }_t{G}_{2,\beta }=A{G}_{2,\beta }\hfill \\ G_{2,\beta }\left(x,0\right)=\mathsf{\delta }\left(x\right)\hfill \end{array}\right.$$

(39)

$$\left\{\begin{array}{l}{\partial }_t{\widehat{G}}_{2,\beta }=-\lambda {\widehat{G}}_{2,\beta }\hfill \\ \lambda \left(k\right)=a^2{\left|k\right|}^2+b^2{\left|k\right|}^{\beta }\hfill \\ {\widehat{G}}_{2,\beta }\left(k,0\right)=1\hfill \end{array}\right.$$

(40)

$${\widehat{G}}_{2,\beta }\left(k,t\right)={\mathrm{e}}^{-\lambda t}={\mathrm{e}}^{-a^2{\left|k\right|}^{2}t}{\mathrm{e}}^{-b^2{\left|k\right|}^{\beta }t}$$

(41)

$$\left\{\begin{array}{l}{G}_{2,\beta }\left(x,t\right)=K\left(x,t\right)\ast G_{\beta }\left(x,t\right)\hfill \\ K\left(x,t\right)=F^{-1}\left\{{\mathrm{e}}^{-a^2{\left|k\right|}^{2}t}\right\}={\displaystyle \frac{1}{\sqrt{\mathsf{\pi }}}}{\displaystyle \frac{1}{2{\left(a^{2}t\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}}{\mathrm{e}}^{-{\displaystyle \frac{{\left|x\right|}^2}{4a^{2}t}}}\hfill \\ G_{\beta }\left(x,t\right)=F^{-1}\left\{{\mathrm{e}}^{-b^2{\left|k\right|}^{\beta }t}\right\}={\displaystyle \frac{1}{{\left(b^{2}t\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}}{S}_{\beta }\left({\displaystyle \frac{\left|x\right|}{{\left(b^{2}t\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}}\right)\hfill \end{array}\right.$$

(42)

$$\left\{\begin{array}{l}{S}_{\beta }\left(z\right)={\displaystyle \frac{1}{\mathsf{\pi }}}{\displaystyle {\int }_0^{+\infty }{\mathrm{e}}^{-u^{\beta }}\cos \left(zu\right)\mathrm{d}u}\hfill \\ z={\displaystyle \frac{\left|x\right|}{{\left(b^{2}t\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}}\hfill \\ u=k{\left(b^{2}t\right)}^{{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}\hfill \end{array}\right.$$

(43)

$$\begin{array}{c}{G}_{\beta }\left(x,t\right)={\displaystyle \frac{1}{{\left(b^{2}t\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}}{S}_{\beta }\left(z\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{1}{\beta {\left(b^{2}t\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}}{W}_{-{\displaystyle \frac{\beta }{2}},1-{\displaystyle \frac{\beta }{2}}}\left(-\left|z\right|\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{1}{\beta {\left(b^{2}t\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}}{\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-\left|z\right|\right)}^n}{n\text{!}\mathsf{\Gamma }\left(-\frac{\beta }{2}n-\frac{\beta }{2}+1\right)}}\hfill \end{array}$$

(44)

where k(·) is the power-law kernel; and Γ(·) is the Gamma function. F⁻¹{·} denotes the inverse Fourier transform; z is the second similarity variable; W_{·}(·) is the Wright function; and ! denotes the factorial operation.

The asymptotic limits across different polarization regimes reveal distinctive dynamical behaviors (starting from Eq. 45; derivations are detailed in Appendix A3):

$$\begin{array}{c}{G}_{2,\beta }\left(x,t\right)=F^{-1}\left\{{\mathrm{e}}^{-a^2{\left|k\right|}^{2}t}{\mathrm{e}}^{-b^2{\left|k\right|}^{\beta }t}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{1}{2\mathsf{\pi }}}{\displaystyle {\int }_{-\infty }^{+\infty }{\mathrm{e}}^{-\left(a^2{\left|k\right|}^2+b^2{\left|k\right|}^{\beta }\right)t}{\mathrm{e}}^{\mathrm{i}kx}\mathrm{d}k}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{1}{2\mathsf{\pi }}}{\displaystyle {\int }_{-\infty }^{+\infty }{\mathrm{e}}^{-\left(a^2{\left|k\right|}^2+b^2{\left|k\right|}^{\beta }\right)t}\left\{\cos \left(kx\right)+\mathrm{isin}\left(kx\right)\right\}\mathrm{d}k}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{1}{\mathsf{\pi }}}{\displaystyle {\int }_0^{+\infty }{\mathrm{e}}^{-\left(a^2{k}^2+b^2{k}^{\beta }\right)t}\cos \left(kx\right)\mathrm{d}k}\hfill \end{array}$$

(45)

(i) Short-time near-field: behavior is dominated by local Gaussian diffusion; long-range correlations contribute only as a minor correction (Eq. 46).

$$G_{2,\beta }\left(x,t\right)\sim {\displaystyle \frac{1}{\mathsf{\pi }}}{\displaystyle \frac{1}{2\sqrt{a^{2}t}}}{\mathrm{e}}^{-{\displaystyle \frac{{\left|x\right|}^2}{4a^{2}t}}}\left(1-{\displaystyle \frac{b^2}{a^{\beta }}}{\displaystyle \frac{1}{\sqrt{\mathsf{\pi }}}}\mathsf{\Gamma }\left({\displaystyle \frac{1+\beta }{2}}\right)t^{1-{\displaystyle \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}+\dots \right)$$

(46)

(ii) Intermediate-time transitional field: the Gaussian kernel and Lévy-type fractional operator exert comparable influence, shaping a relatively smooth spatial gradient (Eq. 47).

$$\begin{array}{c}{G}_{2,\beta }\left(x,t\right)={\displaystyle \frac{1}{\mathsf{\pi }}}{I}_0\left(t\right)+{\displaystyle \frac{x^2}{2\mathsf{\pi }}}{I}_2\left(t\right)+O\left(x^4\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{1}{\mathsf{\pi }}}{\displaystyle {\int }_0^{+\infty }{\mathrm{e}}^{-\left(a^2{k}^2+b^2{k}^{\beta }\right)t}\mathrm{d}k}+{\displaystyle \frac{x^2}{2\mathsf{\pi }}}{\displaystyle {\int }_0^{+\infty }{k}^2{\mathrm{e}}^{-\left(a^2{k}^2+b^2{k}^{\beta }\right)t}\mathrm{d}k}+O\left(x^4\right)\hfill \end{array}$$

(47)

(iii) Intermediate-time far-field and ultimate long-time regime: sub-diffusion dominates, extending from far field to global domain, and the system exhibits pronounced power-law decay (Eqs. 48 and 49). Due to the vanishingly small β, significant memory solidification emerges, causing the thermal field, even after an extremely long duration of repeated heat injection/extraction cycles, to remain implicitly pinned to an approximately constant dissipative plateau (Eq. 50).

(iv) In the extreme limit of very long times and spatial scales far exceeding the fractional scale—though practically unattainable because β is too small—the temperature variation at vast distances retains a power-law form, but its amplitude varies linearly with time (Eq. 51).

$$G_{2,\beta }\left(x,t\right)\sim {\displaystyle \frac{b^{2}t}{\mathsf{\pi }}}\mathsf{\Gamma }\left(1+\beta \right)\sin \left({\displaystyle \frac{\mathsf{\pi }\beta }{2}}\right){\displaystyle \frac{1}{{\left|x\right|}^{1+\beta }}}$$

(48)

$$G_{2,\beta }\left(x,t\right)={\displaystyle \frac{1}{\mathsf{\pi }}}{\displaystyle \frac{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.\mathsf{\Gamma }\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.\right)}{{\left(b^{2}t\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}}$$

(49)

$$G_{2,\beta }\left(x,t\right)\sim G_{\beta }\left(0,t\right)\sim {\displaystyle \frac{1}{\mathsf{\pi }}}{\displaystyle \frac{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.\mathsf{\Gamma }\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.\right)}{{\left(b^{2}t\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}}$$

(50)

$$G_{2,\beta }\left(x,t\right)\sim {\displaystyle \frac{b^{2}t}{\mathsf{\pi }}}\mathsf{\Gamma }\left(1+\beta \right)\sin \left({\displaystyle \frac{\mathsf{\pi }\beta }{2}}\right){\displaystyle \frac{1}{{\left|x\right|}^{1+\beta }}}$$

(51)

where O(·) denotes terms of the same order of infinitesimal; and ο(·) denotes higher-order infinitesimals.

Accordingly, once parameters a², b², and β have been fixed by gradient PCM design—thereby determining the topological structure and dynamical “genetic code” of the thermal field and its conjugate information field—the optimal shaping and modulation of the temporal trajectory of external heat source f become the primary means to extract the system’s maximum possible performance.

3.3. Path-Integral Optimal-Control Framework for Cyclic Heat Recovery

Once the hybrid evolutionary generator has been physically specified, with the dissipative-geometry parameters (a², b², β) determined by the structural action of gradient PCMs, the remaining task is to identify the operating mode of the active heat source f that optimizes cyclic heat-recovery performance. In this setting, the control problem is formulated under the state equation of the thermal field and embedded into the path-integral variational framework established above.

We consider a periodic operating process consisting of heat injection, insulation, heat extraction, and recovery. Let the time interval [0, Τ] denote one complete cycle. The different operating stages and spatial regions are distinguished by indicator functions, as defined in Eq. (52). On this basis, the total injected heat, the total extracted heat, and the heat-recovery efficiency over one cycle are defined in Eq. (53). To characterize the long-time operating history, a sliding-window averaged efficiency can also be introduced. For the derivation of the control framework, however, we first consider a simplified multi-cycle setting in which the injected heat in each cycle is prescribed as a fixed value Eⁱⁿ. The optimization objective is then to maximize the recovered heat, or equivalently the heat-recovery efficiency, under the constraint of the governing evolution equation and the prescribed injection-energy budget.

$$\begin{array}{l}\begin{array}{c}\mathrm{Heat}\mathrm{injection}\mathrm{stage}\text{:}{\chi }^{\mathrm{in}}\left(t\right)=\left\{\begin{array}{l}\begin{array}{cc}1,\hfill & t\in {{\rm T}}^{\mathrm{in}}\hfill \end{array}\hfill \\ \begin{array}{cc}0,\hfill & t\notin {{\rm T}}^{\mathrm{in}}\hfill \end{array}\hfill \end{array}\right.\hfill \end{array}\hfill \\ \begin{array}{c}\mathrm{Heat}\mathrm{extraction}\mathrm{stage}\text{:}{\chi }^{\mathrm{out}}\left(t\right)=\left\{\begin{array}{l}\begin{array}{cc}1,\hfill & t\in {{\rm T}}^{\mathrm{out}}\hfill \end{array}\hfill \\ \begin{array}{cc}0,\hfill & t\notin {{\rm T}}^{\mathrm{out}}\hfill \end{array}\hfill \end{array}\right.\hfill \end{array}\hfill \\ \begin{array}{cccc}\mathrm{Heat}\hfill & \mathrm{injection}\hfill & \mathrm{region}\text{:}\hfill & {\mathrm{R}}^{\mathrm{in}}\in ℝ\hfill \end{array}\hfill \\ \begin{array}{cccc}\mathrm{Heat}\hfill & \mathrm{extraction}\hfill & \mathrm{region}\text{:}\hfill & {\mathrm{R}}^{\mathrm{out}}\in ℝ\hfill \end{array}\hfill \end{array}$$

(52)

$$\left\{\begin{array}{l}{E}^{\mathrm{in}}={\displaystyle {\int }_0^{{\rm T}}{\chi }^{\mathrm{in}}\left(t\right)\mathrm{d}t{\displaystyle {\int }_{{\mathrm{R}}^{\mathrm{in}}}f\left(x,t\right)\mathrm{d}x}}\hfill \\ E^{\mathrm{out}}=-{\displaystyle {\int }_0^{{\rm T}}{\chi }^{\mathrm{out}}\left(t\right)\mathrm{d}t{\displaystyle {\int }_{{\mathrm{R}}^{\mathrm{out}}}f\left(x,t\right)\mathrm{d}x}}\hfill \\ \eta ={\displaystyle \raisebox{1ex}{$E^{\mathrm{out}}$}\!\left/ \!\raisebox{-1ex}{$E^{\mathrm{in}}$}\right.}\hfill \end{array}\right.$$

(53)

where χ and R denote indicator functions; E denotes cumulated heat injection/extraction; η is the heat-recovery efficiency.

To express this problem in a mathematically consistent form, the indicator functions are normalized as shown in Eq. (54), and the corresponding total injected and extracted heat are written in Eq. (55). To avoid unbounded control amplitudes and to account for the actuation cost of the external heat source, a quadratic penalty term is introduced. The resulting objective functional is given in Eq. (56), where ϑ > 0 is the weighting parameter balancing heat recovery and control expenditure. The total injected heat is imposed as an additional equality constraint.

$$\left\{\begin{array}{l}{\upsilon }^{\mathrm{in}}\left(x,t\right)={\chi }^{\mathrm{in}}\left(t\right)\cdot {\mathrm{R}}^{\mathrm{in}}\left(x\right)\hfill \\ {\upsilon }^{\mathrm{out}}\left(x,t\right)={\chi }^{\mathrm{out}}\left(t\right)\cdot {\mathrm{R}}^{\mathrm{out}}\left(x\right)\hfill \end{array}\right.$$

(54)

$$\left\{\begin{array}{l}{E}^{\mathrm{in}}={\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}{\upsilon }^{\mathrm{in}}f\left(x,t\right)\mathrm{d}x}}\hfill \\ E^{\mathrm{out}}=-{\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}{\upsilon }^{\mathrm{out}}f\left(x,t\right)\mathrm{d}x}}\hfill \end{array}\right.$$

(55)

$$\begin{array}{cc}J\left[f\right]=-{\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}{\upsilon }^{\mathrm{out}}f\left(x,t\right)\mathrm{d}x}}-{\displaystyle \frac{\vartheta }{2}}{\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}{f}^2\left(x,t\right)\mathrm{d}x}},\hfill & \vartheta >0\hfill \end{array}$$

(56)

where υ denotes normalized indicator function; J is the objective functional; ϑ is the weighting parameter.

The constrained optimization problem is incorporated into the path-integral framework by introducing the Lagrange multiplier field ψ, which here plays the role of the thermodynamic conjugate information field, together with a scalar multiplier ρ associated with the injection-energy constraint. The augmented Lagrangian functional ℘ is therefore defined by Eq. (57). In this formulation, the optimization is performed over the coupled variables (ϕ, ψ, f), while ρ enforces the cycle-wise energy-input condition.

The optimality system follows from the first variation of ℘. Setting all admissible variations to zero yields a closed forward-backward system. Variation with respect to ψ recovers the state equation for the thermal field ϕ, given in Eq. (58), which governs the forward evolution of the system under the prescribed control. Variation with respect to ϕ yields the adjoint equation for the information field ψ, shown in Eq. (59), together with the corresponding terminal condition. Variation with respect to f gives the pointwise stationarity condition for the optimal control, Eq. (60), and variation with respect to ρ restores the injection-energy constraint, Eq. (61). Eqs. (58)–(61) thus form the complete optimality system for cyclic heat recovery.

$$\begin{array}{c}\wp \left[\varphi ,f,\psi ,\rho \right]=J\left[f\right]\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\rho \left(E^{\mathrm{in}}-{\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}{\upsilon }^{\mathrm{in}}f\left(x,t\right)\mathrm{d}x}}\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\psi \left({\partial }_t\varphi -A\varphi -f\right)\left(x,t\right)\mathrm{d}x}}\hfill \end{array}$$

(57)

$$\begin{array}{c}\because {\mathsf{\delta }}_{\psi }\wp =\left({\partial }_t-A\right)\varphi -f=0\hfill \\ \therefore \begin{array}{cc}\left({\partial }_t-A\right)\varphi \left(x,t\right)=f\left(x,t\right),\hfill & \varphi \left(x,0\right)=0\hfill \end{array}\hfill \end{array}$$

(58)

$$\begin{array}{l}\because {\mathsf{\delta }}_{\varphi }\wp ={\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\psi \left({\partial }_t\left(\mathsf{\delta }\varphi \right)-A\left(\mathsf{\delta }\varphi \right)\right)\left(x,t\right)\mathrm{d}x}}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\mathsf{\delta }\varphi \left(-{\partial }_t-A\right)\psi \left(x,t\right)\mathrm{d}x}}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\begin{array}{cc}0,\hfill & \forall \mathsf{\delta }\varphi \hfill \end{array}\hfill \\ \therefore \begin{array}{cc}\left(-{\partial }_t-A\right)\psi \left(x,t\right)=0,\hfill & \psi \left(x,{\rm T}\right)=0\hfill \end{array}\hfill \end{array}$$

(59)

$$\begin{array}{l}\because {\mathsf{\delta }}_f\wp =-{\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x{\upsilon }^{\mathrm{out}}\mathsf{\delta }f\left(x,t\right)}}-\vartheta {\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\mathrm{d}xf\mathsf{\delta }f\left(x,t\right)}}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\rho \left(-{\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x{\upsilon }^{\mathrm{in}}\mathsf{\delta }f\left(x,t\right)}}\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x\psi \left(-\mathsf{\delta }f\right)\left(x,t\right)}}\hfill \\ \because {\mathsf{\delta }}_f\wp =\begin{array}{cc}0,\hfill & \forall \mathsf{\delta }f\hfill \end{array}\hfill \\ \therefore \left(-{\upsilon }^{\mathrm{out}}-\vartheta f-\rho {\upsilon }^{\mathrm{in}}-\psi \right)\left(x,t\right)=0\hfill \\ \therefore f\left(x,t\right)=-{\displaystyle \frac{1}{\vartheta }}\left({\upsilon }^{\mathrm{out}}+\rho {\upsilon }^{\mathrm{in}}+\psi \right)\left(x,t\right)\hfill \end{array}$$

(60)

$$\begin{array}{c}\because {\mathsf{\delta }}_{\rho }\wp =E^{\mathrm{in}}-{\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}{\upsilon }^{\mathrm{in}}f\left(x,t\right)\mathrm{d}x}}=0\hfill \\ \therefore E^{\mathrm{in}}={\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}{\upsilon }^{\mathrm{in}}f\left(x,t\right)\mathrm{d}x}}\hfill \end{array}$$

(61)

where ℘ is the augmented Lagrangian functional; ρ denotes scalar multiplier.

For numerical implementation, the above system can be solved by a forward-backward sweep procedure. First, an initial control f [0] satisfying Eq. (61) is prescribed; in practice, a nonzero constant distribution can be used. Second, for a given iterate f^[n], the state equation (58) is solved forward in time on [0, Τ] to obtain ϕ^[n], and the adjoint equation (59) is then solved backward in time to obtain ψ^[n]. Substituting ϕ^[n] and ψ^[n] into the stationarity condition (60) gives the updated control expression f^new. Enforcing the constraint (61) then determines the scalar multiplier ρ, whose explicit form is given in Eq. (62). The denominator in Eq. (62) is the space-time measure of the injection domain over one cycle. The control is subsequently updated according to Eq. (63). Third, convergence is checked using the criterion in Eq. (64). If the stopping condition is satisfied, the iteration terminates; otherwise, the relaxation step in Eq. (65) is performed, and the procedure returns to the forward-backward update.

$$\rho =-{\displaystyle \frac{\vartheta E^{\mathrm{in}}+{\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x{\upsilon }^{\mathrm{in}}\left({\upsilon }^{\mathrm{out}}+{\psi }^{\left[n\right]}\right)\left(x,t\right)}}}{{\displaystyle {\int }_0^{{\rm T}}\mathrm{d}t{\displaystyle {\int }_{{ℝ}^n}\mathrm{d}x{\upsilon }^{\mathrm{in}}{\upsilon }^{\mathrm{in}}\left(x,t\right)}}}}$$

(62)

$$f^{\left[n+1\right]}\left(x,t\right)=f^{\mathrm{new}}\left(x,t\right)=-{\displaystyle \frac{1}{\vartheta }}\left({\upsilon }^{\mathrm{out}}+\rho {\upsilon }^{\mathrm{in}}+{\psi }^{\left[n\right]}\right)\left(x,t\right)$$

(63)

$${‖f^{\left[n+1\right]}-f^{\left[n\right]}‖}_{L^2}\le \epsilon$$

(64)

$$n\leftarrow n+1$$

(65)

where ε denotes convergent criterion.

This implementation translates the path-integral optimal-control problem into an iterative solution of coupled state and adjoint field equations under a global cycle-wise energy constraint. More importantly, it preserves the variational structure established in the preceding subsections: the thermal field ϕ describes the physical evolution of the subsurface heat system, the conjugate field ψ represents the corresponding information response, and the control f is determined through their stationarity coupling rather than by an externally imposed heuristic rule. The resulting framework therefore provides a mathematically consistent route for optimizing periodic heat-source operation in underground heat-recovery systems.

For clarity, the above derivation is presented under a simplified fixed-input assumption at the cycle level. In practical applications, the same framework can be extended to include cycle-dependent operating constraints, engineering bounds on the control, and site-specific parameter calibration.

4. Summary and Conclusions

This study derives a hybrid evolutionary generator compatible with both Gaussian diffusion and fractional long-range correlations to characterize medium-to-long-term underground heat recovery. The evolution equation, dictated by this generator and the natural time scale, is fundamentally equivalent to the Euler−Lagrange equation of the Onsager variational principle. The Lagrange density under this architecture highlights the dominant role of the information field: decoding thermal-field information to optimize generator parameters is physically tantamount to synergistically orchestrating the transfer rates and pathways of both thermal energy and its information flow. This framework provides a theoretical basis for achieving intrinsic sustainability in long-term heat recovery. The path-integral control formulation further provides an operational route for maximizing cycle-wise heat-recovery efficiency under prescribed energy-input constraints.

The core findings of this study are summarized as follows:

1. The evolutionary generator, coupled with the natural time scale, defines the underlying spatiotemporal geometry. Dissipation is an intrinsic property of this geometry; it is the natural consequence of the system following the free-energy functional gradient and represents a fundamental systemic attribute that external thermal drivers cannot override.

2. Sustainable underground thermal evolution requires that dissipative geometry possess specific power-law tails. This implies that a balance must be achieved between local Gaussian diffusion and non-local long-range correlation mechanisms.

3. This balance emerges when the action functional exhibits scale invariance, reflecting the conjugate nature between the thermal dissipative system and its information transfer field.

4. Achieving engineering sustainability requires synergistic regulation of thermal energy and its conjugate information flow, forcing the system to approach this critical scale-invariant boundary.

5. The aforementioned regulation—the process of optimizing the generator’s core parameters—is realized in engineering through strategic design and spatial deployment of gradient phase change materials. Consequently, the role of gradient PCMs is reinterpreted: they are no longer only “sensible/latent heat containers,” but serve as physical media intervening directly in information-transport processes.

In conclusion, this study identifies the structural root cause of classical local heat equations failing in cross-cycle predictions. It clarifies the physical basis of “macroscopic sustainability” and illuminates the conjugate relationship between thermal energy and information transfer, providing a pathway to enhanced sustainable efficiency in geothermal engineering.