Isolated Skyrmions in the $CP^2$ nonlinear $\sigma$-model with a Dzyaloshinskii-Moriya type interaction

We study two dimensional soliton solutions in the $CP^2$ nonlinear $\sigma$-model with a Dzyaloshinskii-Moriya type interaction. First, we derive such a model as a continuous limit of the $SU(3)$ tilted ferromagnetic Heisenberg model on a square lattice. Then, introducing an additional potential term to the derived Hamiltonian, we obtain exact soliton solutions for particular sets of parameters of the model. The vacuum of the exact solution can be interpreted as a spin nematic state. For a wider range of coupling constants, we construct numerical solutions, which possess the same type of asymptotic decay as the exact analytical solution, both decaying into a spin nematic state.


Introduction
In the 1960s, Skyrme introduced a (3+1)-dimensional O(4) nonlinear (NL) σ-model [1,2] which is now wellknown as a prototype of a classical field theory that supports topological solitons (See Ref. [3], for example). Historically, the Skyrme model has been proposed as a low-energy effective theory of atomic nuclei. In this framework, the topological charge of the field configuration is identified with the baryon number.
The Skyrme model, apart from being considered a good candidate for the low-energy QCD effective theory, has attracted much attention in various applications, ranging from string theory and cosmology to condensed matter physics. One of the most interesting developments here is related to a planar reduction of the NLσmodel, the so-called baby Skyrme model [4][5][6]. This (2+1)-dimensional simplified theory resembles the basic properties of the original Skyrme model in many aspects. The baby Skyrme model finds a number of physical realizations in different branches of modern physics. Originally, it was proposed as a modification of the Heisenberg model [4,5,7]. Then, it was pointed out that Skyrmion configurations naturally arise in condensed matter systems with intrinsic and induced chirality [8][9][10][11][12]. These baby Skyrmions, often referred to as magnetic Skyrmions, were experimentally observed in non-centrosymmetric or chiral magnets [13][14][15]. This discovery triggered extensive research on Skyrmions in magnetic materials. This direction is a rapidly growing area both theoretically and experimentally [16].
A typical stabilizing mechanism of magnetic skyrmions is the existence of Dzyaloshinskii-Moriya (DM) interaction [17,18], which stems from the spin-orbit coupling. In fact, the magnetic Skyrmions in chiral magnets can be well described by the continuum effective Hamiltonian where m(r) = (m 1 , m 2 , m 3 ) is a three component unit magnetization vector which corresponds to the spin expectation value at position r. The first term in Eq. (1.1) is the continuum limit of the Heisenberg exchange interaction, i.e., the kinetic term of the O(3) NLσ-model, which is often referred to as the Dirichlet term. The second term there is the DM interaction term, the third one is the Zeeman coupling with an external magnetic field B, and the last, symmetry breaking term A |m| 2 + (m 3 ) 2 represents the uniaxial anisotropy.
It is remarkable that in the limiting case A = κ 2 /2J, B = 0, the Hamiltonian (1.1) can be written as the static version of the SU(2) gauged O(3) NLσ-model [19,20] with a background gauge field A 1 = (−κ/J, 0, 0), A 2 = (0, −κ/J, 0). Though the DM term is usually introduced phenomenologically, a mathematical derivation of the Hamiltonian (1.2) with arbitrary A k has been developed recently [19], i.e.; it has been shown that the Hamiltonian can be derived mathematically in a continuum limit of the tilted (quantum) Heisenberg model where the sum ij is taken over the nearest-neighbor sites, S a i denotes the a-th component of spin operators at site i and W i ∈ SU (2). It was reported that the tilting Heisenberg model can be derived from a Hubbard model at half-filling in the presence of spin-orbit coupling [21]. Therefore, the background field A k can still be interpreted as an effect of the spin-orbit coupling.
There are two advantages of utilizing the expression (1.2) for the theoretical study of baby Skyrmions in the presence of the so-called Lifshitz invariant, an interaction term which is linear in a derivative of an order parameter [22,23], like the DM term. The first advantage of the form Eq. (1.2) is that one can study a NLσ-model with various form of Lifshitz invariants which are mathematically derived by choice of the background field A k , although Lifshitz invariants have, in general, a phenomenological origin corresponding to the crystallographic handedness of a given sample. The second advantage of the model (1.2) is that it allows us to employ several analytical techniques developed for the gauged NLσ-model. It has been recently reported in Ref. [20] that the Hamiltonian (1.2) with a specific choice of the potential term exactly satisfies the Bogomol'nyi bound, and the corresponding Bogomol'nyi-Prasad-Sommerfield (BPS) equations have exact closed-form solutions [20,24,25].
Geometrically, the planar Skyrmions are very nicely described in terms of the CP 1 complex field on the compactified domain space S 2 [6]. Further, there are various generalizations of this model; for example, two-dimensional CP 2 Skyrmions have been studied in the pure CP 2 NLσ-model [26][27][28] and in the Faddeev-Skyrme type model [29,30].
Remarkably, the two dimensional CP 2 NLσ-model can be obtained as a continuum limit of the SU(3) ferromagnetic (FM) Heisenberg model [31,32] on a square lattice defined by the Hamiltonian where J is a positive constant, and T m i (m = 1, ..., 8) stand for the SU(3) spin operators of the fundamental representation at site i satisfying the commutation relation Here, the structure constants are given by f lmn = − i 2 Tr (λ l [λ m , λ n ]), where λ m are the usual Gell-Mann matrices.
The SU(3) FM Heisenberg model may play an important role in diverse physical systems ranging from string theory [33] to condensed matter, or quantum optical three-level systems [34]. It can be derived from a spin-1 bilinear-biquadratic model with a specific choice of coupling constants, so-called FM SU(3) point, see, e.g., Ref. [35]. The SU(3) spin operators can be defined in terms of the SU(2) spin operators S a (a = 1, 2, 3) as Using the SU(2) commutation relation S a i , S b i = iε abc S c i where ε abc denotes the anti-symmetric tensor, one can check that the operators (1.6) satisfy the SU(3) commutation relation (1.5).
In the present paper, we study baby Skyrmion solutions of an extended CP 2 NLσ-model composed of the CP 2 Dirichlet term, a DM type interaction term, i.e., the Lifshitz invariant, and a potential term. The Lifshitz invariant, instead of being introduced ad hoc in the continuum Hamiltonian, can be derived in a mathematically well-defined way via consideration of a continuum limit of the SU(3) tilted Heisenberg model. Below we will implement this approach in our derivation of the Lifshitz invariant. In the extended CP 2 NLσ-model, we derive exact soliton solutions for specific combinations of coupling constants called the BPS point and solvable line. For a broader range of coupling constants, we construct solitons by solving the Euler-Lagrange equation numerically.
The organization of this paper is the following: In the next section, we derive an SU(3) gauged CP 2 NLσmodel from the SU(3) tilted Heisenberg model. Similar to the SU(2) case described as Eq. (1.2), the term linear in a background field can be viewed as a Lifshitz invariant term. In Sec. 3, we study exact Skyrmionic solutions of the SU(3) gauged CP 2 NLσ-model in the presence of a potential term for the BPS point and solvable line using the BPS arguments. The numerical construction of baby Skyrmion solutions off the solvable line is given in Sec. 4. Our conclusions are given in Sec. 5.
2 Gauged CP 2 NLσ-model from a spin system To find Lifshitz invariant terms relevant for the CP 2 NLσ-model, we begin to derive an SU(3) gauged CP 2 NLσ-model, a generalization of Eq. (1.2), from a spin system on a square lattice. By analogy with Eq. (1.2), the Lifshitz invariant, in that case, can be introduced as a term linear in a non-dynamical background gauge potential of the gauged CP 2 model.
Following the procedure to obtain a gauged NLσ-model from a spin system, as discussed in Ref. [19], we consider a generalization of the SU(3) Heisenberg model defined by the Hamiltonian whereÛ ij is a background field which can be recognized as a Wilson line operator along with the link from the point i to the point j, which is an element of the SU(3) group in the adjoint representation. As in the SU(2) case [19], the fieldÛ ij may describe effects originated from spin (nematic)-orbital coupling, complicated crystalline structure, and so on. This Hamiltonian can be viewed as the exchange interaction term for the tilted where T stands for the transposition. Let us now find the classical counterpart of the quantum Hamiltonian (2.1). It can be defined as an expectation value of Eq. (2.1) in a state possessing over-completeness, through a path integral representation of the partition function. In order to construct such a state for the spin-1 system, it is convenient to introduce the Cartesian basis where |m = |S = 1, m (m = 0, ±1). In terms of the Cartesian basis, an arbitrary spin-1 state at a site j can be expressed as a linear combination |Z j = Z a (r j ) |x a j where r j stands for the position of the site j, T is a complex vector of unit length [31,36]. Since the state |Z j satisfies an overcompleteness relation, one can obtain the classical Hamiltonian using the state Since Z is normalized and has the gauge degrees of freedom corresponding to the overall phase factor multiplication, it takes values in S 5 /S 1 ≈ CP 2 . In terms of the basis (2.2), the SU(3) spin operators can be defined as whereZ a denotes the complex conjugation of Z a . In the context of QCD, the field n m is usually termed a color (direction) field [37]. The color field satisfies the constraints where d mpq = 1 4 Tr (λ m {λ p , λ q }). Consequently, the number of degrees of freedom of the color field reduces to four. Note that, combining the constraints (2.6), one can get the Casimir identity d mpq n m n p n q = 8/9.
In terms of the color field, the classical Hamiltonian is given by Let us write the position of a site j next to a site i as r j = r i + aǫe k where e k is the unit vector in the k-th direction, ǫ = ±1, and a stands for the lattice constant. For a ≪ 1, the fieldÛ ij can be approximated by the exponential expansion where 1 is the unit matrix andl m are the generators of SU(3) in the adjoint representation, i.e., (l m ) pq = if mpq . In addition, since the model (2.1) is ferromagnetic, it is natural to assume that nearest-neighbor spins are oriented in the almost same direction, which allows us to use the Taylor expansion Replacing the sum over the lattice sites in Eq. (2.7) by the integral a −2 d 2 x, we obtain a continuum Hamiltonian, except for a constant term, of the form where A k = A m k λ m and n = n m λ m . Similar to its SU(2) counterpart expressed as Eq. (1.2), this Hamiltonian can also be written as the static energy of an SU(3) gauged CP 2 NLσ-model where g ∈ SU(3). Note that, however, since the Hamiltonian (2.11) does not include kinetic terms for the gauge field, like the Yang-Mills term, or the Chern-Simons term, the gauge potential is just a background field, not the dynamical one. We suppose that the gauge field is fixed beforehand by the structure of a sample and give the value by hand, like the SU(2) case. The gauge fixing allows us to recognize the second term in Eq. (2.10) as a Lifshitz invariant term. We would like to emphasize that we do not deal with Eq. (2.11) as a gauge theory. Rather, we deem it the CP 2 NLσ-model with a Lifshitz invariant, and show the existence of the exact and the numerical solutions. For the baby Skyrmion solutions we shall obtain, the color field n approaches to a constant value n ∞ at spatial infinity so that the physical space R 2 can be topologically compactified to S 2 . Therefore, they are characterized by the topological degree of the map n : R 2 ∼ S 2 → CP 2 given by Combining with the assumption that the gauge is fixed, it is reasonable to identify this quantity (2.13) with the topological charge in our model 1 .

Exact solutions of the SU (3) gauged CP 2 NLσ-model
In this section, we derive exact solutions of the model with the Hamiltonian (2.11) supplemented by a potential term. We first remark on the validity of the variational problem. As discussed in Refs. [20,25] for the SU(2) case, a surface term, which appears in the process of variation, cannot be ignored if the physical space is noncompact and the gauge potential A k does not vanish at the spatial infinity like the DM term. This problem can be cured by introducing an appropriate boundary term, like [20] where ρ = J/8. Here the gauge potential A k satisfies where n ∞ is the asymptotic value of n at spatial infinity. Note that Eq. (3.2) corresponds to the asymptotic form of the BPS equation, which we shall discuss in the next subsection. Hence, all field configurations we consider in this paper satisfy this equation automatically. Since (3.1) is a surface term, it does not contribute to the Euler-Lagrange equation, i.e., the classical Heisenberg equation. Note that the solutions derived in the following sections satisfy Derrick's scaling relation with the boundary term, which is obtained by keeping the background field A k intact under the scaling, i. e., E 1 + 2E 0 = 0 where E 1 denotes the energy contribution from the first derivative terms including the boundary term (3.1) and E 0 from no derivative terms.

BPS solutions
Recently, it has been proved that the SU(2) gauged CP 1 NLσ-model (1.2) possesses BPS solutions in the presence of a particular potential term [20,24]. Here, we show that BPS solutions also exist in the SU(3) gauged CP 2 model with a special choice of the potential term, which is given by As we shall see in the next subsection, the potential term can possess a natural physical interpretation for some background gauge field. It follows that the Hamiltonian we study here reads where the double-sign corresponds to that of Eq. (3.1). First, let us show that the lower energy bound of Eq. (3.4) is given by the topological charge (2.13). The first term in Eq. (3.4) can be written as It follows that the equality is satisfied if which reduces to Eq. (3.2) at the spatial infinity. Therefore, one obtains the lower bound of the form where the corresponding BPS equation is given by Eq. (3.6). Note that, unlike the energy bound of the CP N self-dual solutions [7,27], the energy bound (3.7) can be negative, and it is not proportional to the absolute value of the topological charge.
As is often the case in two-dimensional BPS equations [7,20], solutions can be best described in terms of the complex coordinates z ± = x 1 ± ix 2 . Further, we make use of the associated differential operator and background field defined as ∂ ± = 1 2 (∂ 1 ∓ i∂ 2 ) and A ± = 1 2 (A 1 ∓ iA 2 ). Then, the BPS equation ( where g ∈ SL(3, C). Note that Eq. (3.9) is not necessarily a pure gauge. Similarly, Eq. (3.8) with the minus sign on the right-hand side can be solved if For the background field (3.9), one finds that the BPS equation (3.8) is equivalent to because, under the SL(3, C) gauge transformation, the fields are changed as n →ñ = gng −1 and A + →Ã ± = gA + g −1 + ig∂ ± g −1 = 0. In the following, we only consider Eq. (3.9) to simplify our discussion.
In order to solve the equation (3.10), we introduce a tractable parameterization of the color field where Z is the continuum counter part of the vector Z in Eq. (2.3) and Y 1 , Y 2 are vectors forming an orthonormal basis for C 3 with Z. Up to the gauge degrees of freedom, the components Y i can be written as Therefore, the vector Z fully defines the color field n. Accordingly, we can writẽ C). It follows that the field Z, which is the fundamental field of the model, is given by Z = g −1 W 3 . Substituting the field (3.13) into the equation (3.10), one finds that Eq. (3.10) reduces to the coupled equation where W −1 l = Y † l g −1 (l = 1, 2). Since the three vectors {Y 1 , Y 2 , Z} form an orthonormal basis, Eq. (3.14) implies ∂ + W 3 = βW 3 where the function β is given by β = βW −1 3 W 3 = W −1 3 ∂ + W 3 . Therefore, the equation (3.10) is solved by any configuration satisfying where D + Φ = ∂ + Φ − (Φ −1 ∂ + Φ)Φ for arbitrary non-zero vector Φ. Moreover, we write where w is a three component unit vector, i.e. |w| 2 = w † w = 1. Then, Eq. (3.15) can be reduced to where P has no overall factor, and P a is a polynomial in z − . Therefore, we finally obtain the solution for the Z field where χ is a normalization factor.

Properties of the BPS solutions
As the BPS bound (3.7) indicates, the lowest energy solution among Eq. (3.19) with a given background function g possesses the highest topological charge. In terms of the explicit calculation of the topological charge, we discuss the conditions for the lowest energy solutions. The topological charge (2.13) can be written in terms of Z as We employ the constant background gauge field A + for simplicity. Then, the matrix g in Eq. (3.9) becomes so that the components of g −1 are given by power series in z + . It allows us to write Eq. (3.20) as a line integral along the circle at spatial infinity with C = −iZ † dZ [27,38], since the one-form C becomes globally well-defined. To evaluate the integral in Eq. (3.22), we write explicitly where g −1 ab is the (a, b) component of the inverse matrix g −1 . Let N a (K ab ) be the highest power in P a (g −1 ab ). Note that though g −1 ab are formally represented as power series in z + , the integers K ba are not always infinite; especially, if a positive integer power of A + is zero, all of K ba become finite because g −1 reduces to a polynomial of finite degree in z + . Using the plane polar coordinates {r, θ}, one can write g −1 ba (z + )P a (z − ) ∼ r Na+K ba exp[−i(N a − K ba )θ] at the spatial boundary and find that only the components of the highest power in r contribute to the integral (3.22). Since we are interested in constructing topological solitons, we consider the case when the physical space R 2 can be compactified to the sphere S 2 , i.e., the field Z takes some fixed value on the spatial boundary. Such a compactification is possible if there is only one pair {N a , K ba } giving the largest sum N a + K ba or any pairs {N a , K ba }, sharing the largest sum, have the same value of the difference. For such configurations, the topological charge is given by where the combination {N a , K ba } yields the largest sum among any pairs {N c , K dc }. This equation (3.24) indicates that the highest topological charge configuration is given by the choice N a = 0 for a particular value of a which gives the biggest K ba . We are looking for the lowest energy solutions with an explicit background field. As a particular example, let us consider where κ is a constant. Clearly, this choice yields the potential term which can be interpreted as an easy-axis anisotropy, or quadratic Zeeman term, which naturally appears in condensed matter physics. In this case, the solution (3.19) can be written as Therefore, the solution with the highest topological charge is given by P 1 = α 1 , P 2 = α 2 z − + α 3 with α i ∈ C, and P 3 being a nonzero constant. Choosing P 1 = P 2 = 0, one can obtain the axially-symmetric solution which possesses the topological charge Q = 2. Note that this configuration also satisfies the BPS equation of the pure CP 2 NLσ-model [26,27,31]. Figure 1 shows the distribution of the topological charge (3.20) of this solution (3.28) with κ = 1. We find that the topological charge density has a single peak, although higher charge topological solitons with axial symmetry are likely to possess a volcano structure, see e.g., Ref. [39]. These highest charge solutions give the asymptotic values at spatial infinity of the color field  It indicates that n takes the vacuum value in the Cartan subalgebra of SU(3). Hence, the vacuum of the model corresponds to a spin nematic, i.e., S 1 = S 2 = S 3 = 0 and (S 2 ) 2 = 0, (S 1 ) 2 = (S 3 ) 2 = 1. Unlike the pure CP 2 model, there is no degeneracy between the spin nematic state and ferromagnetic state in our model because the SU(3) global symmetry is broken. As shown in Fig. 2, the spin nematic state is partially broken around the soliton because the expectation values S a become finite. Fig. 3 shows that (S a ) 2 of the solution (3.28) are axially symmetric, although the expectation values S a have angular dependence.

Exact solutions off the BPS point
where A k is a constant background field, as before. Finally, the boundary term H Boundary is defined by Eq. (3.1) with the negative sign in the r.h.s., the same as before. Note that we also introduced constant terms in Eqs. (3.33) and (3.34) in order to guarantee the finiteness of the total energy. Clearly, the Hamiltonian (3.30) is reduced to Eq. (3.4) as we set ν 2 = µ 2 = 1.
The existence of exact solutions of the Hamiltonian (3.30) with ν 2 = µ 2 can be easily shown if we rescale the space coordinates as x → r 0 x, where r 0 is a positive constant, while the background gauge field A k remains intact. By rescaling, the Hamiltonian (3.30) becomes Setting ν 2 = µ 2 and choosing the scale parameter r 0 = ν −2 , one gets Notice that since the solutions (3.19) with P i being arbitrary constants are holomorphic maps from S 2 to CP 2 , they satisfy not only the variational equations δH ν 2 =µ 2 =1 = 0 but also the equations δH D = 0, where δ denotes the variation with respect to n with preserving the constraint (2.6). Therefore, the solutions also satisfy the equations δH r 0 =ν −2 ν 2 =µ 2 = 0. This implies that, in the limit µ 2 = ν 2 , the Hamiltonian (3.30) supports a family of exact solutions of the form Z(ν 2 ) = exp iν 2 A + z + c , (3.37) where c is a three-component complex unit vector. Since the solution (3.37) is a BPS solution of the pure CP 2 model with the positive topological charge Q, one gets H D [Z(ν 2 )] = 16πρQ. In addition, the lower bound at the BPS point (3.7) indicates that H ν 2 =µ 2 =1 [Z(ν 2 = 1)] = −16πρQ. Combining these bounds, we find that the total energy of the solution (3.37) is given by Since the energy becomes negative if ν 2 < 2, we can expect that for small values of the coupling ν 2 , the homogeneous vacuum state becomes unstable, and then separated 2D Skyrmions (or a Skyrmion lattice) emerges as a ground state.

Axial symmetric solutions
In this section, we study baby Skyrmion solutions of the Hamiltonian (3.30) with various combinations of the coupling constants. Apart from the solvable line, no exact solutions could find analytically, and then we have to solve the equations numerically. Here, we restrict ourselves to the case of the background field given by Eq. (3.25). For the background field (3.25), by analogy with the case of the single CP 1 magnetic Skyrmion solution, we can look for a configuration described by the axially symmetric ansatz Z = sin F (r) cos G(r)e iΦ 1 (θ) , sin F (r) sin G(r)e iΦ 2 (θ) , cos F (r) , (4.1) where F and G (Φ 1 and Φ 2 ) are real functions of the plane polar coordinates r (θ). The exact solution on the solvable line ν 2 = µ 2 with axial symmetry can be written in terms of the ansatz with the functions Further, the solution (3.28) is given by Eq. (4.2) with ν 2 = 1. This configuration is a useful reference point in the configuration space as we discuss below some properties of numerical solutions in the extended model (3.30). For our numerical study, it is convenient to introduce the energy unit 8ρ and the length unit κ −1 , in order to scale the coupling constants. Then, the rescaled components of the Hamiltonian with the ansatz (4.1) become where the prime ′ and the dot˙stands for the derivatives with respect to the radial coordinate r and angular coordinate θ, respectively. The system of corresponding Euler-Lagrange equations for Φ i can be solved algebraically for an arbitrary set of the coupling constants, and the solutions are where m is an integer. Without loss of generality, we choose m = 0 by transferring the corresponding multiple windings of the phase Φ 2 to the sign of the profile function G. Then, the system of the Euler-Lagrange equations for the profile functions with the phase factor (4.7) reads +4 sin 2F cos 2 F − sin 2 F cos 2 G 1 + sin 2 G − sin 2F cos 2F (1 + 2 sin 2 G) , (4.11) δH pot δF = 6r sin 2F, (4.12) +r sin 2F F ′ + sin 2 F (1 − 3 cos 2G) + sin 4 F 1 + 3 cos 2G − 2 cos 2 2G , (4.14) δH ani δG = r 8 √ 2 cos F sin 3 F cos G 1 − 3 sin 2 G + 16 sin 4 F cos 3 G sin G − sin 2 2F sin 2G , We solve the equations for ν 2 = µ 2 numerically with the boundary condition which the exact solution (4.2) satisfies. This vacuum corresponds to the spin nematic state (3.29). Let us consider the asymptotic behavior of the solutions of the equations (4.8). Near the origin, the leading terms in the power series expansion are where c F and c G are some constants implicitly depending on the coupling constants of the model. To see the behavior of solutions at large r, we shift the profile functions as   , which is expected to be −2 by the scaling argument. For ν 2 = 1.5, we used the exact solution (4.2) so that the "Derrick" and topological charge for ν 2 = 1.5 are exact values.
Then, one obtains linearized asymptotic equations on the functions F and G of the forms (4.20) Unfortunately, the equations (4.20) may not support an analytical solution. However, these equations imply that the asymptotic behavior of the profile functions is similar to that of the functions (4.2), by a replacement ν 2 κ with (ν 2 + 3µ 2 )/4. Indeed, the asymptotic equations (4.20) depend on such a combination of the coupling constants, and there may exist an exact solution on the solvable line with the same character of asymptotic decay as the localized soliton solution of the equation (4.8).
To implement a numerical integration of the coupled system of ordinary differential equations (4.8), we introduce the normalized compact coordinate X ∈ (0, 1] via The integration was performed by the Newton-Raphson method with the mesh point N MESH = 2000. In Fig. 4, we display some set of numerical solutions for different values of the coupling ν 2 at µ 2 = 1.5 and their topological charge density Q defined through Q = 2π rQdr. The solutions enjoy Derrick's scaling relation and possess a good approximated value of the topological charge, as shown in Table 1. One observes that as the value of the coupling ν 2 becomes relatively small, the function G is delocalizing while the profile function F is approaching its vacuum value everywhere in space except for the origin. This is an indication that any regular non-trivial solution does not exist ν 2 = 0.

Asymptotic behavior
Asymptotic interaction of solitons is related to the overlapping of the tails of the profile functions of wellseparated single solitons [3]. Bounded multi-soliton configurations may exist if there is an attractive force between two isolated solitons.
Considering the above-mentioned soliton solutions of the gauged CP 2 NLσ-model, we have seen that the exact solution (4.2) has the same type of asymptotic decay as any solution of the general system (4.8). Therefore, it is enough to examine the asymptotic force between the solutions on the solvable line (4.2) to understand whether or not the Hamiltonian (3.30) supports multi-soliton solutions of higher topological degrees. Thus, without loss of generality, we can set µ 2 = ν 2 .
Following the approach discussed in Ref. [3], let us consider a superposition of two exact solutions above. This superposition is no longer a solution of the Euler-Lagrange equation, except for in the limit of infinite separation, because there is a force acting on the solitons. The interaction energy of two solitons can be written as where H sp (R) is the energy of two BPS solitons separated by some large but finite distance R from each other, and H exact stands for the static energy of a single exact solution. Notice that the lower bound of the Hamiltonian (3.30) with µ 2 = ν 2 is given where the equality is enjoyed only by holomorphic solutions. Therefore, we immediately conclude where the equality is satisfied only at the limit R → ∞. It follows that the interaction energy is always positive for finite separation, and the interaction is repulsive. Since the exact solution has the topological charge Q = 2, it implies that there are no isolated soliton solutions with the topological charge Q ≥ 4 in this model. Note that, however, as the BPS solution (3.19) suggests, there can exist soliton solutions with an arbitrary negative charge, which are topological excited states on top of the homogeneous vacuum state.

Conclusion
In this paper, we have studied two-dimensional Skyrmions in the CP 2 NLσ-model with a Lifshitz invariant term which is an SU(3) generalization of the DM term. We have shown that the SU(3) tilted FM Heisenberg model turns out to be an SU(3) gauged CP 2 NLσ-model in which the term linear in a background gauge field can be viewed as a Lifshitz invariant. We have found exact BPS-type solutions of the gauged CP 2 model in the presence of a potential term with a specific value of the coupling constant. The least energy configuration among the BPS solutions has been discussed. We have reduced the gauged CP 2 model to the (ungauged) CP 2 model with a Lifshitz invariant by choosing a background gauge field. In the reduced model, we have constructed an exact solution for a special combination of coupling constants called the solvable line and numerical solutions for a wider range of them.
For numerical study, we chose the background field, generating a potential term that can be interpreted as the quadratic Zeeman term or uniaxial anisotropic term. One can also choose a background field generating the Zeeman term; if the background field is chosen as A 1 = −κλ 7 and A 2 = κλ 5 , the associated potential term is proportional to S 3 . The Euler-Lagrange equation for the extended CP 2 model with this background field is not compatible with the axial symmetric ansatz (4.1). Therefore, a two-dimensional full simulation is required to obtain a solution with this background field. This problem, numerical simulation for non-axial symmetric solutions in the CP 2 model with a Lifshitz invariant, is left to future study. In addition, the construction of a CP 2 Skyrmion lattice is a challenging problem. The physical interpretation of the Lifshitz invariants is also an important future task. The microscopic derivation of the SU(3) tilted Heisenberg model [21] may enable us to understand the physical interpretation and physical situation where the Lifshitz invariant appears. Other future work would be the extension of the present study to the SU(3) antiferromagnetic Heisenberg model where soliton/sphaleron solutions can be constructed [41][42][43].
We restricted our analysis on the case that the additional potential term µ 2 H pot is balanced or dominant against the anisotropic potential term ν 2 H ani , i.e., ν 2 ≤ µ 2 . We expect that a classical phase transition occurs outside of the condition, and it causes instability of the solution. At the moment, the phase structure of the model (3.30) is not clear, and we will discuss it in our subsequent work.
Moreover, it has been reported that in some limit of a three-component Ginzburg-Landau model [44,45], and of a three-component Gross-Pitaevskii model [46,47], their vortex solutions can be well-described by planar CP 2 Skyrmions. We believe that our result provides a hint to introduce a Lifshitz invariant to the models, and that our solutions find applications not only in SU(3) spin systems but also in superconductors and Bose-Einstein condensates described by the extended models, including the Lifshitz invariant.