Spin-orbit coupled Hubbard skyrmions

When the material phases exhibit topological quantum numbers, they host defects protected by the nontrivial topology. Magnetic skyrmions are such ``quantized"objects and although many of them are metals they had been most likely treated in theories as purely classical spin states. Here, we show that the electrons described by the Hubbard model with strong spin-orbit couplings can exhibit various nano- or flake-size skyrmions in their ground state. Importantly, the conducting electrons forming Fermi pockets themselves carry textured magnetic moments in both real and momentum space and on top of that, possess Chern numbers in their energy bands. This quantum and conducting skyrmion is related to small skyrmions observed in atomic-layered compounds. We clarify how the effective magnetic interactions and magnetic anisotropies are tied to the spin-orbit coupling, and how they influence the stability of skyrmions beyond the phenomenology.


I. INTRODUCTION
Skyrmions are topologically protected classical defects first discovered in a nonlinear field theory for mesons in high-energy physics 1 .These studies are followed afterward by the intent discussions in the 1990s, regarding the charged excitations of quantum Hall ferromagnets as skyrmions 2,3 .The latest analog appears in magnetic materials 4 particularly, in a family of noncentrosymmetric B20-type transition-metal silicides and germanides [5][6][7][8][9][10] , which have been experimentally established as platforms of skyrmions that extend from nanometers to micrometer scale.In these chiral ferromagnets with broken inversion symmetry, the classical descriptions based on field theories have successfully explained the phenomenological features of their swirling magnetic structures nowadays called Bloch-type 11 ; the underlying helical magnetic structures running in different directions are combined and thermodynamically converted to skyrmions at a very particular range of finite temperature and magnetic field strength.Later, another series called Néel-type skyrmions was found in the cycloidal magnets in non-centrosymmetric polar crystals [12][13][14][15] .Meanwhile, some other classes of skyrmions in centrosymmetric materials 16,17 show that the frustrated classical magnetic interactions 18 or RKKY interactions [19][20][21] can also bear very similar textures to Bloch-type ones.
These "classical skyrmions" are so far believed to be well understood by the classical-spin-based Ginzburg-Landau theories or models 11,22,23 , which seemingly suggests it be buried in oblivion that most of the skyrmionic materials discovered so far are metals.Importantly, for skyrmions found in laboratories, the electrons itinerate as charge carriers which themselves carry spatially textured magnetic moments.This means that although the ordered moments of long wavelength may fit the classical description, their underlying origins should be far from being purely classical.Previously, the magnetic interactions discussed in classical theories were ferromagnetic Heisenberg and Dzyaloshinskii-Moriya (DM) interactions as well as magnetic anisotropies which are treated as handweaving free parameters allowed from the crystal symmetry point of view; it remains unsettled how they stem from both spin-orbit coupling (SOC) and inversionsymmetry breaking.
The current challenge is to truly understand the SOC as the direct origin of skyrmions visibly formed by conducting electrons.It is found that the SOC Hubbard model can afford metallic skyrmions that carry both charges and spin moments while maintaining its topological feature as Chern numbers on top of their magnetic quantum winding numbers.Quite importantly, the observed skyrmions are different from the Bloch or Néel type sizable skyrmions previously reported at finite temperatures.They are either nano-or flake-size and appear in the ground state for a wide range of Coulomb interactions when the SOC is substantially large.It turns out that the DM interactions and the types of spin anisotropies for these insulating magnets are not at all independent parameters, being strictly governed by the SOC.Our conclusions urge the reconsideration of the models and conclusions drawn from classical theories.

II. MODEL AND METHOD
We consider the Hubbard model at half-filling with Rashba SOC on the triangular lattice whose Hamiltonian is given as where ⟨i, j⟩ runs over all pairs of nearest neighboring sites, c † i = (c † i↑ , c † i↓ ) is the creation operator of up and down spin electrons, n i is an electron number operator, and U and B are on-site interaction and an external magnetic field, respectively.The spin-dependent hopping integral λ originates from the SOC 24,25 , and when combined with the t-term we obtain the form, using the SU(2) gauge field, U ij = e i(θ/2)nij •σ , where σ = (σ x , σ y , σ z ) is the Pauli matrix, t eff = √ t 2 + λ 2 and θ = 2arctan(λ/t).The gauge field rotates the spin orientation by θ about n ij -axis when the electron hops to the nearest neighbor sites as shown in Fig. 1(a).
The direction of the rotation axis n ij is determined by the crystal symmetry of the materials.For example, in the B20-type chiral magnets with P 2 1 3 or P 4 1 32 space groups, n ij should be parallel to the bond directions and favor helical structures.Whereas, for polar magnets of C nv group, n ij should be perpendicular to bonds and host cycloidal magnetic structures.In the Hubbard model, the orientation of n ij turns out to be an important factor determining not only the kind of skyrmions but also its stability.In our calculation, we consider C 3v type n ij placed inside the 2D plane (see Fig. 1(a)) and B perpendicular to the plane unless otherwise noted.
The reason why Eq.( 1) remained unexplored despite the importance of studying electronic models for skyrmions is purely because of technical difficulty.Quantum Monte Carlo method (QMC) suffers a sign problem.Other standard numerical solvers that can be applied to Eq.( 1) use the finite-size cluster that often misfits with periods of incommensurate orders.In such cases, the stability of the true ground state is underestimated or overlooked without knowing its periods a priori.For example, dynamical mean-field theory (DMFT) is one of the most successful solvers for Hubbard models 26,27 and give a good description for Mott transition, whereas it is built on a quality of self-energy included in the frequencydependent Green's function which does not take account of the longer-range fluctuations or spatially long-distance correlation effects.When extended to larger-unit clusters known as c-DMFT 28 , the results show substantial dependence on the size and shape of the clusters, and are not suitable for the description of phases with possible magnetic structures or correlations of large spatial periods.In most cases, the benefit of using c-DMFT or other cluster-based methods overwhelms this disadvantage, but for the present Hubbard model with SOC, the artifact hinders the unbiased search of skyrmions.Recently, a theoretical treatment called sine-square-deformed meanfield theory (SSDMF) has been developed 29 and was successfully applied to the same Hamiltonian as Eq.( 1) on a square lattice 30 , whose phase diagram turned out to host a variety of magnetic structures including large-scale incommensurate spin-density-wave states, spiral, stripe, and vortex phases, which seriously compete with each other.There, the reliability of unbiasedness and quality of the SSDMF results are confirmed through the comparison with four different analytical and numerical methods some of which take full account of the correlation effects; For example, c-DMFT predicted nonmagnetic insulating phases that contradict with the QMC studies, whereas SSDMF agrees with QMC.The two-fold periodic commensurate magnetic order reported in the c-DMFT phase diagram was also an artifact; the density matrix embedding theory (DMET) 31 , which takes account of the full correlation and can capture the large spatial period orders 32,33 favor the SDW phase, in agreement with our SSDMF.
Indeed, SSDMF is not simply one of the lowest approximation method 29 ; The two utmost advantages beyond the simple mean field are, it can get rid of the finite size effect, and the ability to avoid the bias to wavevectors commensurate with the cluster size.In Ref. [29], it is even shown that the charge gap of a 1D and square lattice Hubbard model evaluated by SSDMF agrees with those of the exact solutions.This is due to the realspace renormalization effect caused by SSD 34 .We thus apply the SSDMF to Eq.( 1) for the cluster shown in Fig. 1(a) and successfully disclose the basic properties of the ground state.The SSD combined with the quantum many-body solvers have established that it can quantitatively accurately evaluate incommensurate orders or correlated density-wave states in several quantum spin models and Hubbard models 34,35 .Although SSDMF relies on the Hartree-Fock approximation, the predicting powers of SSDMF and SSD-DMRG about the most dominant charge and spin correlations are almost comparable, and the results agree almost perfectly, which is shown in Appendix A.
We briefly outline the SSDMF theory 29 .We prepare a finite-size hexagonal-like cluster and multiply an envelope function to Eq.( 1) that has a sine-squared functional form gradually decreasing from the maximum at the center (r i = 0) toward the cluster edges (see Fig. 1(a)) given as f SSD (r) = (1 + cos(π|r|/R))/2 with a radius R = max i |r i | + 1/2.The cluster size adopted are N = 90 in Fig. 1(a) and N = 198 to confirm the convergence of the finite size effect.
The site-dependent mean fields ⟨n i ⟩ and ⟨S i ⟩ are introduced and are determined self consistently.Our meanfield Hamiltonian with SSD is given by The obtained solutions are site-dependent, and when given the deformed Fourier transform as ⟨S q ⟩ = i f SSD (r i )⟨S i ⟩e iqri / i f SSD (r i ), they yield an unbiased period that almost exactly reproduces that of the bulk limit.This is in sharp contrast to the results obtained by the periodic boundary condition which show significant size dependence even when the cluster period is prepared as compatible with the periodicity of the expected orderings.Previously, for example, the period of (q, q) (q = (1 − 1/ √ 10)π) is found to be exactly detected using the 8 × 8 cluster 29 .We perform the inverse Fourier transformation on Sµ i (k) and ñi (k) with normal periodic boundary conditions of N = 192 clusters to obtain spin and particle densities on a uniform periodic lattice.

A. Phase diagram
We first show the phase diagrams on the plane of θ and B for U/t eff = 8 and 5 in Figs.1(b) and 1(c).The θ = 0 case is the standard Hubbard model at half-filling; when U/t eff = 8, the system is in a paramagnetic phase at B = 0, which is transformed at B/t eff ≥ 0.2 to the state with all the triangles having two-up and one-down spin orientations, which is known as UUD phase in quantum magnetism.This result is consistent with the corresponding previous theories on the half-filled Hubbard model 36 showing the spin liquid Mott insulator at U/t ∼ 8 − 10 (while it is beyond our description) and 120 • Néel ordering at U/t > ∼ 10 [37][38][39] .At θ = π with t = 0 and λ = 1, an almost fully polarized ferromagnet appears, whose origin will be discussed shortly.In the shaded area of the ferromagnetic phase, the metastable solution reminiscent of the Néel-skyrmion appears which is discussed in Appendix B.
At intermediate θ, t and λ compete and generate various phases.The most prominent feature is the emergence of two small-size skyrmion phases at around 0.2 < ∼ θ/π < ∼ 0.5, which we call nano-and flake-skyrmions based on the similarity with the previously found skyrmions in quantum spin models 40,41 ; their real space spin distributions and spin structural factors in the reciprocal space are shown in Figs.1(d) and 1(e).They have site-centered and bond-centered characters with 7-site and 27-site per unit, respectively.The spin structural factors have peaks positioned at q = (0.955π, 0) and (0.955π, 0.064π), respectively.In the lower field region, the AF-skyrmion based on three sublattice structures of a 27-site unit is observed; the structural factors show peaks at different points from the nano-and flake-skyrmions.All of them are characterized as triple-q states.For classical skyrmions, a winding number is used to characterize the slowly varying continuous spin structures, which takes either an integer or half-integer numbers when the gradients of spins are integrated over the skyrmion unit.In quantum systems, the spin moments fluctuate and squeeze on discrete lattice sites and do not show global continuous structures.The replacement of the winding number is given by the skyrmion number W introduced in a Ref. [41] as the total solid angle of triangular spins as where the summation is taken over the skyrmion unit for triangles ⟨ijk⟩ with indices aligned counterclockwise.The value of W is not quantized and is relatively small 42 .However, it can overall capture the nano-and flake-skyr phases as shown in the insets of Fig. 1.
Another notable feature is the 3sub-cluster phase realized at large B and θ ∼ 0.2 − 0.3π whose magnetic structures are shown in Fig. 2(a) and 2(b), corresponding to U/t eff = 8 and 5, respectively.At first sight, their underlying three-sublattice structures are not visible.However, if we separately plot them for each sublattice in the right panels, we find that each sublattice hosts a skyrmion with the 90-sites unit.Similar three-sublattice skyrmion was experimentally found in MnSc 2 S 4 43,44 , which are analyzed theoretically for the J 1 -J 2 -J 3 model 45 and DM model 46 .Interestingly, this skyrmion is a source of the thermal Hall effect based on SU(3) magnons characteristic of a three sublattice structure 47 .Our results point out the possibility of realizing a similar structure at the electronic level for the first time in the non-centrosymmetric Hubbard model.
For other phases at θ/π ∼ 0.5 with increasing the magnetic field, the system hosts stripes and 120 • -or umbrella-like three-sublattice phases both denoted as "3sub AF".The details of these magnetic structures are given in Appendix B.
Because SSDMF is based on the Hartree-Fock approximation, which generally over stabilizes magnetic orderings or metal-insulator transitions, one might suspect that a variety of phases that appear in the phase diagram may be wiped out when a strong fluctuation due to higher order correlation is seriously included 48 .However, in Ref. [30], for the same SOC Hubbard model with small U , it was shown that the SDW phase in the SS-DMF phase diagram agrees with the DMET result near the phase boundary, U ∼ 2, and further agrees with the prediction of random phase approximation (RPA).
To support the present phase diagram, we performed an RPA study for θ = 0.3π and compared it with our SSDMF solution at U = 5, where the magnetic structure may be relatively fragile.As we show in Appendix C, the RPA and SSDMF show good agreement about the magnetic structure factor.Because U = 5 is already sufficiently larger than U c of RPA, the magnetic orderings found in SSDMF can be safely expected.
It is an established knowledge that the standard mean-field calculation overestimates magnetic orderings.Alongside, SSDMF is confirmed to capture accurately the wave numbers of orders without bias.Combining these two facts, although the transition points or the amplitudes of magnetic orderings shall not be quantitative, the SSDMF will not miss the magnetic orderings if really present in the fully correlated case although not being possible to elucidate them by other methods.The states not captured by SSDMF are the quantum nonmagnetic ones like spin liquid or valence bond crystals which are not expressed by the standard Hartree-Fock approach to electrons.Such states may sometimes appear in frustrated systems like the present triangular lattice due to strong quantum fluctuation, while they are clearly out of the scope of the present study.

B. Metallic skyrmions
The phase diagrams have insulating and metallic regions whose boundaries are indicated by the broken lines in the small panels of Fig. 1(b,c).In our framework, the flake-skyrmion at U/t eff = 5 can naturally appear as metals.
In the metallic phases, the magnetic moments shrink, and it becomes much more difficult to have a rigid regular magnetic structure compared to the insulating case.In our calculation particularly at U/t eff = 5, the UUD, skyr-like, and stripe-like phases surrounding the flake-skyrmion phase often show various numerical solutions with domains or irregular winding structures of spins.Contrastingly, the flake-skyrmion metal is rigid (see Fig. 3(a)).
Let us briefly show how we classify the states to be metallic or insulating.Once we obtain the spatial periods of the spin structures by the dominant peak positions S γ (q), we can construct a system with periodic boundary conditions of the size of that period.Based on the corresponding set of mean fields that work as one-body potentials, the energy band structure and the Fermi surface are obtained, which are shown in Figs.3(b) and 3(c) for flake-skyrmion.The two bands cross the Fermi level which gives two hole pockets and electron pockets in different regions of the Brillouin zone.Since not only the valence and conduction bands as well as other bands do not contact with each other, we can calculate the Chern numbers as for the n th band using the fidelity of Bloch wave function ∆U k along the discretized plaquette of the reciprocal space 49 .The results are shown on the r.h.s. of all the energy bands.The conducting and valence bands have C n = 2 and −1, respectively.Since the present system is gapless, the Chern number does not directly lead to the quantization of transport coefficients, but the underlying Berry curvature can bear an anomalous Hall conductivity.This effect can be described by the U(1) gauge field that is generated from the spin texture of the skyrmion.For comparison, we also show the energy bands of the insulating AF-skyrmion and 3sub-cluster phases with their Chern numbers in Appendix B.
Due to SOC, the spins are no longer good quantum numbers.Therefore, the eigenstates carry spin moments pointing in various directions; we show in Fig. 3(d) the spin textures on valence and conduction bands of the flake-skyrmion.Since the spin moments point almost opposite on different sides of the pockets, one may expect the spin current to be observed if this state is realized.eff /U = 0.4.Panel (a) corresponds to the strong coupling limit of our SOC Hubbard model given in Eq.( 7) with in-plane D and n, (b) is the one obtained by taking K = 0, and (c) shows the case where nij in Eq.( 7) is pointing perpendicular to the xy-plane which gives D and n ⊥ xy-plane.In the lower panel of (a), the variation of J/A = cos θ, D/A = sin θ , K/A = (1 − cos θ), and the detailed descriptions of D and n for the three cases are summarized.The direction of vector nij for the DM and anisotropy terms are bond direction (i → j) dependent.

C. Spin models at strong coupling limit
In our calculation, the large Néel-type skyrmion previously found in experiments and classical spin models of the same symmetry class is not observed even for parameters far outside the presented phase diagrams.To understand why some skyrmions appear for certain parameters and others do not, establishing the connection of Hubbard models with classical spin models is useful, since the skyrmions in the latter types of models are intensively studied.
The effective Hamiltonian of the strong coupling region U/t eff ≫ 1 of Eq.( 1) at B = 0 can be obtained by the perturbation about t eff /U .The second-order Hamiltonian includes the Heisenberg interaction, DM interaction, and the Kaplan-Shekhtman-Aharony-Entin-Wohlman (KSAEW) term [50][51][52] as eff /U , and Importantly, the DM vector D ij and spin anisotropy vector both point toward the direction of the rotation axis of the SU(2) gauge, n ij , in Eq.( 1).
Here, S j in H (2) is a quantum spin-1/2 operator.Whereas taking its large-S limit and regarding Eq.( 7) as classical Hamiltonian of vector spins still works to understand the magnetic nature of the original Hamiltonian; it is shown in the square lattice SOC Hubbard model that the magnetic structure at finite U and agrees well with that of the classical ones 30 .To obtain the classical ground state phase diagram, we compare the classical energies Eq.( 7) of representative magnetic orders that were stabilized in the SSDMF calculation as a function of θ and B. ∼ 0.4π and 0.6π < ∼ θ are dominated by the threesublattice (UUD or 120 • ) phase and polarized ferromagnetic phase, respectively.Here, we prepared all representative states that appear in the Hubbard phase diagram as candidates and compared their energies; 120 • , UUD, nano-, flake-(19 site unit), and AF-skyrmions and largeskyrmion (Neél as well as other types) with an energetically optimized periodicity.For the nano/flake phases, the two compete while nano-skyrmion is relatively stable.The large Neél-skyrmion state was energetically subtle and did not appear as the lowest energy state.This diagram is fully consistent with the SSDMF results for Eq.( 1).
The variation of J, |D ij |, K as functions of θ shown in the lower panel of Fig. 4(a) helps to understand the energetics.At θ < 0.5π, the Heisenberg exchange is antiferromagnetic (J > 0) which becomes ferromagnetic (J < 0) at θ > 0.5π, having maximum amplitudes at θ = 0 and π.The DM interaction takes the maximum at θ = 0.5π, while the K term is always positive and monotonously increases; this term dislikes the adjacent spins to point in a similar direction, particularly in parallel to n ij .Since Fig. 4(a) has n ij pointing in-plane, the large DM interaction at θ ∼ 0.5π induces the typical cycloidal variation of moments along the bond directions, while comparably large K favors the spins to vary quite rapidly in space.Both cooperatively stabilize the nano/flake-skyrmion phase when B is not too large.
The large Néel skyrmion, which may appear if present at around θ/π = 0.6 − 0.9 is not favored here, because it suffers the substantial loss of KSAEW energy compared to the fully polarized state.The reason why the previous classical model could easily host large skyrmions was that the anisotropy interactions (different from KSAEW) and D ij are dealt with independently as free parameters, which is not allowed in Hubbard models.To confirm this scenario, we calculated the phase diagram by artificially setting K = 0 as shown in Fig. 4(b).At a small magnetic field next to the ferromagnetic phase the large Néel skyrmion is recovered.For further reference, we show in Appendix D the above-mentioned cases like D ∥ xy and n ij ⊥ xy which host Néel skyrmion, while they are unphysical in considering the SOC Hubbard model as their origin.
Nonetheless, since the large Néel skyrmion is relevant not only to classical models but also to materials, it is natural to expect another channel to realize it.We thus restart from the model Eq.( 1) with n ij pointing perpendicular to the two-dimensional xy-plane for all bonds.This modification does not change the form of Eq.( 7), and its classical phase diagram is shown in Fig. 4(c); the nano/flake-skyrmions no longer exist while the largeskyrmion appears at large θ region instead of the polarized ferromagnets.The large skyrmion here equivalently favors Néel, Bloch, or their mixtures because n ij • S i × S j depends solely on the relative angles of classical spins projected onto the xy plane, which is identical for the two skyrmions.This situation is realized in Mn 1.4 Pt 0.9 Pd 0.1 Sn where the mixture of Néel and Bloch structures is observed 53,54 .
At the fourth order of perturbation, the ring exchange interaction including SOC appears.We have examined the amplitude of these terms in Appendix E and found that they remain one order smaller than the other magnetic interactions.While these higher-order fluctuations will influence the amplitude of spins, their effect is small since the phase diagrams we saw in Figs.1(b) and 1(c) are similar, although the degree of fluctuation is by orders of magnitude larger for the latter.

IV. SUMMARY AND DISCUSSIONS
We have so far clarified that the metallic and nanoor flake-size skyrmions can be realized in an SOC Hubbard model at half-filling in a certain range of magnetic field perpendicular to the triangular lattice plane.Key differences from the former theories are listed as follows; Firstly, the conducting electrons forming Fermi pockets contribute to skyrmions.They may remind of a spin-density wave state but are different in that the spin momentum is not conserved on each band, namely, each k-point at each energy band carries magnetic moment pointing in different directions.Second, each band carries nonzero-integer Chern numbers originating from the U(1) gauge field created by the skyrmion structure which contributes to the Hall conductivity.Finally, the major effective magnetic interactions well-defined in the Mott insulating state, J, D ,and K, are interrelated and are fully controlled by the single parameter θ = 2 arctan(λ/t), originating from the SOC.Accordingly, the large skyrmions that were easily observed in idealized and phenomenological classical spin models are not found.
To make the discussions relevant to laboratory experiments, we present in Table I the representative noncentrosymmetric materials hosting skyrmions.They are classified into separate columns by double-lines as chiral magnets (P 2 1 3, P 4 1 32), atomic layers (Fe/Ir) or het- erostructures (SrIrO 3 /BaTiO 3 interfaces), and bulk layered materials.
The ones based on helical magnets (he) are known as B20-type structures with distorted cubic lattices.The Bloch-skyrmions found there are mostly as large as 10-200 nm, and are treated phenomenologically by the Ginzburg-Landau-based or lattice-based classical spin models; the major interactions are the ferromagnetic Heisenberg and DM interactions with D ij pointing parallel to the bonds along the x, ŷ, ẑ-directions that favor helical or conical magnets.These phenomenologies successfully reproduced the experimental phase diagram, at finite temperatures within a certain field strength 11,22,23 .
However, in reality, most of these magnets are metallic; for example, MnSi is a typical itinerant-electron magnet with a concentration of moments being 0.4µ B per Mn sites, which cannot be simply treated as a localized large-S magnet.Further, the microscopic origin of ferromagnetic interactions in such a situation is not clarified.Whereas, EuPtSi having the same space group is different; the Eu ions have well-localized S = 7/2 that interacts with conducting s-electrons 63,64 , which can be regarded as the RKKY systems.The theories built on classical Kondo-lattice types of models naturally reconcile with this system and exhibit phase diagrams with Bloch skyrmions 21,72 which are more or less the same as the ones commonly observed in B20-magnets.Recently, the DFT calculations have shown that the mixed valence state of MnSi can be separated into a relatively localized spin part and the conducting part 73 , where they propose the Kondo-like model and show that the Kondocoupling strength determines the magnetic concentration and the types of skyrmions realized.It may compromise the phenomenology and the RKKY model for this primary but specific class of skyrmions.There are a series of studies based on the spin-orbit coupled Kondo-lattice model [74][75][76][77][78] , where the skyrmions are induced by the DM interactions in the double-exchange limit.They find a Neél-type skyrmion at t/λ ≈ 0.55/0.45 76, which corresponds to θ ∼ 0.4π of our model.
The layered cycloidal magnets (cy) hosting Néel skyrmions have a different character.The major DM interactions allowed for these space groups have D ij perpendicular to bonds and are mostly pointing in-plane, which corresponds to the types of SOC we applied.Indeed, the SOC of these materials is substantially large 79 ; the atomic SOC for 3d (Fe, Co), 4d (Ru), and 5d (Ir, Pt) ions are 0.04-0.08eV, 0.1-0.2eV, and 0.3-0.6 eV, respectively, and when combined with the crystal field or the electric potentials coming from strong two-dimensionality, the electronic states on the clusters or ions have strong anisotropy and mixed spin-angular momentum.The electronic Hamiltonian in the presence of such effect is represented at the simplest by our Eq.( 1) 24,80 .The value of λ/t or θ is determined by the strength of SOC and the details of the crystal structures and can take substantially large values 25 .
We thus consider that these materials can be the platform of the skyrmions observed in our model.Their skyrmions overall have a relatively smaller size than those of the B20-chiral magnets, which agrees with our results.In particular, the flake-size skyrmions in our model are the only ones that can explain the metallic electron skyrmions observed in C nv crystals.Here, the Mott insulating phases can be partially related to the previously studied spin models that include a certain amount of quantum fluctuation; the nano-like skyrmion is found in a semi-classical spin-S model 42 that aimed to explain the 1 nm sized skyrmions found in Fe/Ir layers 65 .The flakeskyrmion can be compared with those observed in the fully quantum spin-1/2 system 41 .However, their magnetic interactions do not conform to the ones naturally obtained from Hubbard models (see Appendix D).
We stress that although previous theories assume the ferromagnetic Heisenberg interactions, it is elusive for electronic models without SOC.The superexchange interactions are so far the only natural realization for ferromagnetic exchange, which shall be the case applied to GaV 4 X 8 or VOSe 2 O 5 where the V-ions on V 4 X 4 or VO 5 clusters carry spin-1/2 and interact through ligands 12,13 .
In the Rashba-type SOC Hubbard model on the square lattice before this study, the double-q spiral phase with noncoplanar spin structures is found at 0.45 ≤ θ/π ≤ 0.65 and U/t eff > ∼ 4 30 , which is identified as an antiferromagnetic skyrmion lattice phase, and shall be relevant to VOSe 2 O 5 15 and SrXO 3 /BaTiO 3 hetrostructures 69,70 .The future perspective will be to apply a similar analysis to Hubbard models or other correlated electron models with different fillings or lattices, which may disclose the landscape of antisymmetric SOC Mott insulators and related SOC metals that host a variety of spin structures based on spin-split bands.

V. ACKNOWLEDGEMENT
The work is supported by JSPS KAKENHI Grants No. JP21K03440 from the Ministry of Education, Science, Sports and Culture of Japan and a Grant-in-Aid for Transformative Research Areas "The Natural Laws of Extreme Universe-A New Paradigm for Spacetime and Matter from Quantum Information" (KAKENHI Grant No. 21H05191) from JSPS of Japan.

Appendix A: Comparison of SSDMF and SSD-DMRG
We show that SSDMF can reliably capture a longrange correlation even though we implement a Hartree-Fock level approximation, which is conventionally considered a primitive approximation that cannot accommodate any higher order or subtle long-range correlations.
Figure 5 shows the comparison of DMRG and SSDMF for the one-dimensional Hubbard model with U = 2 and µ = 0.5, where for SSDMF, we reproduce the benchmark results shown in Fig. 1 of Ref. [29].This parameter gives the state away from half-filling with a charge density of n e ∼ 0.88 per site, where a charge density correlation of Q ∼ 0.24π appears together with a spin density correlation Q ∼ 0.88π.To test the quality of SSDMF, we apply a DMRG combined with SSD to the same state which takes account of the correlation effect almost exactly, and is considered the most reliable among all numerical solvers available.The on-site repulsion induces the misfit of the phase of the charge density of up and down electrons, and makes the total charge density ⟨n i ⟩ and spin density ⟨S z i ⟩ oscillate by π in DMRG; this is a typical Friedel oscillation found for DMRG calculations, an artifact known as a boundary effect 81 .If we connect every two lattice points by lines, we find clear periodicities that are intrinsic to this quantum state.
In SSD-DMRG, the amplitudes of this oscillation does not depend on the total charge number of the system because the SSD is not the method to decide the charge density of the total system but to decide µ and adopt the charge density at the center part of the system 34 .Indeed, the average charge density (center level of oscillation) and the period almost perfectly agree between three system sizes, L = 60, 100, 120.By performing a deformed Fourier transform, almost the same wave number, Q ∼ 0.24π − 0.245π, is captured for DMRG as those of SSDMF.
We note that for a usual periodic boundary, the onebody on-site quantities do not show any periodicity due to translational symmetry and the two-point correlation stores the information of the phase.However, in SSD systems the translational symmetry is broken, and the onsite quantities are modulated by the wave number that captures the dominant correlation of the state.
Appendix B: Several magnetic ground states in the phase diagrams We present in Fig. 6 the spin configuration of each phase in Figs.1(a) and 1(b) not referred to in detail in the main text.At U/t eff = 8 the 3 sub-sky and 3 sub-AF (120 • -like) magnetic states appear at large B whose moments are relatively large (Figs.2(a), 6(a)).
The UUD state at θ = 0 given in Fig. 6(b) is the typical phase observed in quantum Heisenberg or Ising models on a triangular lattice in the magnetic field, often forming plateaus.The stripe phase emerges when the system approaches the ferromagnetic phase at a small B which is shown in Fig. 6(c).
The states at U/t eff = 5 are mostly metallic, so the magnetic moments are much smaller than in the above cases.The flake-skyr phase (Fig. 6(d)) is very stable, which is the one discussed mainly in the main text.The stripy-domains (Fig. 6(e)) that have similarity with Néeltype skyrmion.
The one found in the smaller θ region (shaded) of the ferromagnetic phase in Fig. 1(b,c) in the main text has spin structure shown in Fig. 6(f), and is classified as a polarized ferromagnet, although there is a small but finite modulation of spin moments similar to Néel-type skyr.We have carefully examined the energy differences of this SSDMF solution with those of the fully polarized ferromagnet by choosing the proper cluster with periodic boundary condition and recalculating the solution for the uniform (non-SSD) Hamiltonian, finding that the fully polarized ferromagnet slightly overwhelms the Néel-like skyrmions.The parameter range where this kind of solution appears is shaded in the phase diagram.We expect that the Néel-like skyrmions may appear at smaller θ part of this region at finite temperature as in the previously reported theories.
The phase boundaries are obtained by examining how the energy ⟨H MF ⟩ varies with θ and B, whose examples are shown in Figs.6(g) and 6(h).There is a distinct change in the functional form of energy, which can be analyzed to determine the phase diagram together with the careful examination of the real space and k-space magnetic structures.
We next show the energy band structures of insulating AF-skyrmion and UUD 3 sub-like states in Fig. SreffS 1)) with U = 2 and µ = 0.5.We compare the DMRG combined with SSD using L = 60, 100, 120, and SSDMF using L = 60.In DMRG we performed more than 10 sweeps and kept the truncation error to less than the order of 10 −6 .We take the origin of the real space coordinate at the center of the one-dimensional chain of length L. The deformed Fourier transform of ⟨ni⟩ gives the peak at Q ∼ 0.24π for L = 60 and Q ∼ 0.245π for L = 100, 120 in DMRG, consistent with the SSDMF at Q ∼ 0.24π.Because of the Friedel oscillations and strong correlations, there is an antiferromagnetic correlation with a period of two lattice sites.We plot two lines that connect data points of even and odd indices, respectively, to clarify an incommensurate correlation relevant in this state.For the SSDMF in the lower panel of (a), the broken line is taken from Fig. 1(d) in Ref. [29], which is reproduced in solid lines we newly calculate.
ins.They are realized in the U/t eff = 8 phase diagram Fig. 1(b) at (B/t eff , θ/π) = (0.1, 0.5) and (0.45, 0.15), respectively.The energy dispersion is suppressed compared to Fig. 3(b) and a large charge gap of ∼ t eff opens.The Chern numbers are shown on the r.h.s. of the energy bands.Since the real-space magnetic textures (AFskyrmion is shown in Fig. 1(f) and 3 sub-like state similar to 3 sub-cluster in Fig. 6(a)) do not seem to differ too much from the metallic flake-skyrmion shown in Fig. 3(a) (main text), they have a nonzero winding number, while Chern numbers near the Fermi level are both zero for the latter and the valence band is 1 for the former.

Appendix C: Random phase approximations
To support the reliability of the spin structure found in the phase diagram of Fig. 1(c) we perform a random phase approximation (RPA) to elucidate the instability toward the ordered phases.We start from a magnetic  (C1) where µ, ν = x, y, z.Because the model has a SOC, different spin orientations are mixed in forming energy bands, and the corresponding χ 0 (q) is given in the 3 × 3 matrix form, whose elements are evaluated using Here, f (ε) = 1/(exp(ε/k B T ) + 1) is the Fermi distribution function, δ is an infinitesimal positive number, and ε m (k) and u m (k) are the eigenvalue and corresponding eigenvector (m = ±) of the Bloch Hamiltonian given by where e j denotes the primitive translation vectors for the triangular lattice, n j is the direction of the rotation axis for SOC on the edge corresponding to e j , and ϕ is defined by Imv = |Imv| t (cos ϕ, sin ϕ).The RPA susceptibility is given in the 3 × 3 form as Once the eigenvalues of χ 0 (q) are obtained as λ j (j = 1, 2, 3) , the one with the largest amplitude, λ 3 (q), that satisfies 1 − 2U max q λ 3 (q) = 0, will determine the phase transition point U = U c , and we find the possible ordering wave vector Q = argmax q λ 3 (q), at U > U c . Figure 8(a) shows the density plot of λ 3 (q) for θ = 0.3π and from its highest peak value, the critical point is estimated as U c = 4.19.We find very good agreement between the ordering wave vector obtained by RPA and that obtained by SSDMF for θ = 0.3π and B = 0, shown in Fig. 8 (b).Therefore, SSDMF safely captures the incommensurate wave number, which is difficult for conventional numerical methods based on the periodic boundary condition.
In general, the Hartree-Fock or RPA approaches may overestimate the stability of ordered phases and U c shall be increased when the correlation effects are taken into account.However, in the previous paper for the square lattice Hubbard model with SOC, DMET has supported the emergent magnetic order at U = 2, consistently with SSDMF 30 .This may resolve the concern that the incommensurate magnetic orders at the mean-field level may be wiped out by the fluctuation when we introduce higher-order correlation effects.Considering that the Mott transition in triangular lattice Hubbard model without SOC but with anisotropy occurs for most of the results with quantum many-body numerical solvers at around U c ∼ 4 − 6 82 , it is natural to expect that our phase diagram at U ≥ 5 accommodates magnetically ordered phases, which may not necessarily an insulator.

Appendix D: Comparison with several classical spin models
To understand how the types of magnetic interactions used in the previous studies influence the ground state phase diagram, we consider the analogues of second order Hamiltonian H (2) in Eq.( 7) in the main text, taking S i as classical vector of size-1.In the main text, the results for Eq.( 7) (H (2) ) with n placed in the xy-plane perpendicular to the ij-bonds is shown in Fig. 4(a) where the parameters J, D and K are the functions of θ.This result was compared with its analogs naturally derived by setting K = 0 in Fig. 4(b) or by taking n ∥ z-axis in Fig. 4(c).
However, the types of interactions or anisotropies adopted in previous studies are different; in Ref. [42] they took D ⊥ bond in the xy-plane as in our case but considered easy-axis anisotropy which is valid only for S > 1 as To see this effect we calculate in Fig. 9(a) the classical ground state phase diagram by setting K as constant, while varying D and J as functions of θ.The constant values, K = (4/U )(K m / D 2 m + 1) with K m = 0.518 and D m = 2.16, are set to directly compare with the case of Ref. [42].Their SkX1 phase has a similar structure to our nano/flake-skyrmion and is consistent with the stable skyr region in the phase diagram.
Regarding Ref. [83], they used D ⊥ bond in the xyplane but for spin anisotropies they set n ∥ z (corresponding to the XXZ-model, taking the spin anisotropic exchange interaction), which does not conform to the symmetries provided by the SOC in the Hubbard model.The phase diagram obtained by varying J, D, K in the same manner as our model but taking n ∥ z for the anisotropy term to follow Ref. [83] is shown in Fig. 9(b).The ones taking K = √ J 2 + D 2 /4 as constant about θ is also shown in Fig. 9(c).Broken lines are θ = 0.66π and 0.85π studied in Ref. [83] as J = −2D and J = −0.5D,respectively.If we interpret our nano/flake-skyr as their SKX 3a , the results of ours and theirs are slightly inconsistent in that the nano/flake-skyr appear only in J = −2D but not in J = −0.5D,while they argue that both cases have SKX 3a .In that context, our phase diagram is natural since the θ = 0.85π, where |J| > |D| holds, favors larger size skyrmion.
Although these works stress the quantum mechanical magnetic fluctuation as being key to realize small skyrmions, the systematic comparison given here suggests that the relationships between magnetic interactions governed by SOC seems to play a more important role.Appendix E: Effect of ring exchange in fourth-order perturbation In Fig. 4 in the main text, we examined H (2) including interactions originating from the second-order perturbation terms.Since the role of higher order terms is being discussed recently 19 , we additionally consider the fourthorder spin-dependent ring exchange term.By considering the SOC hopping term λ, the ring exchange term is given in the form,  (E2) This conventional ring exchange term is known to give a substantial influence on the magnetism particularly when U/t eff becomes small near the metal-insulator phase boundary.
At θ ̸ = 0, the coupling constant is given as J 4 (θ) = 80t 4 cos(2θ)/U 3 .To see how much contribution this term gives compared to the second-order terms, we show in Fig. 10 the contour plot of J 4 (θ)/ √ D 2 + J 2 as a function of θ and U where we set t eff = 1.At around U = 8, |J 4 |/ √ D 2 + J 2 < ∼ 0.2, showing that the effect of ring exchange interaction is small for the insulating phase.However, when U = 5 the contribution from the fourth-order terms is similar to the second-order terms, while since the system is in the metallic region at θ < ∼ 0.6, these magnetic interactions no longer make sense.We may interpret the enhancement of J 4 as having an itinerancy in the magnetism, which increases the quantum fluctuation effect.Then, based on this analysis, the comparison between the phase diagrams in Figs.1(b) and 1(c) (main text) can be made; although the magnetic phases become less distinct for the smaller U , the basic features of magnetism do not differ much.Therefore, we can conclude that the fourth-order terms may play only a secondary role for the stabilization of the skyrmion phases.We examine how the band structure of the flake skyrmion at U/t eff = 5 (see Fig. 3) evolves with the field.Here, we fix θ/π = 0.35 and plot the energy band and the Chern numbers for several of B in Fig. 11.As the band varies gradually, the Chern numbers change throughout the band, while the two bands at the Fermi level seem to have ±1.

FIG. 1 .
FIG. 1.(a) Schematic illustration of the rotation axis nij included in the SU(2) gauge field Uij when electrons hop from site -j to site-i.N = 90 finite size cluster with SSD boundary conditions mainly used in the present calculation and the envelope function fSSD(r) (see Model and Method).(b,c) Ground state phase diagram of Eq.(1) at U/t eff = 8 and 5.The two small panels give the density plots of the quantum winding number W and the metal-insulator phase boundaries for the two diagrams.(d,e,f) Realspace spin configuration of the SSDMF solutions for the nano-, flake-and tetrahedral AF-skyrmions.The arrows represent xy-component of each spin and the color and size of the circle represent z-component.The lower two panels are the spin structural factors ⟨Sq⟩ in momentum space for the in-plane xy(left) and z-components(right panel).

)FIG. 3 .
FIG. 3. Metallic flake-skyrmion ground state obtained by SSDMF at U/t eff = 5, B/t eff = 0.45 and θ = 0.35π.(a) Spin configuration inside the 27-site unit.(b) Energy band structure consisting of 54 bands with Chern numbers is shown on the r.h.s.The inset shows the magnified bands near the Fermi surface.(c) Fermi surfaces of valence (hole) and conduction (electron) bands.The bright/dark regions are those below/above the Fermi level.(d) Spin textures in momentum space (z-component as density plots and xy-components as arrows) carried by the valence and conduction bands.The Fermi pockets are given as the guide to the eye.

FIG. 4 .
FIG. 4. Phase diagram of the classical spin model including the Heiseberg interaction J = A cos θ, DM interaction D = A sin θ and spin anisotropy term K = √ J 2 + D 2 − J, taking t eff = 1 and U = 10 i.e., A = 4t 2eff /U = 0.4.Panel (a) corresponds to the strong coupling limit of our SOC Hubbard model given in Eq.(7) with in-plane D and n, (b) is the one obtained by taking K = 0, and (c) shows the case where nij in Eq.(7) is pointing perpendicular to the xy-plane which gives D and n ⊥ xy-plane.In the lower panel of (a), the variation of J/A = cos θ, D/A = sin θ , K/A = (1 − cos θ), and the detailed descriptions of D and n for the three cases are summarized.The direction of vector nij for the DM and anisotropy terms are bond direction (i → j) dependent.

Figure 4 (
a) shows the resultant (θ, B) ground state phase diagram, where we set the units of the interactions as constant, A = 4t 2 eff /U = 0.4 corresponding to t eff = 1 and U = 10.The nano/flake-skyr phase appears at around θ ∼ 0.5π and B < ∼ 0.8, and the regions θ < -

FIG. 5 .
FIG. 5. (a) Charge density ⟨ni⟩ and (b) spin density ⟨S zi ⟩ of a one-dimensional Hubbard model (θ = 0 of Eq.(1)) with U = 2 and µ = 0.5.We compare the DMRG combined with SSD using L = 60, 100, 120, and SSDMF using L = 60.In DMRG we performed more than 10 sweeps and kept the truncation error to less than the order of 10 −6 .We take the origin of the real space coordinate at the center of the one-dimensional chain of length L. The deformed Fourier transform of ⟨ni⟩ gives the peak at Q ∼ 0.24π for L = 60 and Q ∼ 0.245π for L = 100, 120 in DMRG, consistent with the SSDMF at Q ∼ 0.24π.Because of the Friedel oscillations and strong correlations, there is an antiferromagnetic correlation with a period of two lattice sites.We plot two lines that connect data points of even and odd indices, respectively, to clarify an incommensurate correlation relevant in this state.For the SSDMF in the lower panel of (a), the broken line is taken from Fig.1(d) in Ref.[29], which is reproduced in solid lines we newly calculate.

FIG. 7 .
FIG. 7. Energy bands of insulating states at U/t eff = 8 with (a) AF-skyrmions (B/t eff , θ/π) = (0.1, 0.5) shown in Fig. 1(f) in the main text and (b) UUD 3 sub-like state at (0.45, 0.15).Chern numbers of bands are shown on the r.h.s.In (a) it has a nonzero value at the valence band.

FIG. 8 .
FIG. 8. (a) Density plot of the maximum eigenvalue of the bare magnetic susceptibility λ3(q) as a function of q, obtained for θ = 0.3π.(b) SSDMF solution for θ = 0.3π and U = 5.Left panel is the density plot of the amplitude of spin moment |⟨Sq⟩| in the Brillouin zone, and the right panel shows the corresponding spin structure in real space.

FIG. 9 .
FIG. 9. Ground state phase diagram of the classical spin model with different types of anisotropies.The interactions J = cos θ, D = sin θ are set as the same as Fig. 4(main text), with B given in unit of A = 4t 2 eff /U = 0.4.For panel (a) we take −K(S z i ) 2 with constant K given as (4/U )(Km/ √ D 2 m + 1) with Km ≡ K/J = 0.518 and Dm = D/J = 2.16 following Ref.[42].In panels (b) and (c) we considered the XXZ type anisotropies by setting n ∥ z and taking K = 1−cos θ for (b) and K = √ J 2 + D 2 /4 being constant for (c), to compare with Ref.[83] whose J = −2D and J = −0.5Dcases are shown in broken lines θ = 0.66π and 0.85π, respectively.

4 j=0
Rj [S a , S b , S c , S d ] cos j θ 2 sin 4−j θ 2 , (E1)where Rj is a four-body operator about the four spins a, b, c, d on a unit plaquette consisting of two triangles, and its form depends on j.At θ = 0 this term reduces to the normal ring exchange term of the Hubbard model at half-filling84 , which takes the form,J 4 abcd (s a • s b )(s c • s d ) + (s a • s d )(s b • s c ) − (s a • s c )(s b • s d ) .

2 FIG. 10 .
FIG.10.Density (contour) plot of J4/ √ D 2 + J 2 on the plane of θ and U , where we set t eff = 1 and compared the ratio of the major fourth-order ring exchange interaction and the 2nd order magnetic interaction √ D 2 + J 2 = A = 4t 2 eff /U .
Appendix F: Magnetic-field dependence of ground state band structures and Chern numbers

TABLE I .
Details of non-centrosymmetric materials hosting skyrmions based on helical (he) and cycloidal (cy) magnets; space group, metallic/insulating (M/I), and the size of the observed skyrmions.Those belonging to B20-group (P 213 or P 4132) are helical/conical magnets that form Bloch skyrmions, and the corresponding DM interactions are Dij ∥ bonds.Here, EuPtSi is classified in a different column since it is an RKKY system where S = 7/2 localized magnetic moments are coupled to conduction electrons.The Cnv, R3m ,and P 4cc systems are cycloidal magnets that form Néel skyrmions, and the corresponding DM interactions are Dij ⊥ bonds which are pointing mainly in-plane.