Thermodynamics of Small Magnetic Particles

In the present paper, we discuss the interpretation of some of the results of the thermodynamics in the case of very small systems. Most of the usual statistical physics is done for systems with a huge number of elements in what is called the thermodynamic limit, but not all of the approximations done for those conditions can be extended to all properties in the case of objects with less than a thousand elements. The starting point is the Ising model in two dimensions (2D) where an analytic solution exits, which allows validating the numerical techniques used in the present article. From there on, we introduce several variations bearing in mind the small systems such as the nanoscopic or even subnanoscopic particles, which are nowadays produced for several applications. Magnetization is the main property investigated aimed for two singular possible devices. The size of the systems (number of magnetic sites) is decreased so as to appreciate the departure from the results valid in the thermodynamic limit; periodic boundary conditions are eliminated to approach the reality of small particles; 1D, 2D and 3D systems are examined to appreciate the differences established by dimensionality is this small world; upon diluting the lattices, the effect of coordination number (bonding) is also explored; since the 2D Ising model is equivalent to the clock model with q = 2 degrees of freedom, we combine previous results with the supplementary degrees of freedom coming from the variation of q up to q = 20. Most of the previous results are numeric; however, for the case of a very small system, we obtain the exact partition function to compare with the conclusions coming from our numerical results. Conclusions can be summarized in the following way: the laws of thermodynamics remain the same, but the interpretation of the results, averages and numerical treatments need special care for systems with less than about a thousand constituents, and this might need to be adapted for different properties or devices.


Introduction
Recent developments in experimental techniques and refinements of synthesis methods allow for direct measurements of very small systems [1][2][3][4][5][6].Examples of such systems include magnetic domains in ferromagnets, which are typically smaller than 300 nm, quantum dots and solid-like clusters that are important in the relaxation of glassy systems and whose dimensions are a few nanometers.Due to the small size of such systems and the high surface/volume ratio, the usual thermodynamic properties of macroscopic systems do not apply in the usual sense, and the analysis of fundamental principles of statistical thermodynamics needs to be reconsidered [7][8][9][10].Phase transitions, specific heat, thermal conductivity, magnetization and other phenomena need careful examination of the corresponding thermodynamic variables used to characterize them.Ergodicity, fluctuations and error bars are among the aspects that are greatly affected by the small size of the systems.
The increasing interest in the theory of small systems is not limited to physics, but also extends to biological and chemical systems, such as clusters [11,12], thin films [13] and biological molecular machines with sizes under 100 nm [14].This broad picture calls attention to the extension of the mathematical treatments valid in the (usually called) thermodynamic limit (TL) to systems formed by a small number of elements.
In 1962, Hill [15,16] developed a formalism for describing the thermodynamics of small systems, where thermodynamic properties such as enthalpy and Gibbs energy are no longer extensive.The pioneer work of Hill established the fundamental laws that govern finite-size effects in thermodynamics.Later, Hill and Chamberlin [17][18][19][20] introduced the term nanothermodynamics to denote the study of systems far away from the TL.Despite these results, there are still many open questions regarding the thermodynamic behavior and possible applications of very small systems.
Magnetic nanoparticles [21,22] are also examples of small systems whose properties are influenced by finite size effects.The competition between core and surface effects determines the state of the particle, which can go from a paramagnetic state to a ferromagnetic state.The temperature characterizing this transition is the well-known Curie temperature (T C ) in the TL.Experimental studies for different materials show that the T C in nanoparticles is not a unique property of the material (as in the bulk), since it is strongly size dependent and decreases as the characteristic dimension of the particle, V 1/3 (where V is the volume of the particle), is reduced [23][24][25].
In a recent work [26], it has been shown that for tiny systems, not only surface sites are special, but also some inner sites might not be equivalent to each other.This fact can have important consequences on the thermodynamic properties of small nanoparticles (up to sizes of the order of ≈10 nm).In the magnetic systems described below, this effect mainly introduces dispersion in the exchange coupling constants, which in turn may affect the values of the thermodynamical parameters.These considerations tend to lose importance as the nanoparticle size increases.Moreover, they do not affect the overall underlying physics, but impose limitations in the interpretation of the magnitudes defined in the TL.
In the present paper, we would like to answer the following questions: How stable and reliable are devices built of subnanoscopic magnetic particles?How do the representative properties of these devices depend on the size, dimensionality, coordination and internal degrees of freedom among other variables?We will consider two examples offered by very small magnetic particles responding to temperature changes; such particles could be components of devices.
Let us first consider the magnetic switch shown in Figure 1a where a tiny nanomagnet (nm) sticks to a ferromagnetic material (F) under a sensitive temperature T * (a way to avoid calling it the critical temperature).To overcome the difficulties of defining temperature for very small isolated systems we will assume that the devices and everything that is connected to them are immersed in a thermal bath at temperature T. Changes in temperature are smooth and slow, so there is time to equilibrate.As temperature raises, nm can decrease its magnetization under a critical value, which will weaken the attraction to F; then an elastic restoring force (like the spring in the figure) brings nm to its natural position away from F. This magnetic switch could activate an alarm that goes off over a prefixed temperature, which ideally is precisely T * .Let us assume that for this system, mechanical relaxation times are longer than electronic ones; then, such a device is not sensitive to magnetization reversals, so it is basically controlled by the absolute value of the magnetization, namely |m(T)|.
Figure 1b presents a tiny rod-like nanomagnet, which can rotate about its center on the plane of the figure.A stronger fixed magnet (NS) with high T C orientates the tiny magnet.When the temperature goes over a sensitive temperature T * * , a magnetization reversal happens in the small magnet.Such polarity reversion triggers a rotation of the small particle disconnecting the terminals that are connected when the conducting element attached to the "black" part of the particle closes the circuit.This device is strongly dependent on the polarity of the nanoparticle, namely on the sign of m(T).Magnetization reversals are usually quite fast as compared to most mechanical processes in materials.Some orders of magnitude for the switching time (defined as the time in which the magnetization reverses polarity) have been proposed at about 200 ps for small Co systems [27].The switching time could be larger depending on the system size, temperature and magnetic anisotropies, but a few tenths of ns have been reported for other small systems in the absence of an external magnetic field [28].We shall come back to these examples later on since this is the kind of system we have in mind, rather than statistical averages over different samples or instances.At this stage, we point to the gross picture, namely the tendency of properties, and not to material-related or geometrical features controlling a particular device.To focus the discussion, we turn our attention to two general and well-known models whose properties can be reexamined in this new scenario.We refer to the Ising model and to the clock model (related to the better-known Potts model).
The Ising model plays a central role in the study of phase transitions in magnetic systems (transitions from paramagnet to ferromagnet, order-disorder transitions, etc.) [29][30][31][32][33][34][35][36][37].Ising [38] gave an exact solution to the one-dimensional lattice problem in 1925.All other cases are expressed in terms of series solution [39,40], except for the special case of two-dimensional lattices in zero external field, which was exactly solved by Onsager [41] in 1944.Close approximate solutions in dimensions higher than one can be obtained, and the two most important of these are the mean-field approach [42] and the quasi-chemical approximation [42,43].
In the last few years, the Ising model has been adopted as a prototype for studying the magnetic properties of systems far away from the TL [44][45][46].In [44], Velásquez et al. studied the magnetic properties and pseudocritical behavior of ferromagnetic pure and metallic nanoparticles in the presence of atomic disorder, dilution and competing interactions.Using variational and Monte Carlo methods, the authors reported the particle size dependence of the ordering temperature for pure and random diluted particles.
Bertoldi et al. [45] solved analytically the Ising model for small clusters in the microcanonical formalism, using the mean-field theory and considering the presence of surfaces.In agreement with the work of Velásquez et al. [44], the authors found that the critical temperature decreases proportionally to the inverse of the radius of the cluster.Recently [46], the Ising model on a simple-cubic lattice was used to analyze the main concepts of nanothermodynamics applied to Monte Carlo simulations.
Finally, we increase the internal degrees of freedom at each site upon generalizing the Ising model to the discretized XY model in what is called the clock model which is a well-known example of the famous Berezinskii-Kosterlitz-Thouless (BKT) transition [47,48].In this way, we will be able to appreciate what is the role of the increase of the configuration space due to the characteristics of the constituents of small systems and not only the fact that they interact among themselves.
This paper is organized as follows: In Section 2, we present the system and the methodology based both on theoretical formulations and on computer simulations.In Section 3, we present the results stressing those points where the analysis concerning the TL is not applicable for very small particles; the discussion and relations among results are done along with the presentation of the results given by plots and particular numerical values.In Section 4, we list and comment on the main conclusions.

General Definitions
Let us begin by considering the Ising model on a square lattice of dimensions L × L = N, where local magnetic moment or "spin" S i at site i can point either along the positive (+1) or negative direction (−1) of a given axis.The isotropic Hamiltonian for such a system can be written as: where J is the exchange interaction to nearest neighbors; the sum extends to all pairs through the lattice, which is indicated by the symbol i > j under the summation symbol.The magnetization m(T, t) of the system at any temperature T and time t is simply obtained by just adding the spins in an algebraic manner, namely, where we have normalized it over the number of sites N, which is equivalent to normalizing it with respect to the saturation magnetization.Large macroscopic systems are modeled by making use of periodic boundary conditions (PBC) where surface states are ignored; in this approximation, Lars Onsager [41] found an analytic solution for the magnetization of the system in 2 dimensions (2D): where k B is the Boltzmann constant, which we can take as 1.0 from now on (this is equivalent to measuring energy and J in T units).The critical temperature T C can be readily obtained at the point where the magnetization vanishes: Then, the exact numeric solution for this problem for PBC can be easily evaluated, yielding: where J has been omitted, as will be done from now on.The reason to have chosen the 2D Ising model to begin with is precisely this result, namely the existence of an analytic solution for the transition temperature in the TL.We can use it as a reference point to compare with other results and simulations where different conditions will be varied.

Exact Theoretical Approach for a Small System
Next, we consider a very small Ising square lattice, L = 4, where we can calculate everything exactly.The partition function can be expressed as: where the coefficients c n give all of the possible spin configurations compatible with energy E n according to the Ising Hamiltonian; λ is the number of energy levels.
We can write the exact partition functions Z P (T) and Z F (T) for PBC and free boundary conditions (FBC), respectively: where coefficients c Pn and c Fn are to be obtained from Tables 1 and 2, respectively.These tables list the corresponding degeneracies p(n, k) (PBC) and f (n, k) (FBC) for the set of states sharing energy E n given by Equation ( 1) and magnetization M n,k : which runs between −16 and +16 at intervals of 2, as shown by the second row in the tables.
The coefficients in the expansions given by Equations ( 7) and ( 8) can be obtained as: where we take into consideration the symmetry in Tables 1 and 2, namely p(n, It should be noticed that Tables 1 and 2 give only those configurations with nil or positive magnetization.However, the Hamiltonian above is symmetric under the simultaneous reversion of all spins, so both tables should be continued to the left to include the states with negative magnetization. Table 1.Coefficients p(n, k) for periodic boundary condition (PBC) configurations of a given energy E n and magnetization M n,k .The first column gives n; the second column gives the corresponding energy E n ; the following columns give the number of configurations for that energy (row) and magnetization (column).This table is symmetric and should be extended to the left by means of the symmetry properties so as to include the 8 negative values of the magnetization following Equations ( 10) and (12).

Numerical Simulations
In addition to some theoretical calculations, most of the work will deal with numerical calculations based on Monte Carlo (MC) simulations.A square lattice L × L is chosen; PBC or FBC are imposed; a site is randomly visited; and the energy cost δ of flipping the corresponding spin is calculated: if the energy is lowered, the change of orientation is accepted; otherwise, only when exp(−δ/T) ≤ r, the spin flip is accepted, where r is a freshly-generated random number in the range [0,1].This is the usual Metropolis algorithm.A Monte Carlo step (MCS) is reached after N = L × L spin-flip attempts.One of the main goals here is to report the sensitive temperatures (T * or T * * analogues to critical temperatures in the TL) for different systems.
The energy is computed according to the Hamiltonian given by Equation ( 1).The instantaneous magnetization at time t and temperature T is obtained from Equation (2).
For each temperature, 2τ MCSs are performed: the first τ MCS is used to equilibrate at a fixed temperature T, while the second τ MCS is used to measure the magnetization every 20 MCSs, reaching a total of 50,000 measurements.Unless specified in a different way τ = 10 6 in the rest of the paper; this τ value gives stable results and leads to coincidence with the analytic expressions obtained as described above.

The Clock Model
This is the discrete version of the famous 2D XY model, which is probably the most extensively studied example showing the Berezinskii-Kosterlitz-Thouless (BKT) transition [47,48].It is often used as a reference of its peculiar critical behavior at the transition point and universal features [49][50][51][52][53]. Instead of the exclusion of an explicit continuous symmetry essential for the BKT transition, it can also emerge from a system without explicit continuous symmetry.The q-state clock model is the discrete version of the XY model, with discrete symmetry.The Hamiltonian of the clock model can be written as: where J > 0 is the ferromagnetic coupling connecting pairs of nearest neighbors i and j; the discrete angle between the spin orientations is given by θ i,j (η, q) = 2πη/q, for η = {0, 1, .., q − 1}.While the exact XY model is recovered only in the limit of infinite q, it has been found that the BKT characters appear in the clock models when q ≥ 5 [54][55][56][57].The nature of the phase transitions in the general clock model has been widely studied with different theoretical and numerical approaches, which however have given mixed results on the characterization of transitions at the lower bound of q (for instance, see the summary of the related debates in [58]).
In spite of all that, in this work, we want to use this model as a generalization of the Ising model presented above allowing for the increase in the number of internal states at each site.All of this is always aimed at the behavior of small systems.We will focus on a square lattice 8 × 8 with FBC, with the clock model ranging from q = 2 (Ising model) to q = 20 (one discrete version of the 2D finite XY model).

Results and Discussion
Let us begin by considering the theoretical approach for the Ising lattice with L = 4 introduced in the previous section.The magnetization per site can be calculated as: and: both of which give a nil result for any T since negative and positive contributions to the magnetization cancel out as they are weighted with the same probability due to the symmetry properties of Tables 1 and 2.
In macroscopic systems, ergodicity is spontaneously broken favoring any of the two halves of the configuration space.As we will see below, this is not always true for very small systems, but let us start from the ergodicity breaking approach so widely invoked to approach the TL.This can be achieved by either considering only the positive magnetization contributions of Tables 1 and 2 plus half of the coefficients for the M n,k = 0 column of the corresponding table.Equivalently, we can force ergodicity breaking by considering all of the contributions to the magnetization, but using their absolute value, namely, and: where e refers to the forced ergodicity breaking.
We invoke now the numeric simulations defined in the previous section to calculate the magnetization for the same system.Here, the same trick is done to break ergodicity, namely it is always |m(T, t)| (obtained from Equation ( 2)) that is considered as the contribution to magnetization from every configuration visited.In this way, ergodicity breaking is forced to occur in any circumstance.We will call µ Pe (T, t) (PBC) and µ Fe (T, t) (FBC) the magnetizations obtained via numeric simulations for a temperature T and at an instant t.
We compare now previous approaches.Results from Equations ( 16) and ( 17) are shown by up and down solid triangles respectively in Figure 2. The results for the numeric simulations for PBC and FBC are shown here by empty diamonds and empty squares respectively.The values for < µ Pe (T) > and < µ Fe (T) > reported in Figure 2 correspond to positive averages over ν = 50,000 measurements.The corresponding distribution of values allows us to report the error bars for each temperature.17) (down yellow triangles) to compare with the results obtained from numeric simulations using PBC (empty magenta diamonds) and FBC (empty black squares); the Onsager solution is also plotted to mark the expected T C in the thermodynamic limit (TL) (solid red circles).
For comparison and discussion, the Onsager solution is also added to this figure by means of solid circles.
Figure 2 deserves a careful analysis and discussion.The excellent agreement between theory and numeric simulation is very impressive.However, this should not be very surprising since the Metropolis algorithm is based on the same probability laws that appear in the partition function, hence, in the magnetization as well; both approaches consider ergodicity breaking, as already pointed out.The agreement is stressed in the sense that the computer algorithms used here are very reliable, and the equilibration times are enough.It can be noticed that both the theoretical and numeric results for PBC coincide with the Onsager solution up to almost T = 2.0.From there on, both treatments overestimate the magnetization values partly due to a small size effect, but also due to the forced ergodicity breaking, which is assumed in these calculations.We will come back to this point below.
Let us now pay attention to the different results obtained for FBC as compared to PBC regardless of the method employed in the calculation.It is PBC that can be used to compare with the results valid in the TL, which are illustrated in Figure 2, by the Onsager solution yielding T C ≈ 2.27.As can be seen, the theoretical results and magnetic simulations using PBC both agree with the Onsager solution up to T = 2.0 approximately; from there on, both solutions remain at high values, decreasing smoothly without an indication for a sharp phase transition.In the case of very small particles, PBC are unthinkable, since in this limit, surface states become very important.This is an indication that FBC should be used to model the properties of very small particles.
The light and smooth decrease in magnetization due to surface effects will determine that the nanomagnets involved in the simple systems of Figure 1 will have different responses as those made of the same basic material, but operating in the bulk sizes, namely in the TL.The sensitive temperatures T * and T * * for which these devices will activate still depend on other features to be examined below, but this is a first result showing that properties characterized in the bulk need to be reexamined for very tiny particles.
Actually, we can define the corresponding sensitive temperature simply as the temperature for which m(T) or < µ(T) > (with any subscripts) takes the value 0.5.This working definition is consistent with the TL, where a square step magnetization curve is expected at T C (equivalent to our sensitive temperature), jumping from magnetization 1.0 to 0.0.Then, the sensitive temperatures can be read from Figure 2 as 2.4 and 3.4 for FBC and PBC, respectively.Both values are larger than the macroscopic expected result, which requires further discussion.
What is the real meaning of the average values we have just reported?Let us revisit the FBC data and do a different exercise: let us plot the instantaneous value of the magnetization µ F (T, t) (no ergodicity breaking here) after the τ instants of measurement (last measured value for each temperature after 2τ total MCSs).This is plotted by solid squares (black) in Figure 3 using the same data already available from FBC on the 4 × 4 lattice.These results allow for just one interpretation: at least at T = 0.8 (or eventually at even lower T values), a magnetization reversal occurred, and ergodicity breaking is no longer valid.Therefore, the previous approach intended for the TL has to be reconsidered to analyze the results obtained in the small world.This is an indication that the ergodicity breaking assumption is valid only for systems over a minimum size and appropriate for the property under investigation.
All of the previous procedures are usual for TL analysis.However, the magnetic switches presented in the Introduction will respond instantaneously and without having the chance of performing derivatives, averages or even more sophisticated mathematical analyses, like the crossing of Binder cumulants [59] or the crossing of autocorrelation functions [60], to decide which is the appropriate value for the sensitive temperature.In these examples and in all nanosized devices, a careful analysis of the thermodynamics of the object in terms of the desired property should be appropriately done.
In the case of the magnetic switch illustrated by Figure 1a, the nanomagnet (nm) will stick to the ferromagnetic material (F) regardless of its polarity.This would be the case of systems whose magnetization reversals are much faster than any mechanical process (like elastic forces) usually involving atomic masses.Therefore, such a device would be insensitive to magnetization reversals and will respond effectively to the magnitude of the magnetization of the nm switching off the circuit at a certain sensitive temperature T * .
On the other hand, in the case of the magnetic switch illustrated by Figure 1b, the magnetic nanorod will always point with its north pole to the fixed magnet on top, which coincides with the black portion of the rod in the upper part of this figure at low T values.As temperature rises, a magnetization reversal is possible in the softer magnetic material in the nanorod; the north pole is now on the white part of the rod, and the nanomagnet will rotate with respect to its fixed center.Therefore, this device is sensitive to the polarity of the nanomagnet and will disconnect the circuit at a certain sensitive temperature T * * .
Even if the nanomagnets in both switches are made of the same material, T * and T * * will be different in general because they originate from different properties: the former is governed by the value of the absolute magnetization, and the latter responds to the reversal process depending on the instantaneous sign of the magnetization.As we will show below, they are different indeed, and for the systems to be studied, T * > T * * .
Figure 3 also yields the average of the actual instantaneous magnetization values obtained for FBC < µ F (T) > over the ν = 50,000 instants, namely, These results are reported by means of empty circles (red) with the corresponding error bars.If we look at these data, we should estimate T * * F ≤ 0.8.By enforcing ergodicity breaking in Figures 2 and 4, we are biasing the data in several senses: magnetization reversals are inhibited as if there is an infinite anisotropy; the instability of the systems is ignored; the curve is artificially smoothed; the average value of the magnetization is always positive; and it does not converge to 0.0 as T → ∞.This last result is usually confused with the tails of the magnetization curves (and other order parameters) for T > T * due to finite size effects.given in the inset: 2.000 corresponds to PBC; 1.875 to FBC; 1.75 to the dilution of the lattice according to the decoration proposed in Figure 5 where four pieces × with = 4 have been withdrawn; 1.500 corresponds to a similar decoration with = 6; 1.00 corresponds to a similar decoration with = 7. B means P when κ = 2.000, and it represents F for all of the other cases.For clarity, not all curves are shown with error bars.
From the previous discussion, it is clear that a unique definition for a Curie temperature in the nanoworld faces huge difficulties.At least the usual definition of a Curie temperature, where all magnetic properties change at one precise temperature, is not longer sustainable.Actually, in the TL, there would not be a difference between the two switches illustrated in Figure 1.However, we have just shown that for a magnetic systems of just very few interacting elements, different properties may have different sensitive temperatures.
The general message here is that when dealing with small particles, instantaneous responses may have a great importance and could determine completely the performance of the device.The details of this behavior will depend on the properties involved, geometry and other characteristics.
In any case, what is quite clear is that < µ(T) > is not the only property governing the response of the device at very small sizes.The large error bars in Figures 2 and 4 indicate that instantaneous measurements of magnetization can be erratic and at times could disconnect the device at unexpected temperatures under the desired one.As sizes get larger, error bars progressively decrease, and < µ(T) > becomes the right parameter to design devices as those shown in Figure 1.
Another important result shown by Figure 2 is the huge size of the error bars, which are likely associated with the small size of our system.To test this hypothesis, we will do numerical simulations similar to those of Figure 2 using both PBC and FBC for L = 16.Results (along with others to be defined immediately below) are presented in Figure 4 by empty squares and solid circles, respectively.In fact, error bars are now smaller; T * has moved to lower values; the slope in the descent is more vertical (approaching the square step shape expected at the TL); and the remnant absolute magnetization towards high temperatures has decreased.Therefore, it is clear that we are moving in the direction where these simulations designed for the TL are effective.
We can make use of this opportunity to discuss another important consideration for small particle behavior: coordination number (or bonding, to put it in a more general context).Let us define effective coordination number κ as the ratio of the number of bonds B (nearest neighbors in our case) over the number of the N individuals (spins) interconnected by these B bonds.For any square lattice with PBC, B = 2N and κ (2) P = 2.0; the upper index (2) refers to 2D.In the case of FBC, we have B = 2L 2 − 2L; N = L 2 , so: which for L = 16, yields κ (2) F = 1.875, as illustrated in the inset of Figure 4.It is clear that in the TL, κ (2) P , but for small particles, the difference might matter.The effective coordination number defined above can be diminished upon diluting the original lattice.This can be achieved by human decoration on the original lattice (just removing atoms in a prescribed way) or could be due to some natural process that takes away magnetic atoms of our 2D sample (corrosion for instance).For simplicity, we will have the first alternative in mind, but the results for the second one would only add randomness to what we will discuss next.
Let us go back to the ferromagnetic lattice 16 × 16, and we will progressively dilute it upon removing four square sectors × each time in the way shown in Figure 5 for the case = 4.To be more specific, the lattice is divided into four sectors (L/2) × (L/2), and the extraction of the 2 sites is done from the lower right corner of each sector.For the example under consideration, 192 and B = 336, so κ = 1.75.Then, we remove four sectors with = 6 and, finally, four sectors with = 7, obtaining κ = 1.5 and κ = 1.0, respectively.In all of these cases, FBC are imposed.The last case is quite interesting because it is equivalent to the linear chain used by Ernest Ising in his original (and only) paper [38], for which the expected T C is 0.0, namely no phase transition for finite temperatures.Results for < µ(T) > are given in Figure 4, where the values of κ are given in the inset.Error bars were also calculated for all κ values, but they are shown for κ = 2.000, 1.750 and 1.000 only for the simplicity of the figure.Direct comparison shows that error bars for the lattice 16 × 16 are less pronounced than the corresponding ones for the lattice 4 × 4.
The main result of Figure 4 is to show the way T * decreases as the effective coordination number decreases.This is even more so if ergodicity breaking is imposed, an effect that we will discuss again below.
Let us stress another shortcoming of the forced ergodicity breaking method mentioned before, namely the fact that < µ(T) > remains positive at high temperatures.In Figure 2, the last point of the simulation corresponds to < µ(4.0) > ≈ 0.32 with error bars of about 0.23 for FBC.For this lattice, we can actually make use of Equation ( 17) to obtain the exact value in the high T limit: < µ(T → ∞) >→ 0.196, which indicates that temperature has to be increased in the numeric simulations to attain a better comparison.We can follow this tendency simulating higher temperatures including larger lattices.Thus, Figure 6 shows results for L = 4, 8 and 16, where the positive asymptotic behavior is clearly shown; the values for < µ(40.0)> are 0.204, 0.104 and 0.052, respectively.The first result is approaching the already quoted theoretical result 0.196.The set of values indicates that this sort of imposed remnant magnetization tends to be negligible for large lattices in the TL.A different way to put this is that magnetization reversals decrease as L increases, which makes the imposed ergodicity breaking coincide with the real ergodicity breaking.This tendency is shown in the inset of Figure 6 where in the logarithmic scale, we plot the number of magnetization reversals during the 50,000 measurements as a function of 1/L for T = 1.5.For higher temperatures, reversals are easier, but the amplitudes are lower; for lower temperatures, reversals are harder, but the amplitudes are higher (which makes the absolute value of µ(T, t) stay close to unity).When ergodicity breaking is imposed at low temperatures, magnetization always oscillates between values close to +1 or −1.The curve for κ = 1.000 in Figure 4 deserves a special analysis since it is expected that the magnetization should go down as soon at T > 0. However, this figure is for < µ(T) >, namely for forced ergodicity breaking, to compare with Figure 2, related to the analytic expressions given by Equations ( 16) and (17).If we go back to the data for this case, we can study each temperature to find that reversions are occurring at all temperatures.For instance, < µ(0.6) > ≈ 0.7 for κ = 1.000 in Figure 4; however, if we look at the instantaneous data for the 100 leading intervals of the measuring interval (immediately after equilibration), we find µ(0.6,t) oscillating as given in Figure 7.For lower temperatures, oscillations are progressively less and less frequent approaching alternations between +1.0 and −1.0only at T = 0.1, where it is necessary to wait a huge amount of time to get the first reversal.In this way, we approach the result for the 1D Ising lattice: no phase transition is possible at a finite temperature.
Let us now investigate the role of the effective coordination number towards the 3D case with PBC (κ P = 3.0).This is shown in Figure 8, where we begin by a single 4 × 4 layer with FBC (κ F = 1.875); then, we continue the series: two layers 4 × 4 with FBC (κ T * moves to higher temperatures as the effective coordination number increases.Moreover, the dimensionality plays also an important role, which can be appreciated when we compare the second case here having κ (3) F1 = 2.0, with the already seen case κ (2) P = 2.0 for the 4 × 4 PBC lattice (shown in Figure 4).It can be appreciated that for the 3D case, T * is larger than for the 2D case., and part of it was already included in Figure 2; it is given here for completeness and comparison purposes).Now, we show results for the q-states clock model simulated on an 8 × 8 square lattice.The density of magnetic states increases as q increases; therefore, we expect a q dependence of the sensitive temperature characterizing the transit from an ordered ferromagnetic initial regime to a magnetically-disordered phase.
Figure 9 shows the sensitive temperature T * vs. 1/q ranging from q = 2 to q = 20.Here, we follow the same approximation as in Figure 2, namely we force ergodic separation starting always from a ferromagnetic state.T * is defined as the temperature where the normalized magnetization is < µ e (T * ) >= 0.5, which is a practical way to define the transition to a disordered state.In order to sample with the same density the q-possible states for each spin in the 8 × 8 square lattice (there are q 64 states for an 8 × 8 lattice), we increase the number of MCSs as q increases; higher q values need larger numbers of MCS to reach equilibration.Therefore, we use τq/2 MCSs for the simulations of the q-clock model.Sensitive temperature T * from a magnetic ordered to a disordered phase for the q-states clock model as a function of 1/q.A square lattice 8 × 8 with FBC is used.Measurements are obtained after 5 × 10 6 MCSs for averaging.
As we see in Figure 9, there is a monotonous decrease of T * with increasing values of q.Larger values of q mean energy states separated by less energy differences, therefore moving the initial ferromagnetic state to other states with different magnetizations in an easier way: this implies that the disorder for higher q values will be reached at lower temperatures as compared with the starting q = 2 case (Ising model).
A very striking feature is shown by Figure 9: two regimes are clearly shown: for q = 2 to q = 4, and from q = 5 to q = 20 (and possible larger values).In the first regime, T * decreases with a relatively low slope as a function of 1/q, whereas for q > 4, the slope of the T * vs. 1/q curve is larger.This change of slope could be associated with different types of transitions that occur for different q values.Namely, there is an order → disorder transition for q <= 4 and a double transition order → quasi-ordered → disordered for q >= 5 as discussed in [58,61].
It is clear from symmetry arguments that as q increases, lower temperatures are sufficient to originate new ferromagnetic states in any of the q directions irrespective of the starting configuration.Therefore, as q tends to infinity, T * tends to zero.This is known as the Mermin-Wagner theorem [62].This theorem establishes that there is no magnetic state at any finite temperature (T > 0) in the XY model.Moreover the formation of vortex-antivortex states is prohibited in an 8 × 8 lattice because of its small size.Therefore, the BKT transition is also absent in samples of such a small size.

Conclusions
As magnetic particles get very tiny, comprising a low number of elements, especially at the nanoscopic level, the ergodic separation is no longer possible in the way that it is handled in the TL.This means that as temperature progressively increases, very soon, all of the configuration space becomes available, with positive and negative components of the magnetization.While the transit time from a state with one magnetization sign to the opposite one is faster than other characteristic times of the system (elastic forces, vibrations, etc.), the particle will present a magnetization whose magnitude could make it stick to a ferromagnetic material, as is shown in Figure 1a.However, as the temperature continues to increase, such a magnitude will weaken, changing the static magnetic properties of the particle.
For tiny magnetic particles, statistical fluctuations are larger than the needed reliability of the particle, which could be reflected by a given experimental accuracy.This fact is reflected in the error bars found in previous measurements suggesting that at any given temperature, oscillations are large, and a unique value of magnetization cannot be found at any temperature.This behavior requires extreme care when using magnetic particles at the molecular level being part of a device since the performance of this tiny magnet could be under what could be needed.If the dimensions available allow it, it would be safer to use slightly larger particles that already behave as macroscopic ones.The sensitive temperature at which the particle changes according to its magnitude is T * , which has to be thought of as a rather large range centered at a value.Towards the TL, such a range should disappear, and T * → T C .
Such considerations are even more important when the response of the particle is based on the orientation of the nanomagnet, as is proposed in Figure 1b.Here, a magnetization reversal (even if the magnitude remains high) can trigger a reorientation of the particle.Such magnetization reversals could happen at almost any temperature for small particles.The main conditions are that they become less and less frequent as temperature is lowered and size increases.The sensitive temperature at which the particle changes according to its polarity is T * * , which again has to be thought of as a range centered at a value.From previous results, it follows that T * * < T * , since reversals at low T tend to fluctuate between extreme magnetization values.Towards the TL, we can expect that T * * → T * → T C .
Previous results show that a definition for the Curie temperature for small systems, even if a defying and interesting task, does not add any new physical insight for what will be the real response of a single nanosized ferromagnet of any given material.The different characteristics of the tiny system contribute to the fluctuations of its magnetization that in the last run dominate its response.Alternatively, the Curie temperature is a thermodynamic parameter arising from ensemble averages; in this way, T C finds a unique definition in the TL.
However, even some of the theoretical treatments intended for the TL need a critical reformulation.One of them is the fact that the high temperature sector of the magnetization curve is lifted since all magnetizations are taken as positive when ergodic separation is assumed and imposed.Then, critical temperatures tend to be overestimated.In addition, the zero magnetization value is not reachable for finite sizes and finite temperatures as presented in Figure 6 since the tails are forced to be always on the positive side of the magnetization; so not always are the long asymptotic positive tails exclusively due to finite size effects, as is frequently assumed.
In the case of very small particles, simulations cannot impose PBC since this would ignore the surface states, which become more important for these systems.The energy spectrum of surface elements means a more compact and homogeneously-distributed density of states for the same size (see Tables 1 and 2).In terms of the magnetization of the system, simulations with FBC lead to softer magnetism (lower T * ) than simulations for a system of the same size, but with PBC.
Dilution of the lattice upon remotion of magnetic sites (intended or accidental) will lead to an increase of the surface, thus lowering the effective coordination number κ. Materials get softer, namely sensitive temperatures move to lower values, as κ decreases.This phenomenon has still another implication since lower effective coordination numbers also mean higher chances of magnetization reversal.Thus, if magnetic tiny particles lose magnetic mass through corrosion, the sensitive temperatures could move to lower values due to aging.Previous effects were tested through several κ values ranging from 1.0 through 3.0 (namely from one to three dimensions).This allowed one to appreciate that for similar κ values, materials are harder in the higher dimensionality realization.
The instabilities due to magnetization reversal events decrease very rapidly with size.In the case of a 16 × 16 lattice, FBC, at T = 1.5, the ergodicity breaking was violated just once for 50,000 instantaneous measurements.This result seems to indicate that for lattices over 32 × 32 (about 10 3 magnetic sites), magnetization reversals can be progressively ignored for most of the practical applications.
When the density of states gets denser due to more internal degrees of freedom, the sensitive temperatures decrease.This was tested with the clock model varying q, the number of independent orientations of the spin on the 2D lattice.Actually, as q increases, T * tends to approach 0.0 in the expected Mermin-Wagner theorem valid for q → ∞.It is quite interesting that this result obtained for macroscopic systems is also obtained for a lattice as small as an 8 × 8 lattice with FBC.
Finally, we want to stress that in the present study, we have not found any evidence for failure of any of the basic thermodynamic principles, which seem to be very robust.However, we have presented evidence indicating that some approximations done in the TL cannot a priori be applied to magnetic particles with a thousand elements or less.Care has to be paid to the way averages are performed or even if taking averages to describe instant responses of the system makes any sense; actually, such a small magnetic particle could exhibit properties associated with its proper fluctuations.

B
number of bonds in the lattice.c Fn degeneracy of the n-th energy level for FBC, (8).c Pn degeneracy of the n-th energy level for PBC, (7).
δ energy cost due to the possible flip of the visited spin.e subindex to represent that ergodic separation has been enforced.E n energy on the n-th level, (6).f n,k degeneracy of the k-th configuration within the n-th energy level for FBC, (11).F denotes free boundary conditions (FBC).

H
Hamiltonian of a systems of N spins, (1).i or j position of spin S i or S j , (1).J exchange interaction (J = 1), (1).
κ effective coordination number, (19).L side of the square lattice L × L. m(T, t) magnetization for temperature T and at time t, (2).m O (T) magnetization for Onsager solution at temperature T, (4).m F (T) Magnetization at temperature T for FBC from analytic expressions.m P (T, t) Magnetization at temperature T for PBC from analytic expressions.M n,k magnetization for the k-th configuration within the n-th energy level, (9).MCS Monte Carlo step (plural: MCSs).
µ F (T, t) Magnetization at instant t and temperature T for FBC from numeric calculations.
ν number of measurements done within the interval of 10 6 MCCs.p n,k degeneracy of the k-th configuration within the n-th energy level for PBC, (10).Sensitive temperature depending on the magnitude of the magnetization.T * * Sensitive temperature depending on the sign of the magnetization.

Figure 1 .
Figure 1.Magnetic switches.(a) Top: Under T * the magnetic nanoparticle (nm) sticks to the ferromagnetic material (F) even against a restoring elastic force represented by a zig-zag line.Bottom: As temperature T goes over T * , nm loses enough magnetization, so the elastic force dominates switching off the contact; (b) Top: Under T * * , a weak magnet labeled NS, with high enough T C , attracts the north pole (painted black) of a rod-like nanomagnet made of a material with lower T C ; the black part has a conducting element that closes a circuit.Bottom: Over T * * , an internal magnetization reversal occurs in the nanomagnet; the north pole is now in the sector painted white; the rod rotates with respect to an axis going through the center of the rod and perpendicular to the figure; in this new position, the circuit is switched off.

Figure 2 .
Figure 2. Magnetization for a lattice L = 4 obtained using the analytic expressions (16) (up blue triangles) and (17) (down yellow triangles) to compare with the results obtained from numeric simulations using PBC (empty magenta diamonds) and FBC (empty black squares); the Onsager solution is also plotted to mark the expected T C in the thermodynamic limit (TL) (solid red circles).

Figure 4 .
Figure 4. Simulations on 16 × 16 lattices for different effective coordination numbers κ whose values aregiven in the inset: 2.000 corresponds to PBC; 1.875 to FBC; 1.75 to the dilution of the lattice according to the decoration proposed in Figure5where four pieces × with = 4 have been withdrawn; 1.500 corresponds to a similar decoration with = 6; 1.00 corresponds to a similar decoration with = 7. B means P when κ = 2.000, and it represents F for all of the other cases.For clarity, not all curves are shown with error bars.

Figure 5 .
Figure 5. Diluted 16 × 16 lattice upon decoration: four sectors 4 × 4 are removed to leave N = 192 spins with a total of B = 336 bonds; FBC are imposed to get κ = 1.750.Other similar decorations are defined in the text and in the previous figure.

Figure 6 .
Figure 6.Average absolute magnetization for lattices of different sizes showing the tails for the artificial remnant magnetization induced by forced ergodicity breaking.The inset shows the decrease with the size of the number of magnetization reversal within 50,000 measurements after equilibration.For clarity, error bars are given for L = 16 only.

Figure 9 .
Figure 9. Sensitive temperature T * from a magnetic ordered to a disordered phase for the q-states clock model as a function of 1/q.A square lattice 8 × 8 with FBC is used.Measurements are obtained after 5 × 10 6 MCSs for averaging.

Table 2 .
Coefficients (11), k) for free boundary condition (FBC) configurations of a given energy E n and magnetization M k .The first column gives n; the second column gives the corresponding energy E n ; the following columns give the number of configurations for that energy (row) and magnetization (column).This table is symmetric and should be extended to the left so as to include the 8 negative values of the magnetization following Equations (9) and(11).