1 Introduction

Since there is no known fundamental principle requiring spacetime to be \((3 + 1)\)-dimensional [1, 2], it has been suggested that our observable Universe might be a \((3+1)\)-dimensional brane in a higher-dimensional space [3–5]. In most models, there are one or more flat 3-branes embedded discontinuously in the ambient geometry [6]. Moreover, ideas with two 3-branes provide a very elegant description of the large hierarchy between the scales of weak and gravitational forces [6, 7] and contain massless modes which reproduce Newtonian gravity at large distances on the brane [6]. In recent years, particle physics extra-dimensional theories beyond the standard model have become a standard part of the array of phenomenological models [8, 9]. Although there are still no experimental evidence supporting extra dimensions, due to various theoretical motivations, extra-dimensional models continue to be widely considered in the literature [10]. In this context, it would be quite useful to have a set of simple and sufficiently general rules which would allow one to test new models [3]. Most extra-dimensional models require the existence of scalar fields, for instance, to generate a domain wall which localizes matter fields [2]. The scalar fields also serve to stabilize the size of the compact extra dimensions [3, 4], and can also help modify the Randall–Sundrum warped-space [11, 12] to a smoothed-out version [5, 6], or to cut off the extra dimension at a singularity [7, 8]. In order to replace an infinitely thin brane with a thick one, a scalar field with soliton behavior is frequently invoked. The nonlinearity in the scalar field and in particular the existence of discrete vacua in the self-interaction of the scalar field lead to the appearance of a stable localized solution, which is a good motivation for building thick brane models [10]. A general method for determining the lowest energy configuration has been worked out in [13, 14].

Recently, braneworld models have also been considered in higher-order curvature gravity and in modified teleparallel gravity. For instance, five-dimensional modified teleparallel gravity was considered in a brane scenario, where analytic domain walls were found to have a double-kink solution in the aftermath of the torsion of spacetime [15]. Furthermore, this model was extended by using a first-order formalism to find analytical solutions for models that include a scalar field with standard and generalized dynamics. In addition to this, it was found that the brane splits, as a result of the deviation from the standard model by controlling specific parameters [16]. Modified gravity in five-dimensional spacetime has also been analyzed in the Palatini formalism. For instance, a thick Palatini f(R) brane described by an anti-de Sitter warped geometry with a single extra dimension of infinite extent, sourced by a real scalar field was studied in a perturbative scenario [17]. Besides, the model of a domain wall (thick brane) in a non-compact AdS space time with only one extra dimension was further analyzed in [18, 19]. The classical tests of general relativity in thick branes were also studied by studying the motion of test particles in a thick brane scenario and the impact of the brane thickness on the four-dimensional path of massless particles was explored in [20]. More specifically, by applying a confinement mechanism of massive tests particles in the domain wall, for instance, that simulates classically the trapping of the Dirac field in a domain wall, the influence of the brane thickness on the four-dimensional (4D) path of massless particles was analyzed. A generalized version of the Randall–Sundrum II model with different cosmological constants on each side of a brane were also discussed, where specific configurations of a scalar field and its stability as a replacing factor of the singular brane were considered [21]. Models of thick branes in non-compact five-dimensional bulk with different anti-de Sitter geometries on each side of the brane were explored [22, 23], and cosmological applications of soliton-like thick branes have also been studied [24]. On the other hand, the existence of brane solutions were also considered as a result of a real scalar field in the presence of five-dimensional f(R) gravity [25]. In addition to this, asymmetric thick braneworld scenarios were studied, by changing the superpotential of the scalar field [26].

It is widely practiced that in brane world scenarios \(Z_{2}\) symmetry is assumed [10, 27–29], which is originally motivated from the \(Z_{2}\) symmetry considered in M-theory [27]. Under this symmetry the bulk metric on the two sides of the brane should be the same [29]. Moreover, under such a symmetry the empty bulk on either sides of the brane have the same negative cosmological constant and as a result they are AdS [28]. Note that these conditions are satisfied in the Randall–Sandrum model. There are, however, brane models in which there is no \(Z_{2}\) symmetry and the bulk is different on both sides of the brane [27]. In the latter, the Friedmann equation for a positive brane tension situated between two bulk spacetimes that posses the same 5D cosmological constant, but which does not possess a \(Z_2\) symmetry of the metric itself was derived, and the possible effects of dropping the \(Z_2\) symmetry on the expansion of our Universe were examined. In some of these models, the cosmological constant differ on both sides of the brane [30], where the effects of including a Gauss–Bonnet combination of higher-order curvature invariants in the bulk action are taken into account. In fact, by considering braneworld scenarios including the Gauss–Bonnet term, it was found that the cosmological dynamics have the same form as those in Randall–Sundrum scenarios but with time-varying four-dimensional gravitational and cosmological constants [29]. Motivated by such a possibility, we consider in this paper several models, namely, the sine-Gordon (SG), \(\varphi ^{4}\) and \(\varphi ^{6}\) brane models which have broken \(Z_{2}\) symmetry in some cases. What is meant in the paper as \(Z_{2}\) symmetry is the symmetry with respect to the position of the brane (not the actual symmetry in the lagrangian). The brane position is shifted from \(z=0\) for some models. As a result, we find that for the \(\phi ^6\) and SG systems this symmetry is broken and the vacua on the two sides of the brane are not degenerate. In several cases, the \(Z_{2}\) symmetry can be restored by a proper choice of model parameters. The origin of symmetry breaking in our models resides in the fact that the modified scalar field potential may have non-degenerate vacua. These vacua determine the cosmological constant on both sides of the brane.

Relative to the stability issue, topological solitons are known for their non-singular structure and a natural localization mechanism which are highly stable. Zeldovich et al. [14, 31] suggested that the soliton of the \(\varphi ^4\) model is a reasonable source for the formation of domain walls. Vilenkin extended this idea to incorporate the general theory of relativity [14, 32–34]. Since there is a close similarity between domain walls in \(3+1\) dimensions and branes in \(4+1\) dimensions, it is natural to think that the soliton idea might have something to do with the existence and stability of branes. Motivated by this idea, we consider soliton models for thick branes and extend some of the existing works. In particular, in [35] the authors have studied the stability of the \(\varphi ^4\) kink brane model in five dimensions, as well as the \(\varphi ^3\) and the inverted \(\varphi ^4\) potential. Furthermore, the properties of fermions coupled to the sine-Gordon brane model were investigated. In [36], double kink-like solutions were considered and the stability of scalar, vector and tensor perturbations were discussed.

This paper is outlined in the following manner: in Sect. 2, we briefly review the thick brane formalism, by presenting the action and the field equations. Furthermore, we also study the geodesic equations along the fifth dimension, in order to explore the particle motion in the neighborhood of the brane. In Sect. 3, we present new soliton models and discuss their fundamental properties, for instance, by exploring the broken \(Z_2\)-symmetry character of the solutions and the confining effects of the scalar field on the brane. In Sect. 4, we analyze the stability of these brane models, where the metric and the scalar field are perturbed about the static brane and the resulting equations of motion are linearized in the proper gauge. Finally, in Sect. 5, we present our concluding remarks.

2 Thick brane formalism

We consider a thick brane, embedded in a five-dimensional (5D) bulk spacetime, modeled by the following action:

$$\begin{aligned} S=\int \mathrm{d}^{5}x\sqrt{|g^{(5)}|}\left[ \frac{1}{4} R[g^{(5)}]-\frac{1}{2}\partial _{A}\varphi \partial ^{A}\varphi -V(\varphi )\right] , \end{aligned}$$
(1)

where \(g^{(5)}\) is the metric and \(R[g^{(5)}]\) the scalar curvature in the bulk; \(\varphi \) is a dilaton field living on the bulk and \(V(\varphi )\) is a general potential energy. Note that we are using \(\kappa _{5}^{2}=8\pi G_{5}=2\).

The simplest line element of the brane, embedded in the 5D bulk spacetime, can be written as [14]

$$\begin{aligned} \mathrm{d}s^{2}_{5}= & {} g_{AB}\mathrm{d}x^{A}\mathrm{d}x^{B}\nonumber \\= & {} \mathrm{d}w^{2}+\mathrm{e}^{2A}(\mathrm{d}x^{2}+\mathrm{d}y^{2}+\mathrm{d}z^{2}-\mathrm{d}t^{2}), \end{aligned}$$
(2)

where A is the warp factor which depends only on the five-dimensional (5D) coordinate w. For the scalar field \(\varphi \) with the potential \(V(\varphi )\), the 5D energy-momentum tensor is given by

$$\begin{aligned} T_{AB}=\partial _{A}\varphi \partial _{B}\varphi -g_{AB}\left[ \frac{1}{2}\partial _{C}\varphi \;\partial ^{C}\varphi +V(\varphi )\right] , \end{aligned}$$
(3)

where \(g_{AB}\) and \(\varphi \) depend only on w.

The 5D gravitational and scalar field equations take the following forms:

$$\begin{aligned}&3A^{\prime \prime }+6{A^{\prime }}^{2}=-\kappa _{5}^{2}\mathrm{e}^{-2A}T_{00}=-\kappa _{5}^{2}\left[ \frac{1}{2}{\varphi ^{\prime }}^{2}+V(\varphi )\right] ,\end{aligned}$$
(4)
$$\begin{aligned}&6{A^{\prime }}^{2}=\kappa _{5}^{2}T_{44}=\kappa _{5}^{2}\left[ \frac{1}{2}{\varphi ^{\prime }}^{2}-V(\varphi )\right] , \end{aligned}$$
(5)
$$\begin{aligned}&\varphi ^{\prime \prime }+4A^{\prime }\varphi ^{\prime }=\frac{\mathrm{d} V(\varphi )}{\mathrm{d}\varphi }, \end{aligned}$$
(6)

respectively, where the prime denotes a derivative with respect to w. In order to obtain a first-order equation, we introduce an auxiliary function W according to [6, 37–40], which requires

$$\begin{aligned} A^{\prime }= & {} -\frac{1}{3}W(\varphi ), \end{aligned}$$
(7)
$$\begin{aligned} \varphi ^{\prime }= & {} \frac{1}{2}\frac{\partial W(\varphi )}{\partial \varphi }, \end{aligned}$$
(8)

while \(V(\varphi )\) takes the following form [6, 37–40]:

$$\begin{aligned} V(\varphi )=\frac{1}{8}\left( \frac{\partial W(\varphi )}{\partial \varphi }\right) ^{2}-\frac{1}{3}W(\varphi )^{2}. \end{aligned}$$
(9)

The \(T_{00}\) distribution on the bulk, which will be analyzed in detail below, is given by [10]

$$\begin{aligned} T_{00}=\mathrm{e}^{2A}\left[ \frac{1}{2}\left( \frac{\partial \varphi }{\partial w}\right) ^{2}+V(\varphi )\right] . \end{aligned}$$
(10)

It can also be shown that, for models with an infinitely thin brane and Dirac delta distributions, the energy density is equal to the cosmological constant of the bulk plus the energy density on the brane, i.e., \(\varepsilon =\Lambda _{5}^{\pm }+k\delta (w)\). Moreover, it may be instructive to calculate the geodesic equation along the fifth dimension in a thick brane, in order to investigate the particle motion near the brane [41]. As mentioned before, the thick brane models considered in this paper do not have Dirac delta singularities which enable easier direct calculations. To this end, we start with the geodesic equation:

$$\begin{aligned} \frac{\mathrm{d}^{2}x^{0}}{\mathrm{d}\tau ^{2}}+\Gamma ^{0}_{AB}\frac{\mathrm{d}x^{A}}{\mathrm{d}\tau }\frac{\mathrm{d}x^{B}}{\mathrm{d}\tau }= & {} 0 \Rightarrow \frac{\mathrm{d}}{\mathrm{d}\tau }\left( -2\mathrm{e}^{2A}\dot{t}\right) =0, \nonumber \\ \frac{\mathrm{d}^{2}x^{4}}{\mathrm{d}\tau ^{2}}+\Gamma ^{4}_{AB}\frac{\mathrm{d}x^{A}}{\mathrm{d}\tau }\frac{\mathrm{d}x^{B}}{\mathrm{d}\tau }= & {} 0 \Rightarrow \ddot{w}+A'\mathrm{e}^{2A}\dot{t}^{2}=0 , \end{aligned}$$
(11)

which leads to

$$\begin{aligned} \ddot{w}+c_{1}^{2}f(w)=0 , \end{aligned}$$
(12)

where \(c_{1}\) is a constant of integration and the function f(w) is defined as

$$\begin{aligned} f(w)=A'(w)\mathrm{e}^{-2A(w)}. \end{aligned}$$
(13)
Fig. 1
figure 1

Soliton solutions as a function of the five-dimensional coordinate w for the following models: a SG for \(b=1\), b \(\varphi ^4\) for \(\alpha =1\) and c \(\varphi ^6\) for \(\alpha =1\) systems. Dashed, dotted-dashed, and continuous curves correspond to solitons with decreasing brane thickness. In the limit of an infinite \(a/\beta \) parameter, the soliton approaches the step function

Full size image

Equation (12) is a second-order differential equation for w and its solution depends critically on whether f(w) / w is positive or negative. For positive values of f(w) / w one obtains periodic (exponential) solutions, respectively. Note that the periodic (negative) motion indicates particle confinement near the brane, while the exponential solutions implies that the reference point is unstable. However, this may point to the possibility that \(w=0\) is different from the localization of the brane. In a periodic situation, by introducing a new quantity \(F(w)=c_{1}^{2}A'(w)\mathrm{e}^{-2A(w)}\), one can write the geodesic equation in the following form:

$$\begin{aligned} \ddot{w}+F(w)=0. \end{aligned}$$
(14)

The equilibrium point \(w_{0}\) satisfies \(F(w_{0})\)=0. On the other hand, by expanding F(w) around \(w_{0}\), we have

$$\begin{aligned} F(w)=F(w_{0})+F'(w_{0})(w-w_{0})+\cdots , \end{aligned}$$
(15)

and the geodesic equation leads to

$$\begin{aligned} \ddot{w}+F'(w_{0})(w-w_{0})=0. \end{aligned}$$
(16)

Taking into account the change of variable \(\tilde{w}=w-w_{0}\), the geodesic equation reduces to

$$\begin{aligned} \ddot{\tilde{w}}+F'(w_{0})\tilde{w}=0, \quad \mathrm{or} \quad \ddot{\tilde{w}}+\Omega ^{2}\tilde{w}=0 , \end{aligned}$$
(17)

where \(\Omega =\sqrt{F'(w_{0})}\), provided that \(F'(w_{0})\ge 0\).

It is essential to emphasize that in the RS-II brane model, the KK zero mode corresponds to the massless graviton and the massive modes form a continuum which results in a small correction to Newtonian gravity at large distances [42, 43]. Free particles are only affected by the gravitational field and not directly by the scalar field. Any field or particle which has a direct coupling with the scalar field will be further affected by an extra force from the scalar field. Moreover, even an exponential potential like \(\mathrm{e}^{-\alpha w^{2}}\) reduces to a harmonic potential for small amplitude oscillations (\(\mathrm{e}^{-\alpha w^{2}}\approx 1-\alpha w^{2}+O(w^{4})\)).

In the following section, we explore several models for thick branes and we will employ the linearized geodesic Eq. (17) for each model.

3 Soliton models for the brane

3.1 Sine-Gordon-based models

The sine-Gordon (SG) model is a well-known integrable model which has found interesting applications in various disciplines [44, 45]. In fact, single and multiple (topological) soliton solutions of this system are found analytically through different mathematical methods [44, 45]. The self-interaction potential for this model reads

$$\begin{aligned} \tilde{V}(\varphi )=\frac{a}{b}\left[ 1-cos(b\varphi )\right] , \end{aligned}$$
(18)

where a and b are free parameters of the model. When considered as a brane potential, however, this potential should be modified to become consistent with the Einstein equations.

The SG system has the following exact static kink solution [44]:

$$\begin{aligned} \varphi (w)=\frac{4}{b}\arctan \left( \mathrm{e}^{\sqrt{ab}w}\right) , \end{aligned}$$
(19)

which is plotted in Fig. 1a, for various values of parameters a and b, which correspond to branes with different thicknesses. The formalism of our investigation is to keep the soliton solution of the flat space nonlinear equation and modify the scalar field potential in such a way that the soliton solution remains a solution of the full gravitating system. This is why the soliton solution remains the same. The form of the potential, however, changes accordingly. Taking into account the scalar field given by Eq. (19), and plugging it into the field equations, we obtain the following quantities:

$$\begin{aligned} \begin{aligned}&W(\varphi )=-\frac{16\sqrt{ab}}{b^{2}}\left[ \cos ^{2}\left( \frac{b\varphi }{4}\right) \right] ,\\&A= - \frac{8}{3b^{2}}\ln \left( \frac{1+\mathrm{e}^{2\sqrt{ab}w}}{\mathrm{e}^{2\sqrt{ab}w}}\right) ,\\&\exp (2A)=\left( \frac{1+\mathrm{e}^{2\sqrt{ab}w}}{\mathrm{e}^{2\sqrt{ab}w}}\right) ^{-\frac{16}{3b^{2}}}, \end{aligned} \end{aligned}$$
(20)

where the warp factor is plotted in Fig. 2a. The corresponding modified potential for this model is given by

$$\begin{aligned} V(\varphi )=\frac{2a}{b}\sin ^{2}\left( \frac{b\varphi }{2}\right) -\frac{64a}{3b^{3}}\left[ 1+\cos \left( \frac{b\varphi }{2}\right) \right] ^{2}, \end{aligned}$$
(21)

which is depicted in Fig. 3a. Notice that this potential has two series of non-degenerate vacua, as in the double sine-Gordon (DSG) system potential [45]. However, in the limit of \(b\gg a\) these vacua tend to the same value (become degenerate), such as the potentials used in [10, 39].

Fig. 2
figure 2

The plots depict the warp factor as a function of the fifth dimension for the a SG with \(a=100\) and \(b=1\), b \(\varphi ^4\) with \(\alpha =1\) and \(\beta =100\) and c \(\varphi ^6\) with \(\alpha =1\) and \(\beta =100\) systems, respectively. For the SG and \(\varphi ^6\) systems there is an asymmetry between the two sides of the brane and the \(Z_2\) symmetry is broken

Full size image
Fig. 3
figure 3

The plots depict the modified soliton potential as a function of the scalar field for a SG with \(a=b=10\), b \(\varphi ^4\) with \(\alpha =\beta =1\) and c \(\varphi ^6\) with \(\alpha =\beta =1\) systems. The potential of the SG and \(\varphi ^6\) systems have non-degenerate vacua. In contrast, the \(\varphi ^4\) potential has degenerate vacua and this leads to a stable, topological solitonic brane. For \(\varphi ^{4}\) and \(\varphi ^{6}\) systems the number of extrema in the figure is only a result of the range of the plot. In a wider range plot, other extrema appear; these may render the system globally unstable. We have only claimed local stability

Full size image
Fig. 4
figure 4

\(T^{0}_{0}\) as a function of the fifth dimension coordinate w for the following systems: a SG with \(b=1\), b \(\varphi ^4\) with \(\alpha =1\) and c \(\varphi ^6\) with \(\alpha =1\). The dashed, dotted-dashed and continuous curves correspond to increasingly thin branes. The brane becomes infinitely thin (delta function) as the parameter \(a/\beta \) tends to infinity

Full size image

This system leads to the following \(T_{00}\), from which the energy density can be obtained:

$$\begin{aligned} T_{00}=\left( \frac{1+\mathrm{e}^{2\sqrt{ab}w}}{\mathrm{e}^{2\sqrt{ab}w}}\right) ^{\left( \frac{-16}{3b^{2}}\right) }\left[ \frac{16}{3}\frac{a\left( 3b^{2}\mathrm{e}^{2\sqrt{ab}w}-16\right) }{b^{3}\left( 1+\mathrm{e}^{2\sqrt{ab}w}\right) ^{2}}\right] , \end{aligned}$$
(22)

which is plotted in Fig. 4a. The energy-momentum tensor is calculated according to

$$\begin{aligned} T_{\mu \nu }=\partial _{\mu }\varphi \partial _{\nu }\varphi -g_{\mu \nu }\left[ \frac{1}{2}\left( \partial _{\alpha }\varphi \right) \left( \partial ^{\alpha }\varphi \right) +V(\varphi )\right] , \end{aligned}$$
(23)

where the modified potential is used for V. The metric and the energy-momentum tensor are checked to satisfy the Einstein equations. Moreover, note that the potential for any soliton model may be shifted by a constant, without affecting the soliton solutions. The minimum value of the potential (i.e., the classical vacuum), if negative, leads to a negative energy density. This can be avoided by adding a positive constant to the self-interaction potential, when the scalar field is not coupled to gravity (i.e., in flat spacetime). In curved spacetime, however, this constant is non-trivial and plays the role of a cosmological constant of the bulk. Of course, a positive cosmological constant violates the strong energy condition. Note that the energy density is localized at the brane and the thickness of the latter is given by

$$\begin{aligned} \triangle =\frac{1}{2\sqrt{ab}}. \end{aligned}$$
(24)

In this model, the Ricci and Kretschmann scalars are given by

$$\begin{aligned} R=\frac{256}{9}\frac{a\left[ -20+3\left( \mathrm{e}^{2\sqrt{ab}w}\right) b^{2}\right] }{b^{3}\left[ 1+\left( \mathrm{e}^{2\sqrt{ab}w}\right) \right] ^{2}} \end{aligned}$$
(25)

and

$$\begin{aligned} K=\frac{16384}{81}\frac{a^{2}\left[ 160-48\left( \mathrm{e}^{2\sqrt{ab}w}\right) b^{2}{+}9\left( \mathrm{e}^{4\sqrt{ab}w}\right) b^{4}\right] }{b^{6}\left[ 1{+}\left( \mathrm{e}^{2\sqrt{ab}w}\right) \right] ^{4}},\nonumber \\ \end{aligned}$$
(26)

respectively. It can be seen that there is no singularity in the Ricci scalar and/or Kretschmann scalar. In the limits of \(w\rightarrow \pm \infty \), the Ricci scalar becomes

$$\begin{aligned} \lim _{w\rightarrow +\infty } R= & {} 0, \end{aligned}$$
(27)
$$\begin{aligned} \lim _{w\rightarrow -\infty } R= & {} -\frac{5120}{9}\frac{a}{b^{3}}, \end{aligned}$$
(28)

and the limit of \(w\rightarrow 0\) yields

$$\begin{aligned} \lim _{w\rightarrow 0} R=\frac{64}{9}\frac{(-20+3b^{2})a}{b^{3}}. \end{aligned}$$
(29)
Fig. 5
figure 5

Ricci scalar as a function of the fifth dimension coordinate w for the following systems: a SG with \(b=1\), b \(\varphi ^4\) with \(\alpha =1\) and c \(\varphi ^6\) with \(\alpha =1\). The Ricci scalar is different on the two sides of the brane, due to the effect of the warp factor

Full size image
Fig. 6
figure 6

The plots depict the Einstein tensor component \(G_{0}^{0}\) for the a SG with \(a=b=1\), b \(\varphi ^4\) with \(\alpha =\beta =1\) and c \(\varphi ^6\) with \(\alpha =\beta =1\) systems. Note that this quantity approaches different constant values for the SG and \(\varphi ^{6}\) systems, while the \(\varphi ^{4}\) system is \(Z_{2}\)-symmetric

Full size image

Moreover, the mixed Einstein tensor components are given by

$$\begin{aligned} G^{0}_{0}= & {} G^{1}_{1}=G^{2}_{2}=G^{3}_{3} =-\frac{32}{3}\frac{a\left( -16+3\mathrm{e}^{2\sqrt{ab}w}b^{2}\right) }{b^{3}\left( 1+\mathrm{e}^{2\sqrt{ab}w}\right) ^{2}},\nonumber \\ \end{aligned}$$
(30)
$$\begin{aligned} G^{4}_{4}= & {} \frac{512}{3}\frac{a}{b^{3}\left( 1+\mathrm{e}^{2\sqrt{ab}w}\right) ^{2}}, \end{aligned}$$
(31)

respectively. Note that all the components of the Einstein tensor in the limits \(w\rightarrow \pm \infty \) become

$$\begin{aligned} \lim _{w\rightarrow +\infty } G^{A}_{B}= & {} 0,\end{aligned}$$
(32)
$$\begin{aligned} \lim _{w\rightarrow -\infty } G^{A}_{B}= & {} \frac{512}{3}\frac{a}{b^{3}}. \end{aligned}$$
(33)

However, in the limit of \(w\rightarrow 0\) (for \(\mu =\nu =0,1,2,3\)) the Einstein tensor is given by

$$\begin{aligned} \lim _{w\rightarrow 0}G^{\mu }_{\nu }=\frac{8}{3}\frac{(16-3b^{2})a}{b^{3}}\,\delta ^{\mu }_{\nu }, \end{aligned}$$
(34)

and one can interpret it as the cosmological constant on the brane, i.e., \(G^{i}_{j}\propto \Lambda \delta ^{i}_{j}\), with \(\Lambda =\frac{8}{3}\frac{(16-3b^{2})a}{b^{3}}\).

These results can be interpreted in that we have a broken \(Z_{2}\)-symmetry in the bulk, as the two sides of the brane differ completely. On the right (\(w\rightarrow +\infty \)), the Einstein tensor and consequently the cosmological constant of the bulk vanish, so the bulk is asymptotically Minkowski. However, on the other side of the brane, these quantities are nonzero and equal to the constant value \(512 a/(3b^{3})\), and as a result the bulk would be de Sitter. The Ricci scalar and the Einstein tensor component \(G^{0}_{0}\) are plotted Figs. 5a and 6a, respectively.

Furthermore, by calculating the field equations,

$$\begin{aligned} G_{AB}=\kappa _{5}^{2}T_{AB} \end{aligned}$$
(35)

one verifies that \(\kappa _{5}^{2}=2\), which is consistent with the usual normalization notation [37, 39, 46]. As pointed out in DeWolfe et al. [6], in the stiff limit where \(ab\rightarrow \infty \), the wall reduces to the step function and the energy density approaches a \(\delta \)-function (see Figs. 1a, 4a).

The local confining gravitational field of the brane is best observed by looking at the geodesic equation of a test particle moving only in the direction of the extra dimension, as explained in the previous section:

$$\begin{aligned} \ddot{w}+c_{1}^{2}\frac{8}{9}\frac{32^{\frac{1}{b^{2}}}2^{\frac{1}{3b^{2}}}a\left( 3b^{2}+16\right) }{b^{3}}w\approx c_{1}^{2}\frac{8}{3}\frac{32^{\frac{1}{b^{2}}}2^{\frac{1}{3b^{2}}}\sqrt{ab}}{b^{2}}. \end{aligned}$$
(36)

This proves the confining effect of the scalar field, and \(c_{1}\) is an integration constant (see Eq. (12)). For small amplitude oscillations, the relativistic motion reduces to a classical motion in a Newtonian classical potential. Since we have considered the potential up to second order, relativistic effects can be ignored and the corresponding quantum energy levels are therefore those of a non-relativistic quantum particle. So, if interpreted as a quantum oscillator, one can assign an energy to each quantum state given by \(E_n=(n+1/2)\hbar \omega \), where

$$\begin{aligned} \hbar \omega =\Omega =\sqrt{F'(w_{0})}\approx c_{1}\sqrt{\frac{8}{9}\frac{32^{\frac{1}{b^{2}}}2^{\frac{1}{3b^{2}}}a(3b^{2}+16)}{b^{3}}}. \end{aligned}$$
(37)

However, an important point is in order. In the following, we explore the confining effect of the brane, only up to second order in the potential. Even if the classical test particle is confined up to this order, large amplitude oscillations will involve nonlinear effects and this might exploit confinement. Quantum mechanically, the full nonlinear potential might lead to tunneling and thus de-confinement. It is well known that massive KK modes are not confined to the brane.

3.2 \(\varphi ^{4}\)-Based model

The \(\varphi ^4\) models are well known for having simple soliton-like solutions, although it is not strictly integrable like the SG system. This model is also the central ingredient in the Goldstone and Higgs mechanisms. Spontaneous breakdown of the \(Z_2\) symmetry in the complex version of the \(\varphi ^4\) model leads to the appearance of the Goldstone mode and once coupled with a gauge field, it causes the gauge boson to acquire mass [47]. This model is therefore frequently used for building thick branes.

For the \(\varphi ^{4}\)-based model, we have [10]

$$\begin{aligned} \tilde{V}(\varphi )=\frac{\beta ^{2}}{2\alpha ^{2}}(\varphi ^{2}-\alpha ^{2})^{2}, \end{aligned}$$
(38)

where \(\alpha \) and \(\beta \) are constants. The kink solution reads

$$\begin{aligned} \varphi (w)=\alpha \tanh (\beta w), \end{aligned}$$
(39)

which is depicted in Fig. 1b. Proceeding in a similar manner as in the previous case, we have the following solutions:

$$\begin{aligned}&W(\varphi )=\frac{2}{3}\frac{\beta \varphi \left( 3\alpha ^{2}-\varphi ^{2}\right) }{\alpha },\end{aligned}$$
(40)
$$\begin{aligned}&A=-\frac{4}{9}\alpha ^{2}\ln \left[ \cosh (\beta w)\right] +\frac{1}{9}\alpha ^{2}\frac{1}{\cosh ^{2}(\beta w)},\end{aligned}$$
(41)
$$\begin{aligned}&\exp (2A)=\left[ \cosh (\beta w)\right] ^{-\frac{8}{9}\alpha ^{2}}\exp \left[ \frac{2}{9}\alpha ^{2}\frac{1}{\cosh ^{2}(\beta w)}\right] , \end{aligned}$$
(42)

respectively, and the potential is obtained:

$$\begin{aligned} V(\varphi )=\frac{1}{2}\alpha ^{2}\beta ^{2}\left( 1-\frac{\varphi ^{2}}{\alpha ^{2}}\right) ^{2}-\frac{4}{27}\varphi ^{2}\alpha ^{2}\beta ^{2}\left( 3-\frac{\varphi ^{2}}{\alpha ^{2}}\right) ^{2}, \end{aligned}$$
(43)

which is plotted in Fig. 3b. It is seen that while \(\tilde{V}(\varphi )\) was \(O(\varphi ^{4})\), \(V(\varphi )\) is \(O(\varphi ^{6})\).

The corresponding \(T_{00}\) is given by

$$\begin{aligned} T_{00}= & {} -\frac{1}{27}\frac{\alpha ^{2}\beta ^{2}\exp \left[ -\frac{2}{9}\alpha ^{2}\tanh ^{2}(\beta w)\right] }{\cosh ^{\left( 6-\frac{8}{9}\alpha ^{2}\right) }(\beta w)} \nonumber \\&[ -4\alpha ^{2} -27\cosh ^{2}(\beta w) \nonumber \\&+16\alpha ^{2}\cosh ^{6}(\beta w)-12\alpha ^{2}\cosh ^{2}(\beta w)], \end{aligned}$$
(44)

which is depicted in Fig. 4b, where one verifies that the energy density is localized at the brane. The brane thickness becomes \(\triangle =\beta ^{-1}\).

Moreover, one can show that the Ricci and Kretschmann scalars are given by

$$\begin{aligned} R= & {} -\frac{16}{81}\frac{\alpha ^{2}\beta ^{2}}{\cosh ^{6}(\beta w)} [-15\alpha ^{2}\cosh ^{2}(\beta w)-5\alpha ^{2}\nonumber \\&+20\alpha ^{2}\cosh ^{6}(\beta w)-27\cosh ^{2}(\beta w)], \end{aligned}$$
(45)

which is depicted in Fig. 5b, and

$$\begin{aligned} K= & {} \frac{64}{6561}\frac{\alpha ^{4}\beta ^{4}}{\cosh ^{12}(\beta w)}[90\alpha ^{4}\cosh ^{4}(\beta w) -80\alpha ^{4}\cosh ^{6}(\beta w)\nonumber \\&-240\alpha ^{4}\cosh ^{8}(\beta w)+160\alpha ^{4}\cosh ^{12}(\beta w)\nonumber \\&+60\alpha ^{4}\cosh ^{2}(\beta w) +10\alpha ^{4}+324\alpha ^{2}\cosh ^{4}(\beta w)\nonumber \\&+108\alpha ^{2}\cosh ^{2}(\beta w)-432\alpha ^{2}\cosh ^{8}(\beta w) \nonumber \\&+729\cosh ^{4}(\beta w)], \end{aligned}$$
(46)

respectively. The mixed Einstein tensor components take the following form:

$$\begin{aligned} G^{0}_{0}= & {} G^{1}_{1}=G^{2}_{2}= G^{3}_{3} =\frac{2}{27} \frac{\alpha ^{2}\beta ^{2}}{\cosh ^{6}(\beta w)} \nonumber \\&[-4\alpha ^{2}-12\alpha ^{2}\cosh ^{2}(\beta w)\nonumber \\&+16\alpha ^{2}\cosh ^{6}(\beta w)-27\cosh ^{2}(\beta w)], \end{aligned}$$
(47)
$$\begin{aligned} G^{4}_{4}= & {} \frac{8}{27}\frac{\alpha ^{4}\beta ^{2}[-3\cosh ^{2}(\beta w)-1+4\cosh ^{6}(\beta w)]}{\cosh ^{6}(\beta w)} \end{aligned}$$
(48)

where the \(G^0_0\) component is depicted in Fig. 6b. As for the previous SG model, we determine the limits \(w\rightarrow \pm \infty \) for all the components of the Einstein tensor components, which are given by

$$\begin{aligned} \lim _{w\rightarrow +\infty } G^{A}_{B}= & {} \frac{32}{27}\alpha ^{4}\beta ^{2},\end{aligned}$$
(49)
$$\begin{aligned} \lim _{w\rightarrow -\infty }G^{A}_{B}= & {} \frac{32}{27}\alpha ^{4}\beta ^{2}, \end{aligned}$$
(50)

and in the limit of \(w\rightarrow 0\) the Einstein tensor components (for \(\mu =\nu =0,1,2,3\)) take the form

$$\begin{aligned} \lim _{w\rightarrow 0} G^{\mu }_{\nu }=-2\alpha ^{2}\beta ^{2}\delta ^{\mu }_{\nu }. \end{aligned}$$
(51)

Note that in this model the cosmological constant on the brane \(\Lambda \) would be \(-2\alpha ^{2}\beta ^{2}\). Taking into account all of the above considerations, we verify that the Einstein equations are given consistently by \(G_{AB}=\kappa _{5}^{2}T_{AB}\) where \(\kappa _{5}^{2}=2\). For the geodesic equation for a test particle moving in the direction of the fifth dimension one obtains

$$\begin{aligned} \ddot{w}+c_{1}^{2}\frac{2}{3}\alpha ^{2}\beta ^{2}w=0, \end{aligned}$$
(52)

which corresponds to a linearized quantum mode of energy \(\hbar \omega =\Omega =\sqrt{F'(w_{0})}=\sqrt{\frac{2}{3}}c_{1}\alpha \beta \).

It is seen that the resulting potential for the \(\varphi ^{4}\) model is an odd function with respect to w, while the A function is even. This property is not verified in the SG and \(\varphi ^{6}\) models, where the latter is discussed below. This results in an asymmetry in the corresponding properties (such as the energy density). However, by an appropriate selection of the model parameters, one can restore the \(Z_2\) symmetry in the \(\varphi ^4\) and \(\varphi ^6\) cases, which is commonly used in brane models.

3.3 \(\varphi ^{6}\)-Based model

For this model, we have the following potential:

$$\begin{aligned} \tilde{V}(\varphi )=\frac{\beta ^{2}}{4\alpha ^{2}}\varphi ^{2}(\varphi ^{2}-\alpha ^{2})^{2}, \end{aligned}$$
(53)

and as a result, the kink solution is given by [48]

$$\begin{aligned} \varphi (w)=\frac{\alpha }{\sqrt{1+\mathrm{e}^{(-\sqrt{2}\alpha \beta w)}}}, \end{aligned}$$
(54)

where \(\alpha \) and \(\beta \) are constant (as in the \(\varphi ^{4}\) model). The kink solution is depicted in Fig. 1c. For this model the potential W, and the warp factor A are given by the following expressions:

$$\begin{aligned}&W(\varphi )=\frac{\sqrt{2}}{4}\frac{\beta \varphi ^{2}(2\alpha ^{2}-\varphi ^{2})}{\alpha },\end{aligned}$$
(55)
$$\begin{aligned}&A=-\frac{\alpha ^{2}}{12}\left[ \frac{1}{1+\mathrm{e}^{-\sqrt{2}\alpha \beta w}} +\ln \left( \frac{1+\mathrm{e}^{-\sqrt{2}\alpha \beta w}}{\mathrm{e}^{-\sqrt{2}\alpha \beta w}}\right) \right] ,\end{aligned}$$
(56)
$$\begin{aligned}&\exp (2A)=\left( \frac{\mathrm{e}^{-w\sqrt{2}\beta \alpha }}{1+\mathrm{e}^{-w\sqrt{2}\beta \alpha }}\right) ^{\frac{\alpha ^{2}}{6}}\mathrm{e}^{-\frac{1}{6}\frac{\alpha ^{2}}{1+\mathrm{e}^{-w\sqrt{2}\beta \alpha }}}, \end{aligned}$$
(57)

respectively. Thus, for the \(\varphi ^{6}\) system the self-interaction potential takes the form

$$\begin{aligned} V(\varphi )=\frac{1}{4}\frac{\beta ^{2}\varphi ^{2}(\alpha ^{2}-\varphi ^{2})^{2}}{\alpha ^{2}}-\frac{1}{24}\frac{\beta ^{2}\varphi ^{4}(2\alpha ^{2}-\varphi ^{2})^{2}}{\alpha ^{2}}, \end{aligned}$$
(58)

which is depicted in Fig. 3c. Note that \(V(\varphi )\) has been raised to \(O(\varphi ^{8})\). The energy-momentum tensor component \(T_{00}\) is given by

$$\begin{aligned} T_{00}= & {} -\frac{1}{24}\left( \frac{\mathrm{e}^{-w\sqrt{2}\beta \alpha }}{1+\mathrm{e}^{-w\sqrt{2}\beta \alpha }}\right) ^{\frac{\alpha ^{2}}{6}} \frac{\mathrm{e}^{-\frac{1}{6}\frac{\alpha ^{2}}{1+\mathrm{e}^{-w\sqrt{2}\beta \alpha }}}}{\left( 1+\mathrm{e}^{-w\sqrt{2}\beta \alpha }\right) ^{4}} \nonumber \\&\times [\alpha ^{4}\beta ^{2}(-12\mathrm{e}^{-2w\sqrt{2}\beta \alpha } \nonumber \\&-12\mathrm{e}^{-3w\sqrt{2}\beta \alpha }+4\alpha ^{2}\mathrm{e}^{-w\sqrt{2}\beta \alpha }+\alpha ^{2}\nonumber \\&+\,4\alpha ^{2}\mathrm{e}^{-2w\sqrt{2}\beta \alpha })], \end{aligned}$$
(59)

which is plotted in Fig. 4c. The thickness of this brane is given by \(\triangle =(\sqrt{2}\alpha \beta )^{-1}\).

For this system the Ricci scalar is given by

$$\begin{aligned} R= & {} -\frac{1}{18} \frac{\alpha ^{4}\beta ^{2}\mathrm{e}^{w\sqrt{2}\beta \alpha }}{\left( \mathrm{e}^{w\sqrt{2}\beta \alpha }+1\right) ^{4}} \Big [20\alpha ^{2}\mathrm{e}^{w\sqrt{2}\beta \alpha }-48\mathrm{e}^{w\sqrt{2}\beta \alpha }\nonumber \\&+\, 5\alpha ^{2}\left( \mathrm{e}^{3w\sqrt{2}\beta \alpha }\right) +20\alpha ^{2}\left( \mathrm{e}^{2w\sqrt{2}\beta \alpha }\right) -48 \Big ], \end{aligned}$$
(60)

which is plotted in Fig. 5c, and the Kretschmann scalar by

$$\begin{aligned} K= & {} \frac{1}{648}\alpha ^{8}\beta ^{4}\frac{\mathrm{e}^{2w\sqrt{2}\beta \alpha }}{\left( \mathrm{e}^{w\sqrt{2}\beta \alpha }+1\right) ^{8}} \Bigg [1152+5\alpha ^{4}\left( \mathrm{e}^{6w\sqrt{2}\beta \alpha }\right) \nonumber \\&+\,40\alpha ^{4}\left( \mathrm{e}^{5w\sqrt{2}\beta \alpha }\right) ^{5}+120\alpha ^{4}\left( \mathrm{e}^{4w\sqrt{2}\beta \alpha }\right) \nonumber \\&+\,160\alpha ^{4}\left( \mathrm{e}^{3w\sqrt{2}\beta \alpha }\right) +80\alpha ^{4}\left( \mathrm{e}^{2w\sqrt{2}\beta \alpha }\right) \nonumber \\&-\,480\alpha ^{2}\left( \mathrm{e}^{3w\sqrt{2}\beta \alpha }\right) +1152\left( \mathrm{e}^{2w\sqrt{2}\beta \alpha }\right) \nonumber \\&-\,768\alpha ^{2}\left( \mathrm{e}^{2w\sqrt{2}\beta \alpha }\right) -384\alpha ^{2}\mathrm{e}^{w\sqrt{2}\beta \alpha }\nonumber \\&-\,96\alpha ^{2}\left( \mathrm{e}^{4w\sqrt{2}\beta \alpha }\right) +2304\mathrm{e}^{w\sqrt{2}\beta \alpha }\Big ]. \end{aligned}$$
(61)

The mixed Einstein tensor components are given by

$$\begin{aligned} G^{0}_{0}= & {} G^{1}_{1}=G^{2}_{2}=G^{3}_{3} =\frac{\alpha ^{4}\beta ^{2}\mathrm{e}^{w\sqrt{2}\beta \alpha }}{12\left( \mathrm{e}^{w\sqrt{2}\beta \alpha }+1\right) ^{4}} \Big (\alpha ^{2}\mathrm{e}^{3w\sqrt{2}\beta \alpha } \nonumber \\&+\,4\alpha ^{2}\mathrm{e}^{2w\sqrt{2}\beta \alpha }+4\alpha ^{2}\mathrm{e}^{w\sqrt{2}\beta \alpha }-12\mathrm{e}^{w\sqrt{2}\beta \alpha }-12\Big ),\nonumber \\ \end{aligned}$$
(62)
$$\begin{aligned} G^{4}_{4}= & {} \frac{1}{12}\frac{\alpha ^{6}\beta ^{2}\mathrm{e}^{2w\sqrt{2}\beta \alpha }\left( \mathrm{e}^{2w\sqrt{2}\beta \alpha }+4\mathrm{e}^{w\sqrt{2}\beta \alpha }+4\right) }{\left( \mathrm{e}^{w\sqrt{2}\beta \alpha }+1\right) ^{4}}, \end{aligned}$$
(63)

where the \(G^0_0\) component is depicted in Fig. 6c. These tensor components reduce to the following in the limit of \(w\longrightarrow \pm \infty \):

$$\begin{aligned} \lim _{w\rightarrow +\infty } G^{A}_{B}= & {} \frac{1}{12}\alpha ^{6}\beta ^{2},\nonumber \\ \lim _{w\rightarrow -\infty }G^{A}_{B}= & {} 0, \end{aligned}$$
(64)

and in the limit of \(w\rightarrow 0\), the Einstein tensor is (for \(\mu =\nu =0,1,2,3\))

$$\begin{aligned} \lim _{w\rightarrow 0} G^{\mu }_{\nu }= \frac{1}{8}\alpha ^{4}\beta ^{2} \left( \frac{3}{8}\alpha ^{2}-1\right) \, \delta ^{\mu }_{\nu }. \end{aligned}$$
(65)

As in the previous cases, the Einstein equation in the bulk \(G_{AB}=\kappa _{5}^{2}T_{AB}\) is found self-consistently with \(\kappa _{5}^{2}=2\), and which leads to the following geodesic equation:

$$\begin{aligned}&\ddot{w}+c_{1}^{2}\frac{1}{192}\exp \left[ \frac{\alpha ^{2}}{12}\left( 2\ln 2+1\right) \right] \alpha ^{4}\beta ^{2}\left( 8+3\alpha ^{2}\right) w \nonumber \\&\quad \approx \frac{1}{16}\exp \left[ \frac{1}{12}\alpha ^{2}\left( 2\ln 2+1\right) \right] \alpha ^{3}\beta \sqrt{2}. \end{aligned}$$
(66)

The quantum mode energy is thus given by

$$\begin{aligned} \hbar \omega= & {} \Omega =\sqrt{F'(w_{0})}\nonumber \\\approx & {} c_{1}\frac{\alpha ^{2}\beta }{24}\sqrt{3}\sqrt{\exp \left( \frac{\alpha ^{2}}{12}\left( 2\ln (2)+1\right) \right) \left( 8+3\alpha ^{2}\right) }. \end{aligned}$$
(67)

4 Brane stability

In this section, we examine small perturbations about the soliton branes obtained in the previous sections. To this end, the metric and the scalar field are perturbed about the static brane and the resulting equations of motion are linearized in the proper gauge [6, 37–40]. The linearized equation turns out to be a Schrödinger-like equation with a potential U(z) which determines the linear modes. Unfortunately, for the models considered, this potential is too complicated to be handled analytically, or to be used for finding the corresponding modes. Therefore, we will only examine the potential near its minimum up to second order in z.

In order to check the stability of the brane models, we use the standard scalar–tensor–vector (STV) decompositions [42]. For this purpose, a new variable z is considered in such a way that \(dz=\mathrm{e}^{-A(w)} dw\), so the metric can be considered to be [42]

$$\begin{aligned} g_{MN}=a^2(z)\eta _{MN}, \end{aligned}$$
(68)

where \(a(z)=\mathrm{e}^{A(z)}\). On the other hand, one can linearize the action of the system given by Eq. (1) along with the metric, via \(\delta \varphi \) and \(\delta g_{MN}\equiv a^{2}(z)h_{MN}\), which is given by [42]

$$\begin{aligned} \delta ^{2}S= & {} \frac{1}{2}\int d^{5}x a^{3}\Big [\partial _{M}h_{NP}\partial ^{P}h^{MN}-\partial ^{M}h\partial ^{N}h_{MN} \nonumber \\&- \frac{1}{2}\partial _{P}h_{MN}\partial ^{P}h^{MN}+\frac{1}{2}\partial ^{P}h\partial _{p}h\nonumber \\&+3\frac{a'}{a}\left( h\partial ^{\mu }h_{\mu z}-h_{zz} h'\right) \nonumber \\&+2\kappa _{5}^{2}(a^{2}V_{\varphi \varphi }\left( \delta \varphi \right) ^{2} +2h^{Mz}\varphi '\partial _{M}\delta \varphi \nonumber \\&+\varphi 'h' \delta \varphi -\partial ^{M}\delta \varphi \partial _{M}\delta \varphi )\Big ], \end{aligned}$$
(69)

where \(\partial ^{M}=\eta ^{MN}\partial _{N}\), \(\partial ^{\mu }=\eta ^{\mu \nu }\partial _{\nu }\), \(h=\eta ^{MN}h_{MN}\), and the prime denotes a derivative with respect to z. Moreover, by introducing the following vector and tensor perturbations [42, 49]:

$$\begin{aligned} h_{\mu z}= & {} \partial _{\mu }F+G_{\mu }, \nonumber \\ h_{\mu \nu }= & {} \eta _{\mu \nu }\varphi +\partial _{\mu }\partial _{\nu }B+\partial _{\mu }C_{\nu }+\partial _{\nu }C_{\mu }+D_{\mu \nu }, \end{aligned}$$
(70)

where \(C_{\mu }\), \(G_{\mu }\), and \(D_{\mu \nu }\) are the transverse vector and tensor perturbations, respectively, one obtains

$$\begin{aligned} \partial ^{\mu }C_{\mu }= & {} 0=\partial ^{\mu }G_{\mu }, \nonumber \\ \partial ^{\mu }D_{\mu \nu }= & {} 0=D_{\mu }^{\mu }. \end{aligned}$$
(71)

By using this STV decomposition method, one can decompose \(\delta ^{(2)}S\) into the following decoupled parts [42, 49]:

$$\begin{aligned} \delta ^{(2)}S_{\mathrm{vector}}= & {} \frac{1}{2}\int \mathrm{d}^{5}x\hat{v}^{\mu }\Box ^{(4)}\hat{v}_{\mu },\nonumber \\ \delta ^{(2)}S_{\mathrm{tensor}}= & {} \frac{1}{4}\int \mathrm{d}^{5}x\hat{D}^{\mu \nu }\left[ \Box ^{(4)}\hat{D}_{\mu \nu }\right. \nonumber \\&\left. +\hat{D}_{\mu \nu }''-\frac{(a^{\frac{3}{2}})''}{a^{\frac{3}{2}}}\hat{D}^{\mu \nu }\right] , \end{aligned}$$
(72)

where

$$\begin{aligned} \hat{v}^{\mu }=a^{\frac{3}{2}}(G_{\mu }-C_{\mu }'), \quad \hat{D}^{\mu \nu }=a^{\frac{3}{2}}D^{\mu \nu }, \quad \Box ^{(4)}=\partial ^{\mu }\partial _{\mu }. \end{aligned}$$
(73)

The scalar perturbations of the action, in turn, lead to two parts [42, 49]:

$$\begin{aligned} \delta ^{(2)}S_{\mathrm{scalar}-1}=\int \mathrm{d}^{5}x a^{3}\left\{ 3\frac{a'}{a}h_{zz}{-}2\varphi '{-}2\kappa _{5}^{2}\mathcal{L}_{X}\varphi '\delta \varphi \right\} \Box ^{(4)}\psi , \end{aligned}$$
(74)

with \(\psi =F-\frac{1}{2}B'\) and

$$\begin{aligned} \delta ^{(2)}S_{\mathrm{scalar}-2}= & {} \frac{1}{2}\int \mathrm{d}^{5}x a^{3} [-3\varphi \Box ^{(4)}\varphi -3h_{zz}\Box ^{(4)}\varphi +6\varphi '\varphi '\nonumber \\&-3\frac{a'}{a}h_{zz}(h_{zz}'+4\varphi ')+2\kappa _{5}^{2}(\delta \varphi \Box ^{(4)}\delta \varphi \nonumber \\&+a^{2}\mathcal{L}_{\varphi \varphi }(\delta \varphi )^{2}+2 h_{zz}\varphi '\delta \varphi '\nonumber \\&+\varphi '(h_{zz}'+4\varphi ')\delta \varphi -(\delta \varphi ')^{2})]. \end{aligned}$$
(75)

After eliminating \(h_{zz}\) in action (75), by taking into account Eq. (74), and doing some simplifications, one obtains [42, 49]:

$$\begin{aligned} \delta ^{(2)}S_{\mathrm{scalar}-2}=\int \mathrm{d}^{5}x \hat{\mathcal {G}}\left\{ \Box ^{(4)}\hat{\mathcal {G}}+\hat{\mathcal {G}}''-\frac{\theta ''}{\theta }\hat{\mathcal {G}}\right\} , \end{aligned}$$
(76)

where \(\mathcal {G}\) and \(\theta \) are gauge invariant variables and specific functions of a, and they are given by [42, 49]

$$\begin{aligned} \mathcal {G}= & {} \frac{\kappa _{5}^{2}}{2}a^{\frac{3}{2}}\left( 2\delta \varphi -\frac{\varphi ' a}{a'\varphi }\right) ,\nonumber \\ \theta= & {} a^{\frac{3}{2}}\frac{\varphi ' a}{a'}, \end{aligned}$$
(77)

respectively. Finally, the equation of normal modes of the linear perturbations are

$$\begin{aligned} \text {vector:}&\quad \Box ^{(4)}\hat{v}_{\mu }=0, \nonumber \\ \text {tensor:}&\quad \Box ^{(4)}\hat{D}_{\mu \nu }+\hat{D}_{\mu \nu }''-\frac{(a^{\frac{3}{2}})''}{a^{\frac{3}{2}}}\hat{D}_{\mu \nu }=0, \nonumber \\ \text {scalar:}&\quad \Box ^{(4)}\hat{\mathcal {G}}+\hat{\mathcal {G}}''-\frac{\theta ''}{\theta }\hat{\mathcal {G}}=0. \end{aligned}$$
(78)

Since the vector perturbations only have zero modes, solutions are stable against vector perturbations [42]. However, for tensor and scalar modes another decomposition should be introduced [42], namely

$$\begin{aligned}&\hat{D}_{\mu \nu }(x^{\lambda },z)=\epsilon _{\mu \nu }\mathrm{e}^{ip_{\lambda }x^{\lambda }}\rho _{p}(z), \nonumber \\&\hat{\mathcal {G}}(x^{\lambda },z)=\mathrm{e}^{iq_{\lambda }x^{\lambda }}\Phi _{q}(z), \end{aligned}$$
(79)

where \(\epsilon _{\mu \nu }\) is the TT polarization tensor and \(\rho _{p}(z)\) and \(\Phi _{q}(z)\) satisfy the following equations:

$$\begin{aligned} \mathcal {A}_{t}\mathcal {A}_{t}^{\dag }\rho _{p}= & {} m_{p}^{2}\rho _{p},\end{aligned}$$
(80)
$$\begin{aligned} \mathcal {A}_{s}\mathcal {A}_{s}^{\dag }\Phi _{q}= & {} M_{q}^{2}\Phi _{q}, \end{aligned}$$
(81)

with \(m_{p}^{2}=-p^{\mu }p_{\mu }\), \(M_{q}^{2}=-q^{\mu }q_{\mu }\), and

$$\begin{aligned} \begin{aligned} \mathcal {A}_{t}&=\frac{\mathrm{d}}{\mathrm{d}z}+\frac{(a^{\frac{3}{2}})'}{a^{\frac{3}{2}}}, \\ \mathcal {A}_{s}&=\frac{\mathrm{d}}{\mathrm{d}z}+\frac{\theta ''}{\theta }. \end{aligned} \end{aligned}$$
(82)

Equation (80) is a Schrödinger-like equation, and it takes the form

$$\begin{aligned} -\rho _{p}''+\mathcal {W}(z)\rho _{p}=m_{p}^{2}\rho _{p}, \end{aligned}$$
(83)

where the effective potential \(\mathcal {W}(z)\) is given

$$\begin{aligned} \mathcal {W}(z)= & {} \frac{(a^{\frac{3}{2}})''}{a^{\frac{3}{2}}}=\frac{\left( \mathrm{e}^{\frac{3}{2}A(z)}\right) ''}{\mathrm{e}^{\frac{3}{2}A(z)}},\nonumber \\= & {} \frac{3}{2}A''+\frac{9}{4}A'^{2}. \end{aligned}$$
(84)

The effective four-dimensional gravity is determined by the spectrum of the tensor KK modes \(\rho _{p}\) [42]. For instance, for zero mode \(\rho _{0}\) with \(m_{0}=0\) and a normalizable \(\rho _{0}\) leads to the four-dimensional Newton’s law [12, 42, 50].

On the other hand, in order to study the stability of the branes, some authors choose an “axial gauge” where the metric is perturbed as [6, 37–40]:

$$\begin{aligned} \mathrm{d}s^{2}=\mathrm{e}^{ 2A(w)}(g_{\mu \nu } + \varepsilon h_{\mu \nu } )\mathrm{d}x^{\mu } dx^{\nu } - \mathrm{d}w^{2}. \end{aligned}$$
(85)

Here \(g_{\mu \nu }\) represents the background metric, \(h_{\mu \nu }\) denotes the metric perturbations, and \(\varepsilon \) is a small parameter [38]. They also consider the transformation \(\varphi \longrightarrow \varphi +\varepsilon \tilde{\varphi }\), where \(\tilde{\varphi }=\varphi (x,w)\) [40]. Moreover, in order to render the metric conformally flat, one can choose \(\mathrm{d}z=\mathrm{e}^{-A(w)} \mathrm{d}w\). In this case, the corresponding Schrödinger equation takes the form [6, 25, 37–40, 51]

$$\begin{aligned} -\frac{\mathrm{d}^{2}\psi (z)}{\mathrm{d}z^{2}}+U(z)\psi (z)=k^{2}\psi (z), \end{aligned}$$
(86)

where the potential is given by

$$\begin{aligned} U(z)=-\frac{9}{4}\Lambda +\frac{9}{4}A'^{2}+\frac{3}{2}A''. \end{aligned}$$
(87)

One can check that this potential and the \(\psi \) function are the same as \(\mathcal {W}(z)\) and \(\rho \) in Eq. (83). Note that \(\Lambda \) is a cosmological constant on the brane, which could be positive, negative or zero corresponding to the 4D spacetime being de Sitter (\(dS_{4}\)), anti-de Sitter (\(AdS_{4}\)) or Minkowski (\(M_{4}\)) [37, 38]. Besides, it is notable that the Hamiltonian corresponding to Eq. (86) can be written in the form [25, 37, 51]

$$\begin{aligned} H=\left( \frac{\mathrm{d}}{\mathrm{d}z}+\frac{3}{2}A'(z)\right) \left( -\frac{\mathrm{d}}{\mathrm{d}z}+\frac{3}{2}A'(z)\right) , \end{aligned}$$
(88)

which is obviously Hermitian and therefore leads to real k (\(k^{2}\ge 0\)). Accordingly, there are no unstable tachyonic excitations in the system [40, 51].

The solution for the zero modes (\(k = 0\)) is [40, 42, 51]:

$$\begin{aligned} \psi (z)=N\mathrm{e}^{\frac{3A(z)}{2}} \end{aligned}$$
(89)

where N is a normalization factor[40, 42, 51] and satisfies

$$\begin{aligned} 1= & {} \int _{-\infty }^{+\infty }\mathrm{d}z|\psi _{0}(z)|^{2}=N^{2}\int _{-\infty }^{+\infty }\mathrm{d}z\mathrm{e}^{3A(z)}\nonumber \\= & {} \frac{N^{2}}{l}\int _{-\infty }^{+\infty }\mathrm{d}y\mathrm{e}^{2A(y)}, \end{aligned}$$
(90)

where \(y=lw\) is a dimensionless variable. The asymptotic behavior of the solutions at large y are checked and the result is that only the \(\varphi ^{4}\) system has a normalizable zero mode and thus stable. On the other hand, as in quantum mechanical systems, we may check for the stability of the system via the existence of a real frequency, in bound states. Since the potentials for the three systems considered in this paper, are too complicated to be solved analytically, we found the corresponding ground state eigenvalues via expansions in terms of the fifth coordinate w. One can deduce the stability up to \(O(w^2)\) by looking at the sign of the \(w^2\) term. It is seen in Figs. 7, 8 and 9 that this coefficient is everywhere positive (stability) for the \(\varphi ^4\) system, while there are regions of the parameter space where the coefficient is negative (instability) for the SG and \(\varphi ^6\) systems. However, as noted in Sect. 3.1, this conclusion is not decisive, since higher-order effects might have drastic effects. The difference between the results of these two approaches is probably caused by the inevitable approximations used in the analysis.

Fig. 7
figure 7

The coefficient of the \(z^2\) term in the potential of the linearized Schrödinger equation as a function of the free parameters a and b for the SG system. Negative values correspond to a first-order instability

Full size image
Fig. 8
figure 8

The coefficient of the \(z^2\) term in the potential of the linearized Schrödinger equation as a function of the free parameters \(\alpha \) and \(\beta \) for the \(\phi ^4\) system. The coefficient is everywhere positive, signaling linear stability

Full size image
Fig. 9
figure 9

The coefficient of the \(z^2\) term in the potential of the linearized Schrödinger equation as a function of the free parameters \(\alpha \) and \(\beta \) for the \(\phi ^6\) system. Negative values correspond to a first-order instability. There are also vast patches in the parameter space which have almost neutral stability

Full size image

5 Conclusion

In this work, we obtained exact thick brane models inspired by well-known nonlinear systems, namely, the sine-Gordon (SG), \(\varphi ^{4}\) and \(\varphi ^{6}\) models. The confining effects of the scalar field in all these three models were confirmed by examining the geodesic equation for a test particle moving normal to the brane. In particular, it turns out that the modified potential for the SG system resembles that of the double sine-Gordon (DSG) system, while those of \(\varphi ^{4}\) and \(\varphi ^{6}\) became \(\varphi ^{6}\) and \(\varphi ^{8}\), respectively. We have extended previous brane models [10, 14, 39] based on SG and \(\varphi ^4\) solitons taking into account different parametrizations. The similarity of the \(\varphi ^4\) model with the generic Higgs model makes this choice particularly interesting, especially as the resulting potential is an odd function of the fifth coordinate and the \(Z_2\) symmetry is respected. We have studied the \(\varphi ^6\) model for the first time. This model is interesting by its own right, since unlike the \(\varphi ^4\) model, we have two pairs of solitons and anti-solitons, which live in different topological sectors.

In the case of the SG model, we have used a more general form of the potential compared to the one used in [10, 39]. The resulting brane does not have \(Z_2\) symmetry, in general, where the center of the brane may be displaced from \(w=0\) and the potential will not be an odd function of w in general. However, by a suitable choice of the model parameters it is possible to make the vacua of the effective potential degenerate, in which case the \(Z_2\) symmetry is restored. In the case of the \(\varphi ^6\) model, however, we could not restore this symmetry via re-parametrization. Finally, using standard procedures, we examined the stability of the thick branes, by determining the sign of the \(w^2\) term in the expansion of the potential for the resulting Schrödinger-like equation. It turns out that the \(\varphi ^4\) brane is stable, while there are unstable modes for certain ranges of the model parameters in the SG and \(\varphi ^6\) branes.

We considered the limiting case in which the brane tends to zero thickness and approaches a thin brane. It should be noted that the topological stability of the soliton brane remains valid even in this limit (at least at the classical level). An interesting question would be whether the thick brane continuous metric develops a discontinuity and whether the Israel junction conditions will be satisfied in this limit. Although one would expect intuitively that this is the case, we have not worked out the detailed calculations. This issue will be explored in separate paper.