Thermostatistical Description of Small Systems in Nonequilibrium Conditions: Energy Conversion and Harvesting

Hysteresis cycles are very important features of energy conversion and harvesting devices, such as batteries. The efficiency of these may be strongly affected by the physical size of the system. Here, we show that in systems which are small enough, the existence of physical boundaries which produce nonhomogeneities of the interaction potential gives rise to inflections and barriers in the associated free energy. This in turn brings on irreversible processes which can be triggered under suitable external conditions imposed by a heat bath. As an example, by controlling the temperature, the state of a small system may be impelled to oscillate between two different structural configurations or aggregation states avoiding equilibrium coexistence and therefore dissipating energy. This cyclical behavior associated with a hysteresis cycle may be prototypical of energy conversion, storage, or generating nanodevices, as exemplified by Li-ion insertion batteries.


I. INTRODUCTION
The concept of stability is essential when describing the behavior and properties of matter [1].At the microscopic level, the stability of atoms becomes manifest through the existence of chemical elements, from which quantization can be deduced.Macroscopic thermodynamic systems show stability through the existence of distinct phases which can coexist under the same physical conditions.
Unlike macroscopic systems, the existence of boundaries in small systems confers on them peculiar characteristics.Finite size induces nonhomogeneities of the interaction energy which may give rise to free-energy barriers separating two different structural configurations or even thermodynamic phases [2,3].Transformations between these structural configurations or phases is a matter of thermodynamic transformation theory, a well-understood problem in classical thermodynamics [4,5].When a macroscopic system is subjected to destabilizing conditions, it separates into two or more phases that may coexist in equilibrium [4,5].This partitioning involves the formation of new free-energy barriers associated with interfaces and finally, from the thermodynamic point of view, to the emergence of new systems with their own free energy determining their physical properties: compressibility, specific heats, etc. [5].
However, when the system is finite and small enough, the formation of an interface could imply an excessively high energetic cost and therefore becomes energetically unfavorable.This energetic restriction has been observed, for instance, in the formation of magnetic domains in ferromagnetic materials [6] or in the hysteresis cycles of Li-ion batteries and other energy storage systems [7][8][9].Even more exotic systems such as atomic nanoclusters with magical numbers show peculiar effects on their behavior and properties similar to those of the small systems we will consider here [10,11].
In this article we will analyze the implications of the finite size of the system on its thermodynamic behavior, and then show that inflections and barriers in the thermodynamic free energy are related to irreversible processes.We show this by calculating the entropy produced in a transformation that implies the transition from one structural configuration or thermodynamic phase to another by overcoming the freeenergy barrier.This is possible due to the existence of a thermodynamic affinity originated by the external conditions.Since these processes may be cyclical, they can be used to generate or convert energy at the nanoscale if the external conditions (or constraints) imposed on the system are appropriate, that is, when the external constraints force the state of the system to reside in an interval of the order parameter domain corresponding to the unstable region of the free energy.Under these circumstances, we show that an oscillation between both two mentioned states separated by the free-energy barrier may be triggered and strongly affected by thermal fluctuations.This cyclical motion persists provided the system is maintained near the top of the free-energy barrier.
The results obtained here are of great interest since they underlie the physics of the energy generation and conversion nanodevices [12].Potential applications have been recently reported for energy nanogenerators [13], systems based on the pyroelectric effect [14], and also for electric energy storage systems such as batteries [7].Coupling of several small systems may lead to interesting effects that are in the thermodynamic origin of the electric hysteresis [8,9].Oscillating behaviors are also useful in energy-converting nanodevices whose operation depends on the pressure conditions imposed by the heat bath [15].This article is organized as follows.In Sec.II we analyze the stability conditions and their relation with small systems.After analyzing the implications of these conditions on the bases of statistical mechanics, we show how the free-energy barrier leads to an irreversible oscillating behavior in the case of small systems.We also obtain the expression for the entropy produced during these oscillations.Section III is devoted to analyzing the oscillations of an order parameter associated with the temperature or average kinetic energy of a toy model of a small system as a function of time for its bistable free energy in the deterministic and stochastic cases.Section IV analyzes the problem from the perspective of a Fokker-Planck equation for the probability density of the order parameter.A numerical solution of the equations shows the oscillation of the average value of this order parameter and is used to determine an apparent stationary, time-averaged distribution function.Section V is devoted to analyzing the case of insertion batteries.Finally, in Sec.VI we present our main conclusions.

II. THERMODYNAMIC STABILITY AND THE ENTROPY PRODUCED WHEN OVERCOMING A FREE-ENERGY BARRIER
In this section we will use thermodynamic stability theory and classical statistical mechanics to analyze the behavior of small systems in contact with a heat bath when the external conditions impel the system to oscillate between two different structural configurations or aggregation states avoiding equilibrium coexistence and therefore dissipating energy.
For the sake of clarity, it is convenient first to introduce the concepts of thermodynamic stability theory and apply them to classical small systems such as droplets or bubbles in which the distinction between bulk and surface can be performed [16].After this, we start our general analysis of small systems which avoids the distinction between bulk and surface, since it is not clear in the general case.

A. Thermodynamic stability
In the framework of thermodynamics [5] the stability condition of a large homogeneous system with respect to thermal fluctuations establishes that the heat capacity at constant volume and particle number must be positive definite, C V ,N > 0. In the canonical ensemble, this condition can be expressed in terms of the Helmholtz free energy F (T ,V ,N) which must exhibit an absolute extremum minimum as a function of the extensive variables or an extremum maximum as a function of the intensive variables.Accordingly which is the concavity condition.
In the case of finite systems the free energy has two contributions, one associated with the bulk and the other one with the boundary separating the system from the surroundings.Therefore the free energy can be written as where F b is the bulk free energy, ρ is the local mass density related to the radial distribution function, and φ(r) is an effective interaction potential taking into account the mentioned boundary effects [17].Here, V ρ(r,T ,N )φ(r)dV = φ T ,N is the average energy which plays the role of an order parameter or reaction coordinate similar to the corresponding one introduced previously [3,10].Hence, the stability condition now reads where C b is the contribution of the bulk to the heat capacity and where for simplicity's sake we have assumed φ T ,N ≡ η(N )h N (T ).Equation (3) shows that the stability of a small system depends dramatically on the curvature of h N (T ) in such a way that if this function has inflection points, provided η(N )∂ 2 h N /∂T 2 > C b /T , the stability condition is violated.Therefore, in this pathological situation a convexity (barrier) in the free energy will appear.Additionally, from Eq. ( 3) it follows that a critical size of the system N c (T ) (stability curve) separating the stable (N > N c ) from the unstable (N < N c ) regions of the (T ,N) diagram exists and can be obtained as a solution of The discussion in the previous paragraph can be illustrated by considering the following expression of the free energy [16]: , where in the righthand side the first term constitutes the bulk contribution with μ being the supersaturation and the second term constitutes a surface free energy.Here, we can identify h N (T ) = A(T ) and η(N ) = N 2/3 , and therefore the necessary condition for the appearance of the barrier is ∂ 2 A/∂T 2 > 0. Note that in the more general case, when the system has a fractal boundary, the function η(N ) will be of the form η(N ) = N ν/3 with ν < 2 and with a similar behavior as described in the previous example.

B. Thermodynamic stability approach to small systems
In this work, we will base our description on the assumption that small systems can be characterized as an effective homogeneous system whose intensive parameters differ from those of the corresponding macroscopic system [16,[18][19][20].Other approaches to the description of processes in systems not in equilibrium with a bath have been recently reported [21,22].
The partition function of the (small) system is defined, in classical terms, through the well-known relation: Z(β) ∼ exp[−βH (p,q)]dpdq, where p and q are the generalized momenta and coordinates, and β ≡ 1/k B T with k B the Boltzmann constant and T the temperature.As follows from its definition, Z(β) does not depend on the phase-space coordinates, but only on thermodynamic variables.These thermodynamic variables can in general be calculated by means of statistical averages; for example, the temperature T is related to the average kinetic energy [23,24].More precisely, the temperature of a small system is defined here as the configurational average of the total mean kinetic energy of the system [24][25][26][27].The difference between the partition function of macroscopic and small systems relies on the fact that the Hamiltonian H (p,q) of the latter contains long-range interactions [28,29], related to boundary effects and not necessarily to external bodies.
In order to show that the free-energy barrier gives rise to irreversible processes, let us assume that the stability condition of a small system described in a global way by its partition function Z(β) is still given by the relation where C V ,N refers to the heat capacity at constant volume of the whole system.
After performing straightforward calculations we deduce that from Eq. ( 4) the following relation for the heat capacity follows: where Z = (∂Z/∂β) V ,N and Z = (∂ 2 Z/∂β 2 ) V ,N .Hence, for a stable system at constant volume and particle number Eq. ( 5) implies Z Z/Z 2 > 1 while for a unstable system Z Z/Z 2 < 1.Thus, let us write with α(β) = 0, such that the system is stable for α(β) > 0 and unstable for α(β) < 0. After successive elementary integrations and arrangements, Eq. ( 6) leads to where More insight can be achieved by rewriting Eq. ( 7) using exp [−R(β)] ≡ exp (−R 0 ) [1 + χ (β,β 0 )] which defines χ (β,β 0 ) and taking also into account that U 0 − F 0 = T 0 S 0 , with S 0 = S(β 0 ); we obtain where F = F − F 0 .In Fig. 1 we show a schematic representation of Eq. ( 8) when a free-energy barrier exists.In light of Fig. 1, due to the presence of the barrier there are two possibilities: (a) For a macroscopic system, its state can be represented along the straight line, i.e., the coexistence line, indicating that this system suffers an equilibrium phase transformation and splits into two different phases having the same temperature.In this case, the integral on the right-hand side of Eq. ( 8) is zero.(b) For a small system, its state can be represented along the solid line indicating that this system undergoes a nonequilibrium transformation.Since β F − T 0 S 0 (β − β 0 ) < 0 (as Fig. 1 shows), then β β 0 χ (β,β 0 )dβ < 0 regardless of the direction of the evolution in β.From our analysis, it follows that the integral on the right-hand side of Eq. ( 8) clearly is related to the bending of the free-energy difference in Fig. 1.
From Eq. ( 8) and by using the thermodynamic relation F = U − (T S), we can obtain the total entropy change which in general can be considered the sum of two contributions: Here, e S is the entropy change due to exchange of matter and energy with the surroundings, as stated by Gibbs' equation: T e S = U − W .The term i S is the entropy change due to the "uncompensated heat transformation," i.e., the entropy produced by irreversible processes within the system which according to the second law of thermodynamics must be nonnegative, i S 0 [30].The sign of the entropy production i S indicates the direction of evolution of the nonequilibrium system in a natural process.In view of this, Eq. ( 9) can be rewritten in the form where W rev is the reversible work.This result states that the term on the left-hand side of Eq. ( 10) is the net work performed in the presence of irreversible processes, which is less than the reversible work due to the existence of the dissipated work W dis = T i S. Thus, one can write For processes occurring at constant volume, W rev = 0, Eq. ( 11) reduces to Thus, from Eqs. ( 8) and ( 12) we deduce the following important relation for the free-energy change: The following important result must be stressed: Eq. ( 13) can be reinterpreted in terms of the entropy produced by the system during the transformation starting from β 0 where now F is a Landau-like free energy and we used the relation β = 1/k B T .Notice that if the transition occurs from β 1 then the previous equation takes the form The superindexes indicate that the process takes place from T 0 or T 1 , respectively.
According to our discussion after Eq. ( 8), i S becomes zero for equilibrium phase transitions and thus it follows that F = S 0 (T 0 − T ).The free-energy change follows the Maxwell construction between potential wells; that is, the system performs an equilibrium phase transformation [5].This occurs for macroscopic systems by virtue of their extensive character.
On the other hand, when F < S 0 (T 0 − T ) then i S > 0; the system performs an irreversible process.These processes may appear when the energetic cost of partitioning the system is too high, as occurs for small systems or systems where finitesize effects are important as those showing hysteresis [7][8][9], or many nanodevices designed and used for energy conversion and generation [7,[12][13][14].
To summarize this section, we may conclude that forcing the system to reside in the region corresponding to the barrier (spinodal region) by controlling the temperature of the heat bath, i.e., maintaining the bath at a critical T c , leads to two different outcomes.On the one hand, if the system is large enough-extensive-it breaks into two different macroscopic systems experiencing an equilibrium phase transformation in which two different phases coexist in equilibrium.On the other hand, if the system is small, partitioning into coexisting phases could imply a too high energetic cost, becoming then an energetically unfavorable process.Thus, the finite system will remain in a stationary state strongly sensitive to fluctuations that impel the system to relax to one of the nearest metastable states, either T 0 or T 1 .However, since the temperature of both states differs from the bath temperature, a thermodynamic affinity appears driving the state of the system towards the unstable state T c where again it cannot remain.The process just described will start again and so on if the temperature of the heat bath is maintained constant at T c or near to it.Notice that, in practice, since the oscillating process is irreversible and therefore an uncompensated heat is produced constantly, an alternative inflow and outflow of energy is required to compensate for this wasted heat.We call critical oscillations to this irreversible cyclical process between the states I and II corresponding to T 0 and T 1 .This is an operation pathway followed in, for instance, energy nanogenerating devices based on the pyroelectric effect [14].For a finite system Eq.( 14) and Eq. ( 15) represent the entropy produced during the critical oscillations between the states I and II at T 0 and T 1 , respectively, and passing through T c .
It is worth mentioning that a similar oscillatory behavior between different metastable states has been observed for simulations of atomic nanoclusters in Refs.[3,11,31].In the literature on this subject, this oscillatory behavior has been described through a single and stationary double-peaked probability distribution from which equilibrium thermodynamic properties have been inferred [11].However, these stationary double-peaked probability distributions seem to not be representative of an equilibrium thermodynamic formalism as we show in the present work.

III. CRITICAL OSCILLATIONS OF THE ORDER PARAMETER
Let us consider now that the heat bath imposes on the system the temperature T c = 1/k B β c , that is, the temperature at the critical point of the free-energy barrier.Thus, we may write the system's temperature as T = T c − T .Substituting this into Eq.( 14) and Eq. ( 15) we obtain (16) where the j = 0,1 distinguishes the processes described by Eq. ( 14) and Eq. ( 15), respectively.Now, expanding in series the free energy around T c we have that since by hypothesis the free energy is a bistable function.
For simplicity we have assumed that the standard state is defined at T c .Substituting now Eq. ( 17) into Eq.( 16) and adding the entropy produced in each irreversible process we obtain where the last term is the irreversible heat exchanged between the system and the bath when passing between states T 0 and T 1 through T c : Q irr ex = 1 j =0 S j (T j − T c ).Note that the term proportional to c = 1 j =0 S j introduces an asymmetry on the wells of the free energy landscape proportional to the average of the entropies associated with states T 0 and T 1 .
Following Onsager nonequilibrium thermodynamics, the dynamics of the system can be derived by taking the time derivative of Eq. (18).We obtain To satisfy the second law statement d i S/dt > 0 we may adopt linear relationships between fluxes dT /dt and thermodynamic forces X = −∂F /∂T .This yields the dynamic equation where τ is a characteristic time (an Onsager coefficient) and the time dependent heat exchange rate, dQ irr ex /dt ∝ −dT /dt, can be modeled through an oscillating function.In the present case we can assume that it is determined by the square wave function j irr q (t): where d is the amplitude of the square wave.This term implies that the external forcing coming from the heat exchange with the bath tilts the bistable potential periodically in two directions, thus alternatively favoring the states which correspond to T 0 or T 1 , acting as an energy source to maintain the oscillating behavior of the system.
Figure 2 shows the time behavior of the state of the system for different values of parameters controlling the tilting and the amplitude of the oscillation of the free-energy landscape due to the heat exchange that alternatively favors states I and II.The red (solid) lines correspond to the numerical solution of the deterministic equation (20), whereas the blue (irregular) points correspond to the stochastic equation (22), that we will discuss in detail in the following subsection.Figures 2(a) and 2(b) show an asymmetric case in which the heat exchanged is large enough to produce deterministic oscillations in the state of the system.The different cases correspond to opposite tilts of the free energy as indicated by the initial behavior of the deterministic solution.Figure 2(c) corresponds to a case with no stationary tilting (c = 0) and a strong amplitude of the energy exchange d = 0.6.
Figures 2(d) and 2(e) correspond to cases in which the heat exchange is weak and produces a small oscillation within the same well.In this case, however, fluctuations in the rate of heat exchange may induce transitions between states I and II.Finally, Fig. 2(f) corresponds to the case when no tilting and no oscillating force is taken into account.In this case, if the noise is sufficiently large, it may spontaneously produce transitions between potential wells, as is well known.The difference between Figs. 2(a)-2(c) and 2(d)-2(f) is that the first ones lead to a nonequilibrium distribution function whose average value oscillates between states I and II.The other cases show other possibilities in which small external oscillations together with noise may also induce random transitions between the two states.

IV. CRITICAL OSCILLATIONS OF THE PROBABILITY DENSITY
In view of the previous analysis, it is important to consider now that the variables characterizing the state of small finite systems may fluctuate.In our case, this implies that the dynamic equation for the temperature of the system should be perturbed by a noise term coming from thermal fluctuations.This term has its origin in the energy fluctuations that a system in contact with a heat bath may experience.Thus, Eq. ( 20) becomes where ξ ( t) is a random Gaussian noise term with zero mean and second moment (fluctuation-dissipation theorem): ξ ( t)ξ ( t ) = 2Dδ( t − t ), with D = T 2 bath /τ the intensity of noise.
The presence of a noise term is interesting because it makes more drastic the oscillations of the system between states I and II.This is shown in Fig. 2, where different cases are shown for the intensity of the tilting term c in Eq. ( 22).In Fig. 2(d) we show the case when the random heat exchange induces a permanent oscillation between states I and II even in the case of negligible deterministic oscillations.The effect of noise is to change the frequency of the oscillation.
A interesting consequence of Eq. ( 22) is the fact that it is equivalent to a Fokker-Planck equation for the nonequilibrium distribution function P (T , t) or equivalently for the microscopic averaged potential.Since the temperature of the system can be defined as the microscopically averaged kinetic energy of the system k B T = K and in view of the virial theorem 2 K = − φ T ,N , one has 2k B T = − φ T ,N = − ( t).The corresponding Fokker-Planck equation may be determined by comparing the moments of (22) with those coming from a Kramers-Moyal expansion; the result is where we have introduced the diffusion coefficient of the energy D = (k B T bath ) 2 /τ , and the dimensionless time t = t/τ .The state of the system spends different times in states I and I I because of the asymmetry of the potential, favoring state I .The time spent at the middle of the oscillation 1.5 is much lower than the time spent at the extrema.Thus, an apparent stationary bimodal distribution can be obtained by averaging P ( ,t) over a window of time larger than the period of the oscillations.Panel (d) shows three different reduced time averages P ( ,t) t of the oscillating probability distribution P ( ,t).The red dashed line represents an average performed on the time interval t = 20, the blue dotted line in the time interval t = 60, and the black solid line in the time interval t = 100 all for an initial condition of t = 2.The influence of the initial condition is evident in the first average.For averages taken over longer intervals the result converges to an apparent stationary bimodal distribution.
The numerical solution of this equation has been calculated by using a BDF implicit method implemented in Mathematica (Wolfram Research) assuming that the initial condition is a Gaussian distribution located at 0 = 1.3 with variance σ 2 = 0.06.The results obtained are shown in Figs.3(a)-3(c).The oscillation of the distribution function is shown by the fact that its maximum oscillates around = 1.5.The oscillation of the energy distribution P ( ,t) is shown in Fig. 4 to make clearer the dynamic ensemble.This result is fully consistent with the oscillations reported in Fig. 1 of Ref. [3] for the short time (microscopic) average of the interaction energy.A similar result is also commented on below Fig. 2 of Ref. [11], where the authors mention that it must be described as a special case of "dynamic equilibrium" (dynamic coexistence); we conclude that the distribution of energy P ( ,t) can also be periodic.
The bimodal distribution represented by the solid line in Fig. 2(b) of Ref. [11] is the result of MC simulations whereas the dashed lines are Gaussian fits.Likewise, the bimodal distribution represented in Figs. 3 and 4 of Ref. [32] can be interpreted as the cumulative effect of two distributions characterizing the two possible states of the clusters, just as is represented in Fig. 2(b) of Ref. [11].Since the local minimum of this distribution of macroscopic states has a finite value, the population of clusters within the unstable region is not negligible.This fact indicates the presence of a macroscopic number of clusters which would suggest that the system is at chemical equilibrium.In fact, it can be probed that in an equilibrated chemical reaction the population of the intermediate state is negligible (strictly zero) whereas for a nonequilibrated reaction the population of the local minimum is finite.As we mentioned above, this rather suggests that the bimodal distributions represented in Figs. 3 and 4 of Ref. [32] does not correspond to an equilibrium distribution of energies.It seems instead a time average of a nonequilibrium distribution over a large enough window of time.
In Figs.3(a), 3(b), and 3(c) we show three snapshots of the probability distribution P ( ,t) in terms of the order parameter for ten different reduced times t = t/τ as obtained by numerically solving Eq. ( 23).For dimensionless time t > 10, P ( ,t) returns periodically to states I 1.2 and I I 1.7.Averaging this behavior over a sufficiently large period of time leads to an apparent stationary double peaked distribution, P ( ,t) t .This fact is shown in Fig. 3(d), where the red dashed line represents a time average performed in the time interval t = 20, the blue dotted line in the time interval t = 60, and the black solid line in the time interval t = 100.As expected, the initial condition deforms the averaged distribution in the first case.Once the average becomes independent from the initial condition, the result takes the form of an apparent stationary bimodal distribution.Unlike the stationary distributions of energy reported in the literature [11,32], this time average cannot be interpreted as representative of chemical equilibrium and cannot be used to extract the equilibrium thermodynamic properties of the system [33].Figure 4 shows the three-dimensional representation of the normalized probability density as a function of for different reduced times.

V. SMALL PARTICLE INDUCED VOLTAGE OSCILLATION IN INSERTION BATTERIES
During charging and discharging of a lithium-ion battery there is an exchange of Li ions that are respectively released and adsorbed by the interstitial lattice sites of the electrode, that may be made up by FePO 4 crystals [9].Experimental studies have reported that the kinetics of this process strongly depends on the size of the FePO 4 crystals that made up the electrode [7][8][9]34].
For small enough sizes (nanometric) of FePO 4 [34] or TiO 2 particles [7], when the cell voltage takes values in the hysteresis zone, instead of the usual coexisting phase separation characterized with zones with large Li concentration and zones with low concentration, typical of large crystals or bulk matter, an alternating process is observed in which the small electrode crystals become successively filled and unfilled.This behavior is an example of the phenomenology already described and represented, in this case, through the oscillations of the hysteresis cycle, according to the different scenarios of an ensemble of particles that have been characterized in Refs.[8,9].
Using our framework to analyze this system, we will describe the kinetics of a single small FePO 4 or TiO 2 particle in the presence of a heat bath controlling the chemical affinity of the system.To achieve this program, we first note that the kinetic change of crystallographic phases is dominated by the mole fraction y of Li ions occupying holes in lattice sites of the electrodes, and therefore its kinetics is more conveniently described by the Gibbs free energy, G(p,T ,y), instead of the Helmholtz free energy we previously used.This free energy has two contributions, one coming from the bulk and other from surface effects G(y) = G b (y) + G s (y), Refs.[7,9], as we have considered previously.
Let us now assume that temperature and pressure are constants and that in this case the bath is made up by a solution in which the Li ions are free to diffuse.Controlling the number of these ions in solution is a way to control the number of ions within the nanocrystals.
In terms of the Gibbs free energy, the stability condition at constant pressure and temperature can be given by implying that for a stable system the free energy has to be convex as a function of the mole fraction y [0,1].The derivative of the chemical potential with respect to the mole fraction at constant pressure and temperature plays the same role as the derivative of the entropy with respect to temperature at constant volume and mole number.However, for small enough particles, it has been suggested in the literature that the molar Gibbs free energy or chemical potential has the form where L is the positive heat of solution, which is assumed constant, and the last term of Eq. ( 25) comes from the surface contribution.From Fig. 5(a), it is clear that the chemical potential of the adsorbing particles has a change of convexity for intermediate values of y, therefore having an unstable [yellow (light gray)] region.In this figure we also show that the distance between the equilibrium molar Gibbs free energy difference G m,eq or chemical potential and its corresponding nonequilibrium value G m , associated with phase coexistence along the Maxwell construction (blue dashed line), is minus the entropy production −T i S .If the system is originally at point y 0 (as occurs for the charging process), then the molar entropy produced can be written as where the factor (y − y 1 )/y corresponds to the fraction of Li ions in y 0 when the system is forced to be at state y by the bath.In similar form, the molar entropy produced when the system is originally at point y 1 is Adding both contributions the entropy production associated with the complete oscillation is Here, the external forcing is given by E ex = 1 j =0 μ eq y b (y bath − y j ) and the constant c by c ≡ μ eq ( y 0 +y 1 y 2 b ).Both expressions follow after considering the leading contribution of the last term of Eqs. ( 26) and (27).In a way similar to Sec.II, the last term in Eq. ( 28) constitutes an alternating external forcing that comes from the bath, which imposes a given value of ion concentration in the solution y bath .
Taking now the time derivative of Eq. ( 28) we obtain the entropy production per unit time where, according to nonequilibrium thermodynamics [37], we introduced the time dependent energy exchange rate with the bath, dE ex /dt ∝ dy/dt, that can be modeled through an oscillating function.From the last result, we may derive the following dynamical equation for the ion concentration in the electrodes: where τ is a characteristic time and I bath = dE ex /dt.The solution of the corresponding Fokker-Planck equation for the probability density P (y,t) also shows the oscillating behavior in Fig. 5.

VI. DISCUSSION
Here, our goal has been to analyze how energy barriers corresponding to a nonconcave free energy as a function of intensive parameters affect the thermodynamic properties and the dynamic behavior of small systems.We follow a global approach starting from the partition function of the whole system which makes no distinction between bulk and surface, since this characterization may be not clear in all situations.Hence, we assume that the small system can be characterized as an effective homogeneous system whose intensive parameters differ from those of the corresponding macroscopic system.In this sense, we coincide with the idea underlying Hill's thermodynamic theory of small systems.
In addition, this approach allows us to take advantage of Onsager's theory of irreversible processes when analyzing the dynamics.As a result of this description, peculiar behaviors of small systems can be predicted that underly the dynamics of energy generating, storage, or conversion nanodevices.
To go into further detail, we have examined the thermodynamic stability for extensive and small nonextensive systems using thermodynamic stability theory and elementary statistical mechanics.We found that under suitable external conditions, small systems do not break down into thermodynamic phases with different aggregation states as macroscopic extensive systems do, or even states with different structural configurations.Instead, since for small systems this process of division may have an excessively high energetic cost, they can only dissipate energy by performing critical oscillations between two metastable states [3] or structural configurations [2].
This behavior is supported by our own numerical analysis (see Figs. 2 and 3) as well as by the results previously reported in the literature of small systems [2,3,8,9,11,14,15].Similar oscillating dynamics has been observed in singleand many-particle energy storage systems [7][8][9] whose state is characterized by the number of lithium atoms stored in electrode nanoparticles.This dynamics is described by means of an equation analogous to our Fokker-Planck equation (23).Oscillations of the probability density of the stored energy are also obtained and shown in Fig. 5; see also [9].
Our results allows us to conclude that the cyclical dynamics described here constitutes a general nonequilibrium behavior of small systems which may be convenient to take into account in the design of energy exchange and storage nanodevices.Other recent advances in nanotechnology [35,36] suggest more potential applications along this line.Further contributions as to the implementation of our theory in practical situations are currently in progress.
No. FIS2011-22603.We also thank the Academic Mobility Program between the University of Barcelona and the National Autonomous University of Mexico.We would also like to thank Iván Latella and J. López Alamilla for their valuable collaboration.I.S.H. acknowledges the hospitality of Prof. J. M. Rubi and the Department of Física Fonamental of the University of Barcelona where this work was initiated during a sabbatical leave.

FIG. 1 .
FIG. 1. Schematic representation of the reduced free-energy change (solid line) and the free-energy bound (solid straight line) as a function of the inverse temperature β. β 0 and β 1 represent two metastable states separated by the free-energy barrier.The dotted lines indicate the free energies corresponding to the new phases if these appear.

FIG. 2 .
FIG. 2. (Color online) Deterministic (red solid lines) and stochastic (blue irregular points) oscillations of the order parameter φ = T in the presence of a noisy heat exchange with the heat bath.In all the simulations we have used normalized variables with a = 0.24, b = 0.7, and a Gaussian random force with amplitude D 1/2 = 0.44.

FIG. 3 .
FIG.3.(Color online) Snapshots of the occupation probability P ( ,t i ) as a function of the order parameter for different reduced times t i .The maximum of the probability oscillates between states I 1.2 and I I 1.8.Panel (a) shows the initial distribution and its initial evolution for three different times.Panel (b) shows intermediate times, whereas (c) shows the behavior of the probability near the second state.The state of the system spends different times in states I and I I because of the asymmetry of the potential, favoring state I .The time spent at the middle of the oscillation 1.5 is much lower than the time spent at the extrema.Thus, an apparent stationary bimodal distribution can be obtained by averaging P ( ,t) over a window of time larger than the period of the oscillations.Panel (d) shows three different reduced time averages P ( ,t) t of the oscillating probability distribution P ( ,t).The red dashed line represents an average performed on the time interval

FIG. 4 .
FIG.4.Occupation probability P ( ,t i ) as a function of the order parameter and different reduced times t i .The maximum of the probability oscillates between states I 1.1 and I I 1.8.Black solid lines correspond to states with maximum around I whereas the gray solid lines to states with maximum at I I .Light gray dotted lines correspond to crossing states.The diagonal lines correspond to states I and I I .

FIG. 5 .
FIG. 5. (Color online) (a) Chemical potential of the Li x FePO 4 system at the electrode.(b) One realization of the oscillating behavior with initial condition at y = 0.05 for a small oscillation amplitude and one realization of a process affected by noise.The right panel is the three-dimensional solution of the Fokker-Planck equation for the probability density P (y,t).

Figure 5 (
b) and the right panel show the dynamics of the filling fraction y as a function of time for a small amplitude of the oscillation.One realization of the stochastic case (irregular points) is also shown.Noise may induce different dynamics associated with stochastic resonance phenomena.