Quantum Skyrmion Lattices in Heisenberg Ferromagnets

Skyrmions are topological magnetic textures that can arise in non-centrosymmetric ferromagnetic materials. In most systems experimentally investigated to date, skyrmions emerge as classical objects. However, the discovery of skyrmions with nanometer length scales has sparked interest in their quantum properties. Here, we simulate the ground states of two-dimensional spin-$1/2$ Heisenberg lattices with Dzyaloshinskii-Moriya interactions and discover a broad region in the zero-temperature phase diagram which hosts quantum skyrmion lattices. We argue that the quantum skyrmion lattice phase can be detected experimentally in the magnetization profile via local magnetic polarization measurements as well as in the spin structure factor measurable via neutron scattering experiments. Finally, we explore the resulting quantum skyrmion state, analyze its real-space polarization profile and show that it is a non-classical state featuring entanglement between quasiparticle and environment mainly localized near the boundary spins of the skyrmion.


I. INTRODUCTION
Magnetic skyrmions are vortex-like quasiparticles characterized by a nontrivial topological invariant in real space [1][2][3][4]. These states are typically found in noncentrosymmetric ferromagnets in a certain range of external magnetic field and temperature, and are stabilized by an antisymmetric spin exchange energy, termed Dzyaloshinskii-Moriya interaction (DMI) [1]. After their first detection in a magnetic system by a neutron diffraction experiment in 2009 [5], and a full microscopic tomography by electron microscopy in 2010 [6], intense follow-up studies revealed intriguing dynamical properties, rendering skyrmions potentially useful for memory and computing devices [7][8][9][10][11]. Usually, the skyrmions encountered in these systems arise from thermal fluctuations and extend over length scales that are much larger than the interatomic distance and thus behave like classical objects. Other possibilities to create skyrmions are through suitable DC current devices, such as those proposed in Refs. [12,13]. However, smaller skyrmions do exist [14] and have already created interest in possible quantum properties of skyrmions. Several works have predicted the quantum behavior of skyrmions by using classical magnetic textures as a starting point and studying quantum corrections in the semiclassical regime [15][16][17][18][19].
Beyond this semiclassical limit, some works have indicated that quantum analogs of classical skyrmions might exist in spin systems. The authors of Ref. [20] used a multiscale approach to demonstrate that mesoscopic magnetization vortices are stabilized by quantum fluctuations [18,21], which suggests the possibility of inherently quantum-mechanical counterparts of these states at zero temperature. So far, attempts to classify skyrmion excitations with sizes comparable to the interatomic spac- * andreas.haller@uni.lu † thomas.schmidt@uni.lu ing have been made in frustrated quantum lattice systems [22] and ferromagnetic lattices with DMI [23,24]. Several geometries have been studied to understand the quantum analogs of classical skyrmions, and quantitative results have been obtained by numerical diagonalization of the Hamiltonian [22][23][24]. Since the dimension of the quantum Hamiltonian scales exponentially with the number of lattice sites, such exact diagonalization (ED) strategies are limited to small system sizes containing at most ≈ 30 spin-1/2 sites (without exploiting symmetries). Although DMI interactions are among the most popular to investigate the formation of classical skyrmion phases, their quantum analogs are analytically hard to handle and quantitative results beyond system sizes amenable for ED are still lacking. As an alternative route to quantum skyrmions that avoids DMI, frustrated spin lattice systems were studied in Ref. [22]. Using ED for small systems and analytical spin wave theory, the authors identified skyrmions with magnon bound states and developed a phenomenological theory based on a trial wave function.
Here, we use the density matrix renormalization group (DMRG) algorithm to explore ferromagnetic phases of quantum spin-1/2 Heisenberg models with DMI and uniaxial anisotropy. As our main result, we discover a zero-temperature quantum phase with a nontrivial magnetic spin texture that signals an emergent quantum skyrmion lattice. This phase was previously overlooked because it appears only beyond a critical system size which for realistic parameters is larger than the system sizes amenable to ED. We identify three ferromagnetic phases that can be directly observed and distinguished in the space-resolved magnetization profile. Furthermore, we argue that the polarization gives access to the zerotemperature phase diagram of the model under investigation. Contrary to similar quasiparticles found in frustrated lattices [22] or quantum skyrmions embedded in a classical magnet [24], we show that the skyrmion lattice phase reported in this letter emerges from entangled spin-1/2 pairs, which bear witness to a genuine quantum mechanical origin without classical analog, a feature that may pave the way towards a microscopic description of skyrmion qubits used for realizing quantum logic elements based on nanoscale devices [25] II. MODEL We study the zero-temperature phase diagram of a quantum spin-1/2 Heisenberg model with DMI and external magnetic field. The Hamiltonian reads: whereŜ r = σ r /2, with Pauli matricesσ α,r for α ∈ {x, y, z}, denotes a spin-1/2 operator at position r. J < 0 is the ferromagnetic exchange coupling strength, D r −r is the DMI vector, and B = Bê z denotes the applied magnetic field along the z axis. The notation r, r implies a sum over all pairs of nearest-neighbor lattice sites.
In Section VI, we will also consider the impact of a uniaxial magnetic anisotropy with strength K, We solve the above Hamiltonian numerically by means of matrix product state (MPS) simulations on different twodimensional Bravais lattices consisting of lattice sites r = i n i a i spanned by basis vectors a 1,2 with a z,i = 0. For details on the implementation and on the 2D-1D mapping required for using MPS, we refer to the Appendix and Ref. [26]. The DMI vectors read with the positive DMI vector amplitude D > 0. For triangular and square lattices, we depicted the orientation of the DMI vectors in Fig. 1. Note that we apply the external magnetic field B parallel to the lattice plane normalê z . Without loss of generality, we assume negative values B < 0 such that field-polarized spins are eigenstates ofŜ z r with eigenvalue +1/2 and therefore align with the plane normal. For convenience, we will express interatomic distances in units of the lattice constants a i = |a i | and energies in units of D, with = 1.
The Hamiltonian (1) may be seen as the quantum counterpart of typical classical spin models which give rise to magnetic skyrmions [1][2][3][4]. In this work, we discuss the emergence of quantum skyrmions and quantum skyrmion lattices in two-dimensional triangular and square lattices at zero temperature and with different boundary shapes (see Fig. 2).

III. INDIVIDUAL QUANTUM SKYRMIONS
For small system diameters (L 5a), the exact eigenvalues and eigenvectors of the full Hamiltonian can be computed numerically, for instance by the Lanczos or Arnoldi algorithms, and we first cross-checked our own exact diagonalization (ED) codes by reproducing the results of Ref. [23]. To ensure the correctness of our findings, we further compared expectation values computed from the matrix product states (MPS) obtained by DMRG with ED results for all system sizes amenable for ED and found quantitative agreement up to a selfimposed accuracy δ. This shows that finite-size ground states of the spin-1/2 Hamiltonian can be approximated faithfully with MPS (for more details, see Section VII and [26]). In fact, the reliability of DMRG for 2D spin-1/2 quantum Heisenberg models has already been demonstrated in numerous works, with a strong bias towards frustrated antiferromagnets and spin liquids, which are prime examples of the most demanding systems to simulate numerically due to the presence of topological order and long-range entanglement . In this work, we focus on the computationally less demanding scenario of magnetically ordered phases hosting skyrmions and fieldpolarized (FP) states in the regime J < 0. Beyond a critical system diameter of L ≈ 8a, we find values of D and B for which the ground state of the HamiltonianĤ hosts skyrmion-like spin textures confined in the interior (bulk) of the lattice, which we display in Fig. 2. We checked that these skyrmion wave functions correspond to approximate eigenstates of the Hamiltonian by computing the energy variance ε and performed a linear extrapolation towards results without numerical errors (see Section VII for a discussion).
We compute the components of the spin magnetization, which are local expectation values m r = Ŝ r . Since m z,r is parallel to the external field, we call this magnetization component the polarization. For B ≈ J = −D/2, the local spin profiles yield magnetization textures similar to those obtained for classical skyrmion configurations of the Néel (hedgehog) type: the central spin is polarized opposite to the applied magnetic field, and the spins wind radially from the center towards the periphery. In Fig. 2, we depict the polarization along the field (m z,r ) . We find quantum skyrmions for system sizes larger than a critical diameter L ≈ 8a, irrespective of the lattice symmetries and boundary conditions. using a color scale and the in-plane polarization (m x/y,r ) by arrows. More detailed radial and angular distributions of panel Fig. 2(e) are depicted in Fig. 3. In Fig. 4, we show that the average polarization m z = 1 N r m z,r , as a function of the external field, yields three disconnected regions uniquely associated with the three phases of Eq. (1).
We estimate the size of an individual quantum skyrmion as the number of lattice sites over which the polarization changes its orientation once and the components orthogonal to the external field vanish. The radius can then be read out from Figs. 2 and 3 and results in r 0 ≈ 3a for J = −0.5D and K = 0. These quantum skyrmion (SK) ground states occur not only for the fine- tuned parameters presented in Fig. 2 but in a wide range of intermediate values of the magnetic field. Furthermore, the emergence of individual skyrmions for small lattices is largely independent of the lattice geometry. While the effect of boundaries cannot be neglected for the small systems considered here, we verified that the size of an individual skyrmion is neither affected by the boundary conditions (see Fig. 2) nor the system diameter (see Fig. 8). As we show in Fig. 2, for small magnetic fields, the system's ground state is a helical spin spiral (HS) state. This is characterized by a degenerate ground state, in which each possible ground state features an oscillation of the magnetization along a symmetry axis of the lattice. The large ground state degeneracy makes the spin spiral phase notoriously difficult to simulate for tensor network states. In contrast, for large magnetic fields, the system reaches a ferromagnetic state, where all spins are polarized parallel to the external magnetic field. The bulk of the fieldpolarized ferromagnet is devoid of entanglement and can thus be most efficiently approximated by an MPS. For parameters that result in SK ground states, DMRG reliably converges within a few dozen sweeps and yields excellent MPS approximations with maximum truncation error ∆ρ ≈ 10 −6 , even for small bond dimensions M = 32 [26].
We want to stress that the fine-tuned regime J = −0.5D is not necessary to enter the skyrmion lattice phase. By increasing the ferromagnetic exchange coupling to a value J = −2D, we obtain qualitatively similar results, which we present in Fig. 4. Note, however, that the corresponding system diameter for J = −2D is dramatically increased compared to J = −0.5D. We attribute this to the fact that the skyrmion radius is controlled by the ratio J/D.  To probe whether the spins constituting a skyrmion are entangled, a genuine quantum feature of many-body systems, we compute the entanglement entropy of a suitable bipartition of the system. The latter is defined as S = −tr(ρ A ln ρ A ), which can be understood as the von Neumann entropy of the reduced density matrix ρ A = tr B (ρ AB ), obtained by splitting the set of lattice sites N into two disjoint sets A and B, and performing the partial trace over subsystem B. For our purposes, it is sufficient to fix A and B as two patches that are symmetric about the central site (see Section VII).
MPS are constructed to target states of small entropy by truncating the reduced density matrix to a dimension dim(ρ A ) ≤ M , which yields an upper bound for the entanglement entropy of an MPS state,S = ln M , a quantity on the order of 1 − 10 for typical simulations. The bulk of FP states can be approximated by product states of spin-1/2 particles aligned to the axis of the magnetic field, and for those states, one finds S = 0 up to small finite-size corrections. In contrast, for systems hosting a single quantum skyrmion, we obtain values for the entanglement entropy in the range 0.2 < S < 0.7, which demonstrates the presence of significant entanglement in the spin-1/2 quantum skyrmion and indicates that they cannot be expressed as a classical product state. Finally, in the HS phase, we find the strongest entanglement (0.7 < S). While we focus mainly on skyrmions with ferromagnetic exchange interaction (J < 0), we find strong signatures of skyrmions and skyrmion lattices for antiferromagnetic exchange (J > 0) as well. However, this phase differs from the ferromagnetic skyrmions and skyrmion lattices, most apparently in the entanglement, which is significantly larger. We postpone a detailed investigation of quantum skyrmions and skyrmion lattices for antiferromagnetic exchange couplings to future work.
The von Neumann entropy targets quantum correlations between a bipartition of the system but does not provide local information about the entanglement between individual spin-1/2 pairs. To access the spatial distribution of entanglement, it is therefore more convenient to calculate the concurrence C r1r2 , defined for two lattice sites at positions r 1 and r 2 . For a generic state it can be expressed through the root of the spectrum of the non-Hermitian matrix R r1r2 = ρ r1r2ρr1r2 , where ρ r1r2 = tr r / ∈{r1,r2} (ρ) is the reduced density matrix of the two sites r 1 and r 2 , andρ r1r2 = (σ y ⊗ σ y )ρ * r1r2 (σ y ⊗ σ y ) is a rotation of this reduced density matrix. The concurrence is constructed from the square roots λ i (ordered in decreasing order) of the eigenvalues of R, It is related to the entanglement of formation: for separable states C vanishes, and it increases monotonically towards the limit C = 1 for two maximally entangled spin-1/2's [58]. Using the concurrence, we obtain the spaceresolved entanglement distribution between spin pairs in the different phases and present its qualitative distribution by green links in Fig. 5. For a more quantitative analysis, we differentiate between concurrences of different ranges up to 3rd nearest neighbor in Table I.
For helical spin spiral (HS) states, we find the largest concurrences up to C HS = 0.16. Therefore, we use C HS as a measure of reference to quantify the entanglement of the remaining ordered states. For the field-polarized (FP) states (large magnetic field), we find almost vanishing entanglement in the bulk, indicating that the bulk spins are separable. Small nonzero values of C ≈ 0.35C HS occur at the boundary due to the finite system size and strong DMI. Finally, for the quantum skyrmion (SK) states (intermediate magnetic field), we find that the spins inside the skyrmion quasiparticle are only weakly entangled, but we find concurrences C ≈ 0.8C HS at the outer rim spins of the skyrmion, signaling significant entanglement of the quantum skyrmion with the field-polarized environment. Interestingly, both the helical state and quantum skyrmion show long-range concurrences between distant spin pairs beyond next-to-nearest neighbors.

V. STRUCTURE FACTOR AND NEUTRON SCATTERING CROSS SECTION
We compute the Fourier components of the spin-spin correlation function as follows, where α, β ∈ {x, y, z}, q = (q x , q y , q z ) a wave vector which is later associated with the scattering vector, and the expectation values of the product of spin operators are evaluated with the ground state obtained by the MPS simulations of small hexagon flakes (panels (d)-(f) of Fig. 2). From S αβ , the elastic magnetic neutron scattering cross section dσ/dΩ at momentum transfer vector q is given by [59]: whereq = q/q = (q x ,q y ,q z ). Concerning experiments, we emphasize that dσ/dΩ corresponds to a scattering geometry where the externally applied magnetic field B = Bê z is parallel to the wave vector of the incoming neutron beam, and where the detector plane is spanned by the two components q x and q y of the scattering vector. In the limit of the small-angle approximation, one can assume q z ≈ 0. In Fig. 7 we display dσ/dΩ for the helical, skyrmion, and field-polarized states (see Fig. 6 for all components S αβ of the structure factor). The geometry and local polarization of those states is displayed in Fig. 2(d)-(f). Generally speaking, long-range magnetic ordering is signaled by the presence of Bragg peaks in dσ/dΩ at momentum transfers q corresponding to the wave vectors of the ordering. We expect additional diffuse magnetic scattering components in dσ/dΩ rooted in spatial variations of the spin orientation.
As expected, the cross section of the FP state in Fig. 7(c) is isotropic and exhibits a single broad peak centered at q = 0 mainly caused by the component S zz parallel to the external field. The contributions in S αβ (α ∈ {x, y}, β ∈ {x, y, z}) are attributed to finite-size effects, caused by a helical winding of S r near the boundaries due to the strong DMI interaction (see Fig. 2). The azimuthal average of the cross section in Fig. 7(c), defined as (2π) −1 2π 0 dϕ(dσ/dΩ) can be well described by the form factor of a uniformly polarized thin circular disc with a radius R corresponding to the cluster radius, i.e., dσ/dΩ(q) ∝ [2J 1 (qR)/(qR)] 2 , where J 1 (z) denotes the first-order Bessel function. To highlight this point, we display the first two minima of [J 1 (qR)/(qR)] 2 (for R = 4a) by dotted white lines in Fig. 7(c) and find a very good agreement to the numerical data of the discrete system. The HS state in Fig. 7(a) is characterized by a superposition of spin spirals with wave vectors q = 0, resulting in six pronounced Bragg peaks in S zz (q). We observe in the SK phase a superposition of the two extreme limits: in particular, we find a Bragg peak at q = 0, together with an off-diagonal Bragg "ring" caused by the radial polarization winding of the skyrmion (see Fig. 7(b)). The radius of the skyrmion can be estimated as r 0 = q −1 0 ≈ 3a with q 0 ≈ 1/(3a) the momentum modulus of the q = 0 Bragg ring, consistent with the estimate given in Fig. 3. Hence, the predicted quantum skyrmion profile yields a distinct signature in the measurable neutron scattering cross section and allows a determination of its size. We want to stress that a classical spin profile compatible to the (normalized) quantum mechanical expectation values can lead to a cross section in qualitative agreement with the ones reported here. However, quantitative deviations are expected because, for classical systems, the spin-spin correlation functions factorize since the zero temperature configuration is non-degenerate for B = 0. Therefore, a direct measurement of the connected spin-spin correlation function σ α,rσβ,r − σ α,r σ β,r will differentiate classical from quantum skyrmion states. Classical and quantum states also differ locally in the spin norm, which is not necessarily conserved in general | Ŝ r | ≤ 1/2. We find that the domain wall spin norm at the outer rim of the skyrmion (where the concurrence is large) is about 4 % lower compared to the field-polarized environment.

VI. QUANTUM SKYRMION LATTICE PHASE
After having discussed the properties of individual quantum skyrmions in the preceding paragraphs, we now turn to the phase diagram of the system. Our numerical technique makes it possible to reach system sizes much larger than that of individual skyrmions, which in principle allows us to extrapolate towards a phase diagram in the thermodynamic limit. While the HS and the FP phase remain unchanged when increasing the system size, at intermediate magnetic fields, the ground state for large lattices features a regular lattice of quantum skyrmions (SKX). Similar to their classical analogs, the individual skyrmions form a dense packing, and for larger system sizes, we thus find quantum skyrmion chains and lattices, for which we plot examples in Fig. 8.
We elucidate the appearance and robustness of quantum skyrmion lattices as a function of the external magnetic field, the strength of DMI, and perturbations of the form Eq. (2) due to uniaxial anisotropy. For this calculation, we have concentrated on a triangular lattice with disk boundary conditions of diameter L = 9 sites, for which we obtain a single centered skyrmion for B = J = −0.5D in the unperturbed case K = 0 (see Fig. 2(e)). We relax the fine-tuned parameter lines of Fig. 2 by variations of B, J, and consider nonzero uniaxial anisotropies for J < 0 by varying K. Based on our simulations, we predict the existence of three distinct quantum phases of our model: (i) a region hosting helical spin spiral states (HS) for weak field amplitudes, (ii) a valley for field strengths of the order of B ≈ −0.5D, which features a lattice formed by quantum skyrmions of radius r ≈ 3a (SKX 3a ), and (iii) a fieldpolarized phase (FP) where spins align parallel to the external field.
Helical spin spiral states are characterized by a van-ishing average polarization m z , whereas field-polarized states are maximally polarized (up to finite-size effects). As shown in Fig. 2, quantum skyrmions are located in a background of field-polarized spins, and as a consequence, the state in the skyrmion lattice phase will have a finite polarization smaller than a corresponding fieldpolarized state. We numerically confirm this intuitive picture and find disjoint intervals of average polarization m z uniquely linked to each phase (see Fig. 4(a)), which can be summarized in the zero-temperature phase diagram presented in Fig. 9. Note that Fig. 9 is obtained by simulations of a fixed flake system size and is therefore only qualitatively correct in the thermodynamic limit. In order to determine the exact position of quantum critical points or the nature of the quantum phase transition, a finite size extrapolation is necessary, a study which we leave for future work.
Our results about the dependence of the skyrmion phase on the uniaxial anisotropy are in qualitative agreement with corresponding classical systems, where it is known that a weak uniaxial anisotropy tends to stabilize skyrmion configurations at smaller magnetic fields [14,60]. Furthermore, we observe that the skyrmion radius increases with −J/D (compare Fig. 9(e) and (h)), such that quite large spin-1/2 systems might be needed to resolve even individual quantum skyrmions. Based on our results for J = −0.5D we conjecture that the skyrmion lattice phase should also exist for such cases where the individual quasiparticles have a larger radius, but due to limitations dictated by the numerical complexity (which we discuss in Section VII), other numerical techniques must be consulted to make quantitative predictions about the phase diagram for −J/D 0.5. Similar values of the exchange coupling are expected in present thin film experiments [61], thereby making our results of practical relevance.
Besides the ferromagnetic SKX 3a phase, which is the focus of this work, we find signatures of a quantum skyrmion lattice phase in the absence of exchange coupling and even for antiferromagnetic couplings J ≤ 0. However, because of the significantly enhanced entanglement (see Fig. 9(c)), the MPS ansatz for this phase requires an exponential scaling of the bond dimension with the system size, and therefore simulations of large clusters are out of reach for DMRG.

VII. SUMMARY
We have demonstrated that the ground state of the two-dimensional ferromagnetic spin-1/2 Heisenberg model in the presence of DMI hosts quantum skyrmions at intermediate magnetic fields B ≈ J = −D/2. The resulting magnetic textures are characterized by a central spin pointing opposite to the direction of the applied magnetic field and winds radially outwards towards the field-polarized environment, similar to a classical Néel skyrmion. For periodic boundary conditions and in the thermodynamic limit, we expect the ground state of the skyrmion lattice phase to be degenerate, scaling with the area of the individual skyrmion quasiparticles, such that the bulk of a system with open boundary conditions corresponds to a spontaneously symmetry broken state. The existence of quantum skyrmions yields experimental signatures in the position-dependent magnetization, the average polarization, and the structure factor, and we showed that these observables allow a distinction between a spin spiral phase at small magnetic fields, a skyrmion phase at intermediate magnetic fields, and a field-polarized phase at large magnetic fields.
While the spin texture is reminiscent of classical skyrmions, we should point out that in the present case, the skyrmion phase arises as a quantum ground state at zero temperature with open boundary conditions. In contrast, classical skyrmions typically occur at finite temperatures and result from a minimization of the free energy. Moreover, our examination of the resulting quantum state using the entanglement entropy and the concurrence has revealed that the quantum skyrmion state features significant entanglement shared between spin pairs of the skyrmion boundary. We argued that the quantum and classical states can be distinguished by the norm of the polarization | Ŝ r | ≤ 1/2 (conserved for classical states) and by connected correlation functions (vanishing for classical states). We therefore conclude that a semiclassical treatment of the quantum skyrmion based on a classical magnetic texture would not necessarily capture the internal degrees of freedom of a quantum skyrmion.
Towards larger system sizes, we found that the quantum skyrmion phase is characterized by a regular lattice of skyrmions. As the size of individual skyrmions is determined by the system parameters B, J, D, and K, a regular lattice requires commensurability between the lattice size and the skyrmion size. While our numerical simulations cannot reach the limit of infinite system size, our results allow us to extrapolate that the ground state in the thermodynamic limit features a dense packing of quantum skyrmion textures. Each of these quantum skyrmions has entanglement localized near its domain wall, but the entanglement between different skyrmions is small, which suggests that they can be approximated as individual quasiparticles.
We expect that our results may guide the development of an effective analytical field theory of the quantum skyrmion phase. Based on our experience, we conclude that variational tensor networks provide a suitable numerical technique to study these systems. This is not surprising for gapped quantum phases with short-range interactions and bounded entanglement. Nevertheless, using MPS for a two-dimensional system is not without pitfalls, as the necessary mapping on a one-dimensional system causes non-local interactions. We have made sure that our results have fully converged for lattice sizes corresponding to individual skyrmions. However, we have seen that the numerical errors grow for the system sizes required for 4 × 4 skyrmion lattices (29 × 29 spin-1/2 in total). For such large systems, we expect our results to be only qualitatively correct.
Regarding alternative numerical schemes, we have verified that our MPS results agree quantitatively with all available results from exact diagonalization. We have also compared our results to variational methods based on neural-network quantum states. However, we found significant deviations between the exact result and neural-network states even for small system sizes, and the error was already on the order of 10% for the energy eigenvalues. This suggests that neural-network states may not provide an efficient variational ansatz for mesoscopic spin systems with DMI. We expect that quantum Monte-Carlo simulations might be useful to go to larger system sizes. However, the inclusion of DMI together with an external magnetic field brings about a sign problem that hinders convergence. Other promising tensor network states for 2D spin systems with DMI are variational tree tensor network states [62] and projected entangled pair states (PEPS) [63]. We expect finite PEPS to outperform MPS for larger spin-1/2 systems hosting skyrmion lattices, especially since recently a more efficient gradient-based optimization has been developed based on automatic differentiation techniques [57,64,65].
Tensor networks provide an important numerical toolbox in computational physics and have been applied successfully to countless interacting and strongly correlated systems [66][67][68][69]. One of the most established algorithms, called DMRG [70][71][72], is understood as a sequential variational optimization of adjacent MPS tensors until convergence is reached. Our DMRG simulations are mainly based on the Julia package ITensors [73], and we made available a condensed version, reproducing Fig. 2(e), on GitHub [26]. Additionally, we cross-checked ITensor with TeNPy simulations [74], and upload the condensed version for the TeNPy framework alongside the julia implementation [26].
A generic state consisting of N spin-1/2 sites reads  The black line corresponds to a least squares fit, whose energy offset E0 and error is displayed in Table II. FIG. 11. The three distinct ground states of Eq. (1) (spin spiral, skyrmion, field-polarized) obtained by two different 2D→1D mapping strategies: zigzag (top row) vs. spiral (bottom row). The polarization is depicted by colors and arrows, and the von Neumann entropy is encoded by links in gray scale. Whereas the von Neumann entropy is rather small and homogeneously distributed for the zigzag order, it appears not only larger but also inhomogeneous for the spiral ordered 1D chain.
where {|i k } forms a canonical basis of the Hilbert space at site k out of N sites in total. For a finite system with Dirichlet boundary conditions, the objects A n ) is called the bond dimension. If M is fixed to an arbitrary integer, the MPS representation of quantum states can be used as a variational ansatz to approximate the minimum energy eigenstate. The quality of this approximation is controlled by M . This is particularly transparent in the so- . For generic many-body states rewritten as MPS, M is an extensive quantity in the number of sites and diverges in the thermodynamic limit. If the target state of a onedimensional system obeys an area law of the quantum entanglement, the von Neumann entropy is guaranteed to be a finite constant [75]. Consequentially, M remains finite in the thermodynamic limit, and MPS becomes exact, which explains the success of DMRG applied to onedimensional quantum systems. Despite its limitations in two dimensions, DMRG is frequently applied to ladder systems and can even yield reliable results for strongly correlated lattices, especially in the case of quantum spin-1/2 Heisenberg models. We expect MPS to reliably capture the physics of the quantum skyrmion lattice phase because the external field polarizes the environment, and the resulting states carry no entanglement in the paramagnetic regions, but localized entanglement around the domain wall of the skyrmion.
We typically start with random MPS initial states of bond dimensions up to M ≤ 1024, followed by sequential variational optimizations ("sweeping") of two adjacent tensors (two-site DMRG). The two-site DMRG allows us to estimate the truncation error ∆ρ =  which we use in Fig. 10 to extrapolate towards results without numerical errors. To ensure that we display converged results only, we carefully monitor local spin expectation values and stop the simulation if changes in the observables become smaller than δ = 10 −10 . Since we use DMRG in two spatial dimensions, convergence to a spin spiral state may require many sweeps, on the order of 100−1000. The quality of the approximate ground state with energy E(M ) = ψ(M )|Ĥ|ψ(M ) can be estimated by the energy variance By construction, ε = 0 for exact eigenstates of the Hamiltonian. Since MPS approximates the wave function with a finite bond dimension M , we have ε(M ) > 0 and lim M →∞ ε(M ) = 0 in general. Similarly, lim M →∞ E(M ) = E 0 converges to the true eigenstate energy. To estimate the numerical error of our approximate wave functions, we perform linear extrapolations of the energy scaled against the truncation error [28] and the energy variance [43]. The outcomes of this extrapolation are presented in Fig. 10 and Table II. We want to stress that the MPS approximations corresponding to skyrmion and field-polarized states easily reach convergence within a few dozen sweeps and follow the expected linear trend in the approximation errors ∆ρ and ε [26].
Changes in the local spin expectation values beyond M = 128 are invisible to the naked eye when displayed on the scales used in the main text such that a detailed error extrapolation is not needed. A word of caution is due in the case of helical spin spiral states: as we already explained in the main text, these states are difficult to simulate using MPS due to the large degeneracy of the ground state manifold. This leads to some issues in reaching convergence (up to 1000 sweeps are needed) which for too small bond dimensions may even cause DMRG to get stuck in local energy minima corresponding to excited eigenstates. In the helical phase and for large lattices, MPS is thus not always reliably converging to approximations of the global ground state but converges under some circumstances to low-lying excited states with less entanglement -with outcomes roughly comparable with those presented in Ref. [32].

APPENDIX B: MAPPING FROM 2D TO 1D
Before we can apply DMRG to the system at hand, the 2D lattice must be mapped to a 1D chain. The map from a 2D lattice to a 1D chain can be performed by a sequential numbering of the lattice nodes with major ordering along an arbitrary axis (zigzag order). We choose the major axis to be a 2 . This can be achieved by f (r(n 1 , n 2 )) = n 2 + n<n1 l(n), where l(n) is an auxiliary function that encodes the lattice open boundary conditions and r(n 1 , n 2 ) = i n i a i .
As a result, the lattice HamiltonianĤ = r,r Ĥ r,r + rĤ r is mapped to a chain Hamilto-nianĤ = r,r Ĥ f (r),f (r ) + rĤ f (r) . To simplify the remaining discussion, we now assume square or rhomboid boundary conditions (see Fig. 12), and 1 < n i < L i , such that l(n) = L 2 . On-site contributions remain local, nearest-neighbor interactions along the major axis remain short ranged, but the interactions along the a 1 axis now have an extended range |f (r)−f (r ±a 1 )| = L 2 . This results in a growth of the dimension of the Hamiltonian matrix product operator M H ∝ L 2 . To obtain reasonable computation times for large skyrmion lattice systems, one must therefore restrict the bond dimension M to significantly smaller values. For the largest quantum skyrmion lattice system we present in Fig. 8, we plot the converged results of M = 128. Note that the choice of our mapping preserves the locality of the interaction in one direction. In an attempt to remove this bias, we checked the resulting MPS quality for a different mapping, starting at the central spin-1/2 site and ordered radially outward (spiral ordering). Using the spiral ordering, we can confirm using the von Neumann entropy that the outer rim of the skyrmion is strongly entangled with its environment, a conclusion we had also reached based on the concurrence. Compared to the other proposition, the spiral mapping results in higher variational energy, likely the result of the inhomogeneous entanglement distribution (see Fig. 11). Since entanglement can be created by non-local transformations, it is known that certain mappings from 2D to 1D are beneficial compared to others, which can be utilized to obtain a substantial improvement of the overall simulation quality [76]. For the results presented in the main text, we consistently use the zigzag order. The phase diagram presented in Fig. 9 is entirely unaffected by this choice.

APPENDIX C: CLASSICAL VS. QUANTUM
A basic understanding of the classical low-energy configurations can be achieved by performing a variational minimization of the energy functional. In particular, we want to solve for the minimum energy spin configuration which satisfies in which E is the classical energy functional E = 1 2 r,r [JS r · S r + D r −r · (S r × S r )] + r B · S r (10) sharing the same notational conventions with the quantum Hamiltonian, except the use of classical spins. We implemented a standard variational optimization with the Optim julia package [77], which is included in our repository [26]. All variational techniques are prone to being trapped in local minima, which sensibly depends on the initial state. To be sure that for small 61-spin flakes we obtain samples of a global minimum energy configuration, we performed the variational optimization with 1000 different initial states where the azimuth and polar angles are sampled with a uniform distribution. The low energy results of the classical setup are obtained analogously to the quantum case: we analyze the lowest energy configurations as a response to a changing external Zeeman field, for which we fix the norm of the classical spins to 1/2. We present a condensed version of the classical low energy results in Fig. 13. We find that the spectral energy and magnetization lines of the different configurations as a function of B/|D| are continuous over a wide range of parameters in the phase diagram. In particular, the spectral lines show crossings, which we identify with a phase transition: the first excited states become ground states and vice versa, causing the sudden jumps in the magnetization. This simple analysis suggests first-oder Zeeman field induced phase transitions in this model, which are conjectured to be present in the quantum case as well.
If we compare the ranges of the skyrmion phase between quantum and classical, we note that the quantum skyrmions are ground states in the regions of the classical field polarized states, which is in agreement with the results presented in [18], namely that quantum fluctuations stabilize skyrmion textures.