Jean-Claude Krapez
email: jean-claude.krapez@onera.fr

Influence of thermal hysteresis on the heat shuttling effect: the case of VO2

Thermotronics has attracted much attention driven by the promising potentials offered by devices such as thermal diodes, thermal transistors, and thermal memristors. Heat shuttling (or heat ratcheting, or heat pumping) is a phenomenon exhibited by nonlinear materials presenting temperature-dependent thermal conductivity which, when sandwiched between two thermal baths with one bath subjected to a time-varying temperature, show non vanishing net heat flow, although the baths share the same average temperature. Phase-change materials (PCMs) like VO 2 were recently taken for illustration due to a strong change in conductivity over a small temperature range; energy extraction from the thermal variations of the environment was envisioned thereupon. However, up to now, the impact of PCM hysteresis has been either overlooked or roughly approximated. On the basis of a thermal model simulating partial hysteresis loops and non-hysteretic branches, we demonstrate that the presence of hysteresis profoundly modifies the appearance of the heat shuttling effect and can constitute a hindrance to its manifestation. Operating configurations to improve its observation have been proposed.

I. INTRODUCTION

Controlled manipulation of heat flow is of high importance for the development of a wide variety of technologies. Nonlinear heat transfer is essential to tailor heat currents; for this purpose, one can leverage the variations of thermal conductivity or emissivity with temperature. This thermal dependence generates an asymmetric system response that has been used to manipulate heat flows in the same manner as the flow of electrons is controlled in electronic circuits. Heat flow can be carried by phonons (phononics), photons (photonics), or electonic conduction; the word thermotronics has been coined for this field of research devoted to the development of materials and devices for the rectification, amplification, and modulation of heat currents through either photons, electrons, or phonons or a combination thereof. This gave rise to applications both at nanoscale 1 and at macroscopic scale [START_REF] Dai | Designing nonlinear thermal devices and metamaterials under the Fourier law: A route to nonlinear thermotics[END_REF] , among them thermal diodes [START_REF] Li | Thermal diode: Rectification of heat flux[END_REF][START_REF] Terraneo | Controlling the energy flow in nonlinear lattices: a model for a thermal rectifier[END_REF][START_REF] Ordonez-Miranda | Conductive thermal diode based on the thermal hysteresis of VO2 and nitinol[END_REF][START_REF] Schmotz | A thermal diode using phonon rectification[END_REF] , thermal transistors [START_REF] Li | Colloquium: Phononics: Manipulating heat flow with electronic analogs and beyond[END_REF][START_REF] Li | Negative differential thermal resistance and thermal transistor[END_REF][START_REF] Ordonez-Miranda | Transistorlike device for heating and cooling based on the thermal hysteresis of VO2[END_REF][START_REF] Latella | Dynamical response of a radiative thermal transistor based on suspended insulator-metal-transition membranes[END_REF] , thermal logic gates [START_REF] Wang | Thermal logic gates: computation with phonons[END_REF] , thermal memories [START_REF] Wang | Thermal memory: a storage of phononic information[END_REF] , and thermal memristors [START_REF] Ben-Abdallah | Thermal memristor and neuromorphic networks for manipulating heat flow[END_REF][START_REF] Yang | Theoretical framework of the thermal memristor via a solid-state phase change material[END_REF] .

When submitted to temporal modulation, devices supporting nonlinear heat flow may give rise to novel and intriguing phenomena, among them heat shuttling, a net heat current that shows up even in the absence of a mean temperature gradient. The shuttling of heat (or heat ratcheting, or heat pumping) has been discussed a few years ago in the context of heat conduction in nonlinear lattices [START_REF] Li | Ratcheting heat flux against a thermal bias[END_REF][START_REF] Li | Shuttling heat across one-dimensional homogenous nonlinear lattices with a Brownian heat motor[END_REF][START_REF] Ren | Emergence and control of heat current from strict zero thermal bias[END_REF] . A salient requirement is the presence of an induced dynamical symmetry-breaking mechanism in conjunction with nonlinearity [START_REF] Li | Colloquium: Phononics: Manipulating heat flow with electronic analogs and beyond[END_REF] . The necessary symmetry-breaking can be obtained by using two coupled lattices with different thermal properties [START_REF] Li | Ratcheting heat flux against a thermal bias[END_REF][START_REF] Ren | Emergence and control of heat current from strict zero thermal bias[END_REF] or by exploiting the nonlinear response induced by the harmonic mixing mechanism stemming from a time varying two-mode modulation of the bath temperature [START_REF] Li | Shuttling heat across one-dimensional homogenous nonlinear lattices with a Brownian heat motor[END_REF] . Later on, the shuttling effect has been demonstrated for the radiative heat flux exchanged between two bodies separated by a thin vacuum gap [START_REF] Latella | Radiative heat shuttling[END_REF] .

Recently, heat shuttling obtained by heat conduction at the macroscopic scale in phase change materials (PCMs) has been discussed [START_REF] Liu | Energy harvesting from thermal variation with phasechange materials[END_REF][START_REF] Ordonez-Miranda | Net heat current at zero mean temperature gradient[END_REF] . The heat shuttling phenomenon was explained by the existence of a temperature-dependent thermal conductivity: when conductivity increases with temperature, a net heat conduction flow is observed from the bath whose temperature is symmetrically modulated to the static bath, whereas the opposite is observed when conductivity decreases with temperature [START_REF] Ordonez-Miranda | Net heat current at zero mean temperature gradient[END_REF] . The theoretical analysis performed by Ordonez et al. [START_REF] Ordonez-Miranda | Net heat current at zero mean temperature gradient[END_REF] was illustrated by taking as an example the properties of VO 2 , a PCM showing a strong increase of conductivity with temperature. The quantification of the heat shuttling effect was nevertheless made in the static (pseudo-modulated) regime while ignoring the hysteretic behavior of the material. The phase change hysteresis causes that the transformation from one phase to the other and the transformation in the opposite direction do not share the same temperature range. As a result, the curves, as a function of temperature, of a series of state variables and properties (e.g. enthalpy, apparent specific heat, thermal conductivity) are shifted with respect to temperature. The numerical modeling described by Liu et al. [START_REF] Liu | Energy harvesting from thermal variation with phasechange materials[END_REF] focused on nitinol and VO 2 while restricting to the thermallythick regime. In addition, the hysteresis regarding heat conductivity and specific heat was modeled by assuming an instantaneous switch from the heating curve to the cooling curve (and vice versa) after a change of the sign of the temperature evolution, which was recognized as a rough approximation [START_REF] Liu | Energy harvesting from thermal variation with phasechange materials[END_REF] .

In this paper we model the conductive heat shuttling effect in a PCM when considering all following phenomena : nonlinear conductivity, heat storage with a temperature-dependent heat capacity, or, equivalently, through nonlinear enthalpy, hysteresis regarding both conductivity and enthalpy. Incomplete phase change process will be modeled by applying to both conductivity and enthalpy the curve-switch model first described by Bony and Citherlet [START_REF] Bony | Numerical model and experimental validation of heat storage with phase change materials[END_REF] for enthalpy alone. Illustrations for the case of VO 2 will be provided. We show that hysteresis has a strong impact on the net heat current: heat shuttling may totally disappear if the amplitude of the temperature modulation is lower than a threshold depending on the hysteresis width. We describe the optimal configuration for maximizing the heat shuttling effect taking as an example the thermal properties of VO 2 . A parametric analysis regarding the modulation frequency will provide a description of the transitory behavior between the thermally thick and thermally thin regimes. The paper ends with an opening-up on a two-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0147225

bath modulation offering a perspective for improving the heat shuttling effect.

II. THERMAL MODEL

A. Heat shuttling through temperature modulation of one thermal bath

Let us consider a system comprising a nonlinear material of thickness l in thermal contact with two reservoirs (see Fig. 1). One of them (say the one at right) is at constant temperature T c , the other one is at a (sinusoidally) modulated temperature with frequency f , amplitude T a and a mean equal to T c :

T L (t) = T c + T a sin (2π f t) , (1)
T R (t) = T c , (2)
for t > 0. These conditions imply a zero mean gradient.

B. Thermal properties

As a nonlinear material we will consider VO 2 , which, in addition to exhibiting a strong thermal conductivity change in a narrow range of temperature, presents a significant hysteresis. These phenomena are related to a phase change, more specifically a metal-insulator transition (MIT) in which the insulating and metallic phases coexist over a finite temperature range [START_REF] Oh | Thermal conductivity and dynamic heat capacity across the metal-insulator transition in thin film VO2[END_REF] . At low temperature (T < 330 K) VO 2 behaves like a dielectric with low thermal conductivity (k i = 3.6 W m -1 K -1 ), while at high temperature (T > 345 K)

(a) (b)

FIG. 1. Scheme of a PCM in which modulated heat conduction, as prescribed by two surrounding thermal baths, gives rise to heat shuttling. The right bath is at constant temperature T c whereas the left bath supports a periodic temperature variation of amplitude T a between T c + T a (a) and T c -T a (b). The color gradation reflects the temperature variations T (x) while the white line describes the thermal conductivity profile k(x) induced by the nonlinearity and the hysteresis (the absence of biunivocality between the profiles T (x) and k(x) is explained by the hysteresis).

it exhibits a metallic phase with a higher thermal conductivity (k m = 6 W m -1 K -1 ). Actually, the transition region is narrower, about 4 K. Because of hysteresis, a shift in the temperature dependence of the VO 2 properties is observed and has been reported to be about 8 K [START_REF] Qazilbash | Mott transition in VO2 revealed by infrared spectroscopy and nano-imaging[END_REF][START_REF] Qazilbash | Infrared spectroscopy and nano-imaging of the insulator-to-metal transition in vanadium dioxide[END_REF] . The exact temperature bounds and width of the MIT actually depend on a variety of factors: observed changes in transition temperature were attributed to variation in the VO 2 microstructure, mainly the average grain sizes [START_REF] Radue | Effect of a substrate-induced microstructure on the optical properties of the insulator-metal transition temperature in VO2 thin films[END_REF] , to the variation of microcrystal dimensions [START_REF] Yoon | Controlling the temperature and speed of the phase transition of VO2 microcrystals[END_REF] , to the variation in the degree of crystallization by different rapid thermal annealing temperatures or annealing times [START_REF] Lee | Rta and stoichiometry effect on the thermochromism of VO2 thin films[END_REF] , and to substitutional doping of VO 2 , with, e.g. tungsten [START_REF] Rajeswaran | Effect of w addition on the electrical switching of VO2 thin films[END_REF][START_REF] Gomez-Heredia | Measurement of the hysteretic thermal properties of w-doped and undoped nanocrystalline powders of VO2[END_REF] or titanium [START_REF] Chen | Unraveling mechanism on reducing thermal hysteresis width of VO2 by Ti doping: a joint experimental and theoretical study[END_REF] . Works showing an influence of the temperature ramp rate on the transition temperature and hysteresis width are rare [START_REF] Mun | Metal insulator transition characteristics of macro-size single domain VO2 crystals[END_REF] . Yet, the MIT rate of VO 2 is very fast: recent work has shown metallization of VO 2 on just picosecond scales [START_REF] Muskens | Antenna-assisted picosecond control of nanoscale phase transition in vanadium dioxide[END_REF] .

Typical values of the characteristic parameters of the phase transition were evaluated in Ref. 32, leveraging upon experimental data obtained in Ref. 21-23. The model used in Ref. 32 is based on the Bruggeman's symmetric effectivemedium formula and explicit expressions derived for the temperature evolution of the volume fraction of the metallic and insulating domains appearing during the heating and cooling of VO 2 . These expressions involve the Lambert function and the complementary error function. An approximation allows expressing the volume fraction as a sigmoidal function involving a simple exponential function of the weighted temperature difference with the transition temperature [START_REF] Ordonez-Miranda | Conductive thermal diode based on the thermal hysteresis of VO2 and nitinol[END_REF][START_REF] Ordonez-Miranda | Transistorlike device for heating and cooling based on the thermal hysteresis of VO2[END_REF] . The technical drawback of these models is that the temperature range corresponding to a change of the volume fraction from 0 to 1 is theoretically unbounded. Restricting the MIT to a bounded temperature interval presents a clear advantage while introducing an insignificant error. As a consequence and according to Ref. 9, we take for T i and T m , i.e. the temperaturetransition ends in the insulator state, resp. metallic state, the values 341 K and 345 K for the heating process, and 333 K and 337 K for the cooling process. In Ref. 9 the functional form suggested for the volume fraction as a function of temperature is a polynomial of degree 5. This is a sigmoid-like function enabling smooth transition between two levels while preserving zero 1 st and 2 nd order derivatives at both edges [START_REF] Perlin | Improving noise[END_REF] . Here, as an alternative and for reasons explained later we model the volume fraction f of the metallic domains as a smooth step function of temperature given by:

f (T ) = 1 2 1 -cos(π T ) ; T = T -T i T m -T i , (3)
and the conductivity as:

k(T ) = k i + (k m -k i ) f (T ). (4)
Based on this model, the conductivity vs temperature curves are as reported in Fig. 2(a) for the heating (in red) and cooling (in blue) processes. The values of the specific heat in the insulating and metallic phases are quite close and don't vary much in the short temperature range of ∼10 K resp. below and above the MIT This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0147225

that we intend to explore [START_REF] Berglund | Electronic properties of VO2 near the semiconductor-metal transition[END_REF] . We thus considered a common and constant value of C = 0.71 kJ kg -1 K -1 . The enthalpy per unit mass of VO 2 (actually sum of sensible enthalpy and latent heat) can thus be written as:

H(T ) = C(T -T re f ) + L f (T ), (5)
where T re f is an arbitrary reference temperature (here taken as 330 K) and L = 51.45 kJ kg -1 is the latent heat of the transition [START_REF] Berglund | Electronic properties of VO2 near the semiconductor-metal transition[END_REF] . Thereby, the enthalpy values reached at the transition ends of the heating and cooling phases are, respectively, {H i , H m } = {7.81, 62.1} kJ kg -1 and {H i , H m } = {2.13, 56.42} kJ kg -1 . The evolution of enthalpy with the temperature according to the present model is described in Fig. 2(b). During the numerical process, it will be necessary to extract the temperature from the current enthalpy value H, which amounts to apply an inversion to Eq. ( 5). Three cases will be considered depending on how H compares to the boundaries H i and H m , whether in the heating or cooling phase:

T (H) =    T re f + H/C , H ≤ H i T i + (T m -T i ) g(H) , H i ≤ H ≤ H m T re f + (H -L) /C , H ≥ H m (6)
where g(H) is the inverse step function of f (T ) when combined with the linear transformation in Eq. (5). There is no closed-form analytical expression for g(H), nevertheless, the following approximation gives satisfactory results: with

g(H) = c 1 H + (1 -c 1 ) 1 + h(1 -H) h( H) -1 , (7)
h(x) = (1 + c 2 x) c 3 -1, (8)
the normalized enthalpy

H = (H -H i ) / (H m -H i ) , (9)
and the empirical coefficients c 1 = 0.211, c 2 = 5140, c 3 = 0.551. The error on temperature is then less than 1.3 • 10 -3 K and the RMS error is 5 • 10 -4 K. The values of the three coefficients c 1 , c 2 , c 3 only depend on the value taken by the ratio of the sensible heat to latent heat C(T m -T i )/L (which is the same for both heating and cooling processes). Even though the functional form adopted for f in Eq. ( 3) has zero derivatives at both edges only up to order one, it has the advantage that the inverse step function g as defined in Eqs. ( 7)-( 9) leads to an approximation error about 5 times smaller than with the 5 th -order polynomial smooth step function in Ref. 9 and 33.

C. Modeling the hysteresis loops

Thermal cycling with complete phase changes from fully insulating phase to fully metallic phase and vice versa generates a full hysteresis loop which is described by the succession of points [1-2-3-4-5-6-1] in Fig. 2. The question is then: what happens if heating stops before the right end of the heating branch (i.e. 345 K) or, symmetrically, if cooling stops before the left end of the cooling branch (i.e. 333 K) ? This was described, for the electric-resistance versus temperature hysteresis, through the so-called first-order reversal curves [START_REF] Ramírez | Firstorder reversal curve measurements of the metal-insulator transition in VO2: Signatures of persistent metallic domains[END_REF] . Incomplete hysteresis loops were also analysed for the electricresistance and the optical reflectance by Gurvitch et al. [START_REF] Gurvitch | Nonhysteretic behavior inside the hysteresis loop of VO2 and its possible application in infrared imaging[END_REF] . In particular, temperature excursions taken from an attachment point on any side of the major hysteresis loop produce minor loops and, for sufficiently small excursions, these minor loops flatten out, degenerating into nonhysteretic branches (NHBs) that are linear in log(ρ), or reflectance, versus T 36 . The slope of these NHBs evolves from the resistance (or reflectance coefficient) observed in the semiconducting phase to the one in the metallic phase [START_REF] Gurvitch | Nonhysteretic behavior inside the hysteresis loop of VO2 and its possible application in infrared imaging[END_REF] . We will capitalize on these observations regarding NHBs and transpose them qualitatively to thermal conductivity. In addition, to keep our model sufficiently simple, we will consider that NHBs apply to minor loops of large excursions as well (i.e. up to the major loop width). A few example will illustrate the possible trajectories of the representative point (RP) of VO 2 in the conductivity versus temperature map. We will admit that, if, after reaching a point along the heating branch as point 2 in Fig. 2(a) cooling starts, the conductivity stays constant while the RP moves along the horizontal (green) line until point 6 is reached. Then, if the cooling continues, the RP moves along the cooling branch towards point 1. Alternatively, cooling can stop along the green line at point 8 and heating resumes. Temperature and conductivity can then oscillate between any pair of points along the horizontal line 2-6.

The transition between points 2 and 6 is without phase change; the fraction of the metallic domains remains constant. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0147225

Hence, enthalpy evolves according to the specific heat of the present semiconductor/metal mixture. Since we have assumed that the specific heat is the same for both phases in the narrow temperature range considered, the line joining points 2 and 6 (or 3 and 5) in Fig. 2(b) is a straight line with the same slope as that of the enthalpy curve to the left of point 1 and to the right of point 4. Advanced numerical models for heat transfer modeling in phase change materials used for energy storage in buildings apply this approach for evaluating the heat exchanges inside the enthalpy-hysteresis loop [START_REF] Bony | Numerical model and experimental validation of heat storage with phase change materials[END_REF][START_REF] Delcroix | Development and numerical validation of a new model for walls with phase change materials implemented in TRNSYS[END_REF] .

Examples of cyclic paths with incomplete phase changes are given in Fig. 2 by the cycles [1-2-6-1] (after initially heating to point 1), [4-5-3-4] (after initially cooling to point 4), and [2-3-5-6-2] (after initially heating from point 1 to point 2 or initially cooling from point 4 to point 5 and then joining the cycle). Examples of cyclic paths experiencing no phase change are given by the cycles [2-8-2] (after having initially heated from point 1 to point 2), [6-7-6] (after having initially cooled from point 4 to point 6), [2-6-2] (after having initially heated from point 1 to point 2 or cooled from point 4 to point 6), and [8-7-8] (by reaching the cycle through the paths 1-2-8 or 4-6-7).

D. Numerical model

Assuming that the heat exchanges with the environment (radiation, convection) taking place all around the PCM, from one end to another, have a negligible impact, then the heat diffusion process is essentially one dimensional along the longitudinal direction x and the heat equation to solve for the temperature T (x,t) reduces to:

ρ ∂ H (T (x,t)) ∂t = - ∂ ∂ x (ϕ (x,t)) , (10a)
with the heat-conduction flux ϕ(x,t) defined by:

ϕ(x,t) = -k (T (x,t)) ∂ T (x,t) ∂ x , (10b)
where ρ = 4670 kg m -3 is density of VO 2 [START_REF] Leroux | From VO2 (b) to VO2 (r): Theoretical structures of VO2 polymorphs and in situ electron microscopy[END_REF] , together with the Dirichlet boundary conditions in Eqs. ( 1) and ( 2). The resulting partial differential equation for T (x,t) is both nonlinear and with memory since conductivity and enthalpy do not only depend on temperature but on whether it decreases or increases, see Fig. 2.

The classical heat conduction equation in Eq. ( 10a) is based on the Fourier equation in Eq. (10b) which itself is based on the thermodynamic equilibrium and the continuum hypothesis. Accepting the local equilibrium principle is possible only if the rate of change in the system macroparameters due to the external influences is much less than the rate of system relaxation to the local equilibrium. If the characteristic microscale of the system and time of its relaxation to equilibrium are significantly less than the characteristic macroscale and total time of the process, then the differential equations derived based on the local equilibrium principle and continuum hypothesis will be local both in space and in time. Thus, we obtain transport equations that do not contain the relaxation time and the characteristic scale of the microstructure [START_REF] Zhmakin | The zoo of non-fourier heat conduction models[END_REF] . As mentioned before, the MIT rate of VO 2 is very fast, in the range of picoseconds [START_REF] Muskens | Antenna-assisted picosecond control of nanoscale phase transition in vanadium dioxide[END_REF] . The present problem with a VO 2 layer in the millimeter range and temperature modulations in the Hz range is thus well described by the classical heat equation. Only for structures at nanoscale or at very high frequencies (e.g. in non-metallic materials, when the thermal diffusion length is in the range or smaller than the phonon mean free path [START_REF] Regner | Broadband phonon mean free path contributions to thermal conductivity measured using frequency domain thermoreflectance[END_REF] ) the wave nature of the heat propagation or the time lags become dominant which would require considering a non-Fourier heat equation.

The problem described by Eqs. (10a),(10b) was solved numerically by the finite volume method while implementing the enthalpy formalism [START_REF] Swaminathan | On the enthalpy method[END_REF] whose advantages make it a prefered tool for solving phase-change problems [START_REF] Al-Saadi | Modeling phase change materials embedded in building enclosure: A review[END_REF] . An explicit time-marching technique was adopted whose counterpart of simplicity is to constrain the time step to remain below an upper limit to grant numerical stability (Fourier criterion). For the treatment of combined hysteresis and partial phase change, the method of Bony and Citherlet [START_REF] Bony | Numerical model and experimental validation of heat storage with phase change materials[END_REF] for enthalpy has been applied and formally extended to conductivity, which formalizes the aspects discussed previously regarding the partial hysteresis loops (see also [START_REF] Delcroix | Development and numerical validation of a new model for walls with phase change materials implemented in TRNSYS[END_REF][START_REF] Klimeš | Computer modelling and experimental investigation of phase change hysteresis of PCMs: The state-of-theart review[END_REF] ). The details of the numerical model and its validation against analytical results in the case of linear diffusion can be found in the Supplementary Material.

III. RESULTS

A. Temperature and flux modulations

To illustrate the heat shuttling effect in phase change materials with hysteresis, we will consider a 1 mm thick VO 2 layer. We first consider the particular case where the mean temperature of the modulated reservoir is equal to the half-rise temperature along the heating branch, namely T c = 343 K. Because of the memory effect which is a direct consequence of the hysteresis, it is important to clarify what the thermal history of the PCM has been before starting the modulation scheme described by the boundary conditions in Eqs. ( 1) and (2), i.e. for t < 0. In this paper, we considered two different scenarios for the prior history. In the first one, the whole system is assumed initially at a temperature lower than the bottom temperature of the cooling branch (T < 331 K), after which the temperature of both reservoirs is increased to T c = 343 K. As a consequence, through thermal conduction, the PCM reaches (asymptotically) the same uniform equilibrium temperature T c from below. Hence, at the starting point of the modulation, i.e. t = 0, the coordinates {T, k} of the RP in the hysteresis loop of any location inside the PCM are {343, 4.8}, which corresponds to point n2 in Fig. 2(a). In the second scenario, the whole system is assumed initially at a temperature higher than the topmost temperature of the heating branch (T > 345 K), after which the temperature of both reservoirs is lowered to T c = 343 K. Similarly, the PCM reaches (asymptotically) the same uniform equilibrium temperature T c , but from above. Hence, the coordinates {T, k} of the RP of any location inside the PCM are now {343, 6}. In the second scenario, the tem-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0147225

perature modulation of the left reservoir starts while the conductivity of the PCM is uniformly at 6 W m -1 K -1 , whereas it is 4.8 W m -1 K -1 in the first scenario. Transient effects are observed after the onset of the temperature modulation; they die out after a few periods. In this paper we will concentrate on the steady-state periodic response after the transient effects have vanished. The temperature distribution inside the PCM, at four particular times during a period (i.e. after 1, 2, 3 and 4 quarters of a period), is reported in Fig. 3(a) (Multimedia view) for the first scenario, and in Fig. 3(b) (Multimedia view) for the second scenario. In both cases, the modulation amplitude of the left reservoir is T a = 8 K and the frequency is f = 0.3 Hz. At this frequency, the thermal diffusion length (a/(π f )) 1/2 , where a is the thermal diffusivity, is about 1.4 mm in the metallic phase and 1.1 mm in the insulating phase, which are of the order of magnitude of the layer thickness. At this frequency the PCM layer is thus neither thermally thin nor thermally thick. Dots with graded colors were used to represent the position in the PCM of the cells of the numerical model (53 cells in the present case) for these four snapshots in the temperature/conductivity representation.

The representative point {T, k} of the right boundary remains stuck at {343, 4.8} after the first prior scenario, resp. at {343, 6} after the second prior scenario. Meanwhile, the representative point of the left boundary follows the same trajectory in the {T, k} plane for both prior scenarios (partial hysteresis loop). The values reached at the four particular times considered in Fig. 3 are {351, 6}, {343, 6}, {335, 4.8} and {343, 4.8} . This means that the left boundary has lost memory of the prior scenario. Moreover, only the left part of the PCM (somewhat less than one quarter of the thickness) sees its thermal conductivity evolve during the cycling. The RP follows a partial hysteresis loop that progressively flattens out as one moves away from the left reservoir. The remaining of the PCM layer shows temperature changes (of progressively lower amplitude while approaching the right boundary) but no conductivity changes during the cycling (the representative point follows either a flat hysteresis loop or keeps in the full metallic phase). This means that less than 25% of the PCM close to the left reservoir can manifest its nonlinearity and more than 75% of the PCM close to the static reservoir, behaves either like a graded-conductivity linear material after the 1 st prior scenario (see Fig. 1 which qualitatively illustrates the two extreme thermal distributions related to this case) or a homogeneous linear material after the 2 nd prior scenario (the conductivity there remains uniformly equal to 6 W m -1 K -1 ). The feared consequence is that the heat-shuttling effect, which has been previously described for a material exhibiting nonlinearity in the whole thickness [START_REF] Liu | Energy harvesting from thermal variation with phasechange materials[END_REF][START_REF] Ordonez-Miranda | Net heat current at zero mean temperature gradient[END_REF] , is attenuated in the current situation.

Figure 4 shows the cyclic evolution of the left bath temperature, and of the heat flux at both boundaries of the PCM for the same conditions as before while assuming the first prior scenario (the equilibrium temperature T c = 343 K is reached from below before starting the cycling of the left bath temperature). The results when assuming the second prior scenario (the equilibrium temperature T c = 343 K is reached from above), although significantly different, are indiscernible by eye, therefore they are not reproduced here. At the considered frequency, i.e. 0.3 Hz, the PCM layer is thermally "intermediate" (i.e. neither thin nor thick), this is the reason why the heat flux at the right boundary is neither equivalent to the one at the left boundary nor considerably lower. This also explains the phase shifts between the three curves. The heat flux at the right boundary is nearly sinusoidal, probably because the right part of the sample behaves linearly. In contrast, the heat flux at the left boundary shows a bump at the beginning of the first quarter and a depression at the end of the third quarter. Those are coincident with the two time intervals where the left part of the PCM undergoes the strongest changes in conductivity.

The mean heat flux during one period (i.e. the net heat flux) is not zero, contrarily to what is observed with a linear material. It is however very small as compared to the magnitude of the extremal heat flux values on both left and right boundaries. We added in Fig. 4 a horizontal line at the level of the net heat flux; it is hardly discernible with the baseline. Indeed, the net heat flux is only 0.124 kW m -2 whereas the heat flux at the left boundary of the PCM oscillates between -78.2 and +79.8 kW m -2 . The flux modulation is not symmetric, therefore we propose the following definition for the heat-shuttling factor S:

S = ϕ(t) ϕ max -ϕ min , (11)
where ϕ max and ϕ min are the maximum, resp. minimum val- ues of the conduction flux at the left boundary (i.e. where the temperature modulation is imposed) and ϕ(t) is the observed net heat flux (shuttling effect). The heat-shuttling factor is 0.079% with the first prior scenario; it is just slightly higher, namely 0.081%, with the second scenario, see Table I. From the previous observations we can infer that the PCM can manifest its nonlinearity dynamically and hence give a chance to exhibit a heat-shuttling effect only if the left reservoir cools down sufficiently to reach the upper boundary of the cooling branch of the hysteresis loop, namely 337 K. The amplitude of the oscillations should thus be higher than 6 K. For lower amplitudes, the PCM behaves like a graded, yet linear, material (1 st prior scenario) or stays fully in the metallic phase (2 nd prior scenario).

A way to see a heat-shuttling effect for lower amplitude values of the left thermal bath is to move the mean temperature T a closer to the median temperature of the main hysteresis loop, i.e. 339 K. Actually, the previous statement can be generalized by saying that the PCM can exhibit a heat-shuttling effect only if T L (t) reaches the conductivity-variable-part of both the cooling branch and the heating branch, at least partially.

Let us now consider the situation where T a = 339 K, all things equal otherwise. The corresponding temperature distributions inside the PCM after 1, 2, 3 and 4 quarters of a period are reported in Fig. 5(a) (Multimedia view) for the first prior scenario and in Fig. 5(b) (Multimedia view) for the second prior scenario. We can notice that 21 (resp. 18) cells on the left side of the PCM (about 40%, resp. 34%, of total thickness) experience conductivity changes, among them 6 (about 11%) experience a full hysteresis loop. The remaining cells This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0147225 in the right part, which amounts to 60%(resp. 66%) of total thickness, see no conductivity change: 35%( resp. 30%) manifest flat loops anchored on the heating branch (resp. cooling branch) and 25%(resp. 36%) near the right reservoir stay at conductivity level of the insulator phase (resp. metallic phase) in the case of the 1 st prior scenario (resp. 2 nd prior scenario). The temperature T L (t) and the fluxes ϕ L (t) and ϕ R (t) are reported in Fig. 4 for both scenarios. Higher differences between the two scenarios are now observed, especially for the flux at the right boundary, ϕ R (t), mainly because in this area the material stays in two very different states, either in the insulator phase or in the metallic phase. The net heat flux was again added, which can now be distinguished from the baseline (altough with difficulty). As reported in Table I the heatshuttling factor now reaches higher values, 0.38% and 0.52% with the 1 st , resp. 2 nd prior scenario. The observed increase, as compared to when T a was set to 339 K, is due to the fact that a larger portion of the layer is allowed to manifest the intrinsic nonlinearity in conductivity of the PCM material.

B. Limiting value of the heat shuttling factor

Finding the maximum value of the shuttling factor S taking into account all the factors present in the considered problem (i.e. sigmoidal-type evolution of the conductivity and enthalpy, hysteresis), seems to be quite difficult due to the nonlinearity of the problem and the absence of an exact analytical solution. However, an evaluation can be attempted for an asymptotic case. Recall that the heat shuttling effect is closely related to the variability of conductivity with temperature and the fact that high and low conductivity values are encountered successively during a time cycle. To simplify the problem and to search for a configuration that can maximize the heat shuttling effect, let us consider a hypothetical material with two phases of low and high conductivity, resp. k min and k max , an isothermal phase change at T c , negligible latent heat, and no hysteresis. Furthermore, we assume that the temperature of the thermal baths is described by Eqs. ( 12)-( 13) with a very slow modulation (quasi-static regime). Asymptotically, for ω → 0, the flux at the left boundary can be approximated, for each half-cycle, by the solution of the linear, constant-property case related to, alternatively k max and k min (see Eqs. ( 22)- (24) in Suppl. Mat.). The mean heat flux (the numerator in Eq. ( 11)) is thus given by T a (k maxk min )/(πl), whereas the flux difference (the denominator in Eq. ( 11)) is T a (k max + k min )/l. The shuttling factor obtained in this idealized case is thus

(k max -k min )(k max + k min ) -1 π -1
. This is an increasing function of the ratio k max /k min with an asymptotic limit equal to 1/π ≈ 0.32. With the conductivity values of the insulating and metallic phases of VO 2 , the hypothetical value of the shuttling factor would be about four times less, i.e. 8.0%.

C. Parametric analysis

A parametric analysis was performed in order to evaluate the sensitivity of the net heat flux to the PCM properties and operational conditions, and hence to find ways to increase the heat-shuttling factor compared to the results obtained so far and in view of the limiting value described in §III B. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Figure 7 describes the influence of the hysteresis width ∆ h on the net heat flux and the heat-shuttling factor. The conditions are those of the first analysis in Fig. 3,4, namely T c = 343 K, T a = 8 K, and f = 0.3 Hz. The results show that an increase of the hysteresis width induces a rapid, almost linear, decrease of the heat-shuttling effect. The prior scenario has nearly no effect on the results. When comparing the real case with a hysteresis width of 8 K with the hypothetical case where there is no hysteresis, the net heat flux is lower by a factor of 0.024 in the real case and the heat-shuttling effect by a factor of 0.045. Figure 8 describes the influence of the mean temperature of the PCM on, again, the net heat flux and the heat-shuttling factor. The other conditions are T a = 8 K, f = 0.3 Hz, and a hysteresis width of 8 K. The net heat flux and the heatshuttling factor present a maximum which corresponds, or is close to, the median temperature of the hysteresis main loop, i.e. 339 K. The results confirm what was already noticed before (see Table I), namely that the shuttling effect is stronger when applying the 2 nd prior scenario, i.e. by reaching the mean temperature from above before starting the cycling (there is approximately a 40% increase in net heat flux and heat-shuttling factor when applying this scenario). One can also notice that with the 2 nd prior scenario, the optimal T c (338.8 K) is slightly lower than the median temperature of the hysteresis main loop (339 K).

PLEASE CITE THIS ARTICLE AS

Figure 9 describes the influence of the oscillation amplitude of the PCM. The other conditions are T c = 339 K, f = 0.3 Hz, and a hysteresis width of 8 K. The net heat flux and the heatshuttling factor monotonically increase with the amplitude of the temperature modulation of the left reservoir. They are higher and their rate of increase is higher as well when applying the 2 nd prior scenario. One may then ask what happens if the amplitude is increased further and how the shuttling factor evolves with respect to the limiting value described in §III B. In fact, the shuttling factor continues to increase, at first almost linearly, and then slows down until it reaches an asymptotic limit. For an amplitude of 50 K and 100 K, the shuttling factor reaches respectively 3.86% and 5.44% in the first scenario, 4.00% and 5.53% in the second scenario. Then, for a hypothetical amplitude of 300 K, the shuttling factor would reach 6.69% and 6.71% for, respectively, the 1 st and 2 nd scenario, to be compared with the limiting value of 8.0% obtained in §III B. Figure 10 describes the influence of the modulation frequency. The other conditions are T c = 339 K, T a = 8 K, a hysteresis width of 8 K and a VO 2 layer 1 mm thick. Whatever the value of the frequency, the net heat flux is always obtained higher (from 20% to 70%) with the 2 nd prior scenario. The net heat flux is highest at vanishing frequency (thermally thin PCM case). After a slow decrease in the thermally thin range ( f ∈[0, 0.03] Hz), the net heat flux decreases more rapidly in the intermediate range ( f ∈[0.03, 1] Hz), and then, depending on the scenario, either it shows a minimum at ∼0.97 kW m -2 for a frequency of about 3 Hz, and finally reaches a plateau at about 1 kW m -2 (2 nd prior scenario), or it continues to decrease but at a slower rate (1 st prior scenario). Because the amplitude of the heat-flux excursions on the PCM left boundary continuously increases with the frequency, the heat-shuttling factor monotonically decreases with the frequency, whatever the scenario. Starting from values close to 2.2% at 0.01 Hz (nearly static regime) the heat-shuttling factor first slowly decreases to 0.8-1% at 0.1 Hz, then decreases faster as f α with α between -2/3 and -1/2.

D. Heat-shuttling with two-bath modulation

Until now, the heat shuttling was produced by temperature oscillations of one reservoir, the other being maintained at the average temperature of the first. We will now explore the This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. case where both reservoirs experience temperature modulations while maintaining a zero mean thermal bias. The analysis will be focused on the process where the thermal baths share the same temperature modulation (same mean T c , same amplitude T a and same frequency f ) but with a phase lag δ :

PLEASE CITE THIS ARTICLE AS

T L (t) = T c + T a sin (2π f t) , (12)
T R (t) = T c + T a sin (max (0, 2π f t -δ )) ,

with δ in [0, 2π]. According to Eq. ( 12),( 13), both thermal baths are at temperature T c at t = 0 which is assumed to be the initial temperature of the PCM as well. As before, this initial state is considered to have been reached either by heating or by cooling. Subsequently, the temperatures are continuously modulated with the right bath behind the left one when δ ∈ [0, π], or ahead when δ ∈ [π, 2π]. Figure 11 describes the temperature/conductivity distribution inside the PCM at four different times of a cycle when the temperatures of the reservoirs are in quadrature. Some differences can be noticed depending on whether one applies the first prior scenario (Figure 11

FIG. 11. Same as in Fig. 3 when both reservoirs are modulated with the same mean temperature T c = 339 K, same amplitude T a = 8 K, same frequency f = 0.3 Hz, and a phase lag δ = 3π/2 (the right bath is in quadrature ahead). In (a) the first prior scenario was applied, i.e. T c reached from below (Multimedia view), whereas in (b) the second prior scenario was applied, i.e. T c reached from above (Multimedia view).

for the 1 st , resp. 2 nd scenario. The results for different values of the phase lag are summarized in Fig. 13 regarding the net heat flux and the shuttling factor. This factor can be defined in two ways depending on whether the heat flux range at the denominator of Eq. ( 11) is taken on the left or right boundary. Interestingly, the direction of the net heat flux depends on the phase lag. It is always towards the reservoir whose modulation is ahead. In addition it vanishes when the modulations are in phase or out of phase. A maximum for the net heat flux is observed for a phase lag δ close to 3π/2 (a maximum in the opposite direction is observed for a phase lag close to π/2). When compared to the case where only the left reservoir is modulated, the net heat flux has increased by 90% (from 1.00 to 1.87 kW m -2 ) with the 1 st prior scenario, and by 34% (from 1.40 to 1.90 kW m -2 ) with the 2 nd prior scenario. The shuttling factor (when considering the flux range at the left edge) respectively increases from 0.377% to 0.759% and from 0.522% to 0.772%. Despite this relative improvement, the shuttling factor remains lower than 1%.

IV. DISCUSSION AND CONCLUSION

It was demonstrated in Ref. 19 that in a nonlinear material submitted to a quasi-static temperature modulation on one side, a significant heat-shuttling effect can be observed. A theoretical illustration was provided with VO 2 which, because of a MIT, shows a strong variation of conductivity (+67%) when This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. temperature increases from values lower than ∼ 333 K (insulating phase) to values higher than ∼ 345 K. The thermal hysteresis of VO 2 was ignored in Ref. 19 and vanishingly slow temperature modulations (quasi-static regime) were assumed. Anywhere inside the material, the oscillations of the representative point {T, k} were restrained to the heating branch of the main hysteresis loop of VO 2 . As an example, for a temperature modulation of 5 K, a net heat flux of 2.9 kW m -2 and a heat-shuttling factor of about 5.8% were obtained by setting the mean temperature at the VO 2 transition temperature along the heating branch, a condition that was found to yield the strongest heat-shuttling effect. The heat-shuttling effect in PCM materials was also studied in Ref. 18, but introducing rough approximations regarding the thermal hysteresis, both for conductivity and apparent specific heat: when reversing the direction of the temperature evolution, an instantaneous switch between the heating and cooling branches was applied, but no indication was provided regarding the shape of the partial hysteresis loops. It was then shown that increasing the hysteresis width (by keeping it lower than the temperature modulation amplitude) induced a reduction of the net heat flux by a few tens of percent. The considered values for frequency and PCM thickness corresponded to the thermally thick regime.

PLEASE CITE THIS ARTICLE AS

In this paper we explored the influence of the PCM hysteresis and modulation frequency on the heat-shuttling effect by modeling carefully the coupled variations of temperature, conductivity and enthalpy, especially for alternating episodes of cooling or heating of low amplitude. Thereby we transposed the concept of non hysteretic branches (NHBs) observed for electric resistance and optical reflectance of VO 2 [START_REF] Ramírez | Firstorder reversal curve measurements of the metal-insulator transition in VO2: Signatures of persistent metallic domains[END_REF][START_REF] Gurvitch | Nonhysteretic behavior inside the hysteresis loop of VO2 and its possible application in infrared imaging[END_REF] to thermal conductivity and enthalpy, while extrapolating it to wider branches, up to the width of the whole hysteresis loop, which, regarding enthalpy, corresponds to the model described in Ref. 20 and which could be called the gradual curve switching method. More sophisticated models such as the curve scale hysteresis model for incomplete phase transitions [START_REF] Ivshin | A constitutive model for hysteretic phase transition behavior[END_REF][START_REF] Barz | Paraffins as phase change material in a compact plate-fin heat exchanger-part ii: Validation of the "curve scale" hysteresis model for incomplete phase transitions[END_REF] could be used in the future after this particular representation of the partial loops has been validated for VO 2 by a detailed characterization of hysteretical properties.

In summary, we confirm the existence of a shuttling of heat by conduction in a nonlinear material under nonequilibrium conditions. However, we found that hysteresis, when present, contributes to reduce to a large extent the net heat flux, moreover it makes that the heat shuttling effect is dependent on the past thermal history, i.e. before starting the temperature modulation of the reservoir. The reduction of the net heat flux is due to the appearance of NHBs inside the material. In fact, since the material undergoes a temperature modulation of decreasing amplitude as one approaches the static thermal bath, the hysteresis loops gradually thin, then flatten and become NHBs. In the vicinity of the static bath where the conductivity hysteresis loops have flattened and transformed to NHBs, the conductivity remains static (although possibly variable in space), hence the potential nonlinearity k(T ) does not show up locally. As a consequence, the part of the PCM layer close to the static reservoir does not participate to the heat-shuttling effect (this is valid in the thermally thin regime -low to vanishing frequency -and the thickness of the non-participating part increases with frequency). In addition, the temperature amplitude of the modulated bath must exceed a value depending on the global mean temperature and on the hysteresis width for the heat-shuttling effect to appear. As an example, when setting the global mean temperature at the transition temperature along the heating branch (i.e. 343 K), the modulation amplitude must exceed 6 K, which is higher than the values considered in Ref. 19. Even with an amplitude of 8 K, the heat-shuttling factor remains at low values, whatever the prior scenario, namely less than 0.1% in the case of a modulation This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0147225 frequency of 0.3 Hz and a 1 mm thick VO 2 layer. The net heat flux and the shuttling factor can be maximized by setting the global mean temperature close to the median temperature of the VO 2 main hysteresis loop (the exact value depends on the prior scenario). Increasing the modulation amplitude beyond a minimum threshold helps to amplify the heat-shuttling effect. Another mean is to reduce the modulation frequency, down to the thermally thin regime. Among the two prior scenarios considered, i.e. with a preliminary equilibrium temperature reached by heating or cooling, we have shown that the second scenario systematically produces a stronger heatshuttling effect. Compared to the limiting value (e.g. 8% for VO 2 , see §III B), the shuttling effect remains low in intensity when considering all coupled thermal and phase-change phenomena. For example, with a temperature amplitude of 8 K, the shuttling factor is no more than about 2% in the quasi-static regime (thermally thin VO 2 layer) and it drops rapidly with frequency (e.g. ∼0.2% at 1 Hz for a 1 mm thick layer). The net heat flux therefore reaches only a very small fraction of the flux amplitude observed at the modulated reservoir. However, there remains the possibility of leveraging the sensitivity of the reported shuttling effect to the amplitude of the periodic thermal excitation in order to increase the shuttling factor somewhat.

Our work offers a useful tool for predicting and interpreting the response of a dynamical thermal system encompassing PCMs with hysteresis. This paper focused on VO 2 , which can explain that the results were presented in a dimensional form. To expand the generality of the simulations, a non-dimensionalization of the variables should be applied. A first advantage is that it helps reducing the number of parameters in the problem formulation. Afterwards, the results are expected to be easily downscaled to correspond to a wider class of materials and geometrical configurations than the single material and geometrical configuration considered when using the dimensional form. Classically, the space variable is nondimensionalized by the PCM thickness l, and the time variable by the ratio l 2 /a where a is the diffusivity of one of the two phases. The temperature (after subtracting the constant temperature at the right boundary) would be nondimensionalized by the temperature amplitude at left boundary, hence: T ≡ (T -T c )/T a . The remaining nondimensional parameters would be the conductivity ratio k m /k i , the Stefan number ρCT a /L, the mean temperature of the phase change (separately for the heating and cooling phases) ((T m + T i )/2 -T c )/T a , the temperature width of the phase change ((T m -T i )/T a , the hysteresis width ∆ h /T a , and the frequency f l 2 /a. The parametric analysis should then encompass joint scanning of all nondimensional parameters, which is left for another paper.

With respect to thermal management applications, it would be worth performing in the near future a comparison with other models, both numerical and analytical, that describe the steady periodic regime in a finite PCM layer (see e.g. Ref. 46-48), not necessarily in term of heat shuttling, but more generally in terms of spatio-temporal distributions of temperature and flux. For instance, it would be of particular interest to know how the antiresonant regime revealed in Ref. 48 is af-fected by the non-isothermal character of the phase change and by hysteresis.

A potential extension is the inclusion of radiative effects whereby emissivity is nonlinear and presents hysteresis as well, with the objective of an analysis of the interplay between nonlinear conduction, heat storage, thermal emission, all presenting hysteresis, for an improved modeling of thermal diodes, thermal transistors, thermal logic gates, thermal memories, and thermal memristors.

B. Equations for the flux at the boundaries

The relations in 5(b) and (c) giving the flux at cell boundaries are valid except for the two boundaries of the domain (PCM layer). There are two boundary conditions of 1 st type (Dirichlet) on both sides of the domain. The flux at the left boundary is thus expressed based on the Taylor expansions of the temperature at the center of first and second cells with respect to the temperature at the left side of the leftmost cell. A weighted difference of these two expansions allows eliminating the second order terms, which yields an expression for the temperature derivative at the left boundary which is of order O(∆x 3 ). Hence, the flux at the left boundary is expressed as:

ϕ n 1/2 = k n 1/2 9T n 1 -8T n 1/2 -T n 2 3∆x , (7)
where the subscript 1/2 refers to the left side of the leftmost cell. A symmetrical operation is performed to get the flux at the right boundary of the domain as a function of the temperature of cell M -1, cell M, and of the temperature on the right side of cell M:

ϕ n M+1/2 = k n M+1/2 -9T n M + 8T n M+1/2 + T n M-1 3∆x , (8)
where the subscript M + 1/2 refers to the right side of the rightmost cell.

C. Initialization

The initial boundary temperatures, T 0 1/2 and T 0 M+1/2 , are those specified by the boundary conditions in Eq. ( 1)-( 2) or ( 12)-( 13) in the main text, at time t = 0, i.e. T 0 1/2 = T L (t = 0) and T 0 M+1/2 = T R (t = 0). They share a common value, T L (t = 0) = T R (t = 0) = T c . Two scenarios were envisioned in this paper, either the thermal baths have previously reached this value by heating (from below 333 K) or by cooling (from above 345 K). Any point of the PCM then reaches (asymptotically) the same uniform equilibrium temperature T c . This isothermal state is considered to be the one at t = 0. The conductivity is then the same everywhere at t = 0, the value of which is calculated from Eq. ( 3)-( 4) in the main text, with reference to the heating or cooling curve of the main hysteresis loop in Fig. 2(a) depending on whether the common value T c has been reached by heating or cooling the heat baths. Enthalpy as well is uniform at t = 0, the value of which is obtained from Eq. ( 3) and (5) in the main text, by referring to the appropriate heating curve or cooling curve in Fig. 2(b).

D. Temperature determination from the enthalpy

The enthalpy increment H n+1 i -H n i is evaluated at each time step and for each cell i = 1, M by applying Eq. 5(a). Two virtual temperatures are then calculated. The first one, T v,S , is based on the assumption that the enthalpy variation is assimilated to sensible heat:

T v,S = T n i + H n+1 i -H n i C , (9)
where C is the specific heat (assumed equal and constant for both insulating and metallic phases in the small temperature range considered). The second one, T v,L , is obtained by inverting the enthalpy value H n+1 i into temperature according to Eq. ( 6)-( 9) in the main text. This means inverting the enthalpy value into temperature according to the heating or cooling curve in Fig. 2(b). The selection between the heating and cooling curves, and correspondingly the limiting joint values H i and H m in Eq. ( 6)-( 9) in the main text, is made depending on the sign of the enthalpy variation H n+1 i -H n i . If it is positive, we refer to the heating curve, and if it is negative, to the cooling curve. Then, if the enthalpy variation is positive, the new cell temperature is set to:

T n+1 i = min T v,S , T v,L . (10)
Alternatively, if it is negative, the new cell temperature is set to:

T n+1 i = max T v,S , T v,L . (11)
E. Conductivity determination from the temperature Once the temperature field at time step n + 1 has been evaluated, T n+1 i , i = 1, M, together with the boundary values T n+1

1/2 = T L (t = t n+1 ) and T n+1 M+1/2 = T R (t = t n+1 ), one has to calculate the corresponding conductivity values, namely k n+1 i , i = 1, M, in addition to k n+1 1/2 and k n+1 M+1/2 . For each of these M + 2 temperature values, a virtual conductivity k v is first calculated based either on the heating curve or the cooling curve of the main hysteresis loop in Fig. 2(a), depending on whether the temperature has increased or decreased from time step n to time step n + 1. This amounts to use Eq. ( 4) in the main text, with the volume fraction of the metallic domains evaluated with Eq. (3) in the main text, which is fed with the proper values of the temperature-transition limits T i and T m (i.e. by referring to the due heating or cooling process). Then, if the temperature has increased from time step n to time step n + 1, the new conductivity is obtained from:

k n+1 i = max (k n i , k v ) ; i = 1/2, 1, 2...M -1, M, M + 1/2, ( 12
)
whereas if it has decreased, the new conductivity is obtained from:

k n+1 i = min (k n i , k v ) ; i = 1/2, 1, 2...M -1, M, M + 1/2. ( 13
)
F. Time step, space discretization, and convergence criteria

The time step is based on the stability criterion related to the explicit (Euler) scheme when applied to a domain with uniform diffusivity a, which reads:

∆t ≤ ∆t max = ∆x 2 2a . (14)
In the present case of a time-varying diffusivity field bounded by the diffusivity values of the insulating phase, a i = 1.09 mm 2 s -1 , and the metallic phase, a m = 1.81 mm 2 s -1 , we adopted a conservative rule by considering the highest possible diffusivity with, in addition, a safety factor of 1/2:

∆t = 1 2 ∆x 2 2a m . (15)
Actually, the largest integer divider of the time period smaller than the previous value was finally selected for ∆t. The number of time steps in one period is then called N.

The simulations lasted over several periods until reaching the stationary periodic regime with sufficient accuracy. Once the J-th period ended, the temperature distribution T JN i was compared with the one obtained one period before, i.e. T (J-1)N i , to provide the RMS error E T :

E T = 1 M i=M ∑ i=1 (T JN i -T (J-1)N i ) 2 1/2 . ( 16
)
On the other side, the flux was averaged over the last period to give:

ϕ i = 1 N n=JN ∑ n=(J-1)N+1 ϕ n i ; i = 1, ..., M. (17)
Then, the time-mean of the flux was also averaged over space to give:

ϕ = 1 M i=M ∑ i=1 ϕ i , (18)
which corresponds to the net heat flux discussed in the main text and was written there as ϕ(t).

In the stationary periodic regime, the time-mean of the flux, ϕ i , should be uniform, also, the residual (relative) error on the net heat flux, E ϕ(t) was defined as:

E ϕ(t) = max (ϕ i ) -min (ϕ i ) max 1, |ϕ| . ( 19
)
The stopping criterion was the joint satisfaction of E T < 10 -4 K and E ϕ(t) < 10 -3 . Depending on the thermal conditions, this was reached in about four to several tens of periods (up to 100 in some cases). The PLM material was discretized into M cells of equal thickness ∆x. Depending on the modulation frequency, different values for M were tested: 29, 53, 105, 165 and 205. One of these values was retained if, when comparing the net heat flux calculated with this value and with the value before in the former list, the relative difference was less than 4 • 10 -3 . Also, for frequencies up to f = 0.3 Hz (this frequency value was considered for most numerical applications in the main text), a number of M = 53 cells was retained. Then, for frequencies up to f = 10 Hz, 30 Hz, and 100 Hz, M = 105 cells were used, resp. 165 and 205.

III. MODEL VERIFICATION

A simple verification was made against an analytical model describing the case of a linear material submitted to the same boundary conditions, namely a fixed temperature T c at the right boundary and a modulated temperature of amplitude T a at the left boundary and same mean temperature T c . For this purpose, we selected for the input data of the numerical model: T c = 350 K and T a = 1 K. As a consequence, the PCM remains constantly and totally in the metallic phase.

The relation between the Fourier transforms of the temperature (Θ L and Θ R ) and of the heat flux (Φ L and Φ R ) at the left and right boundaries is given by the quadrupole expression:

Θ L Φ L = M Θ R Φ R ; M = A B C D , (20)
where

A =D = cosh ξ √ iω , B = 1 b √ iω sinh ξ √ iω , C =b √ iω sinh ξ √ iω , ( 21
)
where b is the thermal effusivity of the material, b = kρC, and ξ is the square root of the diffusion time through the total thickness l, namely ξ = l/ √ a. The product ξ √ iω can also be expressed as α(1 + i) with α = l/µ where µ is the thermal penetration depth, µ = 2a/ω. Knowing that Θ R = 0, the Fourier transform of the heat flux at the left and right boundaries is then obtained as:

Φ L = b √ iωΘ L cosh ξ √ iω sinh ξ √ iω ; Φ R = b √ iωΘ L 1 sinh ξ √ iω . (22)
The verification focused on the heat flux on the left boundary. After some algebraic manipulations, we easily find that the amplitude |Φ L | is expressed by:

|Φ L | = k √ 2|Θ L | µ 1 + tan 2 (α) tanh 2 (α) tan 2 (α) + tanh 2 (α)
1/2 (23) and the phase, noted P(Φ L ), is expressed by:

P(Φ L ) = π 4 + arctan tan(α) tanh 2 (α) -1 tanh(α) (tan 2 (α) + 1) . ( 24
)
When α → ∞ (thermally thick material), we retrieve asymptotically the solution of a semi-infinite medium; in particular, the heat flux sinusoidal modulation is π/4 ahead of that of temperature.

A comparison between the heat flux values obtained with the analytical solution and the numerical model over one period has shown a relative RMS error of 9 • 10 -6 and a maximum relative error of 1.5 • 10 -5 (the differences were normalized by the maximum value of the theoretical flux), which was deemed satisfactory. Animations are provided that complement the figures 3, 5, and 11 in the main text. These figures represent the combined temperature and conductivity distributions inside the VO 2 layer at four particular times during a period P of the stationary cycling (namely at t = P/4, P/2, 3P/4, P). The animations describe fluently the time evolution of these distributions (the time step of the recording was 0.5°of phase). The correspondance between the animations and the figures is described in the following table. Refer to the captions of the figures in the main text for an explanation of the conditions of the PCM thermal simulation.

FIG. 2 .

 2 FIG. 2. Evolution of the thermal conductivity (a) and enthalpy (b) in relation to temperature during the heating (in red) and cooling (in blue) processes of VO 2 . The green lines describe heating/cooling paths without phase change inside the major hysteresis loop. The numbers are used to describe full and incomplete hysteresis loops (see text).

  3(a) and Fig. 4 c see Fig. 3(b) d see Fig. 5(a) and Fig. 6(a) e see Fig. 5(b) and Fig. 6(b)

FIG. 3 .

 3 FIG.3. Distribution of the PCM discretization nodes in the temperature/conductivity space at four particular times during thermal cycling (i.e. after 1, 2, 3 and 4 quarters of a period). VO 2 sample thickness: 1 mm, frequency: 0.3 Hz, temperature amplitude: T a = 8 K, mean temperature T c = 343 K (i.e. the half-rise temperature along the heating branch). Before the onset of cycling, a thermal equilibrium at T c = 343 K is reached by increasing (a) (Multimedia view) or decreasing (b) (Multimedia view) the temperature of the system. The color dots, from red to blue, are representative of each numerical cell from left to right of the PCM (total of 53 cells). Black dots are used to designate points at x = 0, 0.25l, 0.5l, 0.75l, l. The heating branch of the VO 2 hysteresis loop is in red dashed line; the cooling branch is in blue dashed line (from Fig.2(a)).

2 )FIG. 4 .FIG. 5 .

 245 FIG.4. Cyclic variations of the left thermal bath temperature (in blue), of the heat flux at left boundary of the PCM (continuous red line), and of the heat flux at right boundary (dashed red line), as obtained when the preliminary equilibrium temperature of 343 K is reached from below (case described in Fig.3(a)). A thin horizontal red line depicts the net heat flux, here 0.124 kW m -2 , which yields a shuttling factor of 0.079%.

FIG. 6 .

 6 FIG.6. Same as in Fig.4after changing T c to 339 K. In (a) the preliminary equilibrium temperature of 339 K is reached from below, in (b) it is reached from above. The net heat flux is 1.0 kW m -2 in (a) (thin horizontal red line), and 1.4 kW m -2 in (b), which yields a heat-shuttling factor of 0.38%, resp. 0.52%.

FIG. 7 .

 7 FIG. 7. Net heat flux (a) and shuttling factor (b) as a function of the hysteresis width ∆ h for the following conditions: VO 2 sample thickness: 1 mm, frequency: 0.3 Hz, temperature amplitude: T a = 8 K, mean temperature: T c = 343 K. The preliminary equilibrium temperature T c is initially reached by heating (red curve) or by cooling (blue curve).

FIG. 8 .

 8 FIG. 8. Net heat flux (a) and shuttling factor (b) as a function of the temperature mean value T c for the following conditions: VO 2 sample thickness: 1 mm, frequency: 0.3 Hz, temperature amplitude: T a =8 K. The preliminary equilibrium temperature T c is initially reached by heating (red curve) or by cooling (blue curve).

FIG. 9 .

 9 FIG. 9. Net heat flux (a) and shuttling factor (b) as a function of the temperature modulation amplitude T a for the following conditions: VO 2 sample thickness: 1 mm, frequency: 0.3 Hz, mean temperature: T c = 339 K. The preliminary equilibrium temperature T c is initially reached by heating (red curve) or by cooling (blue curve).

FIG. 10 .

 10 FIG. 10. Net heat flux (a) and shuttling factor (b) as a function of the frequency of the temperature modulation of the left bath for the following conditions: VO 2 sample thickness: 1 mm, temperature amplitude: T a = 8 K, mean temperature: T c = 339 K. The preliminary equilibrium temperature T c is initially reached by heating (red curve) or by cooling (blue curve).

2 331

 2 (a) (Multimedia view)) or second prior scenario (Figure11(b) (Multimedia view), but they essentially concern the core of the material. Temperature and heat flux on the boundaries of the PCM are reported as a function of time in Fig.12for the prior scenario where the mean temperature T c = 339 K is previously reached from above. The other case leads to substantially the same curves; in particular, the net heat flux reaches very close values, 1.90 and 1.87 kW m -

2 )FIG. 12 .FIG. 13 .

 21213 FIG.12. Cyclic evolution of temperature (in blue) and heat flux (in red) at the left reservoir (continuous lines), and right reservoir (dashed lines) in the conditions of Fig.11when applying the second prior scenario (the curves for the first scenario are visually the same). A thin red horizontal line indicates the neat heat flux level.

  Figure Animation 3(a) animation-fig-3a 3(b) animation-fig-3b 5(a) animation-fig-5a 5(b) animation-fig-5b 11(a) animation-fig-11a 11(b) animation-fig-11b

TABLE I .

 I Characteristic data (minimum and maximum heat flux at the left boundary of the PCM, net heat flux and shuttling factor) obtained for two different values of the mean temperature T c and depending on whether the initial equilibrium temperature T c is reached from below of from above. Common parameters are: VO 2 sample thickness l = 1 mm, frequency f = 0.3 Hz, amplitude T a = 8 K.

	T c (K)	Initial direction a	Min	Flux (kW m -2 ) Max	Net	Shuttling factor (%)
	343	↑ b	-78.21	+79.79	+0.124	0.0787
	343	↓ c	-78.65	+79.42	+0.128	0.0812
	339	↑ d	-125.6	+141.2	+1.00	0.377
	339	↓ e	-126.4	+142.0	+1.40	0.522
	a Initial equilibrium temperature reached by heating (↑) or cooling (↓)

b see Fig.

  This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

		332	334	336	338 339 340	342	344	346
		0.6
	Shuttling factor (%)	0.2 0.4
		0
		332	334	336	338 339 340	342	344	346
				Mean temperature T c (K)
	PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0147225

TABLE I .

 I Correspondance between the animations (mp4 format) and the figures in the main text.

V. SUPPLEMENTARY MATERIAL

See the supplementary material for the details of the numerical model and a verification performed against an analytical model for the limiting case of a linear material.

VI. DATA AVAILABILITY

The data supporting the findings of this study are available from the author upon reasonable request.

Supplementary material

Influence of thermal hysteresis on the heat shuttling effect: the case of VO2 Jean-Claude Krapez ONERA, The French Aerospace Lab, DOTA, F-13661 Salon cedex Air, France (*Electronic mail: jean-claude.krapez@onera.fr)

I. PRESENTATION

For the sake of clarity and completeness, here we provide some details on the numerical model that was developed to calculate the temperature evolution in a PCM layer with hysteresis, and on the verification performed against an analytical model for the limiting case of a linear material. We also describe the link between a series of six animations provided in mp4 format showing the dynamic evolution of joint temperature and conductivity distributions and the corresponding figures in the main text.

II. NUMERICAL MODEL

We assume 1D heat diffusion along a layer of thickness l which is in perfect contact with two heat baths, where the left one is at a sinusoidally modulated temperature:

and the right one either at constant temperature:

or equally modulated but with a phase lag:

The heat equation to solve is:

with the heat-conduction flux ϕ(x,t) defined by:

A. General equations for the discretized fields

The numerical solution is obtained after discretizing in space and applying the well known finite volume approach. For simplicity we used a constant space discretization interval ∆x = l/M where M is the number of cells; we also applied the explicit (Euler) time marching method which, from the temperature, enthalpy and conductivity discretized fields known at time step n allows expressing the enthalpy field at time step n + 1 according to:

where the subscripts i, i -1, and i + 1 are affected to temperature and enthalpy to designate the mean temperature, resp. enthalpy, of the cells numbered i, i -1, and i + 1 of common thickness ∆x, whereas the subscripts i-1/2 and i+1/2 are affected to heat flux and conductivity to designate the heat flux and conductivity at the left, resp. right boundary of the cell i. The conductivity at these boundaries is calculated according to the rule of thermal resistances in series, based on resp. the mean conductivity of cells i -1 and i on one side and of the cells i and i + 1 on the other side. :

where k n i-1 , k n i and k n i+1 are the conductivity values inferred from the mean temperature of the cells i -1, i, and i + 1 (see §II E for their determination). The Euler explicit finite-difference approximation leads to a truncation error of order O[∆t, (∆x) 2 ]. It is conditionally stable, meaning that there is a constraint on time step for the solution to be stable (i.e. nondivergent and nonoscillatory), see §II F. The Crank Nicolson is second-order accurate both in time and space, that is O[(∆t) 2 , (∆x) 2 ]. In addition, it is unconditionally stable. However, it is known that despite the unconditional character, oscillations with slow damping may appear. The disadvantage with the Crank Nicolson scheme (or any other partially or totally implicit scheme) is that it requires solving a system of equations and because the transfer equations are here non-linear in many respects, the solution can only be obtained by iterations. This considerably slows down the numerical process; moreover, the accumulation of roundoff-errors is exacerbated by the multiplication of intermediate operations. Both approaches, explicit and Crank Nicolson, have their advantages and disadvantages. For the present work, we found that the explicit approach gives satisfactory results, despite the truncation error of order O[∆t, (∆x) 2 ], see §III.