Extended multi-scalar field theories in (1+1) dimensions

We present the explicit construction of some multi-scalar field theories in (1+1) dimensions supporting BPS (Bogomol'nyi--Prasad--Sommerfield) kink solutions. The construction is based on the ideas of the so-called extension method. In particular, several new interesting two-scalar and three-scalar field theories are explicitly constructed from non-trivial couplings between well-known one-scalar field theories. The BPS solutions of the original one-field systems will be also BPS solutions of the multi-scalar system by construction, and therefore we will analyse their linear stability properties for the constructed models.

analytical solutions, e.g the so-called trial orbit method [33]- [36], which allow us decouple the field equations by introducing very specific orbit equations of the form O(φ 1 , . . . , φ n ) = 0, that is constraints in the target space. Although useful, this method has been shown not so efficient when one is looking for new analytical multi-scalar models.
In that scenario, a simplifier tool in searching for kink solutions is provided by the so-called BPS (Bogomolnyi-Prasad-Sommerfield) method [37,38], which allows to find solutions from firstorder differential (BPS) equations instead of second-order Euler-Lagrange equations. BPS solutions correspond to static configurations of minimal energy. Although simpler, the problem of solving analytically first-order differential coupled equations is still not easy, and then it is necessary to use additional procedures to sort it out.
In this work, we will use the extension method, originally proposed in [1,2], to systematically construct several new multi-scalar field theories in (1+1) dimensions supporting BPS states, starting from a system of several one-scalar models. The basic ingredients in the construction are the socalled deformation functions and its inverses [39]- [42], which provide suitable links between the fields to be coupled. In addition, the method has the nice advantage that the BPS solutions of the one-field systems are also solutions for the multi-scalar system.
We aim that the new models constructed in the present work could improve the knowledge and understanding of the analytical solutions of multi-scalar systems, and believe that they have potential applications to cosmological models, to the study of kink scattering process of multisolitons, and also to analyse integrability and self-duality properties in the multi-scalar-models. For that reason, special attention will be given to the theories with periodic potentials with infinitely degenerate vacua, as is the case of sine-Gordon model, or even more exotic models as the one studied in [43,44].
This paper is organized as follows. In section 2, we briefly review some basics aspects of scalar BPS theories, introducing the superpotential function, or sometimes called the prepotential function [13], as the key ingredient of the whole construction. In section 3, we present the main ideas of the deformation procedure, and then we apply it to obtain several examples of deformed theories. In section 4, we introduce the extension method and construct several new interesting two-scalar field theories. The linear stability of the BPS solutions for these new models will be discussed in section 5. In section 6, we construct some new three-scalar fields extended models by applying a straightforward generalization of the extension method for three-field systems [2], and also analyse the linear stability of their BPS solutions. Final remarks and comments of our work are presented in section 7. In appendix A, we have summarized some basics features of the underlying exactly solvable potentials which appear in the linear stability analysis. Finally, appendix B contains some explicit calculations on the derivation of the superpotentials for the three-field systems.

General settings
Let us start considering theories with n real scalar fields φ a (x, t), a = 1, . . . n, in (1 + 1) dimensions described by the following Lagrangian, where µ = {0, 1}, with metric convention η µν = diag(+1, −1), ∂ µ ≡ ∂ ∂x µ , x 0 = t, x 1 = x, in natural units. The corresponding field equations for φ a (x, t) are given by and for the static configurations (∂ t φ a = 0), we get φ a (x) = ∂V ∂φ a , (2.3) where we are using standard conventions φ a ≡ d dx φ a . These equations can be rewritten as follows, It is worth pointing out that there is no summation assumed in eq. (2.4). The corresponding energy functional for the static configurations reads, Finite energy configurations require existence of the boundary condition φ a (±∞) → 0, and a potential possessing at least one vacuum value, V (φ) = 0, such that φ a (±∞) →φ ± a . When two or more minima exist, then the potential supports topological configurations connecting two adjacent minimaφ − andφ + . Now, by introducing a smooth function of the scalar fields W (n) (φ), sometimes named superpotential or pre-potential [13], the potential V can be written as, where W (n) φa stands for ∂W (n) ∂φa . Then, the field equations can be rewritten as a set of coupled first-order differential equations, with energy given by, The solutions of the first-order differential equations (2.7) with non-zero energy (2.8) are named BPS states. These minimum energy static configurations are also solutions of the second-order differential equations (2.3), which can be understood from the self-duality properties of the BPS theories as claimed in [13]. In fact, for a given field theory, the Bogomolnyi bound (2.8) only depends on the boundary conditions, and not on the field configuration, which means that E BPS is a homotopy invariant, that is invariant under any smooth deformation of the field configurations. These interesting properties makes BPS states so attractive, and it will be the main goal of our work to look for them.

Deforming one-scalar field theories
Recently, it has been proposed an interesting procedure to generate infinite families of one-field theories with topological (kink-like) or non-topological (lump-like) solutions, which is now referred as deformation procedure [39,40]. The main idea is to start from a given "seed" one-scalar field theory possessing static solutions, and then perform a field transformation on the target space to obtain a new one-scalar field theory that also supports static solutions. In particular, we will focus in theories supporting BPS solutions. Let us start from a one-scalar field model described by the Lagrangian, which supports BPS solutions satisfying the first-order differential equation, Now, we introduce an invertible smooth function f on the target space, called the deformation function, such that φ(x) = f (ϕ(x)), (3.3) where ϕ is a new (deformed) scalar field. This function also allows us to introduce a new (deformed) one-field model described by the following Lagrangian, which satisfies the first-order equation, (3.5) providing that the two potentials are related to each other through the deformation function as follows, where f ϕ = df dϕ . This also implies that the two superpotentials are related in the following form, It is worth noting that the static solutions for both scalar fields are related by eq. (3.3), and then by replacing in the first-order differential equation, we find that they also satisfy the following important constraint, . (3.8) This relation between the fields will play a central role in constructing multi-scalar field theories supporting BPS kink-like solutions. In fact, this relation has been already used in [45] for studying systems of two coupled fields in (1+1)-dimensions through orbit equation deformations.
Let us now consider a few interesting examples to illustrate the deformation procedure. First of all, we start with the standard φ 4 model [42], whose potential can be written as where α > 0, is a real dimensionless parameter. This potential satisfies the first-order differential equation (3.10) and supports the following static solution, φ(x) = tanh(αx). (3.11) Now, in order to obtain the deformed model, we consider the following function, After using the deformation function, we obtain that the deformed potential describes the so-called ϕ 6 -like model [41]. The corresponding the first-order differential equation is given by, with the following topological solutions ϕ ± (x) = ±(1 + tanh(αx)), (3.15) which are quite similar to the solutions of the standard φ 6 model [42]. This model possesses three minima at the valuesφ = {0, ±1}, and supports two symmetric BPS sectors [41]. Interestingly, this example shows us that the deformation procedure can change the number of vacua of the seed model, and consequently changing the number of topological sectors. As a second example, let us consider again the φ 4 model as the seed model, with superpotential given by (3.10), and introduce the following periodic deformation φ = f (χ) = sin(βχ), (3.16) where χ is the new deformed field. The corresponding deformed model describes the sine-Gordon model given by the following first-order equation This very well-known model has infinite degenerate vacua at the valuesχ k = k − 1 2 π β , with k ∈ Z, and correspondingly an infinite number of equivalent topological sectors. Connecting the minimā χ 0 andχ 1 , its static solution can be written as follows, As a last example, let us consider the bosonic exotic scalar model (E-model) investigated in [43,44] as our seed model, which is described by the following first-order field equation, This model also has infinitely degenerate trivial vacua at the pointsη k = −1 + e (k− 1 2 )π , with k ∈ Z. However, in this case the infinite number of BPS sector are not equivalent, since the BPS energy depends on the topological sector. A simple kink-like solution for this model connecting the vacuā η 0 andη 1 , can be written as follows, η(x) = exp (arctan(sinh(αx))) − 1. (3.20) Now, by considering the following deformation function, we get the sine-Gordon model, which has been already described in eq. (3.17).

Constructing two-scalar fields models
Let us now describe the method to construct two-scalar field theories from one-scalar field theories.
To do that, we will use the deformation method introduced in the last section. The starting point is the first-order equation for the seed one-scalar field model, which supports static solutions. Now, by introducing a deformation function, i.e. φ = f (ϕ), we can rewrite eq. (4.1) in two different (but equivalent) ways, namely where we have made full use of the function φ → f (ϕ) in the first expression in order to make W φ a function depending only on ϕ, while in the second expression we have made partial use of this function in order to make W (1) φ a function depending on both fields φ and ϕ. Of course, there is an ambiguity in obtaining the last expression since it would depend on how this "lifting" from φ-space to (φ, ϕ)-space is made. Then, there will be an infinite number of resulting models once we chose the form of W (1) φ (φ, ϕ). Some of these models would be trivial, and some of them even do not longer support kink-like solutions. However, our main goal here is to construct models that do support BPS solutions, and that will be reach by choosing carefully that form.
Taking into account that the deformed field ϕ also satisfied a first-order differential equation, the same procedure can be applied in order to rewrite the equation as follows, where again there have been both a full as well as a partial use of the deformation (inverse) function, ϕ = f −1 (φ). Note also that now the eq. (3.8) can be rewritten in several different ways. In fact, in order to proceed with the extension method, we define the new two-fields superpotential through the following ansatz, φ (φ) + p 1 g(ϕ) + p 2 g(φ, ϕ) + p 3 g(φ), (4.5) where for consistency, the parameters a i , b i , p i , and q i with i = 1, 2, 3, must satisfy the following constraints and g andg are arbitrary functions required for the consistency conditions, which can be written as follows, Thus, by substituting eqs. (4.5) and (4.6) in eq.(4.8), we get the following constraint (4.9) The above constraint allows us to obtain the specific form of the functions g andg. After doing that, we can go back to the system given by eqs. (4.5) and (4.6), and perform simple integrations to finally determine the form of W (2) (φ, ϕ). In what follows, we will illustrate the extension method by explicitly constructing new interesting two-scalar fields models.

φ 4 model coupled with ϕ 6 -like model
Now, we will consider a model constructed through the coupling of the standard φ 4 model and the ϕ 6 -like model [41]. Let us start from eq. (3.10), namely with the deformation function given in (3.12), namely, The deformed model is the ϕ 6 -like model, whose superpotential satisfies eq. (3.14), Now, we will use the deformation function to write eq. (4.10) in three equivalent arbitrary forms, where one of these must be a function of only φ, another function of φ and ϕ, and the last function only ϕ, that is Similarly, we can use the inverse deformation function to write eq.(4.12) in the following three different forms, where the constant parameter = ±1 for solutions ϕ ± (3.15), respectively. Now, by substituting these expressions directly into the constraint (4.9), we get where we have chosen p 1 = 0 and q 3 = 0, so p 3 = −p 2 and q 1 = −q 2 , for simplicity. In addition, we can use the deformation function, and its inverse, to write By substituting the above results in eqs. (4.5) and (4.6), we obtain respectively, which upon being integrated results in the following superpotential, where we have just renamed the parameters: a ≡ 2a 1 , b ≡ b 2 , and c ≡ 2a 1 + a 2 − b 2 − 2b 1 + 1. This superpotential describes the coupling between the φ 4 model and the ϕ 6 -like model, and therefore from now on we will name it as the extended (φ 4 + ϕ 6 l ) model. Note that there are several models that can be considered depending on the choice of these parameters. This model contains three minima at the following values: m 1 = (−1, 0), m 2 = (1, 2), and m 3 = (1, −2). It supports then three topological sectors, with only two of them are BPS, they are: the sector connecting m 1 and m 2 , and the one connecting m 1 to m 3 , with the explicit solutions given by eq.(3.11) and (3.15), namely φ(x) = tanh(αx), ϕ ± (x) = ±(1 + tanh(αx)), (4.20) which energy is E BPS = 8α/3. On the other hand, the non-BPS configuration, connecting the minima m 2 and m 3 , do not satisfy the first-order equation. In this case, we can write an explicit solution for the specific values a = 1, b = −1, and c = 0, with energy E = 16α/3, which is twice the energy of the BPS sectors, as it was already expected. It is worth noting that the corresponding anti-kink configurations, ± (x) = ±(1 − tanh(αx)), (4.22) are also in the BPS sectors connecting m 2 and m 3 minima to m 1 , respectively. There are several others topological sectors that appear after chosing the values of the parameters. For instance, for a = b = 0 and c = −1, we recovery the BPS sector associated with the φ 4 model, connecting the two minima (±1, 0), with energy E BPS = 4α/3. On the other hand, we can also verify that the trivial configuration φ = 0 does not belong to the minima space of this potential. Finally, we would like to pointing out that for the values a = c = 1, and b = 0, the superpotential W (2) (φ, ϕ) becomes harmonic, and consequently all the solution will be BPS solutions [46,47].

φ 4 model coupled to sine-Gordon model
Our starting point will be again the φ 4 model. Now, by using the deformation function (3.16) we will write the right-hand side of eq. (3.10) in the following forms, and similarly eq. (3.17) as We see that the last two forms for W (1) χ contain square root functionals, and then it is convenient to consider the following choice of parameters b 2 = b 3 = 0. Also, without loss generality we choose p 1 = p 2 = q 3 = 0. This implies immediately that b 1 = 1, p 3 = 0, and q 2 = −q 1 . Then, by solving the constraint (4.9) to determine theg-function, we get which can be rewritten, by using the deformation function, as follows By substituting the above results in eqs. (4.5) and (4.6), we get   Integrating out these expressions we finally obtain the two-fields superpotential, which is given by the following form This superpotential describes the coupling of the φ 4 and sine-Gordon models, which from now on will be named as the extended (φ 4 +sG) model. The static kink-like solutions, are BPS solutions of this model connecting the minima m 1 = (−1, − π 2β ) and m 2 = (1, π 2β ), with BPS energy given by, In general, we can verify that this potential possesses minima at the points (−1, (2k − 1 2 ) π β ), and (+1, (2k + 1 2 ) π β ), with k ∈ Z. It is also worth highlighting the existence of other BPS solutions. For instance, a particular solution is providing that the parameters a 1 and a 2 are restricted to satisfy a 2 < (1 − a 1 ), and a 2 = 2(1 − a 1 ), otherwise φ (−) (x) becomes an exponential or a constant solution, respectively. By chosing k = 0, we notice that this solution connects the minimum m 1 to a new minimum . We see that the BPS energy of this solution is given by, Another possible solution is, In this case a 2 > (1 − a 1 ), and again a 2 = 2(1 − a 1 ). For k = 0, this solution connects the minimum m 2 to a new minimum m 4 = (− 1−a 1 1−a 1 −a 2 , + π 2β ). We also note that the solution (4.34) possesses the same BPS energy as the solution (4.32), given by eq. (4.33). We have plotted φ (±) (x) in figure 1.

E-model coupled to sine-Gordon model
Let us consider now the sine-Gordon model and the E-model described by the first-order field equations (3.17) and (3.19), respectively, together with the deformation function (3.21). Then, we write the following expressions, and By choosing p 1 = q 3 = 0 in the constraint (4.9), we get After integrating, we find the following solution, As it was done before, we can use the deformation function to write Using these results, we obtain We then can construct the corresponding two-fields superpotential, This two-parameters superpotential leads us to a potential V (η, χ) describing the coupling of the sine-Gordon model and the E-model, which from now on we will named as the extended (sG+E) model. This superpotential supports the static kink-like solutions (3.18) and (3.20), connecting the minima m 1 = e −π/2 − 1, −π/2β and m 2 = e π/2 − 1, π/2β , with BPS energy given by We notice that for the particular values of the parameters a 2 = 0 and b 3 = 1, it is possible to obtain other BPS solutions, at least numerically. In this case, we have that η(x) = −1, and χ(x) has to satisfy It is interesting to see that the associated potential V (−1, χ) represents a modification of the sine-Gordon model. In fact, we have verified that it has infinite minima, and supports BPS solutions. The minima are located approximately at the following points, (4.52) Then, there are at least three type of topological sectors, a small one for k < 0, a medium one for k = 0, and the large one for k > 0. The corresponding BPS energies are given by Unfortunately, we have not been able to obtain the corresponding analytical solutions of the firstorder equation (4.51) for the kink solutions associated to each topological sector. However, we did construct them numerically 3 . For instance, we have plotted in figure 2 the numerical kink solutions connecting the minimaχ 0 toχ +1 , for several values of the parameter β. We can see that for β 0.1, the profile tends to fit the sine-Gordon kinks. For greater values it undergoes a rapid deformation. It is worth noting that this first-order approximation fails when β = e (1−2k)π/3 , for k ≤ 0. However, there is nothing special about those values, but in that case a second-order approximation would be necessary.

Linear stability of the BPS configurations
Let us now discuss the linear stability for the two-scalar fields models we have constructed. The main issue is basically to analyse the spectrum of the corresponding Schrödinger-like operator associated with the normal modes of the classical model. The stability will be ensured when this Schrödinger-like is positive semi-definite, implying that negative eigenvalues will be absent from its spectrum, and the zero mode will correspond to the lowest bound state [48]- [51]. First of all, it is well-known for one-field models that the static configurations of the φ 4 model (3.11), the ϕ 6 -like model (3.15), the sine-Gordon model (3.18), and the E-model (3.20) are all stable [29,41,43,52], with the corresponding Schrödinger-like operators related to the so-called Rosen-Morse II potential (or modified Poschl-Teller potential) for the first three models, and to the so-called Scarf-II (hyperbolic) potential in the latter case [53] (see more details in appendix A). Now, the stability analysis for multi-fields models is in general a highly non-trivial problem. Here, we will follow the line of reasoning introduced in [51], to study the stability of static solutions in the two-scalar field models constructed in section 4. The starting point is to consider a pair of static solutions, say φ s (x) and ϕ s (x), and then introduce small fluctuations around these solutions, given in the following form where ρ k and σ k are the small perturbations, when compared to the static configurations. Now, by substituting the fields φ(x, t) and ϕ(x, t) into the second-order equations (2.3), and considering only first-order terms in the fluctuations, we obtain the Schrödinger-like equation Notice that the derivatives of the potential V (φ, ϕ) are written in terms of the static fields φ s (x) and ϕ s (x). In addition, as it can be seen from eqs. (5.1) and (5.2), linear stability requires that the eigenvalues of H have to be positive semi-definite, i.e. w 2 k ≥ 0, with the zero mode HΨ 0 (x) = 0, being given by where the normalization constant N 0 can be chosen to be the unit. When the potential V (φ, ϕ) supports BPS states the Hamiltonian in (5.3) can be written as follows, where the first-order operators have been introduced. Note that A † ± = A ∓ , implies that the Schrödinger-like operator H is always positive semi-definite for the BPS case, thus ensuring the linear stability of the BPS configurations. In this case, the ground state coincides with zero mode, and can be written as Here, we will have an inherent difficulty regarding the explicit determination of the eigenvalue spectrum of the associated Schrödinger-like operator. As it can be seen, the coupling between static fields results in the coupling of the fluctuations in (5.3). However, the problem turns to be more manageable if we take advantage of the first-order operators (5.6), and diagonalize the matrix W, to obtain where the respective eigenvalues are [51] By substituting (5.8) in (5.5), we obtain two decoupled eigenvalue equations where the quantum mechanical potentials are given by It is worth pointing out that in general this method would require certain simplifications since the square root term appearing in eq.(5.9) brings some complications for the explicit analytical calculations. In what follows, we will try then to simplify this term whenever is possible to perform the analytical analysis of the stability of the BPS configurations for the two-fields models we have constructed. Otherwise, the corresponding spectral problems should be analysed from a numerical point of view.

The extended
Let us first study the stability of the BPS solutions (4.20) of the extended (φ 4 + ϕ 6 l ) model. From the superpotential (4.19), and the BPS solutions eq.(4.20), we obtain where we are assuming that α = 1, and b ≥ |1 + c|, (5.15) in order to simplify the square root term. Using these results, we get the corresponding quantum mechanical potentials (see figure 3), They have again the form of the Rosen-Morse II potentials [53] (see appendix A). In this case the parameters of the potential U + (x) (5.16) are given by, Then, from the stability condition (A.3), we see that the parameters have to satisfy Now, let us choose some interesting values for the parameters. From eq. (5.15), we note that if b = 0, then c = −1, and we get that both potentials are equal, and stability can be guaranteed. On the other hand, when b < 0 the condition (5.15) is not satisfied, and then stability cannot be proven in that case, at least analytically. Finally, by considering b > 0, together with the conditions (5.15) and (5.19), we find that We can also see that, since the potential U − has eigenvalues E 0 = 0 and E 1 = 3, it will have common eigenvalues with U + only if In the table 1, we have chosen some particular values for the parameters in order to illustrate our results. For all these cases, the stability of the solutions is guaranteed.

The extended (φ 4 + sG) model
Now, we will analyze the stability of the BPS solutions (4.30) of the extended (φ 4 + sG) model described by the superpotential (4.29). For sake of simplification, we have chosen a 2 = −2a 1 , and α = 1, to get where sgn(x) is the signum function. Then, we obtain the following quantum-mechanical potentials, In general, they are also associated to the Rosen-Morse II potential (A.1), with B = 0 and α = 1. However, these quantum-mechanical potentials are discontinuous, as it can be seen from figure 4, and that novel feature will require special attention in order to determine the eigenvalues. To do that, we will use the procedure introduced in [54,55] to determine the energy levels of composite potentials, which is based on the so-called Green function factorization theorem [56]. The main idea consists in decomposing the discontinuous potential into two "pieces", namely where θ(x) is the unit step function, and U (L/R) (x) are continuous and symmetric (around the origin) potentials, for which the corresponding energy levels and wave functions for all stationary states are assumed to be known. Then, by considering the Green functions G (L/R) (x, x ; E) associated to each potentials U (L) and U (R) , the allowed eigenvalues E of the composite potential U will be given by the solutions of the following transcendental equation [54,55], In our case, both Green functions are associated to the solvable Rosen-Morse II potential, and its explicit formula can be written as follows [57], where A is the parameter given in (A.1), Γ(z) is the Gamma function, and are the associated Legendre polynomials, which are defined in terms of the hypergeometric function F (a, b; c; z). The corresponding discrete eigenvalue spectrum satisfies Therefore, for the U + (x) potential, we have that and then G (L) From these results, we see that the transcendental equation (5.27) becomes which only allows E = 0 in the spectrum of the discontinuous potential U + (x). From the eq.(5.30), we can see that the zero energy eigenvalue is common to the decomposed potentials, which is consistent with the fact that E = 0 is a pole of the Green function (5.28) for both cases. An identical transcendental equation will be obtained for the potential U − (x), since U , which again will allow only the zero energy eigenvalue. These results lead us to ensure the stability for the BPS solutions of the extended (ϕ 4 +sG) model, at least for our particular choice of parameters.

The extended (sG+E) model
Finally, let us study the stability of the BPS solutions (4.49) of the extended (sG+E) superpotential (4.48). In this case, we find that − β 2 e 2 arctan(sinh(αx)) − 1 b 3 + a 2 β 2 e 2 arctan(sinh(αx)) tanh(αx) In order to simplify the root term in eq.(5.34), and study analytically the associated quantummechanical potentials, we could choose a 2 = b 3 = 0, obtaining From the analytical point of view it is quite complicated to study these quantum-mechanical potentials for general values of the parameters a 2 and b 3 . Instead, we will perform a more qualitative and approximated analysis of the bound states for some values of the parameters. At this point, we can only guarantee that stability does exist at least for some very small values of the parameters, that is, the potentials possess the zero-mode as their fundamental bound state, and there is no negative energy eigenvalues. Of course, a more precise analysis requires a deeper numerical study.
Before doing that, let us take a look of the potential deformations for some small values of the parameters. In the figures 6 and 7 we have plotted the potentials for some configurations with a 2 = 0, and small values of b 3 . While, in figures 8 and 9, we have plotted configurations with b 3 = 0, and small values of a 2 .     Now, for very small values of the parameters ( 10 −2 ), it is clear that the quantum-mechanical potentials converge to exactly solvable problem, see also figure 10. In that case it is possible to apply the time-independent perturbation theory for the calculation of the energy eigenvalues corrections. To do that let us first consider the case when a 2 = 0 and b 3 = λ 10 −2 , so we have a perturbed Hamiltonian which consists of two parts with where U ± are the exactly solvable potentials (5.35) and (5.36), and the first-order corrections U (1) ± are given in this case by Then, the eigenvalues E k of the perturbed problem can be expanded in a power series in the parameter λ as follows where E where ρ k (x) and σ k (x) will be given in this case by the Scarf II with energy E Then, by substituting in eq. (5.42) we find that the first-order correction to the zero mode also vanishes. Therefore, we can ensure that in the weak coupling regime, for a 2 10 −2 and b 3 10 −2 , the stability of the extended (sG+E) BPS solutions. Of course, for greater values of the coupling parameters, we should do a more complete analytical or numerical analysis of the spectral problem. We will leave this specific issue to be explored in other future work.

Extended three-scalar fields models
Now, we will construct some new three-scalar field extended models by applying a generalization of the extension method for three-field systems [2] to the one-field systems studied so far.
6.1 φ 4 model coupled to the ϕ 6 -like and the inverted ζ 4I models Let us start by considering the coupling of the standard φ 4 model with the ϕ 6 -like model, and also with the so-called the inverted ζ 4I -model [39]. The starting point is again the first-order equation together with the deformation functions, from which we have the corresponding the first-order equations for the deformed models, where we have defined ω = sgn(φ). Their corresponding static solutions are Now, the main idea of the method can be straightforwardly generalized to three-fields. First, we write the right-hand side of eq. (6.1) now in seven different and equivalent forms by using the deformation functions and their inverse functions, as follows Similarly, for eq. (6.5) and for the eq. (6.6) we have Now, we will use a generalization of the ansatz used in eqs. (4.5) and (4.6) for the case of three-field systems, in the following form φ (φ, ϕ, ζ) + p 1 g(ϕ) + p 2 g(φ, ϕ) + p 3 g(φ) + p 4 g(ζ) +p 5 g(φ, ζ) + p 6 g(ϕ, ζ) + p 7 g(φ, ϕ, ζ), (6.11) where the parameters must satisfy the following conditions In addition, the g-functions are determined from the following constraints (see in appendix B more details of the full derivation), ζϕ . (6.15) Using the explicit results for the g-functions into eqs. (6.11)-(6.13), yields ζ (φ, ϕ, ζ) = −αζ (φ − c 6 (1 + φ − |ϕ|)) . (6.16) After integrating, we finally obtain the following three-field superpotential From now on we will named this three field model as the extended (φ 4 + ϕ 6 l + ζ 4I ) model, for which possesses the static configurations given in eq.(6.7) are BPS solutions, with energy The issue that arises from these results concerns the linear stability of the solutions for this superpotential. Although the stability analysis for three-field systems follows the same steps that the one presented in section 5, it is actually further more complicated mostly because of the diagonalization of the Schrödinger-type operator, which is the key point in order to find the normal mode fluctuations. In this case, we will have where ρ k , σ k , and ξ k are the fluctuations around the static solutions φ s (x), ϕ s (x), and ζ s (x). Considering the dynamics of these three time-dependent fields up to first-order, we will obtain a corresponding Schrödinger-like equation For the case of BPS potentials, we can write this Hamiltonian in terms of linear operators, namely Our strategy again will be trying to diagonalize the matrix W, and then the Schrödinger-type equation will be split into three equations, which will be analysed separately.
In the case of the BPS solutions (6.7) of the (φ 4 + ϕ 6 l + ζ 4I ) model (6.17), this matrix takes the following form where we have chosen the parameters being a 2 = 0 e c 6 = 1 for simplicity. By computing its corresponding eigenvalues, we find Now, by setting β = 1, we will find that the quantum mechanical potentials are given by, (6.24) which are again Rosen-Morse II potentials (A.1). The U 0 potential has parameters A = 2α and B = 0, and possesses eigenvalues E 0 = 0 and E 1 = 3α 2 . The other two potentials U ± have parameters A = 2α and B = ±2α 2 respectively, and only have the ground state E 0 = 0. Therefore, for these choice of parameters, we will have stability guaranteed. Another possible choice would be a 2 = 0 and c 6 = 0, however we will get essentially the same results.

φ 4 model coupled to sine-Gordon model and the E-model
Let us now construct a model obtained by the coupling of the standard φ 4 model with sine-Gordon, and the E-model. The first-order equation for each one of these models are given by eqs. (3.10), (3.17), and (3.19), namely The deformation functions connecting the three models have the following forms, By using these deformation functions and their inverse functions, we get the following expressions, Similarly, we have and Now, from above parametrizations we can derive explicitly the functions g,g, andĝ (see appendix B for more details of the full derivation). After doing that, we find which after being integrated lead us to the three-field superpotential We see that the first quantum-mechanical potential derived from eq. (6.40) is simply given by U 0 (x) = 4α 2 − 6α 2 sech 2 (αx), (6.42) whose energy eigenvalues are E 0 = 0 and E 1 = 3α 2 , which partially guarantees stability. However, the others two potentials U ± have complicated forms (see figure 11), which as a consequence arises some difficulties in obtaining analytical results, except for the c 1 = 0 case which decouples the sine-Gordon field. Instead of that, we will perform an approximated analysis for those cases. We see from the plots of the potentials in figure 12 that for small values of c 1 these potentials approximate to Rosen-Morse II and Scarf II profiles, respectively. Let us consider small values of the parameter, that is c 1 = λ 10 −2 , so we obtain the approximated potentials up to first-order, (6.43) where the unperturbed potentials are given by (6.44) while the first-order corrections are, We notice that the unperturbed potential U  potentials only possess one bound state, the zero mode E (0) 0 = 0. Therefore, the corresponding first-order correction E (1) 0 to the zero energy will be obtained from On the other hand, it is clear that there is enough room for several different topological sectors depending on the values of the four arbitrary parameters. In particular, if we chose c 1 = c 6 = 0, we get the following two BPS solutions, η (±) (x) = −1 + e ±π/2 , χ(x) = 1 β arctan(sinh(αx)), providing that a 5 = 2(1 − a 4 ), and also that a 5 > (1 − a 4 ) for φ (+) (x), and a 5 < (1 − a 4 ) for φ (−) (x). The kink solutions for the φ 4 field in eq. (6.46) have the very same form as the ones previously found eqs. (4.32) and (4.34) for the extended (φ 4 +sG) model, and also interpolates between the values ±1 and ±(1−a 4 ) a 5 −(1−a 4 ) , respectively. Depending on the values of the parameters, we will have two topological sectors with the corresponding BPS energies given by (6.47) By analysing the stability of the φ (−) solution, we will find that the associated W matrix reads, where we have chosen a 5 = 0 for simplicity, and therefore we have that a 4 < 1. The quantummechanical potentials will be given by from where we immediately see that the potential U 0 only possesses the eigenvalue E 0 = 0. In its turn, we can verify that the potentials U ± can be described by shifted Rosen-Morse II potentials, namely where κ = α(1 − a 4 ), and the parameters are and A − = −2αe π a 4 , B − = −2α 2 e π a 4 (1 + 2e π a 4 ), We find that the potential U + possesses the eigenvalues E 0 = 0 and E 1 = 3α 2 (1 − a 4 ) 2 , whereas the potential U − only possesses the eigenvalue E 0 = 0, if we have that − e −π 2 < a 4 < − e −π 4 . (6.55) Thus, this particular solution is stable only if the parameter a 4 satisfies the constraint (6.55), at least for our choice of parameters. Following an analogous procedure, the stability analysis of the solution φ (+) will lead us to similar results.

φ 4 model coupled to two sine-Gordon models
In this last example, we will construct a three-field system that couples the φ 4 field with two different sine-Gordon fields χ and ψ. The first-order equations are The deformation functions are, φ = f 1 (χ) = sin(βχ), (6.60) φ = f 2 (ψ) = sin(γψ), (6.61) As it was already done in the previous models, we use these functions to write the following equivalent expressions, and and also ψ (φ, χ, ψ) = α γ cos(γψ) cos 2 (βχ) − φ 2 cos(γψ) + 2φ sin(γψ) cos(βχ) . (6.65) As before, we use all of these expressions to obtain the corresponding g-functions (see details in appendix B), and then by substituting the results in eqs. (6.11) -(6.13), we have with the corresponding superpotential given by This new extended three-field superpotential describes the coupling of the φ 4 field with two different sine-Gordon fields, and will be named as the extended (φ 4 +sG 1 +sG 2 ) model. The static solutions (6.59) are BPS solutions of its first-order equations, connecting the minima m 1 = (−1, − π 2β , − π 2γ ) and m 2 = (1, π 2β , π 2γ ), with BPS energy given by Now, in order to analyse linear stability of the BPS solutions (6.59), we find that in this case the corresponding matrix W takes the following form, where we have chosen a 5 = −2a 4 and a 2 = −2a 1 , for simplicity. By diagonalizing this matrix, we will find the following eigenvalues u 0 = −2α tanh(αx), u + = −α tanh(αx), u − = −µ tanh(αx), (6.72) where in this case the parameter µ = α−αb 4 (1+ γ 2 β 2 ). Then, the corresponding quantum-mechanical potentials will be given as follows, U 0 = 4α 2 − 6α 2 sech 2 (αx), U + = α 2 − 2α 2 sech 2 (αx), U − = µ 2 − µ(µ + α)sech 2 (αx). (6.73) which are again Rosen-Morse II potentials. We see that the potential U 0 is the same as the one in eq.(6.24), and has the eigenvalues E 0 = 0 and E 1 = 3α 2 . The parameters for the potential U + are A = α and B = 0, and has only one eigenvalue, E 0 = 0. For the potential U − the parameters are A = µ and B = 0. In this case, the number of eigenvalues will be now constrained by 0 ≤ k < 1 − b 4 (1 + γ 2 β 2 ), which requires that b 4 < β 2 (β 2 +γ 2 ) in order to guarantee stability. Therefore, when 0 < b 4 < β 2 (β 2 +γ 2 ) there exists only one eigenvalue E 0 = 0. For b 4 < 0, we note that the number of bound states increases for decreasing b 4 , and then the potential could have more than one non-negative eigenvalue, guaranteeing in this way the stability of the BPS solutions. There are also several others interesting features that can be mentioned about this new model. For instance, the projection of the corresponding potential in the (χ, ψ) plane gives the following, where we have considered a 2 = a 5 = 0, and a 4 = 1 − a 1 , without loss of generality. It is worth pointing out that this potential is not BPS, and even though its minima are located at the static sine-Gordon kinks are no longer solutions of its field equations. Despite of being an interesting potential (see figure 13), we have not been able to find any explicit analytical solutions for it. It would be interesting to look for at least numerical solutions and also further explore this potential. That could be addressed in more detail in another work. On the other hand, when substituting φ = ±1 directly in (6.69), and setting a 1 = a 2 = a 4 = a 5 = 0, we end up with a different effective two-fields superpotential, and its corresponding Figure 14: Effective two-fields superpotential (on the left) and the associated potential (on the right) for the two coupled sine-Gordon fields. For both, we have plotted the values α = 1, β = 1, γ = 2, and b4 = 0.1.

potential, given by
Note that although this potential is somehow contained within the projection V (0, χ, ψ), they are actually different even if we set a 1 = 0 in eq. (6.74), and in this case the static solutions for the sine-Gordon fields given in (6.59) are BPS solutions of the first-order equation for the effective superpotential (6.76). It is worth also noting that the simple coupling between the two sine-Gordon fields contained in the last term of eq. (6.76) differs from some models previously constructed in the literature. In particular, if we eliminate the coupling term by setting b 4 = 0, then our potential will take the form of the non-integrable two-frequency sine-Gordon model considered in [58], where the authors studied how the particle spectrum of the model changes by considering the second interaction as a perturbation of the original integrable sine-Gordon model. In addition, after proper redefinitions our superpotential (6.76), also with b 4 = 0, can be also seen as a limit case of the FKZ (Ferreira, Klimas, and Zakrewski) pre-potential based on the SU (3) Lie algebra 4 [3]. However, their potential V will be quite different since a constant, real and positive-definite matrix η ab , which is basically a modified version of the associated Cartan matrix, is directly involved in the definition of the FKZ models. Despite of these differences, it would be interesting to analyse if there exist any common points between the two methods of constructing multi-scalar field theories. This issue will be addressed in future investigations.

Final remarks
In this paper, we have presented the explicit construction of several interesting new models described by two and three real scalar fields theories in (1 + 1)-dimensions supporting BPS states. The way of constructing such field theories is called the extension method, which was introduced originally in [1,2]. This method requires considering initially several (not necessarily different) one-field systems which are known to support BPS states, and that are also connected through some mappings called deformation functions. Then, the corresponding first-order equations are rewritten in several different but equivalent non-trivial ways by using such functions and their inverse functions. Doing that, the fields are then coupled by introducing an ansatz for the first-order equations for the resulting two-field model eqs. (4.5) and (4.6), and respectively for three-field model eqs. (6.11)-(6.13). To finish the procedure, some functions, called here as g-functions, are then introduced in order to guarantee smoothness of the superpotential, which are properly derived from consistency constraints (6.15). The constructed theories were obtained by coupling basically some known one-scalar field BPS models, namely φ 4 model, the ϕ 6 -like model, the sine-Gordon model, the E-model, and finally the inverse ζ 4I model. One of the most important advantages of this method of constructing multifields models is that it maintains the BPS solutions of the original one-field systems. However, they are not the only possible BPS solutions for the multi-field superpotential. In fact, in some cases we have been able to find analytically (or numerically) other BPS solutions for the resulting model. We have also studied in some details the linear stability of the BPS states for the resulting multiscalar superpotential. In general, these studies lead us with two very-well known exactly solvable quantum-mechanical problems, the Rosen-Morse II and the Scarf II potentials. For several choices of the potential parameters we have been able to perform analytically such analysis, and have found that they are stable with respect of small perturbations. However, in some cases the problem is somehow complicated and we have only been able to study in a qualitative and approximated way, with no full guarantee of the stability. Of course, such analysis could be improved by performing proper numerical simulations. Those investigations represent the next step in our studies on multiscalar field theories and will be done in future works.
There are several other interesting issues that can be addressed in next investigations from our results. For instance, a more complete numerical study of the solutions and their stability, specially two-solitons solutions, would provide a good scenario for investigating the behaviour during kink collisions. In addition, that kind of analysis also could bring some additional information that allow us identify possible quasi-integrable multi-scalar models [59,60]. In particular, we are interested in the two coupled sine-Gordon model obtained in section 6.3, which are slightly related to the FKZ models. We believe that a more detailed analysis would give some interesting connections between the two methods and probably help us answer some unsolved problems from both sides.
Finally, one more question of interest involves the investigations of possible supersymmetric generalization of the extension method. As it is well-known the interest on the study of supersymmetric kinks has a long history, and essentially concerns with the calculations of quantum corrections to the kink mass, and the central charge [44], [61]- [65]. Therefore, it will be interesting to construct general supersymmetric field theories by using the extension method, especially for the ones that possess intrinsically infinite number of degenerate vacua, as it is the case of sine-Gordon and the E-model. These issues are also currently under investigations.

Acknowledgements
Authors would like to thank to CAPES-Brazil for financial support. Authors are also grateful to the Directorate of Innovation and Research of the Federal University of Itajubá (DIP-UNIFEI) for partial financial support at the very initial stage of this project.

A Associated exactly solvable potentials
The very well-known exactly solvable Rosen-Morse II potential (or modified Poschl-Teller potential) can be written in the following form [53], where α > 0, and A and B are arbitrary real parameters. The bound states have the following eigenvalues, By imposing the stability condition, we find that In addition, the corresponding wave eigenfunctions are given by where P (α,β) k are the Jacobi polynomials, and .
Now, let us consider another very well-known exactly solvable potential, namely the Scarf II potential [53], where α, A, and B are real parameters. Its corresponding bound states possess energy eigenvalues given by Their associated eigenfunctions can be written as follows, ψ k (x) = i k (sech(αx)) s e −u arctan(sinh(αx)) P

B Calculation of the g-functions for three-fields systems
Here, we will present the explicit derivations of the g-function for the three-field model constructed in section 6. In principle they are arbitrary function constructed in a similar way as the superpotential, by using the deformation functions and the corresponding inverse functions. The specific form will come out of the following constraints which are basically consistency conditions for the existence of a well-defined continuous superpotential function given by the ansatz (6.11) -(6.13).