Thermal properties of open-shell metal clusters

Thermal properties of open-shell metallic nanoparticles of elements that are not magnetic in the bulk are studied within the framework of the local spin density approximation at finite temperature. For bulk ferromagnetic materials, these nanosized objects undergo a second-order critical phenomenon at a temperature that strongly depends on electron correlations and the size of the nanoparticle. The magnetic susceptibility is found to obey the Curie–Weiss law, which is in agreement with recent experimental measurements performed on mass-selected platinum clusters.


Introduction
Metallic nanoparticles are interesting condensed matter systems from the point of view of fundamental physics as well as for practical applications. In particular, these finite systems are particularly appealing, for instance as catalysts, as efficient scattering species in bioimaging, or as elementary magnetic recording units which could find applications in information storage and spintronics.
Enhancement of magnetism in nanoparticles of elements that are ferromagnetic as bulk solids such as cobalt or nickel has been discussed at length in many theoretical and experimental works and demonstrated through Stern-Gerlach experiments [1]. However, magnetism in nanoparticles of non-magnetic elements such as noble metals is more unexpected. Recently, magnetization measurements of small, mass-selected platinum clusters embedded in a zeolite matrix have been reported and the high spin polarization of these species has been confirmed [2]. Moreover, it has been shown theoretically that silver atoms become magnetic when grouped into a specific size of clusters [3]. In this work [3], first-principles calculations on small silver clusters have been carried out and it has been demonstrated that these species exhibit magnetic behavior that is not present in bulk silver. Finally, it has also been shown experimentally that gold nanoparticles covered with thiol are ferromagnetic at room temperature [4]. In organic chemistry, a thiol is a compound that contains the functional group composed of a sulfur atom and a hydrogen atom (−SH). The ferromagnetism exhibited by thiol-capped gold nanoparticles is not well understood, in particular from the theoretical side. One possible reason put forward to explain this unexpected magnetism is the fact that localized holes are induced in the 5d band by the thiol ligand (generated through Au-SH bonds) at the surface of the nanoparticle leading to the observed magnetism.
In another context, room temperature ferromagnetism in nanoparticles of non-magnetic metal oxides [5] and superconducting materials [6] has been observed. Let us also mention the recent discovery of magnetism in nanoislands of graphene [7] (which is also a non-magnetic material in the bulk) which seems to be due, according to the authors, to the degeneracies of the cluster's electronic spectrum coming from the edge atoms (i.e. from the finite size of the system). Finally, during the completion of this manuscript, we became aware of the theoretical work by Nielsen and Bhatt [8], who found that non-magnetic n-doped nanoscaled semiconductors at low dopant concentrations may exhibit important ferromagnetism if a small concentration of excess electrons is introduced.
Although the origin of ferromagnetism in these nano-objects is still under debate, it is clear that their surface must play a crucial role. Therefore, the ability to tailor the magnetic properties of non-magnetic materials by simply reducing the lengths of certain critical dimensions constitutes a promising field for investigation.
One of the most familiar aspects of magnetism is the ferromagnetic-paramagnetic phase transition which occurs at the Curie temperature. For non-metallic ferromagnetic materials, this phenomenon is well understood [9]. For metallic ferromagnets, the problem is far more complex 2 . This is due, in particular, to the presence of localized and itinerant electrons. 3 A natural question arises: what is the situation for a finite system such as a metallic nanoparticle? Do they undergo a critical phenomenon similar to the case of bulk ferromagnets? At which temperature? In this work, we investigate this fundamental issue.
For simple metal nanoparticles, the magnetism is caused exclusively by finite size effects and by the exchange-correlation interaction between valence electrons. This exchange interaction originates from the Pauli exclusion principle and favors spin polarization 3 . The Coulomb repulsion among the electrons is minimized when the spins in the electron open shell of the nanoparticle are aligned. This constitutes the well-known Hund's rule of atomic physics [9]. Small-and medium-sized metallic nanoparticles may be viewed as 'superatoms' having similar properties with respect to those of individual atoms [2]. The magnetic properties stem from degeneracies and partially filled orbitals. However, real silver nanoparticles do not present such degeneracies in general, except for a few specific sizes and atomic geometries [3].
In this work, we focus on the description of magnetic properties of metallic nanoparticles at relatively high electron temperature T (typically more than 500 K) whereas the ion temperature, T i , is supposed to be much smaller [12,13]. Such an out-of-equilibrium situation can easily be induced in metals [14] and metallic nanoparticles [12,13] using ultrashort optical pulses. The time evolution T (t) can be described using a two-temperature model in which the electrons are coupled to the lattice via an effective electron-phonon interaction [15]. It is worth mentioning that in this work we are investigating the magnetism of individual open-shell metal nanoparticles and not the one resulting from a collection of interacting magnetized clusters.
Spin density functional theory, in particular the local spin density approximation (LSDA), plays an important role in the theory of itinerant magnetism in metals. This theory was also applied successfully to atoms, molecules, clusters, surfaces and solids with nonzero spin and thus nonzero magnetization. An exhaustive description of the results obtained in these cases can be found in Lundqvist and March [16]. The generalization of this approach to finite temperature ensembles has been performed by Mermin [17]. The results presented in this work have been obtained within this theoretical framework. Non-collinear and orbital magnetism has not been considered.
The paper is organized as follows. In the next section, the theoretical models are described. In section 3, the results are presented and discussed. Finally, the concluding remarks are given. In the following, atomic units (h = |e| = m = 4π 0 = 1) are used unless otherwise specified and k B is Boltzmann's constant, a B the Bohr radius and µ B (= 1 2 in atomic units) the Bohr magneton.

LSDA at finite temperature
The metallic nanoparticles are described in the spherical background jellium model, which is known to be a very good approximation for simple metal nanoparticles. This model consists of replacing the real ionic core potential with the potential of a constant positive background corresponding to a uniformly distributed charge density ρ I [18]. For a nanoparticle of radius than the observed values) and the observed Curie-Weiss law above it. Improvements to this theory have been made that take into account the effect of spin fluctuations in a self-consistent renormalized (SCR) manner [11]. 3 The exchange interactions in metals among the conduction electrons themselves is often referred to as itinerant exchange. R = r s N 1/3 having A singly charged ionic cores and N valence electrons, the electrostatic potential energy is given by For Ag we have used r s = 3.02 a.u. In order to take into account the screening coming from the d electrons, the potential, V C ( r , r ), which appears in (1), corresponds to the Coulomb interaction between two charges situated at r and r in a medium constituted by a sphere of radius R and a static dielectric constant d = 3.62 [19]- [21]. The validity of the spherical jellium approximation (spherical jellium model (SJM)) for the modeling of open-shell metal clusters is discussed in section 3.
The electron system is assumed to be at a thermal equilibrium with a temperature T . In the Kohn-Sham formulation of the density functional theory at finite temperature within the grandcanonical ensemble [16,17], [22]- [25], the ground-state electron density n( r ) = n + ( r ) + n − ( r ) of an N = N + + N − -electron system is written, in terms of single-particle orbitals and energies ϕ σ i and ε σ i , as where f σ i = [1 + exp{(ε σ i − µ)/k B T }] −1 are the Fermi occupation numbers (0 f σ i 1), µ the chemical potential and σ = + (spin ↑) or − (spin ↓). These orbitals and energies obey the Schrödinger equation where V σ KS ( r ) is an effective single-particle energy potential given by where V H ( r ) is the Hartree potential energy and V σ xc ( r ) the exchange-correlation potential energy defined by where xc [n + , n − ] = n( r )ω xc (n + ( r ), n − ( r )) d r is the exchange-correlation thermodynamic potential [26,27] and ω xc the exchange-correlation thermodynamic potential per particle of the homogeneous electron gas calculated at the local density. By noting that (∂/∂n σ ) nω xc (n + , n − ) = (∂/∂n) {nω xc (n, m)} + σ (∂/∂m) {nω xc (n, m)} with m( r ) = n + ( r ) − n − ( r ), the expression (6) can be rewritten as [28] V σ where B xc ( r ) = −(1/µ B ) ∂ ∂m {nω xc (n, m)} n=n( r );m=m( r ) is the exchange-correlation magnetic field acting on spin. This is an internal magnetic field. If an external magnetic field, B ext ( r ) = B ext ( r )ẑ, directed along the z-axis is applied to the system one has instead of (5) The temperature appears in the self-consistent procedure only through the occupation numbers and the exchange-correlation thermodynamic potential. For low temperature (i.e. T T F [n( r )], where T F [n( r )] = (1/2k B )(3π 2 n( r )) 2/3 is the local Fermi temperature), xc [n + , n − ] may be safely replaced by its value at T = 0, that is by E xc [n + , n − ]. This approximation is particularly valid inside the particle where the local density is almost constant and close to the bulk one.
For a neutral cluster, the asymptotic behavior of V σ KS is given by the exchange contribution to V σ xc , which behaves at large distances as n( r ) 1/3 . As a consequence, the effective Kohn-Sham potential V σ KS decreases exponentially to zero, i.e. it does not reproduce the correct 1/r asymptotic behavior. This problem does not appear in the Hartree-Fock (HF) theory, because the HF exchange potential exactly compensates the self-interaction term contained in the Hartree potential. Following Perdew and Zunger [29], we have added a self-interaction correction that restores the correct asymptotic behavior of the effective potential. In the past, this procedure has been successfully applied to the study of both ground-and excited-state properties of atoms [30], molecules and small metal clusters [31,32]. The corrected Kohn-Sham potential V σ SIC,i is then given by where V σ xc,i ( r ) = ∂ ∂n σ nω xc (n + , n − ) It is easy to check that, as in the HF theory, the resulting potential has the correct 1/r behavior at large distances. Note also that this correction leads to KS wave functions which are not orthonormal. In order to be sure that this drawback does not disturb our calculations, a Gram-Schmidt orthonormalization procedure has been employed. The differences observed with and without using this procedure were found to be completely negligible. All the results presented in this paper have been obtained using this scheme.
Since the bulk Fermi temperature of noble metals is usually very high (T F ∼ 6.4 × 10 4 K for Ag), the above approximation is fully justified. However, the local density near the surface decreases by two to three orders of magnitude, leading to a much smaller local Fermi temperature. Therefore, the temperature variation of xc [n + , n − ] close to the surface is important and must be taken into account. In the present study, we have assumed that xc [n + , n − ] = E xc [n + , n − ] and ω xc [n + , n − ] = xc [n + , n − ]. Since the functional form of E xc is generally not known, several approximations have been proposed. In the LSDA, the electron density of the inhomogeneous system is considered to be locally uniform. Thus, the results of the homogeneous interacting electron gas may be used as an input data to investigate inhomogeneous electron systems as metallic nanoparticles. In this work, we have used the form obtained by Perdew and Zunger [29] and Ballone et al [33]. 6 The chemical potential is determined self-consistently by requiring the conservation of the total number of electrons from equation (3), i.e.
The electron specific heat at constant volume is defined as where E is the total electron energy of the system given by and S the grand-canonical entropy (for fermions) [26,27]: In the above expression, We note from (4) and (15) that we have leading to the definitions of the kinetic energy K and the potential energy E p : One also defines the magnetization M as and Finally, the Gibbs thermodynamic potential reads as [26] If B ext ( r ) = B ext , the above formula reduces to the well-known expression = E − ST − µN − M B ext . The quantities T , N and B ext are fixed.
Since we are working in the framework of the grand-canonical ensemble, it should be noted that N , which appears in (11), is the average particle number. The fluctuations of this quantity can be expressed in terms of the occupation numbers. In the case of fermions, we have [27] 7 In order to justify the use of the grand-canonical ensemble for the description of a finite electron system, the quantity N /N should be small. In the thermodynamic limit, when the number of particles in the system approaches infinity, the particle number fluctuations vanish and the grand-canonical ensemble is equivalent to the canonical ensemble. As discussed previously, an explicit temperature dependence of the exchange-correlation thermodynamic potential energy xc should in principle be included in the model. The finitetemperature field formalism seems to be well adapted to undertaking this task [23,34,35]. Compared with the exchange part, which can be solved exactly (see below), the treatment of the correlation part is much more difficult. This issue is out of the scope of the present study and will be addressed in a future work.
For the unpolarized electron gas, within the framework of the grand-canonical ensemble and the local-density approximation, let x [n; T ] be the exchange contribution of the thermodynamic potential of the system ( xc = x + c ). In the Hartree-Fock approximation [34], we have where is the Fermi-Dirac function which depends on p only through the electron energy p = p = p 2 /2m, 4 and µ F is the chemical potential of the Fermi gas. This latter quantity is a solution of the equation [36,37] with the Sommerfeld's function and where µ Bol is the Boltzmann chemical potential and λ T is the thermal de Broglie wavelength The evaluation of (21) can be carried out easily [23] and one finds We can then naturally extend the above result for the spin-polarized electron gas leading to where V σ x ( r ; 0) is given in [29]

Mean molecular field (MMF)
The first thermodynamic model to describe magnetic ordering transitions was developed in 1907, almost 100 years ago, by Weiss [38]. He proposed a phenomenological theory of ferromagnetism in which he assumed that the atomic magnetic moments are localized and interact with one another through a molecular field proportional to the average magnetization. Heisenberg attributed this field to the quantum mechanical exchange interaction between neighboring atoms (direct exchange interaction). Within this model, the Curie-Weiss law where χ is the magnetic susceptibility), which is verified by almost all ferromagnetic materials, is naturally explained. Defining the dimensionless quantitiesM = M/M 0 andT C = T /T C , the fundamental equation of the Weiss mean-field model is whereᾱ C = µ B B ext /k B T C and T C is the critical temperature (also called the Curie temperature for bulk materials) associated with the ferromagnetic-paramagnetic phase transition. Let us recall that the order parameter for this phase transition is the magnetization M. Within the framework of this model, we have the following two limiting cases: These two relations will be used for the interpretation of the LSDA results given in section 3.

Stoner model of the free electron gas
If ε is the energy of an electron in zero magnetic field, the energy of an electron with magnetic moment µ B σ in a field B ext is ε − σ µ B B ext . The magnetic moment per unit volume of an electron gas is therefore f 0 being the Fermi-Dirac function and g(ε) the density of states. For a parabolic band, one has In the above expressions, the chemical potential µ is determined by the condition where γ = µ B B ext , α B and the function S(x, y) are defined in (24) and (23), respectively. If γ = 0, one recovers equation (22) for the non-polarized electron gas.
In the Stoner theory of collective electron ferromagnetism [11,39], the assumption that the energy now includes a term proportional to the spontaneous magnetization is made, and is given by ±(U/2)δ with δ = M/nµ B , the sign depending on whether the electron spin is antiparallel or parallel to M (along the z-axis). The previous two equations are modified by replacing γ by γ + γ , where γ = (U/2)δ.
It can be easily shown that [11] at T = 0 and B ext = 0, a necessary condition for the existence of a spontaneous magnetization (ferromagnetism) is that I ≡g(E F )(U/n) > 1, where I is a dimensionless parameter andg(ε) the density of states per spin. This is usually called the Stoner condition. We see that magnetism is favored when the density of states is high at the Fermi level. For simple metals such as sodium, or noble metals such as silver, it is well known that this condition is not fulfilled [9], leading to nonmagnetic bulk materials. Some remarks have to be made: (a) If γ = 0 one recovers the Pauli paramagnetism for a free electron gas. At T = 0, since |γ | is only of the order of 10 −4 E F ∼ 10 −4 eV, by Taylor expanding (33), we obtain M = χ 0 B ext with χ 0 = µ 2 B g(E F ). This is the Pauli paramagnetic susceptibility. At finite temperature, dx (T and f (x) have been defined after equation (27)). At low temperature, we have the usual behavior, lim t→0 h(t) = 1 − π 2 /12 t 2 + O(t 4 ).
In the LSDA, the electrons within the metallic nanoparticle are treated as itinerant (or delocalized) electrons. Therefore, one may expect that the LSDA magnetic susceptibility as a function of temperature would follow the behavior predicted by the Stoner model. As we will see in the next section, this is not the case. Surprisingly, the behavior is the one obtained from the MMF which has been discussed in the previous subsection and which is based on localized magnetic moments. The reason underlying this effect is that, for a finite electron system, the situation is quite different with respect to a Fermi gas (which is employed in the Stoner model through a parabolic band) since the level density is a discrete energy function g σ (ε) = i δ(ε − ε σ i ). In the next section, many results of simulation are shown which illustrate the crucial role played by finite-size effects on the itinerant magnetism at finite temperature.

Results
Only clusters with partially filled shells (i.e. unpaired spins) can experience a net magnetic moment in the absence of an external field. This fact is a direct consequence of Hund's law in atomic physics. We learn from standard textbooks on magnetism that materials composed of atoms having partially filled d or f shells possess attractive magnetic properties. According to the concept of the 'superatom', the same must be true for partially filled d or f cluster shells. Generally speaking, metallic clusters having high symmetry favor electronic degeneracies and consequently the enhancement of magnetism. In particular, first-principles calculations performed on small silver clusters [3] have demonstrated that Ag 13 -which has a half-filled d-shell according to the spherical shell model-has an icosahedral symmetry and the highest magnetic moment in the cluster size range 2 n 22. In this work, by imposing a spherical geometry on to the ionic background (SJM), we somehow mimic the icosahedral symmetry of the real Ag 13 . Therefore, we expect that our methodology is suitable for the description of this category of clusters and for not very small sizes. This is why we have chosen to show essentially the results for the half-filled shell and medium-sized clusters that should exhibit the highest symmetry. Most of the calculations presented here have been performed for an openshell silver cluster containing 27 electrons. In agreement with Hund's law and the spherical-shell model we have checked that the LSDA ground state of Ag 27 corresponds to seven unpaired f electrons having the same spin projection (see figure 1). It leads to a magnetization M 0 = 7 µ B at T = 0 K or equivalently to 0.26 µ B per atom (for comparison, a platinum cluster having 13 atoms has a magnetic moment as high as 0.65 µ B per atom [2]). For the exchange-correlation term we have employed the functional of Perdew and Zunger [29]. Our numerical procedure has been checked by reproducing the results given in [33] for open-shell sodium clusters at T = 0.
It is well known that closed-shell metal clusters are spherically symmetric and therefore the use of the SJM is fully justified. Indeed, major spherical-shell edges have been experimentally identified for clusters containing 2, 8, 20, 40, 58, 92, etc. atoms and explained theoretically in the framework of the SJM [39,40]. Generally, the description of open-shell systems requires going beyond the SJM by incorporating non-spherical deformations of the ionic background. The first successful attempt is due to Clemenger [41], who used a shell model including axially symmetric distortions. In the latter, the ionic charge distribution is defined as where the quantity R(θ, ϕ) is determined by imposing that In the above expression, is the Heaviside step function. Of course, in the spherical case one has R(θ, ϕ) = R. For Ag 27 , we have performed an LSDA calculation using an ellipsoidal jellium geometry. Among all the possible geometric deformations, the absolute energy minimum is obtained for an ellipse axis ratio of 1.09. This value, being very close to 1, leads one to conclude that the geometry of the ionic background of Ag 27 is only slightly modified with respect to the spherical one. Furthermore, we have found that, in this deformed configuration, the ground-state electronic configuration is 1s 2 1p 6 1d 10 2s 2 1f 7 , in agreement with the ordering of levels obtained in the case of a particle enclosed in a spherical cavity of radius R.

Thermodynamic properties
In figure 2, the chemical potential µ, the grand potential , the entropy S and the specific heat at a constant volume c V of Ag 27 are depicted as a function of temperature. The relation ∂ E/∂ T = T (∂ S/∂ T ) has been checked numerically (see (12)), S being computed from (14) and E from (13).
The heat capacity exhibits a finite discontinuity around T C = 1350 K. This behavior is characteristic of a second-order phase transition in a finite-size system [43]. Let us stress that in  the thermodynamic limit (N → ∞), the heat capacity diverges. As expected for a second-order critical phenomenon, the other thermodynamic quantities µ, and S are continuous at T C but not their first derivative. It is also worth mentioning that the phenomenon of phase transition is only clearly defined at the thermodynamic limit.
Another well-known feature of a second-order phase transition in the bulk is the divergence of the magnetic susceptibility at the transition temperature. This latter quantity is defined as the second derivative of the free energy with respect to the external magnetic field. As it can be seen in figure 3, our finite-size system exhibits the same behavior.
In order to illustrate the critical character of this transition, we also show in figure 4 the number of iterations in the self-consistent Kohn-Sham procedure as a function of temperature. One clearly sees a dramatic increase of the iteration number at the critical temperature indicating the presence of large fluctuations.
As shown in figure 5, where magnetization as a function of temperature is depicted, the critical phenomenon is associated with the full demagnetization of the nanoparticle: above T C the system becomes paramagnetic. Thus, we have a ferromagnetic-paramagnetic 'phase transition' of the cluster spin system and T C is the associated critical temperature. In the inset of figure 5, the LSDA results are presented in normalized units (solid squares) along with the predictions of the molecular mean field (full line) given by equation (30) with B ext = 0. Surprisingly, the two models (LSDA and MMF) give almost exactly the same results. Naturally, the MMF model does not provide either the value of T C or the value of M 0 . It must be emphasized that even if LSDA includes electron correlations, the temperature dependence of the magnetization can be obtained from a pure mean-field model. As already mentioned, the LSDA results and those from the mean-field model are almost exactly the same. This is surprising because the mean-field model is a bulk model and therefore does not incorporate finite-size effects. In order to see some differences between the two approaches, we have to investigate the region located close to the phase transition. Indeed, by definition, in a bulk phase transition the correlation length becomes infinite at the critical point, which is incompatible with the finite size of the nanoparticle. This is one of the main reasons why the study of phase transitions in finite systems is so complex.
In order to better characterize the critical phenomenon, it is necessary to study the behavior of the order parameter M for some limiting cases. In figure 6, the reduced magnetization is shown close to the critical temperature. In the mean-field theory, the so-called magnetization exponent β is defined asM ∼ (T C − T ) β for T < T C and T → T C [43]. From equation (30), we see immediately that in the case of zero external field,M 2 ∼ = 3T C ((T C − T )/T C ) leading to β = 1 2 . In figure 6, we see that the LSDA magnetization starts to deviate from the mean-field result around 0.1 K (log 10 (T C − T ) ≈ −1) below the critical temperature. At this temperature, the mean-field magnetization has already reached its asymptotic behavior (T C − T ) 1/2 , which is  represented by a dotted line in the figure. From the LSDA results, one may be tempted to extract a critical exponent (similar to β for the MMF). Unfortunately, in order to do so and due to the finite size of the studied systems, one should normally use scaling techniques which are well beyond the scope of this work [44].
Near zero temperature, the mean-field model predicts that the magnetization deviates from its saturation value by a term 2 exp(−2/T C ). In figure 7, the behavior of the magnetization near zero temperature is depicted. We see that both models, MMF and LSDA, follow the same law but with different parameters.   (29)).
It is well known that the calculation of magnetic properties requires high numerical precision, in particular for the simulations performed close to the critical temperature T = T C . In this work, we have checked that the convergence of the KS procedure can be achieved if one uses a precision criterion as small as 10 −7 for the KS energies. This point is clearly illustrated in figure 8, where the magnetization, computed with different numerical precision, is shown close to the critical temperature.
Another point of interest concerns the role played by electron correlations on the magnetic properties. In figure 9 we show the magnetization of Ag 27 calculated including or not including (i.e. V σ x instead of V σ xc in equation (5)) the correlation term in the LSDA model. It is clearly seen that the critical temperature is strongly enhanced (more than two times higher) when electron correlations are absent. The reason correlation effects favor the unmagnetized state is simple: if only the exchange term is included (symbol X in figure 9), the exclusion principle keeps electrons of parallel spin separated spatially, while electrons of antiparallel spin are spatially uncorrelated. In passing from the non-magnetic to the ferromagnetic state, the large potential energy associated with these uncorrelated pairs is lost. In the correct (correlated, symbol XC in figure 9) wave function, however, even electrons of antiparallel spins will avoid one another somewhat, and so there will be so much potential energy to be lost. Therefore, correlation should favor unmagnetized states, i.e. should decrease the Curie temperature of the system. This is the trend that is observed in figure 9.
We have also checked the modification of the results due to an explicit temperature dependence of the exchange thermodynamic potential energy given in equation (29). As seen in figure 9, the magnetization is only slightly modified. Owing to its smallness, this correction will be neglected in the following.
In figure 10, the spin densities of states (SDOS) as functions of temperature are depicted. For the sake of clarity the bar energy levels (i.e. KS energies) have been convoluted with a Lorentzian profile of 0.086 eV FWHM. One clearly sees that the thermal demagnetization is due to the variation with increasing temperature of the chemical potential (the solid line moves from top to bottom) of the KS energies (the spin down (up) energies decrease (increase)) and the Fermi occupation numbers. All the occupied orbitals of the ground state at zero temperature (five for this nanoparticle see figure 1) and the lowest unoccupied molecular orbital (LUMO) (the horizontal dashed line in figure 10 shows the energy position of this orbital which is empty at T = 0 K) play a major role in the thermal evolution. More excited states above the LUMO are involved but their contribution is minor compared with the one of the LUMO. It is worth noticing the position of the chemical potential (solid line in the figure) which is situated, according to the laws of quantum statistical mechanics, at an equal distance in energy between the LUMO and the HOMO.

Charged clusters
Owing to the problem of mass selection, most of the experiments on free mass-selected nanoparticles or clusters are carried out with charged species. We have computed the magnetization for two isoelectronic free open-shell silver clusters having 27 valence electrons. Furthermore, it is also of interest to investigate the influence of the cluster charge on the magnetic properties. The results are presented in figure 11. With respect to the neutral case, the magnetization changes are weak, the critical temperature being slightly higher (lower) for the negatively (positively) charged cluster. As seen in the next subsection, this result is consistent with the evolution of the temperature with the mass of the system.

Size evolution
Here we investigate the magnetism of partially filled f-shell nanoparticles. The reduced magnetization of Ag n with n = 21-33 as a function of T C − T is depicted in figure 12. These open-electron systems correspond with the gradual filling of the f-shell. Of course, the two closed-shell clusters Ag 20 and Ag 34 bracketing this size interval are non-magnetic. With the exception of four species (Ag n with n = 21, 27,32,33), all the curves are almost bunched together. This is true despite the fact that clusters having the same magnetization at zero temperature have critical temperatures which differ notably (see the inset in figure 12 where the magnetization of Ag 23 and Ag 31 is shown). Surprisingly, if one represents the critical temperature as a function of the number of atoms (see figure 13), the highest temperature does not correspond to Ag 27 (which has a half-filled f-shell and therefore the highest magnetization at T = 0 K) but is associated with Ag 26 . This is a clear size effect. Indeed, if one uses only the argument based on the number of aligned spins (which leads to associate the highest T C to Ag 27 ), one forgets the mass dependency of the critical temperature (see below) which goes in the opposite direction. Let T ref be the critical temperature for Ag 27 . This system constitutes our reference for the following discussion. In order to be more quantitative, we have performed two complementary numerical simulations: (i) for a nanoparticle of 26 electrons (six unpaired electrons) having the same volume as Ag 27 , we have found a critical temperature smaller than  T ref . This result is logical since the density, with respect to Ag 27 , is smaller; (ii) for a nanoparticle of 27 electrons (seven unpaired electrons) having the same volume as Ag 26 , we have found a larger critical temperature. Once again, this result is logical since the density is now larger compared to that of Ag 27 . Therefore, the fact that Ag 26 has a slightly larger critical temperature compared to Ag 27 is due to a combination of these two effects: size and the number of unpaired electrons. Finally, let us stress that the subtle interplay between them does not allow us to make any predictions in advance.
In the field of cluster physics, a common size-dependent effect is related to the fraction of atoms (electrons) situated at the surface [45]. The fraction of electrons at the surface n s scales with the surface area divided by the volume, which is proportional to n −1/3 . Since silver is All the nanoparticles (except Ag 27 ) have a half-filled d shell and M 0 = 5 µ B . As expected, the temperature range in which the system is magnetic diminishes strongly with increasing nanoparticle size and, in the bulk limit, the magnetism disappears. Since the magnetism studied here is related to the presence of electrons at the surface of the nanoparticle, one may expect that the critical temperature scales with n s and thus with n −1/3 . This behavior is clearly demonstrated in figure 15 where T C is plotted against n −1/3 for the same species as in figure 14.

Density
In figure 16, the magnetization of a nanoparticle containing 27 electrons as a function of temperature is shown for different bulk electron densities. We observe that the critical temperature decreases with the electron density. This trend is similar to the one observed for bulk ferromagnetic materials.
The physical explanation underlying this behavior is the following: exchange-correlation forces, which are responsible for the ferromagnetism, get more and more important when the density increases. Therefore, the thermal energy (with the associated critical temperature) needed to destroy the spin alignment increases proportionally.
It is also worth noticing that we observe the same effect as seen in figure 13: a higher spin density (small value of r s ) leads to a stronger ferromagnetism (higher value of T C ). In figure 13, the volume of the nanoparticle is kept almost constant while the number of unpaired spins is varied, leading to a change in the spin density.

Magnetic field
We end this section by studying the influence of an external magnetic field applied along the z-axis (B ext = 0) to the cluster Ag 27 . Above the critical temperature the system becomes  paramagnetic. In figure 17, the magnetization of Ag 27 as a function of B ext for different temperatures is shown. It is worth noticing that the magnetization of Ag 27 is far from the saturation magnetization at B ext = 15 tesla. The magnetic susceptibility is defined as M = χ B ext . This quantity is extracted from the linear behavior of the LSDA magnetization with respect to the external magnetic field. As we have seen previously, the LSDA results for the magnetization are very close to the predictions of the molecular mean-field model of Weiss (see figure 5). The latter theory predicts that, above the critical temperature, the magnetic susceptibility χ is proportional to 1/(T − T C ) which constitutes the well-known Curie-Weiss law. Therefore, the LSDA model is also expected to verify this law. In figure 18, the inverse magnetic susceptibility 1/χ is plotted against T − T C .  The linear behavior is clearly seen as confirming the expected tendency. Let us stress that this result is in agreement with recent experimental measurements performed on mass-selected platinum clusters [2].
It is also worth noticing that according to our previous discussion in (2.3), the Stoner model using a parabolic band is not able to reproduce the Curie-Weiss law. Here, it is the discrete nature of the electron spectrum and therefore the finite size of the system that plays a crucial role in explaining this behavior. Assuming an energy gap between the HOMO and the LUMO (typically = 0.4 eV and k B T C = 0.12 eV for Ag 27 ), it can be proven analytically using equation (32) that for temperatures T above T C , k B (T − T C ) , the magnetic susceptibility is proportional to 1/(T − T C ). The existence of the HOMO-LUMO gap which is related to the finite size of the system is therefore essential in explaining the Curie-Weiss law. When the size of the nanoparticle increases, tends to zero and we recover the 1/(T − T C ) 2 behavior of the magnetic susceptibility predicted by the bulk Stoner model.

Conclusions
Thermal magnetic properties of open-shell metallic nanoparticles of elements that are not magnetic in the bulk have been studied within the framework of the LSDA at finite temperature and the jellium model. The magnetism of these nanosized objects, which may be viewed as 'superatoms', is only due to finite size effects and the electron spins are oriented in accordance with Hund's law.
A magnetic critical phenomenon that has the characteristics of a second-order ferromagnetic-paramagnetic 'phase transition' has been identified and studied. It constitutes the main result of this work.
The critical temperature vanishes in the bulk limit, scales with N −1/3 and strongly depends on the electron correlations. Surprisingly, the LSDA-reduced magnetization expressed in terms of T /T C is almost identical to the one obtained from the molecular mean-field approach of Weiss. Accordingly, it is found that the magnetic susceptibility above the critical temperature obeys the Curie-Weiss law which is verified by almost all bulk ferromagnetic materials. This result is in agreement with recent experimental measurements performed on small mass-selected platinum clusters [2]. This behavior seems to be due to the discrete character of the electron spectrum resulting from the finite size of the system.
In order to observe this magnetic critical phenomenon, it is important to be in a highly nonequilibrium regime where T T i . As described in [12,13], such a high electron temperature can be easily induced in noble metal nanoparticles using ultrashort optical pulses.