Statistical Thermodynamics of Chiral Skyrmions in a Ferromagnetic Material

Solitons are a challenging topic in condensed matter physics and materials science because of the interplay between their topological and physical properties and for the crucial role they play in topological phase transitions. Among them, chiral skyrmions hosted in ferromagnetic systems are axisymmetric solitonic states attracting a lot of attention for their dazzling physical properties and technological applications. In this paper, the equilibrium statistical thermodynamics of chiral magnetic skyrmions developing in a ferromagnetic material having the shape of an ultrathin cylindrical dot is investigated. This is accomplished by determining via analytical calculations for both Néel and Bloch skyrmions: (1) the internal energy of a single chiral skyrmion; (2) the partition function; (3) the free energy; (4) the pressure; and (5) the equation of state of a skyrmion diameters population. To calculate the thermodynamic functions for points (2)–(5), the derivation of the average internal energy and of the configurational entropy is crucial. Numerical calculations of the thermodynamic functions for points (1)–(5) are applied to Néel skyrmions. These results could advance the field of materials science with special regard to low-dimensional magnetic systems.


Introduction
The thermodynamic description of topological defects and topological phase transitions has been one of most important challenges of the modern condensed matter physics. The leading works on the thermodynamics of topological defects and the relevant underlying physics are the ones dating back to the 70s of Berezinskii and of Kosterlitz and Thouless on defect-mediated phase transitions in two-dimensional (2D) XY superfluid model. In particular, Berezinskii studied the low-temperature state of one-dimensional (1D) and 2D classical and quantum systems such as crystals, isotropic magnetic substances, and superfluids and superconductors having continuous symmetry, showing that in both kinds of systems the long-range order is destroyed due to the increasing fluctuations of the ordering parameter with increasing size [1,2]. Kosterlitz and Thouless deepened this type of investigation, arguing the existence of topological defects having the form of vortices in physical systems described by the XY model, such as superfluids [3][4][5]. They studied the thermodynamic behavior of these systems by means of the calculation of the Helmholtz free energy F showing that, for T→0 K and with increasing size, F can be minimized if no vortices appear, while above a critical temperature F is minimized if there is the formation of couples of unpaired vortices and anti-vortices. On the basis of this analysis, they introduced the concept of defect-mediated topological phase transition in condensed matter physics systems, the so-called infinite-order Berezinskii-Kosterlitz-Thouless phase transition.
After these seminal works, a great deal of attention has been paid to study the interplay between topology and physics characterizing solitons in condensed matter systems and its relevant implications. diameters population was calculated as an expectation value in a way similar to Shannon's information entropy. Because of this definition, the skyrmion diameter dependent Helmholtz free energy, F(D sky ) = E − T S (with E = E(D sky ) the skyrmion energy and D sky the skyrmion diameter), also exhibits a square diameter dependence for each T and the mean square fluctuation of energy of the skyrmion diameters around its average value population is equivalent to the mean square fluctuation of free energy F(D sky ) around its average value.
Recently, in some excellent papers the skyrmion energy has been analytically calculated using either a 1D distribution [40] that was used to describe properties of thick-walled magnetic domains in uniaxial platelets [41] or a two-dimensional domain distributions [25,42], or numerically calculated using an accurate 2D magnetization distribution recently proposed to describe the magnetization texture, namely the orientation of the magnetization vector in the xy plane and with respect to the z direction, of the Néel (hedgehog) chiral skyrmion [35,43]. Additionally, in [42] the magnetostatic contribution of the volume charges has also been included in the skyrmion energy computation.
However, a systematic theoretical study describing the classical statistical thermodynamic behavior of chiral magnetic skyrmions at equilibrium is still lacking in the literature. In this work, a statistical thermodynamic description of chiral magnetic skyrmions that are hosted in ultrathin cylindrical dots based on a classical approach is given, strengthening the analogy with the ideal gas. This is accomplished via the following key results obtained by means of the analytical and numerical calculation of the following physical quantities at a given T: (1) the calculation of the skyrmion energy for both Néel and Bloch skyrmions that is regarded as the internal energy; (2) the determination of the partition function and of the free energy of a skyrmion diameters population within a microcanonical ensemble; (3) the derivation of the skyrmion pressure from the free energy of a skyrmion diameters population and of an equation of state linking the thermodynamic variables P, V, and T.
In addition, it is shown that the configurational entropy calculated for a Nèel skyrmion diameters population takes the same form for a Bloch skyrmions population. This feature also characterizes the partition function and the pressure and results from the fact that the skyrmion internal energy has a quadratic dependence on the skyrmion diameter (or radius) independently of the topological texture and skyrmion number. Therefore, it can be concluded that the statistical thermodynamics qualitative description of axisymmetric solitonic states does not depend on the skyrmion magnetization texture of the chiral skyrmion under study.
In this study the skyrmion energy is analytically computed starting from the 2D magnetization distribution proposed to describe the magnetization texture of the Néel (hedgehog) chiral skyrmion [35,43]. This accurate 2D magnetization distribution recovers, in the isotropic case and for dominating exchange interaction, the Belavin-Polyakov soliton [44] solution and can be easily extended to describe the magnetization texture of the Bloch (vortex-like) chiral skyrmion. In this respect, here it is shown that, for both magnetization textures, the skyrmion energy can be expressed as a combination of elementary transcendental functions and is actually regarded as an internal energy. Owing to some reasonable approximations, from the skyrmion energy a simple analytical form of the equilibrium skyrmion radius in the region of metastability depending on the scaled magnetic parameters is obtained.
The partition function of a skyrmion diameters population is a key physical quantity to describe their statistical thermodynamics. Indeed, it allows to understand the connection between the occupation of microscopic states by the skyrmions population and the macroscopic thermodynamic variables of state, such as the skyrmion free energy and the entropy.
The definition of a skyrmion pressure for a skyrmion diameters population of average volume <V> at a given temperature T is accomplished exploiting the analogy with the pressure exerted by the molecules of a gas on the walls of the container. The derivation of the skyrmion pressure p also allows writing an equation of state linking the thermodynamic variables pressure p, volume V, and T equivalent to the one for an ideal gas. This general relation is valid for any single skyrmion magnetization texture. It is shown that, for D sky = <D sky >, the equation of state reduces to the well-known one for an ideal gas pV = k B T for the number of particles N = 1. However, for diameters around the average diameter, unlike the ideal gas characterized by the universal constant R, it is not possible to define for a skyrmions population a universal constant.
Calculations (1)- (3) are not only important for understanding the thermodynamic properties of chiral skyrmions that form in ultrathin ferromagnetic dots but could enable a further comprehension of the magnetic properties of ferromagnetic materials hosting these topological objects, suggesting new thermodynamic measurements able to confirm the predictions. Especially, calculations (2) and (3) differentiate the analysis and the aim of this work with respect to that of recent studies of the literature that mainly focus on the dynamical properties of magnetic skyrmions treating the equilibrium statistical properties of magnetic skyrmions in the region of metastability within a microcanonical ensemble.
This paper is organized as follows: in Section 2, Methods, the model is presented. First, the skyrmion energy is analytically calculated for both Néel and Bloch magnetic skyrmions according to a simple variable transformation and to some reasonable approximations. Then, the statistical thermodynamics of a skyrmion diameters population is studied within a microcanonical ensemble via the calculation of the partition function, the Helmholtz free energy, and the pressure and by means of the derivation of an equation of state that is similar to the one for an ideal gas. In Section 3, Results and Discussion, the model outlined in Section 2 is applied to a diameters population of Néel skyrmions and the behavior of the corresponding thermodynamic quantities is discussed. The model is also benchmarked via the comparison of the equilibrium diameters as a function of the external magnetic field with available experimental data taken from the literature. In Section 4 Conclusions are drawn.

Methods
In the following calculations the polar coordinates in the dot plane ρ = (ρ,ϕ) are introduced, where ρ is the radial coordinate and ϕ is the azimuthal coordinate defining spherical angles (θ,Φ) of the magnetization vector as functions of ρ. The magnetization unit vector is m = M/M s .
Replacing the above expressions, one gets ε exch = A dθ dρ 2 + sin 2 θ ρ 2 for both Néel and Bloch skyrmions and ε exch does not depend on chirality.
The general form of the bulk DMI energy density is in front of D refers to χ = +1 or counter-clockwise (χ = −1 or clockwise) chirality. Therefore, in the ultrathin film limit, the bulk DMI energy density of a Bloch skyrmion assumes the same expression of the IDMI energy density of a Néel skyrmion [41]. D > 0 (D < 0) corresponds to the case of the heavy metal under (over) the ferromagnetic material. In general, for the proper skyrmion chirality either for a Néel or a Bloch skyrmions, the DMI energy lowers the skyrmion energy.

Anisotropy Energy Density
The anisotropy energy density gets contributions from the perpendicular uniaxial anisotropy and the demagnetization (magnetostatic) energy densities, respectively. The former term takes the compact form ε ani = −K u m 2 z , K u being the perpendicular uniaxial anisotropy constant. In explicit form expressed as a function of θ it reads ε ani = −K u cos 2 θ. Instead, the latter term takes the form ε dem = − 1 2 µ 0 M s m· H dem , with H dem = (0,0,H dem ) being aligned along +z and µ 0 = 4π × 10 −7 H/m being the vacuum permeability. Since skyrmions that are hosted in ultrathin ferromagnetic dots are studied, the ultrathin approximation is applied, according to which the magnetostatic source is given only by the face surface charges of the dot, H dem = −M s m z . Indeed, according to this approximation, in the ultrathin limit the contributions resulting from the side surface charges of the dot and from the volume charges can be safely neglected. This leads to ε dem = 1 2 µ 0 M 2 s cos 2 θ that has the form of a local term. The total anisotropy energy density can be written as an effective anisotropy by defining an effective anisotropy constant K eff = K u − 1 2 µ 0 M 2 s so that ε eff ani = −K eff sin 2 θ. However, for a more realistic description it is convenient to refer the anisotropy energy density to the uniform state (m along +z), ε ani (θ =0) = −K u obtained for θ = 0 expressing it as ε ani = K u sin 2 θ, namely as the difference ε ani (θ) − ε ani (θ =0) recovering the definition of ferromagnetic thin films. Hence, ε ani eff = K u − K eff cos 2 θ and the term proportional to K u leads to a constant upshift of the skyrmion energy from negative to positive values due to the high K u magnitude for typical ferromagnetic materials without affecting its trend vs. the skyrmion radius.

External Field Energy Density
The Zeeman energy density due to the interaction of the static magnetization with the external bias field H ext is ε extfield = −µ 0 M s m·H ext . In terms of the polar angle θ between the magnetization m and the z axis and the angle θ Hext the external bias forms with the z axis it is ε extfield = −µ 0 M s H ext cos(θ-θ Hext ). If H ext is aligned along +z (-z), the Zeeman energy density can be written as ε extfield = −µ 0 M s H ext cos θ, with H ext > 0 (H ext < 0). In both cases, the sign of the Zeeman energy density depends on the sign of cosθ, which is in turn related to the m orientation with respect the z axis.
Instead, if this contribution is expressed as the energy gain with respect to the Zeeman energy density of the perpendicular uniform state (m aligned along + z, θ = 0, with H ext aligned along +z The total energy density is written down using the dimensionless radial coordinate r = ρ/l exch with l exch = √ 2A/(µ 0 M s 2 ) being the exchange length. The total energy density for both skyrmion magnetization textures as explicit function of the polar angle θ reads: with A = A/l 2 exch and D = D/l exch . The total energy density in dimensionless units is calculated at equilibrium using the static magnetization distribution θ(r) = Θ 0 (r) at equilibrium of a chiral Néel skyrmion derived in [35]. This distribution is a trial solution of the nonlinear and transcendental differential equation resulting from the minimization of the skyrmion energy functional complemented by the boundary conditions on the distribution itself and its radial derivative. In the limit of dominating exchange isotropic interaction (quality factor Q = 1 with Q = 2K u /(µ 0 M s 2 ) and DMI parameter D = 0), this solution recovers the well-known Belavin-Polyakov soliton solution [42,44]. Θ 0 (r) describes the radial dependence of the static magnetization of a chiral Néel skyrmion having static magnetization along -z (corresponding to skyrmion number S = −1) hosted in a circular dot. It takes the form: where B = r sky e ξ r sky , with r sky =R sky /l exch the dimensionless skyrmion radius and R sky the skyrmion radius, and ξ = √ Q − 1. According to Equation (3a) it is Θ 0 (r = 0) = π(m z = −1, static magnetization along −z in the core center) and Θ 0 (r → ∞) = 0 (m z = +1, static magnetization along +z at the skyrmion border).
Analogously, the corresponding magnetization distribution at equilibrium for static magnetization along +z with skyrmion number S = +1 reads: According to Equation (3b) it is Θ 0 (r = 0) = 0 (m z = +1, static magnetization along +z in the core center) and Θ 0 (r → ∞) = π (m z = −1, static magnetization along -z at the skyrmion border). Note that the trial magnetization distribution expressed by Equation (3) is the solution of the variational nonlinear differential equation given in [35]. This trial solution is complemented by the boundary conditions Θ 0 (r = 0) =0,π, Θ 0 (r = r sky ) =π/2 and ∂Θ 0 ∂r (r = r d ) = |D|l exch 2A with r d = R d /l exch the dimensionless dot radius (R d is the dot radius) and is valid also for a Bloch skyrmion because the variational differential equation takes the same form (the bulk DMI energy density has indeed the same expression as the IDMI energy density) and the boundary conditions are the same. In principle, these boundary conditions do not take into account the canting of the magnetization at the borders that would become important when the skyrmion radius is comparable to the dot radius. However, in the analytical and numerical calculations presented and discussed in Section 3, this issue does not arise. Indeed, the equilibrium skyrmion radius calculated in the metastability region for vanishing external bias field and in the presence of H ext is at most less than 1 4 of the dot radius for all the cases examined including the comparison of the numerical calculations with the experimental data.
From Equations (3a) and (3b) one gets: where the minus (plus) sign is referred to the skyrmion magnetization texture with the static magnetization along -z (+z) in the core center corresponding to skyrmion number S = −1 (S = +1). Substituting Equation (3) in Equation (2) yields: Here, the upper (lower) signs in the first term of the exchange energy density and in the two terms of the DMI energy density refer to the magnetization texture with the static magnetization along −z (+z) in the core center (r = 0) corresponding to S = −1 (S = +1). Instead, the upper (lower) sign in the Zeeman energy density refers to the static magnetization m z along −z (+z) in the core center (r = 0) and along +z (−z) at the skyrmion border being cos π − 2arctan B e −ξr r = − cos 2arctan B e −ξr r and the external magnetic field H ext applied either along +z (H ext > 0) or along -z (H ext < 0). In other words, the "−" ("+") sign in front of the Zeeman energy density contribution refers to the cases where the static magnetization m z (r = 0) is anti-parallel (parallel) to H ext . By using some trigonometric expansions (see the Appendix A, Equations (A1)-(A4)), one gets: Equation (5) is Equation (A5) of the Appendix A. The expression of the total energy density given by Equation (5) is valid for any Néel and Bloch magnetic skyrmion either for chirality χ = +1 (+ sign in front of the DMI term) or χ = −1 (− sign in front of the DMI term). Here, it is H ext > 0 (H ext < 0) if the external magnetic field is along +z (−z). For D > 0, it is always m z along −z in the core center, r = 0 (S = −1), and χ = +1 (radially outward Néel skyrmion and counter-clockwise Bloch skyrmion).

Skyrmion Energy
In this subsection the calculation of the skyrmion total energy E is outlined. The energy density is integrated over the dot volume E r sky = ε tot r, r sky dV that, in explicit form, is E r sky = ε tot r, r sky dV = (l exch ) 2 ε tot r, r sky r dr so that: where t is the dot thickness. The total energy expressed in Equation (6) is computed starting from the energy density of Equation (5). The total energy consists of the following energy contributions: (1) ferromagnetic exchange; (2) DMI (either IDMI or bulk DMI); (3) effective anisotropy; and (4) contribution dependent on the external bias field. The following dimensionless integrals are evaluated analytically: (1) three I i-exch for the ferromagnetic exchange; (2) four I i-DMI for the IDMI; (3) I ani for the effective perpendicular anisotropy; and (4) I extfield for the contribution dependent on the external bias field. All integrals have a dependence on the skyrmion radius r sky . In compact form: with E = E(r sky ). The explicit expressions of the integrals appearing in Equation (7) are in the Appendix A (Equations (A6)-(A14)). For every contribution of the total skyrmion energy, the radial integration is carried out based on two substitutions. The first substitution is a change of variable: (1) ξ r = s. This leads to e ξ r = e s and to r = 1/ξ s. The second substitution is a real substitution based on the change of variable given in 1): (2) s e s = x so that x = ξ r e ξr . The transcendental equation s e s = x admits the solution s = W(x), allowing to write W(x) e W(x) = x, where W(x) is the Lambert function defined in the real domain as a function of the real variable x and ds = W(x)/(x (1 + W(x))) dx, with W(0) = 0. This allows writing every integral in terms of the variable x, with 0 ≤ x ≤ C and C = ξr d e ξr d a quantity depending on the dot radius. Hence, in these calculations the domain of the Lambert function is Hence, all the integrand functions can be written as the product of a rational fractional function f = f (B,x) and a fractional function g(W(x)) dependent on W(x). While the former function has an explicit dependence on the skyrmion radius r sky via B and on the radial coordinate r, the latter function has an explicit dependence only on the radial coordinate r being x = ξ r e ξ r . The generic integral I for every skyrmion energy contribution is expressed in the form: where K = K(B,ξ) is a coefficient that is dependent on B and ξ and assuming different forms depending on the considered integral and g(W(x)) is either g(W(x)) = W(x)/(1 + W(x)) or g(W(x)) = (W(x)) 2 /(1 + W(x)) depending on the considered integral. The function g(W(x)) appearing in the integral of Equation (8a) is plotted in Figure 1 for 0 ≤ x ≤ C for the two different functions appearing in the skyrmion energy integrals (see the next subsections for the details). In Figure 1a the function W(x)/(1 + W(x)) is shown, while in Figure 1b it is the function (W(x)) 2 /(1 + W(x)). In the skyrmion energy integrals (see Section 2.2.1 for the details) it appears either W(x)/(1 + W(x)) or (W(x)) 2 /(1 + W(x)), depending on the integral. For the special case shown corresponding to the parameters at T = 0 K (see Section 3 for their numerical values) C = 2 × 10 13 . Both functions vanish for x = 0, are monotonically increasing functions but, while W(x)/(1 + W(x)) tends asymptotically to 1 for increasing x, (W(x)) 2 /(1 + W(x)) diverges for x → ∞, exhibiting a trend similar to that of W(x). The smooth behavior of W(x)/(1 + W(x)) for small x that is masked by the large x interval is shown in more detail in the inset to Figure 1a for a reduced interval of x (0 ≤ x ≤ 100).
Because of the rescaled parameters used at higher temperatures, a reduction of the value of C occurs without altering the qualitative behavior of the two functions.
Materials 2019, 12, x FOR PEER REVIEW 9 of 30 smooth behavior of W(x)/(1 + W(x)) for small x that is masked by the large x interval is shown in more detail in the inset to Figure 1a for a reduced interval of x (0 ≤ x ≤ 100). Because of the rescaled parameters used at higher temperatures, a reduction of the value of C occurs without altering the qualitative behavior of the two functions.
The generic integrand function h(x) of Equation (8a) does not have a primitive function for every considered energy contribution, therefore the integral I cannot be solved analytically. However, from a numerical check studying the behavior of the integrand functions, it has been found that, for every energy contribution, the integral I is well approximated if calculated in the form: namely, moving out of the integral the function g(W(x)) depending on the radial coordinate only, and attributing to it a trend depending on the skyrmion radius via the dependence on B (with B = B (rsky)) and writing g(W(Bξ)) in place of g(W(x)). In this way, the radial dependence contained in g(W(x)) has been interchanged with the skyrmion radius dependence. This has been accomplished, first by numerically calculating the integral of Equation (8b) whose integrand f(B,x) has a primitive function for every energy contribution of Equation (7) (see the next subsections), and then by numerically computing the ratio between the numerical integral of Equation (8a) and the integral of Equation (8b), both depending on rsky for the values of the skyrmion radius in the interval 0 ≤ rsky ≤ rd. The calculation was done for every energy contribution of Equation (7). The above-mentioned ratio turned out to be approximately equal to the function g(W(Bξ)) for all skyrmion radii ranging in the interval 0 ≤ rsky ≤ rd. In particular, , with 0 ≤ B ≤ 1/ξ C depending on the integral studied (see the following subsection for the details).

Calculation of Ferromagnetic Exchange Energy
In this subsection, the ferromagnetic exchange energy is computed. Taking into account substitution (1), multiplying numerator and denominator by e 4s , and substitution (2) 1-exch I is rewritten (Equation (A6)) as: Taking into account substitution (1), multiplying numerator and denominator by e 4s , and substitution (2), 2 -exch I (Equation (A7)) is rewritten as: The generic integrand function h(x) of Equation (8a) does not have a primitive function for every considered energy contribution, therefore the integral I cannot be solved analytically. However, from a numerical check studying the behavior of the integrand functions, it has been found that, for every energy contribution, the integral I is well approximated if calculated in the form: namely, moving out of the integral the function g(W(x)) depending on the radial coordinate only, and attributing to it a trend depending on the skyrmion radius via the dependence on B (with B = B (r sky )) and writing g(W(Bξ)) in place of g(W(x)). In this way, the radial dependence contained in g(W(x)) has been interchanged with the skyrmion radius dependence. This has been accomplished, first by numerically calculating the integral of Equation (8b) whose integrand f (B,x) has a primitive function for every energy contribution of Equation (7) (see the next subsections), and then by numerically computing the ratio between the numerical integral of Equation (8a) and the integral of Equation (8b), both depending on r sky for the values of the skyrmion radius in the interval 0 ≤ r sky ≤ r d . The calculation was done for every energy contribution of Equation (7). The above-mentioned ratio turned out to be approximately equal to the function g(W(Bξ)) for all skyrmion radii ranging in the interval 0 ≤ r sky ≤ r d . In particular, g(W(Bξ)) = W(Bξ)/(1 + W(Bξ)), g(W(Bξ)) = [W(Bξ)] 2 /(1 + W(Bξ)), or g(W(Bξ)) =1/(2ξ 3 ) [W(Bξ)] 2 /(1 + W(Bξ)), with 0 ≤ B ≤ 1/ξ C depending on the integral studied (see the following subsection for the details).

Calculation of Ferromagnetic Exchange Energy
In this subsection, the ferromagnetic exchange energy is computed. Taking into account substitution (1), multiplying numerator and denominator by e 4s , and substitution (2) I 1-exch is rewritten (Equation (A6)) as: Taking into account substitution (1), multiplying numerator and denominator by e 4s , and substitution (2), I 2-exch (Equation (A7)) is rewritten as: Now, by considering the sum I 1-2exch of the two integrals of Equations (9a) and (9b) that does not depend on W(x): Performing the integral of Equation (9c) one gets: with P = B/C and, in explicit form, P = r sky ξ r d e ξ(r sky −r d ) . Hence, P does depend on r sky that is a variable quantity.
Finally, taking into account substitution (1), multiplying numerator and denominator by e 4s , and substitution (2), I 3-exch (Equation (A8)) becomes: Owing to the approximation of Equation (8b) the integral can be rewritten in the form: . Solving the integral yields: The corresponding exchange energy E exch = 8πt A I exch with I exch = 3 i=1 I i-exch reads: As r sky →0, B and P→0 so that: Hence, the limit 8π t A is recovered representing the absolute minimum of the exchange energy of a continuous spin structure with integer topological charge [42,44].

Calculation of DMI Energy
This subsection is devoted to the calculation of DMI energy. As shown above this calculation is valid for both IDMI and bulk DMI energy. Taking into account substitution (1), multiplying the numerator and denominator by e 4s , and substitution (2), I 1-DMI (Equation (A9)) becomes: Because of the approximation of Equation (8b) I 1-DMI can be rewritten as: 1+W(C Pξ) . The computation of the integral leads to: Taking into account substitution (1), multiplying numerator and denominator by e 4s , and substitution (2), I 2-DMI (Equation (A10)) becomes: According to the approximation of Equation (8b) the integral can be rewritten in the form: Carrying out the integration yields: Taking into account substitution (1), multiplying numerator and denominator by e 2s , and substitution (2), I 3-DMI (Equation (A11)) reads: The integral is rewritten in the form: Solving the integral yields: Taking into account substitution (1), multiplying numerator and denominator by e 2s , and substitution (2), I 4-IDMI (Equation (A12)) is rewritten as: where the −(+) in front refers to chirality χ = +1 (χ = −1). The integral can be approximated in the form taking into account the approximation of Equation (8b): The calculation of the integral yields: The total DMI energy is expressed as E DMI = ±4π t D l exch I DMI with I DMI = 4 i=1 I i-DMI and reads:

Calculation of the Effective Anisotropy Energy
Taking into account substitution (1), multiplying numerator and denominator by e 2s , and substitution (2), I eff ani (Equation (A13)) turns out to be: Adding and subtracting 4(Bξ) 2 x 2 at the numerator of the integrand function yields: In particular: and: I eff ani1 is solved exactly yielding: Because of the approximation of Equation (8b), I eff ani2 is rewritten as: Solving the integral one gets: with I eff ani = I eff ani1 + I eff ani2 > 0. Substituting Equations (17e) and (17g), the effective anisotropy energy E ani = −2π (l exch ) 2 t K eff I eff ani reads: This term must be added the quantity K u V that leads to a rigid shift of the skyrmion internal energy.

Calculation of the External Field Energy
Taking into account substitution (1), multiplying numerator and denominator by e 2s , and substitution (2), I extfield (Equation (A14)) takes the form: Adding and subtracting 2(Bξ) 2 at the numerator of the integrand function yields: In particular: and: I extfield1 is solved exactly yielding: Owing to the approximation of Equation (8b), singling out the divergence of the integrand function in x = 0 and checking numerically the integral of Equation (19d), I extfield2 is rewritten as proportional to the coefficient b weighted by 1/(2ξ 3 ): with ε→0 and x (ε,C). Solving the integral one gets: The Zeeman energy due to the external bias field reads: Attention must be paid to the computation of E extfield because of the last term strictly depending on the value of ε with ε→0. In the numerical calculations, from the comparison with the exact calculation it is found that ε = 0.001. The total energy E = E exch + E IDMI + E ani + E extfield of the magnetic chiral skyrmion (either Néel or Bloch magnetization texture) in the presence of a perpendicular external magnetic field is the sum of the energy contributions appearing in Equations (11a), (16), (18), and (20). It has been shown that, according to the most accurate magnetization distribution for chiral skyrmions hosted in ultrathin cylindrical dots available, the total energy can be expressed as a combination of transcendental functions.
As magnetic parameters M s , A, D, and K u vary with T according to the following scaling laws [35] also ξ = ξ (T), the exchange l exch = l exch (T), the upper integral limit C = C(T). Moreover, the total skyrmion energy is E = E (T), an internal energy depending on temperature. Note that the thermal effects incorporated in the temperature-dependent magnetic parameters affect also the skyrmion modes. Indeed, they activate not only the skyrmion breathing mode, the only mode that preserves the skyrmion symmetry, but also other skyrmion modes such as, e.g., the skyrmion translational mode that, in principle, breaks the skyrmion symmetry.
The skyrmion radius can be obtained by minimizing the internal energy E, ∂E/∂r sky = 0, and ∂ 2 E/∂r sky 2 > 0. For the sake of convenience, all terms of the skyrmion energy are divided by 8 π t A The energy minimization condition is rewritten as ∂g(P)/∂r sky = 0, with the dimensionless quantity g(P) = E/(8 π t A) having the meaning of a reduced energy given by: where the + (−) sign in front of the coefficient α/ ξ 2 multiplies the − (+) signs inside the round brackets, b = b (CPξ), c = c (CPξ), and the − (+) sign in front of the coefficient γ/ ξ 2 refers to m z = −1 (m z = +1) in the core center (r = 0) and m z = +1 (m z = −1) at the skyrmion border (corresponding to S = −1 (S = +1)).
Since it has been found from the numerical calculations that the equilibrium skyrmion radius, r sky (E = E min ) << r d for the range of temperatures investigated with E min the skyrmion energy minimum, b is expanded to the second order, viz. b ≈ (CPξ) 2 and c to the first order, viz. c ≈ CPξ getting: As P = P(r sky ) it is r sky = 0 so that the minimization condition is rewritten in the form: The minimization condition yields a fifth-degree equation in P. This gives two physical solutions, the smallest one corresponding to the equilibrium skyrmion radius r 0 sky at the minimum skyrmion energy (∂ 2 E/∂r sky 2 > 0) and the larger one to the skyrmion radius at the energy maximum ((∂ 2 E/∂r sky 2 < 0) and three unphysical solution.
To calculate r 0 sky (defined as m z (r 0 sky ) = 0) from Equation (22a), some reasonable approximations have been made in order to get a simple analytical expression of the equilibrium skyrmion radius. First, all terms of the form 1 + powers of (Pξ) were approximated to 1, taking into account that in correspondence of the energy minimum the value of the skyrmion radius is such that P→0 (B << C, B→0). Second, arccot[Pξ] was expanded to the zero-order, viz. arccot[Pξ] ≈ π/2, introducing a coefficient ∆ = π/2 + η (T) with ∆ = ∆ (T). Here, η is a very small coefficient (see Section 3 for the numerical details) depending on T and decreasing approximately with the same rate as ξ (T) that includes the effects of the other approximations and that, for the sake of simplicity, was added, as a small correction, to π/2. For the minimization calculation at H ext 0, the term proportional to ln (P ξ) 2 1+(P ξ) 2 − ln ε C 2 that appears in the expression of the derivative ∂g(P)/∂P has been neglected, assuming that B = C P has the same order of infinitesimal of ε (ε = 0.001) for the values of B crucial for searching the equilibrium radius ranging between 0 and those corresponding to the equilibrium radius itself that are still very small if compared to C. Owing to these approximations, it is possible to extract the equilibrium skyrmion radius from the following first degree equation in B via B = P C: The equilibrium skyrmion radius in dimensionless units is calculated from B 0 = r 0 sky exp (ξ r 0sky ) .
In particular, r 0sky in the absence of an external magnetic bias field, γ = 0, reads: Instead, r 0sky when the skyrmion is subjected to a perpendicular external magnetic bias field, γ 0, takes the form: with α = α (T), β = β (T), γ = γ (T), and ∆ = ∆ (T). The skyrmion radius is proportional to the Lambert function W(y) depending on the magnetic parameters that are, in turn, scaled as a function of T and the parameter ξ. Here, ±2γ refers to m z = −1 (m z = +1) in the core center (r = 0) and m z = +1 (m z = −1) at the skyrmion border (corresponding to S = −1 (S = +1)). Moreover, 8β −2αξ ∆ ± 2γ/ξ 3 +ξ 2 > 0 for the scaled parameters used at every T and for any external magnetic field investigated ensuring the positivity of the argument y of W(y). This result is general and valid for both magnetization textures. Note that the determination of the skyrmion size as a function of the magnetic parameters has been the subject of several investigations [45][46][47]. The expressions of r 0 sky of Equations (22c) and (22d) have some similarities with other expressions of the equilibrium skyrmion radius recently derived in the literature from the minimization of the skyrmion energy, where r 0 sky is a function of the magnetic parameters appearing in a fractional form [45,46]. However, note that, in the present approach, the equilibrium magnetization distribution is an accurate 2D distribution that recovers for Q = 1 and D = 0 the Belavin-Polyakov soliton solution, while in those studies the magnetization distribution had the form of an approximated cosine-like variation or resulted from a 360 • domain-wall distribution.
According to Equations (22c) and (22d), as T increases, ξ decreases for the scaled parameters used and W(x) increases leading to an increase of r 0sky for any fixed external field amplitude and for any magnetization texture. The equilibrium skyrmion radius in dimensional units reads R 0 sky = r 0 sky l exch at any T and for any external magnetic field and the equilibrium skyrmion diameter is D 0 sky = 2 R 0 sky .
For the case considered in the numerical calculations (Nèel skyrmion with S = −1 and radially outward magnetization, χ = +1), it has been found that this behavior is confirmed. This trend characterizing the equilibrium skyrmion diameter as a function of T is indeed consistent with the one obtained calculating the skyrmion energy via the numerical evaluation of the integrals and of the equilibrium skyrmion diameter (for the comparison of the equilibrium skyrmion diameter calculated by means of Equations (22c) and (22d) and numerically, see Section 3.1). Because of the invariance of the DMI energy, this trend occurs also for a Nèel skyrmion with S = +1 and radially inward magnetization, χ = −1. However, due to the generality of Equations (22c) and (22d), this behavior characterizes r 0sky also for a Bloch skyrmion (either S = +1 and counter-clockwise chirality, χ = +1, or S = −1 and clockwise chirality, χ = −1).
Moreover, at fixed T, r 0sky has different trends with increasing H ext (either positively or negatively). In particular, for a Néel or a Bloch skyrmion with S = −1 and χ = +1 and H ext > 0 (H ext < 0), r 0sky reduces (increases) with increasing the positive (negative) amplitude of the external magnetic field leading to a positive (negative) increase of γ, γ > 0 (γ < 0), confirming recent micromagnetic and numerical calculations [43].
Finally, according to Equations (22c) and (22d), at fixed magnetic parameters for a given T, the equilibrium skyrmion radius does not have an explicit dependence on the dot radius. By means of a comparison with the equilibrium radius calculated from the numerical minimization of the skyrmion energy it has been found that, for the magnetic parameters used, this is approximately valid for a dot radius R d > 50 nm. For R d ≤ 50 nm R 0 sky depends on the dot radius and reduces with decreasing R d as found by Tejo et al. [44] so that, in this regime Equations (22c) and (22d) are not anymore valid.

Statistical Thermodynamic Properties of Skyrmion Diameters Population
In this section the statistical thermodynamics aspects of chiral magnetic skyrmions at equilibrium are investigated taking into account the strict analogy with the particles behavior in an ideal gas. In particular, key quantities to understand their thermodynamic behavior such as the partition function and the free energy are calculated within a microcanonical ensemble. Finally, the pressure is determined and an equation of state linking pressure, volume, and temperature is proposed.

Partition Function and Free Energy of a Skyrmion Diameters Population within a Microcanonical Ensemble
Let us calculate the partition function and the free energy of a skyrmion diameters population within a microcanonical ensemble. Very recently, it has been shown that, starting from a canonical ensemble, a skyrmion diameters population can be approximately described within a microcanonical ensemble due to the small fluctuations of the energy around the average energy. To calculate the partition function, the expression of the configurational entropy at thermodynamic equilibrium resulting from an average over the 3D equilibrium Maxwell-Boltzmann distribution of skyrmion diameters is recalled [39]: with a = a(T) a parameter proportional to the parabola curvature, <D sky > = <D sky (T)> the average skyrmion diameter, and S 0 = k B (2 ln 2+1/2 ln π). <D sky > was calculated as a spatial integral of D sky over the 3D Maxwell-Boltzmann distribution [39]. The reason for performing the analysis in a 3D space instead of a 2D (the magnetization distribution is indeed planar) results from the observation, via micromagnetic simulations, that the thickness of the nanodot (even though less than 1 nm) can be important to establish the behavior of the diameters distribution. Strictly speaking, Equation (23) was derived for a Nèel skyrmion with S = −1. However, as already outlined in [39], this expression can be considered valid also for a Nèel skyrmion with skyrmion number S = +1 and applies also to the Bloch skyrmion with S = ±1. Indeed, the energy takes the same form for any magnetization texture and its trend in the vicinity of the minimum can be approximated by a quadratic dependence on the skyrmion diameter. Note that, as for the ideal gas, configurational entropy expressed by Equation (23) fulfills the requirement of additivity (entropy is indeed an extensive thermodynamic quantity), is consistent with the definition of temperature and with the second principle of thermodynamics, and with the adiabatic invariance for slow changes occurring in the system. On the other hand, due to its classical derivation, S does not fulfil the third principle of thermodynamics or Nernst principle according to which S = 0 J/K for T = 0 K, but S → −∞ as T→0 K, showing a similar behavior to that of the Sackur-Tetrode equation expressing the entropy for an ideal gas as a function of T [39].
Equating the calculated configurational entropy to S = k B lnW representing the general entropy definition valid for a microcanonical ensemble, W can be determined. In statistical mechanics, the quantity W represents the number of microscopic possibilities to realize the macroscopic state. More specifically, W can be regarded as the statistical multiplicity of the energy level having value <E> or as the degeneracy of the ground state of the skyrmions population having average energy <E>. One gets: where K = 4 √ π e 3/2 /t is a constant at fixed dot thickness t. The statistical multiplicity of a skyrmion diameters population depends not only on T but also on the geometric and magnetic parameters of the magnetic system that, in turn, have an intrinsic dependence on T.
It can be noted that W = 0 for T = 0 K: at the absolute zero the ground state is no longer degenerate. With increasing T there is an enhancement of the degeneracy. This means that the higher the temperature the greater the probability to populate the ground state levels of skyrmion diameters population.
In [39], it has been also shown that the partition function of the skyrmion diameters population within a canonical ensemble can be approximated by the one of the microcanonical ensemble that in the discrete limit reads Z ≈ W e − <E>/(k B T) . Hence, replacing W given in Equation (24), the explicit expression of the partition function within a microcanonical ensemble is obtained: where <E> (5 k B T+2a(D 0 sky ) 2 )/2 is the average energy calculated in [39] as <E> a <D sky 2 > and D 0 sky = D 0 sky (T) is the equilibrium skyrmion diameter with D 0 sky <D sky >, especially at low T. Looking at Equation (25) it can be observed that the partition function of a skyrmion diameters population depends on fractional powers of T. However, there is a combined dependence on T 1/2 and T 3/2 if compared to the partition function for an ideal gas whose particles follow a 3D MB distribution where there is only a T 3/2 dependence. Moreover, the fractional dependence on T appearing in the rational fraction is weighted by an exponential term that also depends on T. This exponential dependence differentiates the partition function Z of a skyrmions population from that of the particles of a 3D ideal gas.
Similarly to the configurational entropy, the partition function has a combined dependence on geometric and magnetic parameters. Taking into account the relation between the Helmholtz free energy F and the partition function of the microcanonical ensemble expressed as F = −k B T lnZ, the relation F <E> − T S can be easily found, with F = F(T) for every T (the symbol results from the approximated form of Z), the well-known thermodynamic relation generally used for ferromagnets. For this case, because of the definition of S as an expectation value, F, while depending on T, is regarded as an average free energy. Note that the Zeeman energy is included in <E> and, for this reason, still one deals with the Helmholtz free energy F. However, if this contribution is thought as separated from the other skyrmion energy contributions to the average energy one would deal with the Gibbs free energy G. Within this framework, the two approaches are equivalent.

Skyrmion Diameters Population Pressure and Equation of State
From the free energy, it is possible to derive, in complete analogy with the pressure exerted by the particles of a gas on the walls of a container, the pressure p of the skyrmion diameters population.
Taking into account the infinitesimal relation dF = − S dT − p dV it is p = − (∂F/∂V) T . In explicit form: Equation (26) gives the pressure generated by skyrmion diameters population at a fixed temperature T. If a container filled with an ideal gas is expanded instantaneously, the temperature of the gas does not change at all. Analogously, if the skyrmion size is allowed to fluctuate instantaneously, the skyrmion diameter deviates from the average value and, as a result, the skyrmion volume fluctuates around its average value. This process occurs without affecting the temperature, as occurs for particles in an ideal gas.
To calculate the pressure, one could substitute either F <E> − T S or F = −k B T ln Z in Equation (26). However, unlike for the case of an ideal gas, the thermodynamic variables <E>, S, and Z within this framework have a dependence not on the generic volume V but on the average skyrmion volume <V> with in addition <V> = <V(T)>, that is the dependence on V is mixed with that of T. Therefore, it is convenient to argue in terms of densities of thermodynamic variables. Let us introduce a free energy density f = f (D sky ) dependent on skyrmion diameter with f = e − T s, where e = e(D sky ) is the internal energy density (internal energy per unit volume e = E/V) and s = s(D sky ) the configurational entropy density (configurational entropy per unit volume). These quantities are expressed taking into account the approximation on the skyrmion energy in the vicinity of the minimum via a quadratic dependence [39] on D sky , viz. E = a (D sky − D 0sky ) 2 used to calculate <D sky > and the skyrmion configurational entropy expressed in Equation (23). Note that, for the total computation of the energy, an energy shift of b = E (D = D 0 sky ) would be added to this expression that gives the energy in the minimum. This harmonic approximation is reasonable because the main contribution to the configurational entropy results from skyrmion diameters around the energy minimum. In particular: e a (D sky − <D sky >) 2 /V (27) where V = 1 4 π D sky 2 is the skyrmion volume for a given D sky and the symbol was introduced because, for consistency with the definition of entropy density (see below), D 0 sky has been replaced with <D sky >. For the scaled parameter used, this approximation holds in the region of metastability being the difference between <D sky > and D 0 sky less than 10 %. The volume V is at most equal to the Instead, the entropy density can be derived from the definition of the Boltzmann order function at equilibrium equivalent (apart from the sign) to the Shannon information entropy. First, the Gaussian diameters distribution or probability density appearing in the Boltzmann order function at equilibrium , with ∆D sky = D sky − <D sky > expressing the deviation from the average skyrmion diameter at a given T and N the normalization constant. The normalization constant has the dimension of an inverse of a volume and is obtained by normalizing to one skyrmion having different sizes and diameters in different instants of time, giving rise to a skyrmion diameters population at each T in the range of temperatures corresponding to the metastable state, as shown by micromagnetic simulations [39]. Here, f 0 has the meaning of a Gaussian probability density. In principle, one should normalize the Gaussian distribution in the usual way, performing a volume integration of f 0 . This results in a value of N that has an explicit dependence on T, on <D sky >, and on magnetic parameters through a, and is crucial for the calculation of the configurational entropy. However, to have a strict connection with the magnetic skyrmion volume V = V(D sky ) that strictly depends on the skyrmion diameter D sky for each T, one takes, without loss of generality, C =1/V. Hence, . Note that f 0 is not anymore, strictly speaking, a Gaussian distribution because V∝ (D sky ) 2 leading to a singularity of f 0 for vanishing D sky . However, for diameters different from zero its trend is almost superimposable to that of the corresponding Gaussian, with a slight downshift of the maximum (with respect to the maximum for D sky = <D sky > of the Gaussian) that increases with increasing temperature and a slight variation of the standard deviation and of the full width at half maximum.
In view of the above arguments the entropy density s = −k B f 0 ln f 0 can be written in the form: This expression can be derived from the integrand of the Boltzmann order function H 0 defined in the continuum limit as an integral over spatial coordinates and evaluated at equilibrium replacing the normalization constant C obtained integrating the Gaussian distribution over the spatial coordinates [39] with the volume V of the skyrmion and rescaling lnf 0 to ln(f 0 <V>) for dimensional reasons. According to these definitions, the pressure is rewritten in the form: where <V> 1 4 π <D sky > 2 is the average volume according to the approximation <D sky 2 > ≈ <D sky > 2 that holds for the range of temperatures considered. Indicating the source of pressure having energy nature with p E = −<V> (∂e/∂V) T , one obtains substituting Equation (27): One notes that this pressure contribution is always positive and, at fixed a, it decreases with increasing V.
Labeling with p S = <V> T (∂s/∂V) T the source of pressure having entropic nature, via Equation (28), yields: with <V>/ V = <D sky > 2 /D sky 2 . This contribution is negative for low volumes, positive for intermediate volumes close to the average volume, and again positive for high volumes bigger than the average volume and is thus responsible for the oscillatory behavior of the pressure curve, especially the one present at high volumes (see Section 3.4). Moreover, p S is about one order of magnitude less than p E . Pressure p is obtained by summing the two contributions of Equations (29b) and (29c), viz. p = p S + p E : where the symbol results from the approximated energy contribution to the pressure. Straightforwardly, p written in terms of dimensionless variables reads: dimensionless variable and f c (∆D sky >) = 1 − e -a (∆D sky )2/k B T . At fixed T, pressure reduces with increasing V but exhibits a superimposed fluctuating behavior for V > <V> for each T (for the details of the numerical calculations, see Section 3) due to the statistical dependence on the skyrmion diameters contained in the coefficient of k B T and in the second term on the second member. This statistical aspect differentiates the skyrmion population behavior from that of pressure exerted on the walls of a box by the particles of an ideal gas that has a monotonic behavior as a function of V for fixed T. By inspecting Equations (29b) and (29c), one notes that both sources of pressure depend on the confining skyrmion potential or skyrmion energy whose curvature in the vicinity of the absolute minimum is determined by the coefficient a, in turn depending on the magnetostatic field and the exchange interactions. From Equation (29e) one gets: Equation (30) is the equation of state for a skyrmion diameters population and is the main result of this study. The second member contains a term proportional to T weighted by a coefficient in turn depending on f (∆D sky >) and a term proportional to the square deviation of D sky from <D sky > in turn weighted by f c (∆D sky >). This equation describes any single chiral skyrmion diameters population independently of its magnetization texture and ferromagnetic material considered.
Straightforwardly, for D sky = <D sky > it is p =k B T/V, the pressure exerted by a single particle (N = 1) on the walls of the container of volume <V> = V that, in turn, yields pV = k B T, the equation of state for an ideal gas for N = 1.
It is interesting to derive the asymptotic behavior of the equation of state in the limit for a(D sky − <D sky >) 2 << k B T for small fluctuations of diameters around the average value <D sky >. This regime is named the small fluctuations regime. This occurs at low temperatures, especially in the presence of an applied bias field characterized by sharp Gaussian distribution [35].
For a(∆D sky ) 2 << k B T, the exponential e −a (∆D sky )2/k B T in Equation (30) is expanded to the first-order in (∆D sky ) 2 and the fourth-order term a 2 (∆<D sky >) 4 /k B T is neglected being a 2 (∆<D sky >) 4 << k B T, yielding: Equation (31) is the equation of state of a skyrmion diameters population in the small fluctuations regime and is valid for values of D sky very close to <D sky >. For D sky = <D sky > one gets g(v) = 1 and v = 1 so that again the ideal gas law pV = k B T for N = 1 is obtained.

Results and Discussion
In this section the main numerical results are discussed. The numerical calculations were carried out for an outwardly Neel skyrmion with S = −1 (magnetization in the core center along −z) hosted in a Co dot of R d = 200 nm and t = 0.8 nm using the following magnetic parameters at T = 0 K: saturation magnetization M s = 600 KA/m, exchange stiffness constant A = 20 pJ/m, IDMI parameter D = 2.0 mJ/m 2 , K u = 0.6 MJ/m 3 (for example Co), and the scaled values at T 0 according to the scaling laws expressed in Section 2 [35]. However, note that the model can be numerically applied to other magnetization textures present in ultrathin cylindrical dots. In Table 1 the values of a and <D sky > used in the numerical calculations presented in this Section are summarized [39].  The values of the equilibrium diameters for every temperature and external bias field are listed in Table 2.

Numerical and Analytical Equilibrium Diameters: A Comparison
In this subsection the equilibrium diameters calculated by means of Equations (22c) and (22d) via D 0sky = 2R 0sky are compared with the ones determined numerically. The numerical calculations were performed minimizing the skyrmion energy of Equation (7) using Equations (A6)-(A14) and determining the equilibrium radius as a function of temperature and for different bias external fields. The results of this comparison are shown in Figure 2. The overall agreement is very good for the whole range of temperatures investigated, with a more pronounced discrepancy, especially at high temperatures.  The values of the equilibrium diameters for every temperature and external bias field are listed in Table 2.

Numerical and Analytical Equilibrium Diameters: A Comparison
In this subsection the equilibrium diameters calculated by means of Equations (22c) and (22d) via D0sky = 2R0sky are compared with the ones determined numerically. The numerical calculations were performed minimizing the skyrmion energy of Equation (7) using Equations (A6)-(A14) and determining the equilibrium radius as a function of temperature and for different bias external fields. The results of this comparison are shown in Figure 2. The overall agreement is very good for the whole range of temperatures investigated, with a more pronounced discrepancy, especially at high temperatures.  In the numerical calculations, the values of the coefficient η appearing in ∆ = π/2 + η of Equations (22b), (22c) and (22d) summarized in Table 3 were used. The coefficient η was fitted at T = 0 K to the value of D 0 sky evaluated numerically for three different amplitudes of the external field and linearly reduces with T with the same rate of decrease of ξ with T. Table 3. Calculated values of η used in the numerical calculations. These results are easily generalizable to the case with skyrmion number S = +1 and other magnetization skyrmion textures (e.g., to a Bloch skyrmion) in the region of metastability.

Energy of the Skyrmion State and of Perpendicular Uniform State: A Comparison
In Figure 3a, the skyrmion energy evaluated at the energy minimum as a function of temperature in the skyrmion state in the region of metastability is displayed and compared to the skyrmion energy of the ideal uniform state in the presence of an external magnetic field aligned along +z (magnetization distribution along +z in the whole dot and parallel to the external magnetic field).

Microcanonical Partition Function and Free Energy of the Skyrmions Population
The partition function and the Helmholtz free energy of the skyrmion diameters population studied within a microcanonical ensemble are shown in Figure 4. The partition function was calculated according to Equation (25) by varying the temperature T at fixed external magnetic field. The Helmholtz free energy is computed via the relation F = −kBT lnZ. The values of the average skyrmion diameters entering in the expressions of Z and F are the ones listed in Table 1 for the three different amplitudes of the external bias field, μ0Hext = 0 mT, μ0Hext = 25 mT, and μ0Hext = 50 mT.
Looking at Figure 4a, ln Z increases with increasing T, reflecting the increase of Z with T, a number ranging between 0 and 1. Instead, one notes that the effect of the external magnetic field is to The skyrmion energy was calculated by numerically integrating Equation (7) and expressing the total skyrmion energy using Equations (A6)-(A14). Instead, the energy of the ideal uniform state along +z is expressed as: As expected, the energy of the uniform state is downshifted with respect to that of the skyrmion state at each T [44]. Between the two minima, there is an energy barrier, mainly due to the ferromagnetic exchange and the i-DMI contributions, and the energy minimum of the uniform state is lower with respect to that of the skyrmion state. This behavior marks the metastability of the skyrmion state. In particular, with increasing temperature, the two minima tend to merge. To benchmark the model, the equilibrium diameters computed as a function of the external bias field (continuous black line) have been compared with the ones obtained with spin-polarized scanning tunnel microscopy (red circles) and the result of the comparison is shown in Figure 3b. Also in this case, the equilibrium diameters were calculated by numerically minimizing the skyrmion energy of Equation (7) using Equations (A6)-(A14). The geometric and magnetic parameters used in the calculations are: R d = 25 nm, t = 0.408 nm, A = 2.0 pJ/m, D = 3.9 mJ/m 2 , K = 2.5 MJ/m 3 , and M s = 1.1 MA/m. The magnetic parameters are typical values for thin-film ferromagnetic systems, such as the bilayer of PdFe on Ir (111) single crystal substrate [48]. Note that in that case the external field is anti-parallel to the magnetization at the center of the skyrmion core but with m z (r = 0) = +1 and H ext along −z. This configuration is equivalent to the one studied in this work with m z (r = 0) = −1 and H ext along +z. The agreement is very good for the whole interval of external bias fields studied.

Microcanonical Partition Function and Free Energy of the Skyrmions Population
The partition function and the Helmholtz free energy of the skyrmion diameters population studied within a microcanonical ensemble are shown in Figure 4. The partition function was calculated according to Equation (25) by varying the temperature T at fixed external magnetic field. The Helmholtz free energy is computed via the relation F = −k B T lnZ. The values of the average skyrmion diameters entering in the expressions of Z and F are the ones listed in Table 1 for the three different amplitudes of the external bias field, µ 0 H ext = 0 mT, µ 0 H ext = 25 mT, and µ 0 H ext = 50 mT.  The effect of the external magnetic field leading to a downshift of Z and F would be similar considering the case with Hext applied along -z and mz(r = 0) = +1 (S = +1). Instead, an opposite behavior leading to an upshift of the same thermodynamic quantities in the cases of an applied field parallel to mz (either along −z (mz(r = 0) = −1 and S = −1) or along +z (mz(r = 0) = +1 and S = +1) should be expected.

Skyrmions Population Pressure and Equation of State
The pressure of the skyrmions population as a function of the volume V of the skyrmion diameters population is shown in Figure 5 at fixed T. For each temperature T at fixed external bias field, volumes considered correspond to the diameters belonging to the range of about ±3σ<Dsky>, Looking at Figure 4a, ln Z increases with increasing T, reflecting the increase of Z with T, a number ranging between 0 and 1. Instead, one notes that the effect of the external magnetic field is to reduce Z. This can be understood taking into account the physical meaning of the partition function in statistical thermodynamics. Z encodes how the probabilities are partitioned among the different microstates, based on their individual energies. The partition function Z for T→0 K vanishes and, for T 0 K, increases with T. Due to the disordering effect introduced by the temperature on the thermodynamic system, with increasing T the partitioning probability among different microstates of the skyrmion diameters population increases. The ordering effect of the external bias field as in the case of configurational entropy [39] partially contrasts with the disordering effect of T attenuating the increase of Z with T.
On the other hand, the trend of the Helmholtz free energy (with F > 0) as a function of T in the presence of an external bias field is almost constant while, in its absence, F decreases with T. F merges asymptotically towards a value of about 6 × 10 −21 J as T→0 K (see Figure 4b), masking almost completely the effect of H ext . The behavior of the Helmholtz free energy F as a function of T is different with respect to the one of <E> that increases with T.
The effect of the external magnetic field leading to a downshift of Z and F would be similar considering the case with H ext applied along -z and m z (r = 0) = +1 (S = +1). Instead, an opposite behavior leading to an upshift of the same thermodynamic quantities in the cases of an applied field parallel to m z (either along −z (m z (r = 0) = −1 and S = −1) or along +z (m z (r = 0) = +1 and S = +1) should be expected.

Skyrmions Population Pressure and Equation of State
The pressure of the skyrmions population as a function of the volume V of the skyrmion diameters population is shown in Figure 5 at fixed T. For each temperature T at fixed external bias field, volumes considered correspond to the diameters belonging to the range of about ±3σ <Dsky> , where σ <Dsky> is the standard deviation of the Gaussian distribution centered at the average diameter <D sky >. In Figure 5a, pressure p vs. volume V in the absence of an external bias field at different temperatures (isotherms) is displayed. p was determined by means of Equation (29e) and its calculation takes into account the Gaussian distribution of skyrmion diameters fluctuations from <D sky > at each temperature and external bias field, with <D sky > given in Table 1 and the volume approximately computed as V ≈ 1 4 π < D sky > 2 .
The general trend at all temperatures is the dramatic increase of p with decreasing V below a given volume, depending on the temperature T. This means that, for a skyrmion of reduced size, the pressure exerted on the region of the ferromagnet just outside the skyrmion core is very high if compared to a skyrmion of intermediate size. Moreover, at fixed volume, pressure reduces with reducing the skyrmion temperature in a way similar to that of the particles of an ideal gas with decreasing T.
Looking at Figure 5b, it can be observed that the pressure curve (depicted for T = 100 K) in the presence of an external magnetic field along −z that is anti-parallel to the static magnetization at the core center (m z (r =0) = +1) at fixed volume shifts towards lower values so that, at a given value of V, the effect of a bias field is to reduce pressure in a way similar to the reduction of the configurational entropy [39]. This trend is also typical for other temperatures and higher external magnetic fields. The skyrmion pressure curves intersect the curve p =k B T/V of the ideal gas at V = <V>: for µ 0 H ext = 0 mT it is V = 7.78 × 10 −25 m 3 , for µ 0 H ext = 25 mT it is V = 4.70 × 10 −25 m 3 and for µ 0 H ext = 50 mT it is V = 3.28 × 10 −25 m 3 . One notes that the intersection with the ideal gas curve occurs at lower volumes in the presence of an external magnetic field. At intermediate and high volumes (V > <V>) there is an oscillatory behavior of the pressure curve superimposed to the monotonic behavior that is enhanced with increasing the amplitude of the external bias field (see Figure 5c). Similar conclusions are drawn for other values of temperature in the region of skyrmion metastability. Because of this oscillatory behavior introduced by the diameters distributions f and f c , unlike for the universal law for an ideal gas characterized by the universal constant R = PV/T at any T, the classical thermodynamics of a skyrmions population cannot be strictly characterized by a universal constant. external magnetic field applied along −z and a Néel skyrmion with the magnetization at the centre of the core mz(r = 0) = +1 (S = +1). Instead, for the cases where the the effect is to increase the skyrmions population pressure at a given volume V.
These results can be easily generalized to other cases (e.g. external magnetic field is parallel to mz(r = 0) either along +z or along −z) and to Bloch skyrmions in the region of metastability, leading to conclusions about the qualitative trends of p obtained similar to the ones drawn for the case studied.  For values of V close to <V>, Equation (31) expressing pressure P in the small fluctuation regime can be applied. A similar shift towards lower values of the pressure curve would occur for the external magnetic field applied along −z and a Néel skyrmion with the magnetization at the centre of the core m z (r = 0) = +1 (S = +1). Instead, for the cases where the the effect is to increase the skyrmions population pressure at a given volume V.
These results can be easily generalized to other cases (e.g., external magnetic field is parallel to m z (r = 0) either along +z or along −z) and to Bloch skyrmions in the region of metastability, leading to conclusions about the qualitative trends of p obtained similar to the ones drawn for the case studied.

Conclusions
In this paper, the equilibrium statistical mechanics and thermodynamic properties of a skyrmion diameters population in an ultrathin cylindrical dot were investigated. An analytical expression of the internal skyrmion energy valid for any type of magnetization texture, consisting of a combination of transcendental functions, was obtained. From the minimization of the internal skyrmion energy it was found that the equilibrium diameter of a single chiral skyrmion, both in the hedgehog and vortex-like configurations, can be expressed in terms of a Lambert function whose argument depends explicitly on the magnetic parameters characterizing the magnetic skyrmion at each temperature. Unlike other formalisms based on 1D domain wall distributions, this was accomplished starting from a recently proposed 2D equilibrium magnetization that reduces to the Belavin-Poliakov solitonic solution in the isotropic case and for dominating exchange interaction. Exploiting the analogy between the behavior of particles in a gas and a population of skyrmions diameters, the Helmholtz free energy and the partition function were calculated treating the diameters population within the framework of a microcanonical ensemble. It was shown that, like for particles in an ideal gas, the skyrmions population is also characterized by a pressure strictly depending on the confining potential at each temperature, and in turn related to the exchange and magnetostatic interaction characterizing magnetic skyrmions. It was found that skyrmions pressure decreases monotonically with the volume but for volumes higher than the average volume an oscillatory behavior is superimposed onto the monotonic behavior because of the entropic contribution to the pressure that is one order of magnitude less than the average energy contribution. This trend occurs at each temperature and external field amplitude. An equation of state for a skyrmions population, regarded as the equivalent of the equation of state for an ideal gas, was derived and it was shown that it is valid for any type of magnetization texture. The equation of state reduces to that of an ideal gas when the skyrmion volume equals the average skyrmion volume at each temperature and external bias field.
These theoretical results on the thermodynamic properties of chiral magnetic skyrmions could pave the way for calorimetric measurements that allow determining the heat released or absorbed by the ferromagnetic system hosting the skyrmion and could enable the measurement of the free energy, exploiting the relations between heat exchanged and entropy, and heat exchanged, internal energy, and work done by the system. In this way, a direct measurement of the specific heat characterizing a population of skyrmions in a ferromagnetic system could also be carried out. These measurements would also enable verifying the predicted equation of state characterizing the statistical thermodynamic behavior of a skyrmions population and to experimentally determine the corresponding pressure. Finally, the analysis extended to very low temperatures would allow exploring quantum effects on the thermodynamic variables.
Funding: This research received no external funding.