The skyrmion-bubble transition in a ferromagnetic thin film

Magnetic skyrmions and bubbles, observed in ferromagnetic thin films with perpendicular magnetic anisotropy, are topological solitons which differ by their characteristic size and the balance in the energies at the origin of their stabilisation. However, these two spin textures have the same topology and a continuous transformation between them is allowed. In the present work, we derive an analytical model to explore the skyrmion-bubble transition. We evidence a region in the parameter space where both topological soliton solutions coexist and close to which transformations between skyrmion and bubbles are observed as a function of the magnetic field. Above a critical point, at which the energy barrier separating both solutions vanishes, only one topological soliton solution remains, which size can be continuously tuned from micrometer to nanometer with applied magnetic field.


Introduction
Skyrmions, are topological solitons which present particle-like properties: they have quantized topological charges, interact via attractive and repulsive forces, and can condense into ordered phases. The concept of skyrmions has spread over various branches of physics [1] including condensed matter, as for example in the case of liquid crystals [2], quantum Hall magnets [3,4] and Bose-Einstein condensates [5]. In ferromagnets, skyrmions were originally called two-dimensional (2D) topological solitons or magnetic vortices and their existence has been predicted in isotropic ferromagnets [6], uniaxial ferromagnets [7,8,9] and non-centrosymetric magnets [10,8]. Some early indirect experimental evidences of their existence have been obtained in quasi-2D antiferromagnets [11,12,13]. More recently indirect [14] and direct observations of skyrmion have been reported in chiral magnets [15] and ultrathin layers of conventional transition-metal-based ferromagnets in contact with heavy metals [16,17]. The nanoscale size and non-trivial topology of skyrmions make them particularly attractive for information technologies. The idea of a skyrmion-based memory was already discussed in the 80's [18] at the time when the research on magnetic bubbles was at its apogee. Magnetic bubbles are cylindrical magnetic domains which appear in ferromagnetic films with perpendicular magnetic anisotropy when the demagnetising energy compensates the domain wall (DW) energy [19]. Memories based on magnetic bubbles displaced by a rotating magnetic field were commercialised in the 90's but were progressively replaced by transistor based memories and hard disk drives due to their higher capacity and lower cost. The recent observations of skyrmions at room temperature (RT) and their fast displacement with low electrical currents [20,21,22] has triggered a revival of the quest for a memory based on topological solitons, taking the form of a skyrmion racetrack memory [23,24,25]. Magnetic bubbles and skyrmions are close relatives as they can share the same topology [26]. However, their characteristic sizes differ and while classical bubbles present a long lifetime at RT, much shorter lifetimes were found for nanometer-sized skyrmions in recent experimental [27,28] and theoretical works [29,30,31]. The case of intermediate-size solitons is more favorable, as stable RT topological solitons with sizes of a few hundred to a few tens of nanometers have been reported in multilayers [32] and even in a single ferromagnetic layer [33,34]. These topological solitons are sometimes called skyrmion bubbles when the demagnetising energy plays a role in their stabilization. In this context, the necessity to clarify whether a fundamental difference exist between skyrmions and bubbles appears. In previous works, the difference between magnetic bubbles with a large number of collinear spins in their center and skyrmions with a compact core has been described [35,36]. The dynamic transformation between a magnetic bubble and a skyrmion using a magnetic tip has also been demonstrated [37]. As minimising the total energy of the soliton is numerically expensive, we derive in the present work an analytical topological soliton model in order to build a skyrmion-bubble phase diagram and obtain a better physical insight of the differences between skyrmions and bubbles. This allows us to calculate skyrmions and bubbles equilibrium solutions out of a single model from nanometer to micrometer radius and demonstrate the existence of transitions between them as a function of magnetic field.

Topological soliton model
The model is developed in view of describing isolated skyrmions and bubbles spin textures in an infinite ferromagnetic thin film. As both spin textures have the same topology, we will use the more general name topological soliton to discuss the model solutions in the first place and later specify what type of topological soliton solution presents the characteristic of a magnetic skyrmion or bubble. The system we consider is a 1 nm-thick ferromagnetic thin film with an easy axis perpendicular to the plane. The material is described by its thickness t, its spontaneous magnetisation M s , its exchange constant A ex. , its volume magnetocrystalline anisotropy K u and its micromagnetic Dzyaloshinskii-Moriya interaction (DMI) parameter D. The magnetocrystalline anisotropy and the DMI are written as volume-related quantities K u and D. The soliton is assumed to have a cylindrical symmetry and cylindrical coordinates are used, where ρ is the radial distance and the z−axis is perpendicular to the film. We restrict ourselves to D values sufficiently large so that DWs will have a Néel nature (based on Ref. [39]) and consequently we also assume a Néel-type for the topological soliton (represented in Figure 1b). In these conditions, the magnetisation direction is defined only by the angle θ between M and the z−axis (θ = 0 is up, θ = π is down). A magnetic field µ 0 H may be applied, for which the positive direction is in the +z direction i.e. antiparallel to the magnetisation in the soliton center (see Figure 1a and b). The film is thinner than the exchange length l ex. = 2A ex. /µ 0 M 2 s so we assume that magnetisation does not depend on z. In the following, we will present the energy functional containing all the energy terms required to describe the topological soliton energy (Section 2.1). Then, in a first step, we will use the local approximation for the demagnetising energy, generally used for skyrmions [40,41,26,42], to calculate the topological soliton profile and energy (Section 2.2). In a second step, we will show that analytical expressions of the energy terms can be used to calculate the (2) (black squares) and Eq. (9) (black line). The planar DW energy density σ w is indicated by a dashed line. (e) Energies versus topological soliton radius calculated from Eq. (2) (symbols) and Eq. (9) (lines) with the parameters given in (c) and normalised by the RT thermal energy E RT = k B T 293K . The normalised zero radius energy is indicated as ε 0 = E 0 /E RT ∼ 31. All the programs used to create the data for the figures can be found in a repository [38]. skyrmion solution (Section 2.4) and we will extend the model by taking into account the long range surface demagnetising effect which is at the origin of the bubble solution stabilisation (Section 2.4.4).

Energy functional
The topological soliton energy E s [θ(ρ)] is the sum of 5 terms: exchange energy E exch. , anisotropy energy E anis. , DMI energy E DMI , demagnetising energy E dem. and Zeeman energy E zee. : The demagnetising energy can be decomposed into two terms E dem. = E vol. dem. + E surf. dem. corresponding to the contributions from volume and surface magnetic charges. The volume charges appear only inside the Néel domain wall. The surface charges are present inside the soliton but also in the infinite region with uniform magnetisation. Due to the long range nature of the demagnetising effects, it is analytically impossible and numerically difficult to find the function θ(ρ) which minimizes the E s functional. This has been done in the early 70's, in absence of DMI, in the case of magnetic bubbles stabilised by the surface demagnetising energy [43,44] and by Kiselev and co-workers in presence of DMI [35]. In other works, the demagnetising energy is often taken into account only as a local energy contribution [40,41,26,42] as described in the next Section.

Local approximation for the demagnetising field
We restrict ourselves to the ultrathin film limit where the DW width is much larger than the film thickness and where E vol.
dem. is negligible compared to E surf. dem. . Indeed, in a thin film, the demagnetising contribution from the volume charges inside a Néel DW is proportional to the ratio between the film thickness and the DW width [45]. This consideration is also valid for a Néel skyrmion which possesses a 2π DW cross section. In addition, the film thickness is much smaller than the exchange length and we make a short range uniform magnetisation assumption. In these conditions the surface charge contribution to the demagnetising field is The local surface demagnetising energy becomes equivalent to an anisotropy and the anisotropy constant in Eq. (1) is replaced by an effective anisotropy constant Considering now the energy difference between the topological soliton and the uniform state we obtain the following Euler Equation [26]:

Soliton energies versus radius
We solved the Euler-Lagrange Equation (2) using the shooting method starting from an inverse tangent try function and the two boundaries values : θ(0) = π and θ(∞) = 0. The energetically minimised soliton profile θ(ρ) is shown in Figure 1c. The topological soliton radius r is defined as the ρ value at which θ(ρ) = π/2 (M z = 0). The exchange (blue triangles), anisotropy (red triangles) and DMI energies (black circles) are plotted versus topological soliton radius r in Figure 1e. The Bloch DW width defined as π∆ = π A ex. /K eff. is indicated. The DMI and anisotropy energies present a linear variation versus r down to r < π∆. On the contrary, in the same r range, the exchange energy deviates from linearity and tends toward a non-zero constant value. This is due to the curvature 1/ρ exchange term in Eq. (1) which is related to the rotation of the spins in the (x,y)-plane: when the topological soliton radius is decreased, the angle between two adjacent spins in this plane increases and the exchange energy increases, leading to a non zero exchange energy limit when r → 0. Consequently the topological soliton wall energy density σ s defined as the sum of the exchange, anisotropy and DMI energy densities deviates from the planar DW energy density σ w = 4 √ A ex. K eff. − πD in the low r range as illustrated in Figure 1d. This non-linearity in the exchange term is at the origin of the presence of a local minimum in the total energy in presence of a sufficiently large DMI as we will discuss in Section 3.1. As E anis. , E DMI and E exch. show relatively simple r dependences, we will in the following derive analytical expressions to reproduce these dependences.

Analytical model
We now introduce analytical expressions for the 5 different energy terms which constitute the total energy of the topological soliton. The total energy, as well as each energy term is defined as the energy difference between an isolated topological soliton with a down magnetisation in its core and the uniform ferromagnetic state in the up direction.

Anisotropy and exchange energies
For r π∆ the sum of exchange and anisotropy energies is proportional to 4 √ A ex. K eff. and both energies contribute equally : For r < π∆, the exchange energy deviates from linearity as discussed in Section 2.3 and we modify this expression by adding a low r correction: which gives the exact E 0 = 8πA ex. t zero-radius limit for the exchange energy, (Belavin and Polyakov [6], see also [26,46]) and a zero (r = 0) dE exch. /dr derivative. The analytical expressions of E anis. and E exch. as well as their sum appear as lines in Figure 1.e. This expression reproduces the exchange energy obtained in Section 2.3 with a 3% maximum error in the full radius range.

DMI energy
The DMI energy is proportional to the total π rotation of the spins (from the center to the periphery of the topological soliton) and varies linearly with r as observed in Figure 1.e. It expresses as: We have chosen the rotation chirality which lowers the energy and D > 0. E DMI is negative, thus it favours the expansion of topological solitons.

Zeeman energy
We use an approximate expression for the Zeeman energy of the topological soliton which represents the Zeeman energy difference between a film with a uniform +M s magnetisation, containing a magnetic cylinder of radius r with an opposite uniform −M s magnetisation, and the Zeeman energy of the uniform +M s state.
The error associated with this approximation will be discussed in Section 4.3.

Demagnetising energy
As discussed in Section 2.1 the demagnetising energy, cannot be expressed analytically and approximations have to be used [44]. The local effect of the demagnetising energy in the region where the spin rotates is taken into account using a local approximation and replacing K u by K eff. in Eq. (1) as discussed in Section 2.2. This local approximation neglects the long range demagnetising effect which becomes non negligible as the skyrmion radius grows. This long range demagnetising energy contribution is at the origin of classical magnetic bubbles formation [19] and its role in the stabilisation of skyrmions in thin films has been recently shown [34,46]. We will use a classical expression [47] for the surface demagnetising energy which represents the demagnetising energy difference between a magnetic cylinder with uniform magnetisation pointing in one direction and the uniform ferromagnetic state with magnetisation pointing in the other direction: where and where K(u) and E(u) are the complete elliptic integrals of the first and second kind. This formula is a very good approximation to obtain the surface demagnetising energy when the topological soliton radius is much larger than the DW width r π∆. For r ∼ π∆ it leads to an overestimation of this energy as it assumes an abrupt variation of the surface charge density instead of a progressive variation along the DW. The impact of this overestimation is discussed in Section 4.3 and Section 4.4.

Total energy
We obtain the following analytical expression for the soliton energy with respect to the homogeneous state: where σ s is the topological soliton wall energy density defined in Section 2.3. When the radius decreases, it deviates from σ w , as shown in Fig. 1d, and the radius dependent correction coming from the curvature of the wall is becoming comparable to σ w itself.

Topological soliton solutions
Topological soliton solutions are minima in the soliton energy E s as a function of the topological soliton radius r. We call the equilibrium topological soliton radius r s . In order to study the conditions giving rise to these minima, we fix the parameters A ex. , M s , K u and t and vary D. We have checked that changes in the fixed parameters do not modify qualitatively the results presented here but rather shift the main features to different D values. The different energy terms (except the Zeeman energy) are plotted as a function of r (up to 1 µm) in Figure 2.a, using the same parameters as in Figure 1, and three different DMI values. The resulting topological soliton energies E s are shown in Figure 2b. As the different energy terms compensate, a small variation of D is enough to modify the slope of E (1) s , E (2) s and E s . In the following sections we will describe three type of topological solition solutions, observed for increasing D values.

Skyrmion solutions
As discussed in Section 2.3 and illustrated in Figure 1e, dE exch. /dr decreases with r in the r < π∆ range. This non linearity may lead to the formation of an energy minimum as observed in Figure 2c. For r < r s (indicated with an α) , the soliton energy variation is dominated by the dE DMI /dr variation while for r > r s (indicated with a β), the d(E exch. +E anis. )/dr variation is taking over. This soliton solution corresponds to what is usually referred to as a skyrmion: its radius is of the order of a few ∆ or less and it exists down to H = 0. The skyrmion radius increases with D and depends weakly on the applied magnetic field as we will see in Section 4.

Bubble solutions and coexistence of skyrmions and bubbles
The second situation occurs when the positive and negative energy terms nearly compensate over a wide r range (see E with r, before E zee. with a r 2 variation takes over for larger r (δ part). Consequently, in presence of magnetic field, two soliton solutions are observed, separated by a local maximum of energy. The lower radius solution is the skyrmion solution described in Section 3.1 and shown in Fig. 2c. The second solution presents the characteristics of what is usually named a magnetic bubble: it collapses when increasing the magnetic field (see Fig. 2g) and its size diverges at H = 0. This coexistence of a skyrmion and a bubble solution was evoked in a pioneering work on skyrmions [8] and described in theoretical works from Kiselev at al. [35] and Büttner et al. [46]. The local maximum of energy creates an energy barrier which is at the origin of hysteretic behaviours in the M (H) loops of magnetic bubble materials [19].

Solutions above a critical D cs value
The third situation occurs when the DMI reaches a critical D cs value above which the local maximum of energy, as observed in Fig. 2d, disappears (see Fig. 2e). In this case, the total energy E s presents a negative slope at all r in absence of applied magnetic field (Figure 2b). In presence of magnetic field, an energy minimum is restored in E (3) s as the positive Zeeman energy variation dominates at sufficiently large radius due to its r 2 variation (δ part in Fig.  2e). For increasing magnetic field, this solution can be compressed to very low radius, as it is the case for skyrmions, without encountering a collapse field. When the magnetic field is decreased, the topological soliton radius will increase and diverge at H = 0 as it is the case for bubbles. The critical D cs value above which these solutions appears will be further discussed in Section 4.

Topopogical solitons phase diagram
In Figure 3, we present the evolution of the topological soliton radius r s , calculated with the same A ex. , M s , t, and K u parameters as in Figure 1 and  (Figure 3f). The asymmetry of the topological soliton phase diagram with respect to the magnetic field comes from the fact that we are describing metastable states compared to the state uniformy magnetised in the positive z magnetic field direction.

Description of the topological soliton phase diagram
We have divided the topological soliton phase diagram in different zones appearing in Figure  3b and f and described in the following.

Single topological soliton: zones 1 and 2
The skyrmion solutions described in Section 3.1 appears in zone 1. These solutions persist down to H = 0 and for negative applied magnetic fields (see Figure 3.a) and disappear along a blue line visible in Fig. 3b which was defined as the skyrmion bursting line in previous works [48,35]. This suppression of the skyrmion solution for decreasing magnetic field is illustrated in Figure 2f. The skyrmion bursting line ends at a critical point (D cs ,H cs ) indicated by a blue dot in Fig. 3b at which no local maximum is observed any more in E s (r). In the works where the long range demagnetising effect is neglected [48], this critical point appears at H = 0 and D cw = 4 √ A ex. K eff. /π. In our case, the critical topological soliton D value D cs is lowered compared to D cw , because we take into account the long range contribution in the demagnetising energy, which stabilizes the topological soliton similarly to the DMI energy term. The skyrmion radius in zone 1 is always smaller than the radius at the critical point r cs ( ∼100 nm in our case). In addition the skyrmion radius shows a small susceptibility dr s /dH, except close to the critical skyrmion bursting line (see Fig. 3c). This is due to he fact that the Zeeman energy term is a second order contribution at low r due to its r 2 variation. Zone 2 contains the solutions described in Section 3.3. As observed in Fig. 3c the topological solitons in this zone present a skyrmionic behaviour at high positive magnetic field: low dr s /dH susceptibility and no collapse field. At low magnetic field their susceptibility is large and increases with decreasing field, as it is the case for bubbles. The dashed line in Figs. 3a,b and d,f is the line at which the isolated topological soliton energy E s (r s ) becomes negative indicating that the uniform ferromagnetic state is not any more the ground state. Above this line the low energy cost of topological solitons and DWs favours the formation of a topological soliton lattice or a stripe/helical phase as described in the works from Bogdanov et al. [48] and Kiselev et al. [35]. We point out that zone 2 is extending above the critical D cw value. In this zone, for high magnetic fields, the magnetic field compresses r s down to r ∼ π∆ and the non-linearity in the exchange term causes the soliton energy to increase and become positive again: the isolated soliton solution is restored despite a negative σ w (see Figure2e). The D values are ranging from D/D cw = 0.59 to 1.13 for (c) and from D/D cw = 0.87 to 0.93 for (g). The green line in (g) correspond to D/D cw = 0.896, close to D cs . All the programs used to create the data and the figures can be found in a repository [38].
Experimental observations of such metastable isolated topological soliton at high magnetic field can be found in the work from Romming et al. [17] who reported skyrmions in systems with D > D cw under high magnetic fields of a few Tesla.

No topological soliton: zone 0 and 3
Zone 1 starts for D values larger than D min. = 2 √ A ex. K eff. /π. Below this value, in zone 0, the DMI energy is too low to compensate the anisotropy energy, the α slope in Figure 2.c becomes positive, and there is no energy minimum in E s (r). This defines a lower DMI value for the formation of a static topological soliton, however, dynamical solitons [49,50], not discussed in the present work, can still be created for example using spin polarised currents [51,52]. In addition, higher order exchange energy terms, not taken into account in the present model can also lead to the stabilisation of skyrmions in absence of DMI [8,9,53,54]. Zone 3 is delimited by the blue skyrmion bursting line at which the skyrmion solution disappear and the H = 0 red line. In this zone the negative magnetic field suppresses the isolated metastable state at low D. For sufficiently large D, at low magnetic field, a topological soliton lattice or a stripe/helical phase is predicted in this zone [48,35,26]).

Bubbles solutions and skyrmion and bubble coexistence: zone 4
The region where bubbles appear is defined as zone 4. This zone overlaps with zone 1: a second energy minimum corresponding to the bubble solution appears in E(r) without modifying the skyrmion radius. In Figs. 3d and e we show the topological soliton solution at low magnetic field and for a narrower DMI range compared to Figure. 3a. When two solutions coexist, as discussed in Section 3.2, the larger/smaller solution is shown respectively in Figs. 3d and e. The coexistence zone is delimited by the H = 0 line, the bubble collapse red line and the skyrmion bursting blue line. The bubble collapse line and the skyrmion bursting line meet at the critical (D cs ,H cs ) point. The coexistence is also visible in the vertical cross Section of Figs. 3d and e shown in Fig. 3g. When D is decreased, the collapse field decreases and zone 4 vanishes because the topological soliton wall energy becomes to large to be compensated by the long range demagnetising energy.

Transitions close to the critical (D cs ,H cs ) point
In Figure 3g we show the r s evolution with magnetic field very close to D cs . Below D cs the energy barrier which separates the two solutions, visible in Fig.2d, causes hysteretic behaviors in the bubble-skyrmion transformations. A bubble transform into a skyrmion via collapse for increasing magnetic field and skyrmions abruptly expand when the magnetic field is increased in the negative direction up to the bursting field (see Fig. 2f and g and Fig. 3g). At D = D cs and H = H cs (blue dot in Fig. 3g), the hysteretic behaviour is suppressed (see green line in Fig. 3g). This is due to the suppression of the energy barriers separating the two solutions which leads to a remarkably flat E s (r) energy profile close to r cs as visible in Fig. 2h. This particularity, which is only observed close to the critical point, is due to an almost perfect compensation, when the radius is varying, of the topological soliton wall energy cost by the surface demagnetising energy gain. The skyrmion-bubble transition observed here is reminiscent of the critical phenomena observed in the liquid-gaz second order phase transition. Firstly, the transitions occur along lines that terminate at a critical point (D cs ,H cs ). Secondly, for D < D cs an interval where both solitons coexist is observed similarly to the gas and liquid mixture observed in the temperature versus density plane of the liquid/gas phase diagram. Thirdly, we observe numerically a divergence in the topological soliton compressibility dr s /dH at the critical point (Figure 3.g). In analogy to the opalescence phenomena observed in liquids at the critical point, topological soliton at this point should present remarkable behaviours due to strong thermally activated fluctuations of its radius. shifts the critical D value at which the compensation between positive and negative energy terms occurs. However, as discussed in Section 2.4.4 the expression we use to estimate E long range dem.
leads to an overestimation of this energy. Consequently, while the D cs calculated with the parameters used in the present work is equal to 0.88 · D cw , the critical D cs in a real system may be closer to D cw . This is confirmed by our micromagnetic simulations presented in Section 4.4 where the D cs value is found to be 8% larger leading to D cs = 0.95 · D cw .

Zeeman energy
The error associated with the Zeeman energy analytical expression that we use comes from a non compensation of the up and down M z component in the region where the spins rotate. When r π∆ this error is negligible (< 1%). For r = 2π∆ = 20 nm the error in E Zee is of the order of 5%. For r = 5 nm it reaches 30% of E Zee . However the impact of this underestimation on the topological soliton radius is limited by the fact that E Zee decreases quadratically with r and that the Zeeman energy and its variations are always one to several orders of magnitude smaller than the total energy value and total energy variations in the r < ∆ range.

Comparison with micromagnetic simulations
We have carried out micromagnetic simulations in order to confirm the predictions obtained using our analytical model and to estimate the impact of our approximations. The simulations have been performed using the Mumax3 open source software from Ghent University [55]. The skyrmion stability under an applied external magnetic field has been computed after discretising the system into orthorhombic cells (finite difference approach). Periodic boundary conditions has been used to limit size effects and the damping coefficient is set to 0.5 to speed up convergence (no magnetization dynamics). The system is a 2048 x 2048 nm 2 square box with a mesh size of 0.5 x 0.5 x 0.5 nm 3 , in which a nanometer sized Néel-type skyrmion is first relaxed. Once the skyrmion is initialized, a 5 mT perpendicular magnetic field is applied and reduced in steps of 0.5 mT. A series of simulations is run for variable D values, ranging from 3.60 to 3.90 mJ/m 2 . We present the relaxed topological soliton radius as function of applied magnetic field (Figure 4a) and for comparison the equilibrium topological soliton radius obtained analytically using the same parameters except the D values, ranging from 3.40 to 3.60 mJ/m 2 (Figure 4b). The results support qualitatively the predictions from the analytical model. In zone 1, the skyrmion solution shows a low dr s /dH susceptibility. In zone 2, the topological soliton is larger and this susceptibility is increased. In zone 4, the susceptibility varies strongly with decreasing magnetic field. Close to the critical D value D cs =3.8 mJ/m 2 for Fig. 4a and D cs =3.5 mJ/m 2 for Fig. 4b we observe a skyrmion burst when the magnetic field is decreased to the critical H cs value. Close to these D cs and H cs values we observe hystretic behaviours in the topological soliton radius variation versus magnetic field in the micromagnetic simulation (not shown here). This behaviour, which reveals a bi-stability, is very similar to both what was reported in a recent micromagnetic study very similar to the one carried out here [56] and to what is predicted by our analytical model (see Figure  3g). To finish we would like to emphasize that such micromagnetic results must be considered with care, especially for bubbles which have a diameter of a fraction of the simulation box. Indeed, when D becomes large (D > 3.86 mJ/m 2 typically) or when the applied field is small the skyrmion starts to develop a squareness [57]. In addition, the skyrmion starts to feel the edges of the simulation box and confinement effects cannot be neglected any more, explaining the fact that the measured radius is not diverging at very low field.

The topological soliton skyrmionic factor
In order to estimate the role of the long range surface demagnetising energy in the stabilization of a given topological soliton with radius r s , we introduce the skyrmionic factor S which represent the ratio between the topological soliton wall energy cost and long range demagnetising energy gain: The S values corresponding to solutions from Figure 5.a are shown in Figure 5.b. For r s π∆ the skyrmionic factor is large: S 1 as E long range dem.
is negligible in this range compared to the topological soliton wall energy cost which tends toward E 0 (Figure 1.e). On the contrary, large skyrmions with r s π∆ are only formed when the topological soliton wall energy cost is lower than the energy gain due to long range demagnetising effect which implies S < 1.
To check this correlation between the size and the S factor of a given topological soliton, we plot the soliton radius as a function of S in Figure 5.d. Bubble and skyrmion solutions from zone 4 and 1 are shown respectively in red and black while topological soliton above the critical point (zone 2) appear in yellow in Figure 5.d. Skyrmion and bubble solutions appear in two distinct r ranges, above and below r cs . However above the critical point, the radius of topological solitons as well as their S factor can be tuned continuously across r cs . Topological solitons with S close to 1 shows a large r s distribution while the r s of skyrmions and compact topological solitons at high magnetic field are strongly correlated with S. We conclude that the size criteria is relevant to distinguish between skyrmions and bubbles only below the critical (D cs ,H cs ) point.  Figure 3f. All the programs used to create the data and the figures can be found in a repository [38].

Topological solitons stability
The analytical topological soliton model derived here allows us to calculate the energy barrier protecting the solitons from collapse. This collapse energy represents the energy necessary for the topological soliton to annihilate via compression. This gives an estimation of the stability of a topological soliton, keeping in mind that our continuous model may become irrelevant at the atomic scale and that other annihilation mechanism with lower energy paths may exist, in particular in the presence of defects or edges [30]. We define the topological soliton collapse energy barrier as E c = E 0 − E min s where E min s is the local minimum in E s (r). The bubble collapse energy is defined as E c = E max s −E min s where E max s is the local maximum. In figure 5.c we plot the collapse energies corresponding to solutions from figure 5.a divided by E RT = k B T 293K . In figure 5.e we show the topological soliton equilibrium radius r s as a function of their collapse barrier E c /E RT where the topological solitons from zone 2 are represented in yellow while skyrmion and bubbles solutions appear respectively in black and red. The segregation of skyrmion and bubbles in two different r s ranges, due to the presence of an energy barrier between them in E s (r), appears again clearly. In addition, bubbles have a strong dispersion in their collapse energies while the skyrmion stability is strongly correlated with its size with a power-law dependence. Topological solitons above the critical point in yellow show large collapse energy barriers and their energy at a given size is tunable. When the magnetic field increases, these solutions show the same power-law collapse energy versus radius dependence as skyrmions. Topological solitons and skyrmions of 20 nm and above present collapse barriers larger than 21k B T 293K , which corresponds to lifetimes longer than 1s (Arrhenius-Néel law with a try rate 1/τ 0 = 10 9 Hz) and their stability increases with size and can further be increased by parameter engineering (see also [46]).

Conclusion
We have developed an analytical topological soliton model containing expressions of the long range demagnetising and exchange curvature energies, two key ingredients to stabilize bubbles and skyrmions in ferromagnetic thin films. This allowed us to study systematically topological soliton solutions over a wide range of parameters and explore quantitatively the possible transitions between small and large topological solitons. The observed skyrmion-bubble transition present similarities with the liquid-gas transition, in particular a critical point is present above which the transformation between both spin textures becomes continuous. While distinct characteristics of skyrmions and bubbles remain, their common nature as topological solitons is emphasised. Above the critical (D cs ,H cs ) point, the topological soliton can not be strictly named a skyrmion or a bubble, as it possesses some characteristics of both spin textures. This hybrid between a bubble and a skyrmion may be referred to as a supercritical skyrmion.
on the manuscript. Special thanks to André Thiaville, Stanislas Rohart, Olivier Fruchart for discussion on the demagnetising field. Thanks to Albert Fert, Vincent Cros and Henrik Rønnow for their support.