text:
# Mode-coupling approach to non-Newtonian Hele-Shaw flow

## I Introduction

The celebrated Saffman-Taylor instability MW\_3  arises when a less viscous fluid pushes a more viscous one in the thin gap of a Hele-Shaw cell. The less viscous fluid can be either injected at an end of a channel-shaped cell (rectangular geometry) MW\_3 , or from the center of the cell (radial geometry) Pat . In both geometries the interface separating the fluids may deform, leading to the formation of fingerlike patterns, commonly known as viscous fingers. During the last four decades, the viscous fingering instability has been extensively studied, both theoretically MW\_1  and experimentally MW\_1b . Much of the research in this area has examined the case in which the fluids involved are Newtonian. For Newtonian fluids, it is observed that the fingers grow and compete dynamically, resulting in a single stable finger in the rectangular geometry, and in patterned structures markedly characterized by the spreading, and subsequent splitting of the finger tips in the radial setup.

A whole different class of interfacial patterns arise when the Saffman-Taylor instability is studied by taking the displaced fluid as non-Newtonian MW\_1b ; Van . In contrast to most Newtonian fluids, non-Newtonian fluids differ widely in their hydrodynamic properties, with different fluids exhibiting a range of effects from elasticity and plasticity to shear thinning and shear thickening. Experiments using non-Newtonian fluids in radial and rectangular Hele-Shaw cells have revealed a wide variety of new patterns, showing snowflake-like shapes MW\_1b ; Buka  and fracturelike structures Lemaire ; Zhao . Instead of the traditional, tip-splitting-dominated Newtonian patterns, these experiments MW\_1b ; Buka ; Lemaire ; Zhao  exhibit dendritic fingers and side branching. This morphological diversity and rich dynamical behavior motivated a number of theoretical studies of the problem Sader ; Bonn ; Kondic\_96 ; Kondic\_98 ; Kondic\_long . One major difficulty faced by researchers is that, as opposed to the Newtonian case, the pressure field is no longer Laplacian. This implies a serious theoretical challenge, since a priori one would not be allowed to directly apply a Darcy’s law approach to attack the problem.

The instability of radial Hele-Shaw flows involving non-Newtonian fluids has been studied theoretically by Sader et al. Sader . They considered power law fluids, and performed a linear stability analysis without making use of Darcy’s law. Essentially, they showed that decreasing the power law index dramatically increases the growth rates leading to a more rapid development of the fingering patterns. A Darcy’s law-type approach has been proposed by Bonn and collaborators Bonn . Bonn et al. suggested a modified Darcy’s law including a shear rate dependent viscosity, and showed that, within their approach, the pressure field remains Laplacian. In a series of interesting papers Kondic et al. Kondic\_96 ; Kondic\_98 , and subsequently Fast et al. Kondic\_long  extended reference Bonn  ideas and derived a generalized Darcy’s law from first principles, where viscosity depends upon the squared pressure gradient. It turns out that the Darcy’s law formula proposed by Bonn et al. Bonn  follows from the more basic version rigorously derived by Kondic et al. Kondic\_96 ; Kondic\_98 . Efficient numerical simulations performed in references Kondic\_98 ; Kondic\_long  have shown that shear thinning can suppress tip-splitting and leads to the formation of dendritic structures, presenting a clear side branching behavior.

Theoretical studies of the fully nonlinear stages of Hele-Shaw flow with non-Newtonian fluids rely heavily on intensive numerical simulations Kondic\_98 ; Kondic\_long . On the analytical side, the structure of the fingering dynamics in such complex fluids is largely restricted to linear stability investigations Sader ; Kondic\_96 . Much less attention has been paid to the analytical investigation of the dynamics that bridges the initial (purely linear) and final (fully nonlinear) time regimes. In addition, theoretical as well as experimental analyses of flow of shear thickening fluids in Hele-Shaw cells still need to be addressed. In this paper we carry out the analytical weakly nonlinear analysis for the intermediate stages of evolution, and examine both shear-thinning and shear thickening cases. We adapt a weakly nonlinear approach originally developed to study Newtonian Hele-Shaw flows MW\_rectang ; MW\_rad , to the non-Newtonian situation. We focus on the onset of the nonlinear effects, and try to understand how mode-coupling dynamics leads to basic morphological features and behaviors observed in non-Newtonian Hele-Shaw flows.

The article is organized as follows: section II formulates our theoretical approach. We perform a Fourier decomposition of the interface shape, and from an alternative form of Darcy’s law study the influence of weak shear effects on the development of interfacial patterns. In contrast to the analysis for Newtonian fluids, conventionally based on a scalar velocity potential, we employ a vector velocity potential capable of describing vorticity arising from the non-Newtonian fluid flow. Coupled, nonlinear, ordinary differential equations governing the time evolution of Fourier amplitudes are derived in detail. Section III discusses both linear and weakly nonlinear dynamics. It concentrates on the effect of shear thinning and shear thickening on finger tip-splitting and side branching. Section III.1 briefly discusses our linear stability results. Linear results are useful and instructive, but do not allow accurate predictions of important interfacial features. In section III.2 we show that some of such features can indeed be predicted and more quantitatively explained by our analytical mode-coupling approach. At second order we find a mechanism for finger tip-splitting in non-Newtonian Hele-Shaw flow: it is suppressed (favored) for shear thinning (thickening) fluids. Our results indicate absence of side branching in the weak shear limit and early flow stages, but suggest that it could be enhanced (inhibited) for shear thinning (thickening) fluids. Section III.3 discusses mode-coupling in rectangular flow geometry. Our chief conclusions are summarized in Sec. IV.

## II The mode-coupling differential equation

The Hele-Shaw cell (Fig.1) consists of two parallel plates separated by a small distance. The cell thickness $`b`$ is considered to be much smaller that any other length scale in the problem, so that the system is essentially two-dimensional. Consider the displacement of a viscous, non-Newtonian fluid by a Newtonian fluid of negligible viscosity in such confined geometry. The surface tension between the fluids is denoted by $`\sigma `$. The Newtonian fluid is injected at a constant areal flow rate $`Q`$ at the center of the cell.

The initial circular fluid-fluid interface is slightly perturbed, $`=R(t)+\zeta (\theta ,t)`$ ($`\zeta /R1`$), where the time dependent unperturbed radius is given by

$$R(t)=\sqrt{R_0^2+\frac{Qt}{\pi }},$$
(1)
$`R_0`$ being the unperturbed radius at $`t=0`$. The interface perturbation is written in the form of a Fourier expansion

$$\zeta (\theta ,t)=\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}\zeta _n(t)\mathrm{exp}\left(in\theta \right),$$
(2)
where $`\zeta _n(t)=(1/2\pi )_0^{2\pi }\zeta (\theta ,t)\mathrm{exp}(in\theta )d\theta `$ denotes the complex Fourier mode amplitudes and $`n`$=0, $`\pm 1`$, $`\pm 2`$, $`\mathrm{}`$ is the discrete polar wave number. In our Fourier expansion (2) we include the $`n=0`$ mode to keep the area of the perturbed shape independent of the perturbation $`\zeta `$. Mass conservation imposes that the zeroth mode is written in terms of the other modes as $`\zeta _0=(1/2R)\underset{n0}{}|\zeta _n(t)|^2`$.

The relevant hydrodynamic equation for Newtonian Hele-Shaw flows, Darcy’s law MW\_3 ; MW\_1 , states that the gradient of the pressure is proportional to the fluid velocity, and oriented in the opposite direction with respect to the fluid flow

$$p=\frac{12\mu }{b^2}𝐯,$$
(3)
where $`𝐯=𝐯(r,\theta )`$ and $`p=p(r,\theta )`$ are, respectively, the gap-averaged velocity and pressure in the fluid. The viscosity of the fluid is represented by $`\mu `$. To model non-Newtonian Hele-Shaw flows, we use a suitable form of the Darcy’s law in the weak shear thinning or thickening limit. We follow Bonn et al. Bonn  and consider their shear rate dependent viscosity in the weak non-Newtonian limit $`(\omega \tau )^21`$

$$\mu (\omega ^2)\mu _0\left[1\left(1\alpha \right)\omega ^2\tau ^2\right],$$
(4)
where $`\omega v/b`$ is the shear rate, $`v=|𝐯|`$, $`\tau `$ denotes the characteristic relaxation time of the fluid, and $`\mu _0`$ is a constant, zero-shear viscosity. The parameter $`\alpha `$ measures the shear dependence: $`\alpha =1`$ corresponds to the Newtonian fluids, and $`\alpha <1`$ ($`\alpha >1`$) gives the shear thinning (thickening) case.

By substituing Eq. (4) into Eq. (3) we obtain an alternative form of Darcy’s law ideally suited to describe weak non-Newtonian effects

$$P=𝐯+\delta v^2𝐯.$$
(5)
Here $`P=[b^2/(12\mu _0)]p`$ represents a generalized pressure field, and $`\delta =(1\alpha )(\tau /b)^2`$ is a small parameter that expresses the non-Newtonian nature of the displaced fluid: $`\delta =0`$ corresponds to the Newtonian case, while $`\delta >0`$ ($`\delta <0`$) describes the shear thinning (thickening) case. Taking the divergence of Eq. (5) and using incompressibility ($`𝐯=0`$) the pressure is seen to be anharmonic (nonvanishing Laplacian). Our subsequent analysis incorporates this effect, in contrast to the treatment in Ref. Bonn . Indeed, an alternate route to our Eq. (5) is to start with viscosity depending on the square of the pressure gradient (as in Kondic, et al. Kondic\_96 ), then approximate the pressure gradient with the velocity for small $`\delta `$.

Our pertubative approach keeps terms up to the second order in $`\zeta `$ and up to first order in $`\delta `$. Considering the fact that the viscosity of the Newtonian fluid is negligible, the generalized pressure jump condition at the interface can be written as MW\_1

$$P|_{}=\gamma \kappa _{}|_{},$$
(6)
where $`\gamma =b^2\sigma /(12\mu _0)`$ and $`\kappa _{}`$ is the curvature in the direction parallel to the plates.

The weakly nonlinear approach to radial, Newtonian Hele-Shaw flow developed in reference MW\_rad , related the fluid velocity to a scalar velocity potential $`𝐯=\varphi `$, this replacement made possible by the irrotational nature of the flow for Newtonian fluids. For non-Newtonian fluids, in contrast, flows governed by the modified Darcy’s law (5) exhibit vorticity. Hence we perform our calculations using a vector potential $`𝐯=\times 𝐀`$. The most general form of the vector potential can be written as

$$𝐀=\left[\frac{Q}{2\pi }\theta +\underset{m,n0}{}A_{mn}\left(\frac{R}{r}\right)^m\mathrm{exp}(in\theta )\right]\widehat{𝐳},$$
(7)
where $`A_{mn}`$ are the Fourier coefficients of the velocity vector potential and $`\widehat{𝐳}`$ is the outward unit-normal to the upper cell plate. The radial and polar components of the fluids velocities are

$$v_r=\frac{Q}{2\pi r}+\underset{m,n0}{}inA_{mn}\left(\frac{R^m}{r^{m+1}}\right)\mathrm{exp}(in\theta ),$$
(8)
and

$$v_\theta =\underset{m,n0}{}mA_{mn}\left(\frac{R^m}{r^{m+1}}\right)\mathrm{exp}(in\theta ).$$
(9)
Note that the vector potential reduces to the unperturbed steady flow with a circular interface ($`v_r=Q/2\pi r`$, $`v_\theta =0`$) when $`R0`$ and also when $`r\mathrm{}`$.

We exploit the fact that $`P`$ must be curl free, and impose the so-called solvability condition $`\times P=0`$. It simplifies the general form of the vector potential expansion given in Eq. (7). The solvability condition reveals that, without loss of generality, one can rewrite the vector potential as

$$𝐀=\left\{\frac{Q}{2\pi }\theta +\underset{n0}{}A_n\left(\frac{R}{r}\right)^{|n|}\mathrm{exp}\left(in\theta \right)+\delta \left[\underset{n0}{}B_n\left(\frac{R}{r}\right)^{|n|}\frac{1}{r^2}\mathrm{exp}(in\theta )\right]\right\}\widehat{𝐳},$$
(10)
replacing the array of coefficients $`A_{mn}`$ with the simpler set of $`A_n`$ and $`B_n`$. Observe that the vector potential (10) is simply a superposition of a purely Newtonian term ($`\delta ^0`$, coefficients $`A_n`$) and a non-Newtonian contribution ($`\delta ^1`$, coefficients $`B_n`$)

$$𝐀=𝐀_N+𝐀_{NN}.$$
(11)
The flow described by $`𝐀_N`$ is irrotational, while $`𝐀_{NN}`$ has a curl.

Similarly, we express the pressure of the outer fluid as a sum of Newtonian and non-Newtonian pressures, and propose a general form for their Fourier expansion

$$P=P_N+P_{NN},$$
(12)
where

$$P_N=\frac{Q}{2\pi }\mathrm{log}\left(\frac{r}{R}\right)+\underset{n0}{}p_n\left(\frac{R}{r}\right)^{|n|}\mathrm{exp}(in\theta ),$$
(13)
and

$$P_{NN}=\delta \left[\left(\frac{Q}{2\pi }\right)^3\frac{1}{2r^2}+\underset{n0}{}q_n\left(\frac{R}{r}\right)^{|n|}\frac{1}{r^2}\mathrm{exp}(in\theta )\right].$$
(14)
The gradient of the complex pressure field (12) must satisfy the non-Newtonian Darcy’s law given by Eq. (5). By inspecting the $`r`$ and $`\theta `$ components of (5), and by examining the Newtonian and non-Newtonian components of it, we can express the Fourier coefficients of $`P_N`$, $`P_{NN}`$ and $`𝐀_{NN}`$ in terms of the Fourier coefficients of $`𝐀_N`$,

$$p_n=\mathrm{sgn}(n)\left(iA_n\right),$$
(15)
$$q_n=iA_n\beta (n)\frac{n}{|n|}\left(\frac{Q}{2\pi }\right)^2+\frac{Q}{4\pi \left(|n|+1\right)}\underset{m0,mn}{}mA_m\left(nm\right)A_{nm}k(n,m),$$
(16)
$$B_n=A_n\alpha (n)\left(\frac{Q}{2\pi }\right)^2+\frac{Q}{4\pi \left(|n|+1\right)}\underset{m0,mn}{}m\left(iA_m\right)\left(nm\right)A_{nm}h(n,m),$$
(17)
where in order to keep the results in a compact form, we introduced the coefficients

$$\alpha \left(n\right)=\frac{|n|\left(|n|1\right)}{2\left(|n|+1\right)},$$
(18)
$$\beta \left(n\right)=\frac{|n|\left(|n|+3\right)}{2\left(|n|+1\right)},$$
(19)
$$h(n,m)=\left(|n|+2\right)\mathrm{sgn}(nm)\frac{n}{2}\left\{3\mathrm{sgn}\left[m\left(nm\right)\right]\right\},$$
(20)
and

$$k(n,m)=n\mathrm{sgn}(nm)+\frac{1}{2}\left(|n|+2\right)\left\{3\mathrm{sgn}\left[m\left(nm\right)\right]\right\}.$$
(21)
Note that $`\mathrm{sgn}(n)=1`$ if $`n>0`$ and $`\mathrm{sgn}(n)=1`$ if $`n<0`$.

Using Eqs. (15)-(17), which are consistent with the solvability condition and Darcy’s law (5), we can derive the general expression of the vector potential Fourier coefficients in terms of the perturbation amplitudes. To fulfill this goal, consider the kinematic boundary condition MW\_1 ; MW\_1b

$$\frac{}{t}=\left[\frac{1}{r}\frac{}{\theta }\left(v_\theta \right)+v_r\right]_|_{},$$
(22)
which states that the normal components of each fluid’s velocity at the interface equals the velocity of the interface itself. By expanding Eq. (22) up to the second order in $`\zeta `$ and up to first order in $`\delta `$ we find the coefficient of the vector potential corresponding to the $`n`$-th evolution mode, $`A_n^{(k)}`$, in terms of $`\delta `$ and the $`k`$-th order in $`\zeta `$ ($`k`$ = 1, 2)

$$iA_n^{\left(1\right)}(t)=\left[\frac{R}{n}\dot{\zeta }_n+\frac{\dot{R}}{n}\zeta _n\right]\left[1+\delta \dot{R}^2\alpha \left(n\right)\right],$$

$`iA_n^{\left(2\right)}(t)`$ $`=`$ $`{\displaystyle \frac{\dot{R}}{R}}{\displaystyle \underset{m0,n}{}}\left[{\displaystyle \frac{|m|}{m}}+\delta \dot{R}^2u(n,m)\right]\zeta _m\zeta _{nm}+{\displaystyle \underset{m0,n}{}}\left[{\displaystyle \frac{1}{n}}+{\displaystyle \frac{|m|}{m}}+\delta \dot{R}^2v(n,m)\right]\dot{\zeta }_m\zeta _{nm}`$
$``$ $`{\displaystyle \frac{\delta \dot{R}R}{2\left(|n|+1\right)}}{\displaystyle \underset{m0,n}{}}h(n,m)\dot{\zeta }_m\dot{\zeta }_{nm}{\displaystyle \frac{\delta \dot{R}^2}{2\left(|n|+1\right)}}{\displaystyle \underset{m0,n}{}}h(n,m)\zeta _m\dot{\zeta }_{nm},`$
where the overdot denotes total time derivative, $`\dot{R}=Q/(2\pi R)`$ represents the unperturbed interface velocity, and the coefficients

$$u(n,m)=\alpha (n)\frac{|m|}{m}2\frac{\alpha (m)}{m}\frac{h(n,m)}{2\left(|n|+1\right)},$$
(24)
$$v(n,m)=\alpha (n)\left(\frac{1}{n}+\frac{|m|}{m}\right)2\frac{\alpha (m)}{m}\frac{h(n,m)}{2\left(|n|+1\right)}.$$
(25)
To conclude our derivation we need one more step. The vector potential coefficients can be introduced into the pressure jump condition (Eq.(6)) and using Darcy’s law (Eq.(5)) one can finally find the equation of motion for perturbation amplitudes $`\zeta _n`$. We present the evolution of the perturbation amplitudes in terms of $`\delta `$ and the $`k`$-th order in the perturbation amplitude $`\zeta `$

$$\dot{\zeta }_n=\dot{\zeta }_n^{\left(1\right)}+\dot{\zeta }_n^{\left(2\right)},$$
(26)
where

$$\dot{\zeta }_n^{\left(1\right)}=\lambda \left(n\right)\zeta _n,$$
(27)
$$\lambda (n)=\frac{\dot{R}}{R}\left(|n|1\right)\frac{\gamma }{R^3}|n|\left(n^21\right)+\delta \dot{R}^2\frac{|n|}{R}\left[\dot{R}\frac{|n|1}{|n|+1}\frac{2\gamma }{R^2}\frac{|n|\left(n^21\right)}{|n|+1}\right],$$
(28)
is the non-Newtonian linear growth rate first reported in Kondic, et al. Kondic\_96 , and

$`\dot{\zeta }_n^{\left(2\right)}`$ $`=`$ $`{\displaystyle \underset{mn,0}{}}\left[F_N(n,m)+\delta F_{NN}(n,m)\right]\zeta _m\zeta _{nm}+{\displaystyle \underset{mn,0}{}}\left[G_N(n,m)+\delta G_{NN}(n,m)\right]\dot{\zeta }_m\zeta _{nm}`$ (29)
$`+`$ $`\delta {\displaystyle \underset{mn,0}{}}H_{NN}(n,m)\zeta _m\dot{\zeta }_{nm}+\delta {\displaystyle \underset{mn,0}{}}J_{NN}(n,m)\dot{\zeta }_m\dot{\zeta }_{nm}.`$
In Eq. (29) the coefficients $`F_N,F_{NN},G_N,G_{NN},H_{NN}`$, and $`J_{NN}`$ represent the second order Newtonian (N) and non-Newtonian (NN) terms. Their detailed functional form are presented in the appendix. An important feature of the second order coefficients is that they present special reflection symmetries

$`𝒞(n,m)=𝒞(n,m)`$
$`𝒞(n,m)=𝒞(n,m)`$ (30)
where $`𝒞=F_N,F_{NN},G_N,G_{NN}`$, $`H_{NN}`$, and $`J_{NN}`$, respectively. Equation (26) is the mode-coupling equation of the non-Newtonian Saffman-Taylor problem in radial geometry. It gives us the time evolution of the perturbation amplitudes $`\zeta _n`$, accurate to second order, in the weak shear limit. The rest of the paper uses Eq. (26) to study the development of interfacial instabilities, and to examine how the non-Newtonian parameter $`\delta `$ affects pattern morphology.

## III Discussion

In the next three subsections we use our mode-coupling approach to investigate the interface evolution at first and second order. To simplify our discussion it is convenient to rewrite the net perturbation (2) in terms of cosine and sine modes

$$\zeta (\theta ,t)=\zeta _0+\underset{n=1}{\overset{\mathrm{}}{}}\left[a_n(t)\mathrm{cos}(n\theta )+b_n(t)\mathrm{sin}(n\theta )\right],$$
(31)
where $`a_n=\zeta _n+\zeta _n`$ and $`b_n=i\left(\zeta _n\zeta _n\right)`$ are real-valued. Without loss of generality we may choose the phase of the fundamental mode so that $`a_n>0`$ and $`b_n=0`$.

### III.1 First order

Although at the level of linear analysis we do not expect to detect or rigorously predict important nonlinear effects such as tip-splitting and side branching, some useful information may still be extracted. A nice example of how purely linear results can help to understand complicated morphological features appearing in non-Newtonian radial Hele-Shaw flows is found in references Kondic\_96 ; Kondic\_98 . By studying their linear growth rate, Kondic et al. Kondic\_96 ; Kondic\_98  found that shear thinning decreases the wave number of maximal growth, increases the maximum growth rate, and tightens the band of unstable modes. Based on this increased selectivity of wavelengths, they postulated that shear thinning can lead to suppression of tip-splitting. Their speculations have not been further investigated analytically, but instead have been supported by their own intensive numerical simulations Kondic\_98 ; Kondic\_long .

Our first order results, namely our linear growth rate expression (28), agrees with Kondic et al. Kondic\_96  equivalent formula, if we set their constant $`C`$ equal to $`\dot{R}^2`$. We use growth rate (28) to gain insight into shear thinning/thickening behavior. Figure 2 plots $`\lambda (n)`$ as a function of mode number $`n`$ for three different values of the non-Newtoniam parameter: a) $`\delta =0`$; b) $`\delta =0.05`$; and c) $`\delta =0.05`$. By inspecting Fig. 2, we notice that, unlike the shear thinning case ($`\delta >0`$) discussed in references Kondic\_96 ; Kondic\_98 , shear thickening ($`\delta <0`$) widens the band of unstable modes, and decreases growth rate for the wave number of maximal growth. So, by applying similar arguments as those used by Kondic et al. Kondic\_96 ; Kondic\_98  for the weak shear thinning case, we postulate that shear thickening would lead to enhanced tip-splitting. Unfortunately, for the shear thickening radial Hele-Shaw flow, both numerical simulations and experiments are not available in the literature to confirm this claim. Of course, it is also of interest to study such a possibility analytically. This is one of the topics we examine in Sec. III.2.

### III.2 Second order

On the analytical side, one must go beyond linear analysis in order to investigate in more detail the main morphological features found in Hele-Shaw flows with both shear thinning and shear thickening fluids. To this day, besides a few first order linear growth studies Kondic\_96 ; Kondic\_98 ; Kondic\_long , there are no other analytical results concerning non-Newtonian fluid evolution in radial Hele-Shaw geometry. If on one hand the non-Newtonian problem is vastly richer, it is also true it is much less amenable to analytic attack than its Newtonian counterpart. Our mode-coupling approach intends to provide a better and more complete analytical understanding of the complex non-Newtonian time evolution dynamics in Hele-Shaw cells.

#### III.2.1 Action of the parameter $`\delta `$ on finger tip-splitting

We use the mode-coupling equation (26) to investigate the coupling of a small number of modes. At second order the most noteworthy effect refers to the action of the non-Newtonian parameter $`\delta `$ on finger tip-splitting. Tip-splitting is related to the influence of a fundamental mode $`n`$ on the growth of its harmonic $`2n`$ MW\_rad . For consistent second order expressions, we replace the time derivative terms $`\dot{a}_n`$ and $`\dot{b}_n`$ by $`\lambda (n)a_n`$ and $`\lambda (n)b_n`$, respectively. Under these circumstances the equation of motion for the cosine mode $`2n`$ is

$$\dot{a}_{2n}=\lambda (2n)a_{2n}+\frac{1}{2}T(2n,n)a_n^2,$$
(32)
where the function that multiplies $`a_n^2`$, $`T(2n,n)`$, is called the tip-splitting function and its general expression is

$`T(n,m)`$ $`=`$ $`F_N(n,m)+\lambda (m)G_N(n,m)+\delta [F_{NN}(n,m)+\lambda (m)G_{NN}(n,m)`$ (33)
$`+`$ $`\lambda (nm)H_{NN}(n,m)+\lambda (m)\lambda (nm)J_{NN}(n,m)].`$
Equation (32) shows that the presence of the fundamental mode $`n`$ forces growth of the harmonic mode $`2n`$. The function $`T(2n,n)`$ acts like a driving force and its sign dictates if finger tip-splitting is favored or not by the dynamics. If $`T(2n,n)<0`$, $`a_{2n}`$ is driven negative, precisely the sign that leads to finger tip-widening and finger tip-splitting. If $`T(2n,n)>0`$ growth of $`a_{2n}>0`$ would be favored, leading to outwards-pointing finger tip-narrowing.

To investigate the influence of the non-Newtonian parameter $`\delta `$ on tip-splitting behavior at second order, we plot in Fig. 3 the behavior of $`T(2n,n)`$ as a function of $`\delta `$, for a few Fourier modes ($`n=4,5,6`$). To simplify our analysis we adopt an instantaneous approach: we consider a particular $`\dot{R}`$ and $`R`$ combination, using the identity $`\dot{R}R=Q/2\pi `$, at the onset of growth of mode $`2n`$ \[using the condition $`\lambda (2n)=0`$\] in the Newtonian limit $`\delta =0`$, where we know $`T(2n,n)`$ is negative MW\_rad . We see from Fig. 3 that the curve for mode $`n=4`$ lies below the curves associated to $`n=5,6`$. This is an expected behavior since smaller values of $`n`$ would mean more room for the existing fingers to split. We also observe that the curves associated to smaller $`n`$ show a stronger angular inclination with respect to the horizontal $`\delta `$-axis. Therefore, lower Fourier modes would be more sensitive to variations in $`\delta `$.

Further inspection of Fig. 3 reveals that, for $`\delta >0`$, $`T(2n,n)`$ becomes less negative as $`\delta `$ increases, meaning that the interface has less tendency toward tip-splitting. In contrast, for $`\delta <0`$ we observe that $`T(2n,n)`$ becomes more negative as the magnitude of $`\delta <0`$ increases, indicating an enhanced tendency of the fingers to split at their tips. We conclude that tip-splitting is suppressed for shear thinning fluids, and enhanced for shear thickening ones. This fact is more clearly illustrated in Fig. 4. In Fig. 4 we plot the fluid-fluid interface for a certain time ($`t=30\mathrm{s}`$), considering the interaction of two cosine modes (a fundamental $`n=4`$ and its harmonic $`2n=8`$), for three diferrent values of $`\delta `$: (a) $`\delta =0`$ (solid curve); (b) $`\delta >0`$ (black dashed curve); and (c) $`\delta <0`$ (gray dashed curve). From Fig. 4 it is evident that there is a stronger splitting in the shear thickening case.

These second order effects regarding tip-splitting are consistent with the first order effects described in section III.1. With respect to shear thinning behavior, our analytical results agree with numerical simulations Kondic\_98 ; Kondic\_long  and experiments MW\_1b ; Buka  of fully nonlinear stages of interface evolution. In addition, we detect favored tip-splitting in the shear thickening case, a relevant nonlinear behavior not previously reported in the literature. We have also verified that, for a given mode $`n`$, smaller values of the surface tension parameter $`\gamma `$ lead to enhancement (suppression) of splitting events for $`\delta <0`$ ($`\delta >0`$).

A physical mechanism that seems to be at work in shear thinning case has been proposed in references Kondic\_96 ; Kondic\_98 ; Kondic\_long : viscosity is lowered in high-shear regions, which are in front of the fingers. Hence the tip of a finger experiences less resistance than the part of the interface near the tip. This suggests that tip-splitting would be suppressed if the fluid shear-thins. Our second order mode-coupling results indicate that a similar mechanism is at work in the shear thickening case: the resistance would be increased at the finger tips, suggesting that tip-splitting would be increased.

#### III.2.2 Action of the parameter $`\delta `$ on side branching

Another relevant non-Newtonian effect that can be studied at second order refers to the side branching phenomenon. In the framework of a mode-coupling theory, side branching requires the presence of mode $`3n`$. If the harmonic mode $`a_{3n}`$ is positive and sufficiently large, it can produce interfacial lobes branching out sidewards which we interpret as side branching.

Consider the influence of the fundamental mode $`n`$, and its harmonic $`2n`$, on the growth of mode $`3n`$. The equation of motion for the cosine $`3n`$ mode is

$$\dot{a}_{3n}=\lambda (3n)a_{3n}+\frac{1}{2}S(3n)a_na_{2n},$$
(34)
where the side branching function $`S(3n)=\left[T(3n,n)+T(3n,2n)\right]`$ can be easily obtained from Eq. (33). By analyzing Eq. (34) we observe that mode $`3n`$ can be spontaneously generated due to the driving term proportional to $`a_na_{2n}`$, such that it enters through the dynamics even when it is missing from the initial conditions. The existence and phase of mode $`3n`$ depends on the interplay of the modes $`n`$ and $`2n`$. Side branching would be favored if $`a_{3n}>0`$.

To study the growth of mode $`3n`$ as the non-Newtonian parameter $`\delta `$ is varied, we plot $`S(3n)`$ as a function of $`\delta `$ in Fig. 5. We consider a particular $`\dot{R}`$ and $`R`$ combination that corresponds to the onset of growth of mode $`3n`$ (i.e. obeying the condition $`\lambda (3n)=0`$) in the Newtonian limit $`\delta =0`$. From Fig. 5 one can verify that in the weak shear limit we consider in this work, $`S(3n)`$ is negative for all values of $`\delta `$. As shown in section III.2.1, starting with a fundamental mode $`a_n`$, the harmonic mode $`a_{2n}`$ is driven negative. Hence the product $`S(3n)a_na_{2n}`$ in Eq. (34) is positive, driving $`a_{3n}>0`$.

Whether side branching actually occurs depends on the magnitude of $`a_{3n}`$. Indeed, our discussion shows the presence of $`a_{3n}>0`$ even in the case of Newtonian fluids, where side branching does not occur. Why does side branching occur for shear thinning fluids? We speculate its origin may lie in the cubic form of our modified Darcy’s law, Eq. (5), which will add a term of the form $`\delta a_n^3`$ to $`\dot{a}_{3n}`$. Unfortunately, we have not calculated this term since we stopped our derivation at second order. Third order calculations for non-Newtonian Hele-Shaw flows lead to several new and complex mode-coupling terms, whose analysis and interpretation are not at all obvious, and go beyond the scope of our current work

Furthermore, there is a delicate interplay of mode $`3n`$ and mode $`2n`$, described by the evolution equation (34) and a similar expression for the growth of mode $`2n`$,

$$\dot{a}_{2n}=\lambda (2n)a_{2n}+\frac{1}{2}[T(2n,n)a_n^2+S(2n)a_na_{3n}],$$
(35)
where $`S(2n)=[T(2n,n)+T(2n,3n)]>0`$. Hence, side branching via a positive $`a_{3n}`$ will tend to drive $`a_{2n}`$ positive (or at least less negative), inhibiting tip splitting but also reducing the growth rate of $`a_{3n}`$ itself.

### III.3 Rectangular geometry limit

It is interesting to study how the main interfacial features examined in the radial geometry behave if we consider flow of non-Newtonian fluids in rectangular Hele-Shaw cells. A mode-coupling equation describing the system can be obtained if we take the “rectangular geometry limit”: $`R\mathrm{}`$ and $`Q\mathrm{}`$, such that $`Q/(2\pi R)v_{\mathrm{}}`$ and $`n/Rk`$ remain constant, where $`v_{\mathrm{}}`$ is the flow velocity at infinity and $`k`$ denotes the wave number of the disturbance. In this limit the interface evolution reverts to the evolution of the rectangular flow geometry, with linear growth rate given by

$$\lambda (k)=|k|\left[v_{\mathrm{}}(1+\delta v_{\mathrm{}}^2)\gamma |k|^2(1+2\delta v_{\mathrm{}}^2)\right].$$
(36)
In addition, the only nonzero second order mode-coupling terms are

$$G_N(k,k^{})=|k|\left[1\mathrm{sgn}(kk^{^{}})\right],$$
(37)
$$G_{NN}(k,k^{})=|k|v_{\mathrm{}}^2\left[1\mathrm{sgn}(kk^{^{}})\right],$$
(38)
and

$$J_{NN}(k,k^{})=\frac{v_{\mathrm{}}}{2}\left\{3\mathrm{sgn}[k^{^{}}(kk^{^{}})]+2\mathrm{sgn}[k(kk^{^{}})]\right\}.$$
(39)
Note that all rectangular limit expressions (36)-(39) are time-independent, and preserve the reflection symmetries shown in Eq. (II). Based on these findings, and using Eq. (26) we obtain a more compact mode-coupling equation for non-Newtonian flow in rectangular Hele-Shaw cells

$`\dot{\zeta }_k`$ $`=`$ $`\lambda (k)\zeta _k+{\displaystyle \underset{k^{^{}}k,0}{}}\left[G_N(k,k^{^{}})+\delta G_{NN}(k,k^{^{}})\right]\dot{\zeta }_k^{^{}}\zeta _{kk^{^{}}}`$ (40)
$`+`$ $`\delta {\displaystyle \underset{k^{^{}}k,0}{}}J_{NN}(k,k^{^{}})\dot{\zeta }_k^{^{}}\dot{\zeta }_{kk^{^{}}}.`$
The rectangular mode-coupling Eq. (40) is useful to study the influence of weak shear effects on the development of the fluid-fluid interface in rectangular cells. First, we examine finger tip-splitting related issues (finger narrowing/widening). As discussed in section III.2.1, we analyze the influence of the fundamental $`k`$ on the growth of its harmonic $`2k`$. The equation of motion for the cosine mode $`2k`$ is

$$\dot{a}_{2k}=\lambda (2k)a_{2k}+\frac{1}{2}T(2k,k)a_k^2,$$
(41)
where

$$T(2k,k)=\delta \left[2k^2v_{\mathrm{}}(v_{\mathrm{}}\gamma k^2)^2\right].$$
(42)
Note that the driving force term $`T(2k,k)`$ vanishes in the Newtonian limit $`\delta =0`$. This agrees with the results obtained in Ref. MW\_rectang  for Newtonian rectangular flow. Notice further, that for the non-Newtonian case ($`\delta 0`$) the harmonic mode $`2k`$ grows spontaneously even if missing from the initial conditions. The selected sign of $`a_{2k}`$ is given by the sign of $`T(2k,k)`$ and thus is dictated by $`\delta `$: if $`\delta >0`$ ($`\delta <0`$), $`T(2k,k)`$ is positive (negative), and hence $`a_{2k}>0`$ ($`a_{2k}<0`$).

Considering the case $`\delta >0`$, this means that the fingers become narrower as the shear thinning behavior becomes more pronounced. Moreover, by inspecting Eq. (42) we observe that the width of the finger tips decreases for increasingly larger values of the flow velocity $`v_{\mathrm{}}`$. Of course, these effects are minimized for larger values of the surface tension parameter $`\gamma `$. All these results are in agreement with recent experimental Lindner  and theoretical investigations Ben1 ; Ben2  for weak shear-thinning flows in rectangular geometry, in which shear-thinning narrows the width of steady state fingers.

If we take $`k=k^{}\sqrt{v_{\mathrm{}}/3\gamma }`$ as the fastest growing mode, then $`\lambda (2k)<0`$ so exponential growth of $`a_{2k}`$ is prevented. Unlike the case of radial flow where any mode $`n`$ eventually goes unstable for sufficiently large $`R`$, the growth rates $`\lambda (k)`$ in the rectangular geometry are time independent. Instead of exponential growth of $`a_{2k}`$, we instead expect its amplitude to saturate at the value $`T(2k,k)a_k^2/(2\lambda (2k))`$ obtained by setting $`\dot{a}_{2k}=0`$ in Eq. (41). Alternatively, one could specially prepare an initial condition with an initial perturbation of wavevector $`k`$ sufficiently small that $`2k`$ still lies within the band of unstable modes (up to about $`\sqrt{v_{\mathrm{}}/\gamma }`$), in which case mode $`2k`$ will be spontaneously generated and able to grow to a considerable magnitude.

For completeness, we also discuss the shear thickening case: for $`\delta <0`$ the function $`T(2k,k)`$ is negative, favoring finger tip-widening. Even though this shear thickening behavior has not been studied experimentally, the broadening of the fingers in shear thickening flow in rectangular cells has been theoretically predicted by numerical simulations Ben2 . Our analytical results reinforce the correctness of such numerical predictions.

Finally, we briefly discuss side branching behavior in rectangular Hele-Shaw cells. We consider the influence of modes $`k`$ and $`2k`$ on the evolution of mode $`3k`$. In this case, the relevant equation of motion has the form

$$\dot{a}_{3k}=\lambda (3k)a_{3k}+\frac{1}{2}S(3k)a_ka_{2k},$$
(43)
where, as it was in the radial flow case, the coefficient $`S(3k)`$ is the sum of $`T(3k,k)`$ and $`T(3k,2k)`$

$$S(3k)=\delta \left[8k^2v_{\mathrm{}}(v_{\mathrm{}}\gamma k^2)(v_{\mathrm{}}4\gamma k^2)\right].$$
(44)
The favored sign of $`a_{3k}`$ depends on the sign of the product $`a_ka_{2k}`$ and also on the sign of $`S(3k)`$. If $`\delta `$ is positive (shear thinning) then $`S(3k)`$ is negative in the band of wavevectors

$$\sqrt{v_{\mathrm{}}/4\gamma }<k<\sqrt{v_{\mathrm{}}/\gamma }$$
(45)
and positive outside this band. As we saw previously in the case of mode $`2k`$, if the fundamental mode is taken as the fastest growing mode $`k=k^{}`$, then exponential growth of the harmonic $`3k`$ is inhibited because $`\lambda (3k)<0`$. Instead it will saturate at a magnitude of $`S(3k)a_ka_{2k}/(2\lambda (3k))`$.

However, for a specially prepared initial condition perturbed at wavevector $`kk^{}/3`$, the harmonics $`2k`$ and $`3k`$ will both be able to grow. By Eq. (41) a positive harmonic mode $`a_{2k}>0`$ will appear spontaneously. Then, by Eq. (43), modes $`a_k`$ and $`a_{2k}`$ will conspire to create mode $`a_{3k}`$. Because $`kk^{}/3`$ lies outside the wavevector band (45), the value of $`S(3k)`$ is positive. Hence the mode $`a_{3k}`$ is driven positive, the sign needed to create side branches. While this conclusion is only tentative (due to our neglect of the complete set of third order terms) it would be interesting to test via numerical simulation or experiment.

## IV Concluding remarks

Visually striking patterns arise when a less viscous Newtonian fluid displaces a more viscous non-Newtonian fluid in the confined geometry of a Hele-Shaw cell. These complex patterned structures are the result of finger tip-splitting and dendritic interfacial instabilities. Due to the complicated dynamics of the system, and also from the limitations imposed by purely linear analysis, the majority of the theoretical studies in this area of research rely heavily on sophisticated numerical methods. In this work, we developed an alternative theoretical approach to the problem which allowed us to address important nonlinear issues analytically. A key point in our derivation was the introduction of velocity vector potential (as opposed to scalar potential) which was demanded by the non-Newtonian fluid flow.

We started our investigation by analyzing the role of fluid viscosity anisotropy on the development of the Saffman-Taylor instability in radial Hele-Shaw cells. To approach the problem analytically we considered the weak shear limit, and focused on the onset of nonlinear effects. In order to examine the influence of shear thinning and shear thickening on the shape of the emerging patterns, we derived a mode-coupling equation which is ideally suited to describe the weakly nonlinear interface evolution. Our analytical results show that finger tip-splitting is enhanced (diminished) in the case of shear thickening (thinning) fluids. We applied the rectangular geometry limit to our radial mode-coupling equations of motion, and also studied finger-tip and side branching behavior in rectangular cells. Our main findings predict finger tip-widening (narrowing) for shear thickening (thinning) fluids.

In neither the radial nor the rectangular geometry limit were we able to show the spontaneous appearance of side branching, except for some specially prepared initial conditions in rectangular geometry. We speculate that extending the current approach to third order could shed more light on this problem.

In summary, our analytical mode-coupling approach detects and predicts several features of the patterns formed in non-Newtonian Hele-Shaw flows in both radial and rectangular geometries. Furthermore, it predicts some behaviors not yet investigated in the literature, such as the evolution mechanisms involving tendency towards advanced finger tip-splitting and reduction of side branching for shear thickening fluids. We hope the main results presented in this paper will prompt further theoretical and experimental work on non-Newtonian Hele-Shaw flows, especially in the area of early stage transients in rectangular geometry where our theory makes specific predictions that have not been subject to experimental or numerical check.

###### Acknowledgements.

J.A.M. thanks the Brazilian Research Council - CNPq (through its PRONEX Program) for financial support. Work of M.W. was supported in part by the National Science Foundation grant No. DMR-0111198. M. C. acknowledges support from Carnegie Mellon University and University of Maryland (National Science Foundation grant No. DMR-00-80008).

*

## Appendix A Second order mode-coupling terms

This appendix presents the expressions for the second order Newtonian (N), and non-Newtonian (NN) mode-coupling coefficients which appear in Eq. (29)

$$F_N(n,m)=\frac{|n|}{R}\left\{\frac{\dot{R}}{R}\left[\frac{1}{2}\mathrm{sgn}(nm)\right]\frac{\gamma }{R^3}\left(1\frac{nm}{2}\frac{3m^2}{2}\right)\right\},$$
(46)
$`F_{NN}(n,m)`$ $`=`$ $`{\displaystyle \frac{|n|}{R}}\dot{R}^2\{{\displaystyle \frac{\dot{R}}{R}}[{\displaystyle \frac{2\alpha (m)}{m}}{\displaystyle \frac{|n|}{n}}{\displaystyle \frac{2\beta (n)}{|n|}}+{\displaystyle \frac{f(n,m)}{2(|n|+1)}}]`$ (47)
$``$ $`{\displaystyle \frac{\gamma }{R^3}}{\displaystyle \frac{2|n|}{(|n|+1)}}(1{\displaystyle \frac{nm}{2}}{\displaystyle \frac{3m^2}{2}})\},`$
$$G_N(n,m)=\frac{|n|}{R}\left[1\mathrm{sgn}(nm)\frac{1}{|n|}\right],$$
(48)
$$G_{NN}(n,m)=\frac{|n|}{R}\dot{R}^2\left[\frac{2\alpha (m)}{m}\frac{|n|}{n}\frac{2\beta (n)}{|n|}+\frac{f(n,m)}{2(|n|+1)}\right],$$
(49)
$$H_{NN}(n,m)=\frac{|n|}{2(|n|+1)}\frac{\dot{R}^2}{R}f(n,m),$$
(50)
and

$$J_{NN}(n,m)=\frac{|n|}{2(|n|+1)}\dot{R}f(n,m),$$
(51)
where

$$f(n,m)=\left[\frac{|n|}{n}h(n,m)+k(n,m)\right].$$
(52)


format: text

text:
# I Introduction.

## I Introduction.

The Skyrme modelSkyrme  was proposed in the sixties as a model for the strong interactions of hadrons and was very successful in describing nucleons as quantum states of the chiral soliton in original and generalized Skyrme model NikNov . The Skyrme model gives us very unusual instrument to study new physics especially in the light nuclei region. In this region traditional one nucleon degrees of freedom are possibly not so important as finite sizes of nucleon, comparable to the nuclear radiuses Nik , hep-ph0109192 . There is no analytic solutions for the Skyrme model equations of motion. We still have to use variational approaches. The most popular in between them is the so called rational map ansatz Mant  leading to the a number of the solution with discrete space symmetries and topological charges corresponding to light nuclei atomic numbers up to 22. They are very like to fullerene structures more usual for the larger molecular scale Battye . In any way such solutions are like pure numerical solutions obtained in Braaten0  for topological charges 2,3,4,5 and 6.

Here we try to search soliton with axially symmetric baryon charge distribution. Quantization procedure for the states with baryon number equal to 2,3 and 4 was worked out in Weigel  -Kozhevnikov  without vibrations have been taken into account and including the breathing modeNikolaev1  -Nikolaev2 . The Nik.Sov.23  and NikJaf  describe quantization rules for axially symmetric soliton we are considering here.

The variational ansatz we use here was proposed independently in Nikolaev3 ,Kurihara  and Sorace ). The ansatz being very simple, gives the possibility to do analytical analysis of a part of the nuclear problem.

In this paper we present the results of our variational calculations of the classical soliton structure with baryon charge $`B=12`$ in the framework of the original $`SU(2)`$ Skyrme model. After the quantization procedure some of these solitons could be identified with shape isomers of $`C^{12}`$.

## II Ansatz for the Static Solutions.

We follow our papers Nikolaev5  and Nik24  with some modifications. In variational form of the chiral field $`U`$:

$`U(\stackrel{}{r})=cosF(r)+i(\stackrel{}{\tau }\stackrel{}{N})sinF(r).`$ (1)
we use the next general assumption about the configuration of the isotopic vector field $`\stackrel{}{N}`$ for axially symmetric soliton:

$`\stackrel{}{N}=\{cos(\mathrm{\Phi }(\varphi ,\theta ))sin(T(\theta )),sin(\mathrm{\Phi }(\varphi ,\theta ))sin(T(\theta )),cos(T(\theta ))\}.`$ (2)
In eq.(2) $`\mathrm{\Phi }(\varphi ),T(\theta )`$ are some arbitrary functions of angles $`(\theta ,\varphi )`$ of the vector $`\stackrel{}{r}`$ in the spherical coordinate system.

## III Mass Functional and Solutions for Static Equations.

After some algebra (1), (2) and the Lagrangian density $``$ for the stationary solution

$`={\displaystyle \frac{F_\pi ^2}{16}}Tr(L_kL_k)+{\displaystyle \frac{1}{32e^2}}Tr[L_k,L_i]^2,`$ (3)
expressed through the left currants $`L_k=U^+_kU`$ lead to the expression

$`=_2+_4,`$ (4)
where

$$_2=\frac{F_\pi ^2}{8}\left\{\left(\frac{F}{x}\right)^2+\left[\left(\frac{sinT}{sin\theta }\frac{\mathrm{\Phi }}{\varphi }\right)^2+\left(\frac{T}{\theta }\right)^2+sin^2T\left(\frac{\mathrm{\Phi }}{\theta }\right)^2\right]\frac{sin^2F}{r^2}\right\}$$

and

$`_4`$ $`=`$ $`{\displaystyle \frac{1}{2e^2}}{\displaystyle \frac{sin^2F}{r^2}}\{{\displaystyle \frac{sin^2T}{sin^2\theta }}\left({\displaystyle \frac{T}{\theta }}\right)^2\left({\displaystyle \frac{\mathrm{\Phi }}{\varphi }}\right)^2{\displaystyle \frac{sin^2F}{r^2}}`$ (5)
$`+`$ $`[{\displaystyle \frac{sin^2T}{sin^2\theta }}\left({\displaystyle \frac{\mathrm{\Phi }}{\varphi }}\right)^2+\left({\displaystyle \frac{T}{\theta }}\right)^2+sin^2T\left({\displaystyle \frac{\mathrm{\Phi }}{\theta }}\right)^2]\left({\displaystyle \frac{F}{x}}\right)^2\}`$
The variation of the functional $`L=𝑑\stackrel{}{r}`$ with respect to $`\mathrm{\Phi }`$ leads to an equation which has a solution of the type

$`\mathrm{\Phi }(\varphi )=k(\theta )\varphi +Const`$ (6)
with a constrain:

$`{\displaystyle \frac{}{\theta }}\left[sin^2T(\theta )sin\theta {\displaystyle \frac{k(\theta )}{\theta }}\right]=0.`$ (7)
It is easily seen from eq.Nikolaev5  that function $`k(\theta )`$ may be piecewise constant function (step function) in general case:

$`\mathrm{\Phi }(\theta ,\varphi )=\{\begin{array}{c}k^{(1)}\varphi +\rho ^{(1)},\mathrm{for}0\theta <\theta _1,\\ \\ k^{(2)}\varphi +\rho ^{(2)},\mathrm{for}\theta _1\theta <\theta _2,\\ \\ ...\\ \\ k^{(l)}\varphi +\rho ^{(l)},\mathrm{for}\theta _{l1}\theta <\pi .\end{array}`$
Moreover $`k(\theta )`$ must be integer in any region $`\theta _m\theta \theta _{m+1}`$, where $`\theta _m`$, $`\theta _{m+1}`$ are successive points of discontinuity. The positions of these points are determined by the condition

$`T(\theta _m)=m\pi ,`$ (8)
with integer $`m`$, as follows from eq.(7).

Now we have the following expression for the mass of the soliton

$`M=\gamma [aA+bB+C],`$ (9)
where $`\gamma =\pi F_\pi /e`$ and $`x=F_\pi er`$ and the $`a,b`$ and $`A,B,C`$ are the following integrals:

$`a={\displaystyle \underset{0}{\overset{\pi }{}}}\left[k^2{\displaystyle \frac{sin^2T}{sin^2\theta }}+(T^{})^2\right]sin\theta 𝑑\theta ,b=k^2{\displaystyle \underset{0}{\overset{\pi }{}}}{\displaystyle \frac{sin^2T}{sin^2\theta }}(T^{})^2sin\theta 𝑑\theta ,`$ (10)
$`A={\displaystyle \underset{0}{\overset{\mathrm{}}{}}}sin^2F\left[{\displaystyle \frac{1}{4}}+(F^{})^2\right]𝑑x,B={\displaystyle \underset{0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{sin^4F}{x^2}}𝑑x,C={\displaystyle \frac{1}{2}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}(F^{}x)^2𝑑x.`$ (11)
Here we use the symbol prime to denote the following derivatives

$`\mathrm{\Phi }^{}={\displaystyle \frac{\mathrm{\Phi }}{\varphi }};T^{}={\displaystyle \frac{T}{\theta }};F^{}={\displaystyle \frac{F}{r}}`$ (12)
We consider the configurations with finite masses. The only configurations which obey the finiteness of mass condition are the configurations with $`F(0)=n\pi `$ where $`n`$-is some integer number. Without loss of generality we take $`F(\mathrm{})=0`$. As it was shown inNikolaev5  $`T(\theta )`$ has the following behaviour near the boundary of the domain of its definition

$`T(\theta )\theta ^k,for\theta 0;T(\theta )\pi l(\pi \theta )^k,for\theta \pi .`$ (13)
Here $`l`$ is an integer number. Thus we have the following estimation for the number of discontinuity points $`d`$:

$`0dl1.`$ (14)
Now all solutions $`U_{l\{k_i,n_i\}}`$ are classified by a set of integer numbers $`l`$, $`k_0,\mathrm{},k_{l1}`$ and $`n_0,\mathrm{},n_{l1}`$. The functions $`F(x)`$ and $`T(\theta )`$ have to obey the equations (14,15) fromNikolaev5  in arbitrary space region with given number $`k`$.

## IV Baryon Charge Distribution and the Soliton Structure.

Now consider more carefully the structure of solitons. For that purpose let us calculate the baryon charge density

$`J_0^B(\stackrel{}{r})={\displaystyle \frac{1}{24\pi ^2}}ϵ_{0\mu \nu \rho }Tr(L_\mu L_\nu L_\rho ).`$ (15)
The straightforward calculation gives

$`J_0^B(r,\theta )={\displaystyle \frac{1}{2\pi ^2}}{\displaystyle \frac{sin^2F}{r^2}}{\displaystyle \frac{dF}{dr}}{\displaystyle \frac{sinT}{sin\theta }}{\displaystyle \frac{dT}{d\theta }}{\displaystyle \frac{d\mathrm{\Phi }}{d\varphi }}`$ (16)
Equation (16) immediately results in the expression for the corresponding topological charge

$`B={\displaystyle \underset{m=0}{\overset{l1}{}}}(1)^mn_mk_m`$ (17)
InNikolaev5  we have investigated toroidal multiskyrmion configurations with baryon numbers $`B=1,2,3,4,5`$ and more complicated nontoroidal (including antiskyrmions ($`\overline{S}`$)) configurations.

It is obvious that setting $`k_m<0`$ for even m, $`k_m>0`$ for odd m and $`n_m>0`$ for all m, we obtain configuration with positive baryon charge. In the general case for obtaining configuration with positive baryon charge we must require that

$$n_mk_m>0foroddl,$$

$$n_mk_m<0forevenl.$$

## V The Masses and baryon charge distributions.

Here we reproduce the mass values and baryon charge distributions corresponding to the obtained local minima of the energy functional for the Skyrme field. We have to point out that we discuss multiskyrmion configurations we search for not only classically stable configurations (The decay in two or more skyrmions is forbidden energetically). Nonstable configurations are also in our attention because they may become stable after the quantization procedure Nikolaev5  or pion field Casimir energy would taken into account in full quantum description of the considered solitons.

We restrict ourself to configuration with a symmetric distribution of energy (mass) density in the ($`x,y`$) plane. This mean that from the class of all solutions considered, characterized by the numbers $`l,\{k_m,n_m\}|_1^l`$, we choose only the solution satisfying the condition

$`k_i=k_{l+1i},n_i=n_{l+1i}\mathrm{for}\mathrm{odd}l,`$ (18)
$`k_i=k_{l+1i},n_i=n_{l+1i}\mathrm{for}\mathrm{even}l.`$ (19)
In Table1 we present the masses of shape isomers of $`C^{12}`$. The calculated soliton masses are given in $`(\pi F_\pi /e)`$ units.

From Table 1 one see that in calculated part of spectrum all configurations have very different structures. Presence of such isomers could probably be seen in high energy ion-ion scattering experiments.

The reason we are looking for possible not pure toroidal solitons is that there is a number of solutions with smaller masses among them. For example a configuration composed from three toroidal multibaryons 3{4.1-4.1-4.1} has a smaller mass. Now this state can not decay into 12 classical skyrmions with $`B=1`$ ($`M_{1\{1.1\}}=11.60608\pi F_\pi /e`$) or into three toroidal skyrmion with $`B=4`$ ($`M_{1\{4.1\}}=47.67478\pi F_\pi /e`$). We point such a skyrmion as classically stable configuration. The configuration 3{2.2-2.2-2.2} do not obey the condition of classical stability.

For calculated axially symmetric configurations we present baryon density distributions integrated on $`d\mathrm{\Omega }=\mathrm{sin}(\theta )d\varphi d\theta `$ at Figure 1. Here we use dimensionless coordinates $`x=F_\pi er`$.

In according to our calculations the solitons from the same topological sector can have strongly different masses if they have different structure. We also have to point out that a number of the states have shell like structure. There are states which has $`n_i1`$. Such shape isomer can give different specific contribution to physics processes in light nuclei.

## VI Conclusions.

The axially symmetric solitons with baryon number $`B=12`$ have been investigated in the framework of the very general assumption about the form of the solution of the Skyrme model equations. The obtained solitons could be seen in nuclear reactions as isomer contributions in reactions involving $`C^{12}`$. Such isomers correspond to different form of baryon density distribution. We have to point out that used ansatz leads to stable solitons with $`B=12`$ and shell like structure of the baryon density distribution.

## VII Acknowledgements.

We thank A.N. Antonov and V.K. Lukyanov for useful discussions of obtained results.
This work is supported in part by the programme ”University of Russia” UR.02.01.020


text:
# Effective theory for close limit of two branes

## I Introduction

The purpose of this paper is to write down the effective theory for two branes system at the closed limit. This kind of system is a situation just before the collision of two branes. Since the brane is a key object in superstring theory and braneworld scenarios, the collision process of branes definitely could be an important and fundamental process. In a braneworld the brane-brane collision could provide us a new picture of the creation of the big-bang universe. A brane-bubble collision also will be interesting because it might give us a new scenario of the inflation and reheating process. Thus it is worth deriving an effective theory on the brane for close limit of two branes.

To derive an effective theory on the brane, we must solve the bulk spacetime. In the ordinary Kaluza-Klein theory it was a trivial task. In the braneworld, it is quite non-trivial because the gravity on the brane is non-trivially coupled to the gravitational fields in the bulk. Recently it has been reported that we can do so for a two branes system if we employ the derivative expansion scheme, that is, long wave approximations. The two junction conditions on the branes gives us sufficient boundary conditions to determine the bulk geometry. Thus, the bulk geometry could be uniquely determined. Finally, a low energy effective theory can be derived (See also Refs.  and  for earlier works and Refs.  for the other issues.). In this paper, on the other hand, we point out that we can apply the Taylor expansion to the closed limit of two branes. In this method, the high energy effects can be easily seen, which are difficult to deal with in low energy approximations. The price of it is that the separation between two branes should be small.

In this paper we will use a toy model to concentrate on the dealing of the gravitational effect and dynamics for two branes system of the Randall-Sundrum model where the bulk stress energy is a cosmological constant. The brane will be treated as the thin domain wall. In this approximation, we will be able to use the Israel junction formalism. However a problem comes out here. From the derivation it is obvious that the Israel junction condition cannot be applied to just a moment of the collision of two branes with $`Z_2`$-symmetry. That is, under $`Z_2`$-symmetry, the size of the extra dimension shrinks to zero at the collision. So it is fair to say that our present analysis cannot treat the collision in Randall-Sundrum type models. Our formalism can be applied just before the collision of two branes. In order to treat the collision, we need the interaction between branes, that is, the brane field action which will be reduced to something such as the Born-Infeld action at the thin wall limit. This is, however, beyond the scope of present study.

In this paper we will derive the effective theory describing the close limit of two branes and discuss the radion dynamics, which expresses the evolution of the brane separation, and its relation to the bulk geometry. We also give an insight into the braneworld black holes. Our present study can be regarded as the generalization of that by Khoury and Zhang who considered the Friedman brane cosmology in an empty bulk and discussed the meaning of dark radiation. We will work in the covariant curvature formalism which is most useful to keep the four dimensional covariance and obtain the full view.

The rest of this paper is organised as follows. In Sec. II, we review the covariant curvature formalism. In Sec. III, we solve the bulk geometry by Taylor expansion and then derive the four dimensional effective theory at the close limit. In Sec. III the radion dynamics and the braneworld black holes will be addressed. Finally we will give the summary and discussion in Sec. IV.

## II Covariant curvature formalism

In this section we review the covariant curvature formalism developed in Refs. . Note that the central issue is how to solve the bulk seriously. Without this step we cannot obtain a correct effective theory on the brane.

We employ the metric form

$`ds^2=e^{2\varphi (x)}dy^2+q_{\mu \nu }(y,x)dx^\mu dx^\nu .`$ (1)
In the above it is supposed that the positive and negative tension branes are located at $`y=0`$ and $`y=y_0`$, respectively. The proper distance between two branes is given by $`d_0(x)=_0^{y_0}e^{\varphi (x)}𝑑y=y_0e^{\varphi (x)}`$. $`q_{\mu \nu }(y,x)`$ is the induced metric of $`y=`$const hypersurfaces. At the close limit of two branes we expect that the coordinate used in the metric of Eq. (1) does not break down for a wide class of the braneworld spacetimes.

From the Gauss-Codacci equations, we have two key equations

$`{}_{}{}^{(4)}G_{\nu }^{\mu }`$ $`=`$ $`{\displaystyle \frac{3}{\mathrm{}^2}}\delta _\nu ^\mu +KK_\nu ^\mu K_\alpha ^\mu K_\nu ^\alpha `$ (3)
$`{\displaystyle \frac{1}{2}}\delta _\nu ^\mu (K^2K_\beta ^\alpha K_\alpha ^\beta )E_\nu ^\mu `$
and

$`D_\mu K_\nu ^\mu D_\nu K=0,`$ (4)
where $`D_\mu `$ is the covariant derivative with respect to $`q_{\mu \nu }`$. $`\mathrm{}`$ is the bulk curvature radius. $`{}_{}{}^{(4)}G_{\mu \nu }^{}`$ is the four-dimensional Einstein tensor with respect to $`q_{\mu \nu }`$ and $`K_{\mu \nu }`$ is the extrinsic curvature of $`y=`$constant hypersurfaces defined by

$`K_{\mu \nu }={\displaystyle \frac{1}{2}}\text{£}_nq_{\mu \nu }=_\mu n_\nu +n_\mu D_\nu \varphi ,`$ (5)
where $`n=e^\varphi _y`$. Here note that $`a^\mu =n^\nu _\nu n^\mu =D^\mu \varphi (x)`$. $`E_{\mu \nu }`$ is a part of the projected Weyl tensor defined by

$`E_\nu ^\mu `$ $`=`$ $`{}_{}{}^{(5)}C_{\mu \alpha \nu \beta }^{}n^\alpha n^\beta `$ (6)
$`=`$ $`D^\mu D_\nu \varphi D^\mu \varphi D_\nu \varphi \text{£}_nK_\nu ^\mu K_\alpha ^\mu K_\nu ^\alpha `$ (8)
$`+{\displaystyle \frac{1}{\mathrm{}^2}}\delta _\nu ^\mu ,`$
where $`{}_{}{}^{(5)}C_{\mu \nu \alpha \beta }^{}`$ is the five-dimensional Weyl tensor.

The junction conditions on the branes give us the relation between the extrinsic curvatures and the energy-momentum tensor on the branes:

$`[K_\nu ^\mu \delta _\nu ^\mu K]_{y=0}={\displaystyle \frac{\kappa ^2}{2}}\left(\sigma _1\delta _\nu ^\mu +T_{1\nu }^\mu \right)`$ (9)
and

$`[K_\nu ^\mu \delta _\nu ^\mu K]_{y=y_0}={\displaystyle \frac{\kappa ^2}{2}}\left(\sigma _2\delta _\nu ^\mu +T_{2\nu }^\mu \right).`$ (10)
$`T_{1\nu }^\mu `$ and $`T_{2\nu }^\mu `$ are the energy-momentum tensor localised on the positive and negative branes. $`\sigma _1`$ and $`\sigma _2`$ are the brane tensions. If one substitutes the above conditions to Eq. (3), we might be able to derive the Einstein equation on the brane. However, it turns out that $`E_{\mu \nu }`$ is not still unknown, because it is a five dimensional quantity. For the single brane $`E_{\mu \nu }`$ is identical to the contribution from the Kaluza-Klein modes in the linealised theory and it will be higher order correction at large distance, that is, $`E_{\mu \nu }`$ is negligible at low energy in the linearlised theory. For the two branes system, on the other hand, $`E_{\mu \nu }`$ is not negligible as seen later. Anycase, we need the evolutional equation for $`E_{\mu \nu }`$ to the bulk.

To evaluate $`E_{\mu \nu }`$ in the bulk, we derive its evolutional equation. The result is given by

$`\text{£}_nE_{\alpha \beta }=D^\mu B_{\mu (\alpha \beta )}+K^{\mu \nu }{}_{}{}^{(4)}C_{\mu \alpha \nu \beta }^{}+4K_{(\alpha }^\mu E_{\beta )\mu }`$ (11)
$`{\displaystyle \frac{3}{2}}KE_{\alpha \beta }{\displaystyle \frac{1}{2}}q_{\alpha \beta }K^{\mu \nu }E_{\mu \nu }+2D^\beta \varphi B_{\beta (\mu \nu )}`$ (12)
$`+2\stackrel{~}{K}_\alpha ^\mu \stackrel{~}{K}_{\mu \nu }\stackrel{~}{K}_\beta ^\nu {\displaystyle \frac{7}{6}}\stackrel{~}{K}_{\mu \nu }\stackrel{~}{K}^{\mu \nu }\stackrel{~}{K}_{\alpha \beta }`$ (13)
$`{\displaystyle \frac{1}{2}}q_{\mu \nu }\stackrel{~}{K}_{\alpha \beta }\stackrel{~}{K}_\rho ^\alpha \stackrel{~}{K}^{\rho \beta },`$ (14)
where $`B_{\mu \nu \alpha }=q_\mu ^\rho q_\nu ^\sigma {}_{}{}^{(5)}C_{\rho \sigma \alpha \beta }^{}n^\beta `$ and $`\stackrel{~}{K}_{\mu \nu }=K_{\mu \nu }\frac{1}{4}q_{\mu \nu }K`$. Since the right-hand side contains $`B_{\mu \nu \alpha }`$, $`K_{\mu \nu }`$, and $`{}_{}{}^{(4)}C_{\mu \nu \alpha \beta }^{}`$, we need their evolutional equations too. After a long calculation we obtain

$`\text{£}_n{}_{}{}^{(4)}R_{\mu \nu \alpha \beta }^{}+2{}_{}{}^{(4)}R_{\mu \nu \rho [\alpha }^{}K_{\beta ]}^\rho +2D_{[\mu }B_{|\alpha \beta |\nu ]}`$ (15)
$`+2(D_\mu D_{[\alpha }\varphi +D_\mu \varphi D_{[\alpha }\varphi )K_{\beta ]\nu }`$ (16)
$`2(D_\nu D_{[\alpha }\varphi D_\nu \varphi D_{[\alpha }\varphi )K_{\beta ]\mu }`$ (17)
$`2B_{\alpha \beta [\mu }D_{\nu ]}\varphi 2B_{\mu \nu [\alpha }D_{\beta ]}\varphi =0,`$ (18)
$`\text{£}_nB_{\mu \nu \alpha }+2D_{[\mu }E_{\nu ]\alpha }+2D_{[\mu }\varphi E_{\nu ]\alpha }`$ (19)
$`B_{\mu \nu \beta }K_\alpha ^\beta +2B_{\alpha \beta [\mu }K_{\nu ]}^\beta `$ (20)
$`+({}_{}{}^{(4)}R_{\mu \nu \alpha \beta }^{}K_{\mu \alpha }K_{\nu \beta }+K_{\mu \beta }K_{\nu \alpha })D^\beta \varphi =0,`$ (21)
and

$`e^\varphi _yK_\nu ^\mu `$ $`=`$ $`D^\mu D_\nu \varphi D^\mu \varphi D_\nu \varphi K_\alpha ^\mu K_\nu ^\alpha +{\displaystyle \frac{1}{\mathrm{}^2}}\delta _\nu ^\mu `$ (23)
$`E_\nu ^\mu .`$
Equation (23) is a rearrangement of Eq. (8).

The junction condition directly implies the boundary condition on the branes for $`K_{\mu \nu }`$ and $`B_{\mu \nu \alpha }`$ because of

$`B_{\mu \nu \alpha }=2D_{[\mu }K_{\nu ]\alpha }.`$ (24)
For the discussion in the next section, we write down the second derivative of $`K_\nu ^\mu `$:

$`_y^2K_\nu ^\mu `$ $`=`$ $`e^{2\varphi }[K_\nu ^\alpha D^\mu \varphi D_\alpha \varphi +K_\alpha ^\mu D_\nu \varphi D^\alpha \varphi `$ (29)
$`+3K_\alpha ^\mu (D^\alpha D_\nu \varphi +D^\alpha \varphi D_\nu \varphi )`$
$`+K_\nu ^\alpha (D^\mu D_\alpha \varphi +D^\mu \varphi D_\alpha \varphi )K_\nu ^\alpha E_\alpha ^\mu +K_\alpha ^\mu E_\nu ^\alpha `$
$`+{\displaystyle \frac{3}{2}}KE_\nu ^\mu +{\displaystyle \frac{1}{2}}\delta _\nu ^\mu K_\beta ^\alpha E_\alpha ^\beta +D_\alpha \varphi D^\mu K_\nu ^\alpha K_\nu ^\mu (D\varphi )^2`$
$`+B_{\nu \alpha }^\mu D^\alpha \varphi D^\alpha \varphi (B_{\alpha \nu }^\mu +B_{\alpha \nu }^\mu )+\chi _\nu ^\mu ]`$
where $`\chi _\nu ^\mu `$ is defined by

$`\chi _\nu ^\mu `$ $`=`$ $`{\displaystyle \frac{1}{2}}(D^\alpha B_{\alpha \nu }^\mu +D^\alpha B_{\alpha \nu }^\mu )K^{\alpha \beta }{}_{}{}^{(4)}C_{\alpha \nu \beta }^{\mu }`$ (32)
$`+2K_\alpha ^\mu K_\beta ^\alpha K_\nu ^\beta {\displaystyle \frac{2}{\mathrm{}^2}}K_\nu ^\mu 2\stackrel{~}{K}_\alpha ^\mu \stackrel{~}{K}_\beta ^\alpha \stackrel{~}{K}_\nu ^\beta `$
$`+{\displaystyle \frac{7}{6}}\stackrel{~}{K}_\nu ^\mu \stackrel{~}{K}_\beta ^\alpha \stackrel{~}{K}_\alpha ^\beta +{\displaystyle \frac{1}{2}}\delta _\nu ^\mu \stackrel{~}{K}_\beta ^\alpha \stackrel{~}{K}_\alpha ^\rho \stackrel{~}{K}_\rho ^\beta .`$
See Ref.  for the derivation of the low energy effective theory by using the present covariant formalism. In Ref.  we employed the derivative expansion scheme which corresponds to the long-wave approximation, that is, the typical scale on the brane is much larger than the bulk curvature scale. In other words, the energy scale of the energy momentum tensor localized on the branes should be sufficiently lower than the tension of the branes. As said before, in this paper, we will focus on the close limit and employ the Taylor expansion assuming that the distance between two branes $`d_0`$ is much shorter than the bulk curvature scale $`\mathrm{}`$. In this method, it is not necessary to take the Randall-Sundrum tuning of the tensions and low energy approximation $`\kappa ^2T_\nu ^\mu /\sigma 1`$.

## III Effective theory at closed limit

In this section, we solve the bulk geometry using the Taylor expansion in the covariant curvature formalism. What we will do is the evaluation of $`E_{\mu \nu }`$, that is, the writing down of $`E_{\mu \nu }`$ in terms of the four dimensional quantities on the branes. Then we derive the effective equation on the branes. This theory can be applicable to the non-linear regime as long as we consider the close limit.

The point of the derivation is simple. We consider the Taylor expansion of the extrinsic curvature, $`K_\nu ^\mu (y_0,x)=K_\nu ^\mu (0,x)+_yK_\nu ^\mu (0,x)y_0+\mathrm{}`$. The derivative of the extrinsic curvature $`_yK_{\mu \nu }`$ contains the bulk informations($`E_{\mu \nu }`$), and then the junction condition determines them in terms of the four dimensional quantity. This is the specialty in two brane systems. If one considers a single-brane system, we need a boundary condition in the bulk somewhere such as Cauchy horizon.

In this paper we will think of the positive tension brane. The extension of the present discussion to the negative tension brane is very easy. We begin with the Taylor expansion of the extrinsic curvature $`K_\nu ^\mu `$ around the positive tension brane

$`K_\nu ^\mu (y_0,x)`$ $`=`$ $`K_\nu ^\mu (0,x)+_yK_\nu ^\mu (0,x)y_0`$ (34)
$`+{\displaystyle \frac{1}{2}}_y^2K_\nu ^\mu (0,x)y_0^2+O(y^3).`$
Using Eqs (23) and (29), Eq. (34) becomes

$`K_\nu ^\mu (y_0,x)`$ $`=`$ $`K_\nu ^\mu (0,x)d_0[{\displaystyle \frac{D^\mu D_\nu d_0}{d_0}}+K_\alpha ^\mu K_\nu ^\alpha +E_\nu ^\mu `$ (41)
$`{\displaystyle \frac{1}{\mathrm{}^2}}\delta _\nu ^\mu ]+{\displaystyle \frac{1}{2}}d_0^2[\chi _\nu ^\mu +3{\displaystyle \frac{D^\alpha D_\nu d_0}{d_0}}K_\alpha ^\mu `$
$`+{\displaystyle \frac{D_\alpha D^\mu d_0}{d_0}}K_\nu ^\alpha +K_\nu ^\alpha {\displaystyle \frac{D^\mu d_0D_\alpha d_0}{d_0^2}}`$
$`+K_\alpha ^\mu {\displaystyle \frac{D_\nu d_0D^\alpha d_0}{d_0^2}}K_\nu ^\mu \left({\displaystyle \frac{Dd_0}{d_0}}\right)^2`$
$`+{\displaystyle \frac{D^\alpha d_0}{d_0}}\left[D^\mu K_{\alpha \nu }+B_{\alpha \nu }^\mu +2B_{\nu \alpha }^\mu \right]`$
$`K_\nu ^\alpha E_\alpha ^\mu +K_\alpha ^\mu E_\nu ^\alpha +{\displaystyle \frac{3}{2}}KE_\nu ^\mu +{\displaystyle \frac{1}{2}}\delta _\nu ^\mu E_\beta ^\alpha K_\alpha ^\beta ]`$
$`+O(d_0^3),`$
where $`K_\nu ^\mu `$ and $`E_\nu ^\mu `$ in the right-hand side is evaluated at $`y=0`$.

Let us regard Eq. (41) as the equation for $`E_\nu ^\mu `$. To solve the equation, we assume that two branes are close enough to satisfy the conditions

$$d_0/\mathrm{},\kappa ^2T_1d_0,\kappa ^2T_2d_0<1.$$
(42)
Then we expand $`E_\nu ^\mu `$ as

$`E_\nu ^\mu =\stackrel{\left(1\right)}{E_\nu ^\mu }{\displaystyle \frac{1}{d_0}}+\stackrel{\left(0\right)}{E_\nu ^\mu }+O(d_0).`$ (43)
Substituting the above into Eq. (41) and solving iteratively we see that $`\stackrel{\left(1\right)}{E_\nu ^\mu }`$ and $`\stackrel{\left(0\right)}{E_\nu ^\mu }`$ has the following form:

$`\stackrel{\left(1\right)}{E_\nu ^\mu }=K_\nu ^\mu (y_0,x)K_\nu ^\mu +D^\mu D_\nu d_0`$ (44)
and

$`\stackrel{\left(0\right)}{E_\nu ^\mu }`$ $`=`$ $`K_\alpha ^\mu K_\nu ^\alpha {\displaystyle \frac{1}{\mathrm{}^2}}\delta _\nu ^\mu {\displaystyle \frac{1}{2}}(K_\nu ^\alpha \stackrel{\left(1\right)}{E_\alpha ^\mu }+K_\alpha ^\mu \stackrel{\left(1\right)}{E_\nu ^\alpha }`$ (50)
$`+{\displaystyle \frac{3}{2}}K\stackrel{\left(1\right)}{E_\nu ^\mu }+{\displaystyle \frac{1}{2}}\delta _\nu ^\mu K_\beta ^\alpha \stackrel{\left(1\right)}{E_\alpha ^\beta })`$
$`{\displaystyle \frac{1}{2}}(3K_\alpha ^\mu D^\alpha D_\nu d_0+K_\nu ^\alpha D_\alpha D^\mu d_0)`$
$`{\displaystyle \frac{D^\mu d_0D_\alpha d_0}{2d_0}}K_\nu ^\alpha {\displaystyle \frac{D_\mu d_0D^\alpha d_0}{2d_0}}K_\alpha ^\mu `$
$`+{\displaystyle \frac{1}{2}}K_\nu ^\mu {\displaystyle \frac{(Dd_0)^2}{d_0}}{\displaystyle \frac{1}{2}}D_\alpha d_0D^\mu K_\nu ^\alpha `$
$`+{\displaystyle \frac{1}{2}}D^\alpha d_0(2B_{\alpha \nu }^\mu B_{\alpha \nu }^\mu ).`$
In the above we supposed that $`D^2d_0`$ is the same order as $`K`$. We should notice that it is not necessary to impose low energy conditions $`T/\sigma \kappa ^2T\mathrm{}<1`$. What we imposed here is $`\kappa ^2Td_0<1`$.

From the fact that $`E_\nu ^\mu `$ is traceless we can have the equation for $`d_0`$

$`{\displaystyle \frac{1}{d_0}}\left(1+{\displaystyle \frac{d_0}{\mathrm{}_1}}\right)D^2d_0{\displaystyle \frac{4}{d_0}}\left({\displaystyle \frac{1}{\mathrm{}_2}}{\displaystyle \frac{1}{\mathrm{}_1}}\right)+4\left({\displaystyle \frac{1}{\mathrm{}_1\mathrm{}_2}}{\displaystyle \frac{1}{\mathrm{}^2}}\right)`$ (51)
$`{\displaystyle \frac{\kappa ^2}{6d_0}}\left[\left(1+{\displaystyle \frac{d_0}{\mathrm{}_2}}\right)T_1+\left(1{\displaystyle \frac{d_0}{\mathrm{}_1}}\right)T_2\right]`$ (52)
$`+{\displaystyle \frac{\kappa ^2}{2}}\left(T_{1\beta }^\alpha {\displaystyle \frac{1}{3}}\delta _\beta ^\alpha T_1\right)D^\beta D_\alpha d_0`$ (53)
$`+{\displaystyle \frac{\kappa ^2}{2}}\left(T_{1\beta }^\alpha {\displaystyle \frac{1}{6}}\delta _\beta ^\alpha T_1\right){\displaystyle \frac{D^\beta d_0D_\alpha d_0}{d_0}}`$ (54)
$`{\displaystyle \frac{(Dd_0)^2}{\mathrm{}_1d_0}}{\displaystyle \frac{\kappa ^4}{4}}\left(T_{1\beta }^\alpha T_{2\alpha }^\beta {\displaystyle \frac{2}{9}}T_1T_2\right)`$ (55)
$`{\displaystyle \frac{\kappa ^2}{12}}D^\alpha d_0D_\alpha T_1=0,`$ (56)
where we defined

$`{\displaystyle \frac{1}{\mathrm{}_1}}={\displaystyle \frac{1}{6}}\kappa ^2\sigma _1\mathrm{and}{\displaystyle \frac{1}{\mathrm{}_2}}={\displaystyle \frac{1}{6}}\kappa ^2\sigma _2.`$ (57)
Since we write down $`E_\nu ^\mu `$ in terms of four dimensional quantities, we can derive the effective gravitational equation on the brane

$`{}_{}{}^{(4)}G_{\nu }^{\mu }`$ $`=`$ $`6\left({\displaystyle \frac{1}{\mathrm{}^2}}{\displaystyle \frac{1}{\mathrm{}_1^2}}\right)\delta _\nu ^\mu {\displaystyle \frac{9}{\mathrm{}_1}}\left({\displaystyle \frac{1}{\mathrm{}_2}}{\displaystyle \frac{1}{\mathrm{}_1}}\right)\delta _\nu ^\mu `$ (74)
$`+{\displaystyle \frac{3}{d_0}}\left({\displaystyle \frac{1}{\mathrm{}_2}}{\displaystyle \frac{1}{\mathrm{}_1}}\right)\delta _\nu ^\mu +{\displaystyle \frac{\kappa ^2}{2d_0}}(T_{1\nu }^\mu +T_{2\nu }^\mu )`$
$`+{\displaystyle \frac{3\kappa ^2}{2}}\left({\displaystyle \frac{1}{\mathrm{}_2}}{\displaystyle \frac{1}{\mathrm{}_1}}\right)T_1\delta _\nu ^\mu +{\displaystyle \frac{\kappa ^2}{2\mathrm{}_1}}(T_{1\nu }^\mu 3T_{2\nu }^\mu )`$
$`+\left(1{\displaystyle \frac{d_0}{\mathrm{}_1}}\right){\displaystyle \frac{D^\mu D_\nu d_0\delta _\nu ^\mu D^2d_0}{d_0}}`$
$`+{\displaystyle \frac{2D^\mu d_0D_\nu d_0+\delta _\nu ^\mu (Dd_0)^2}{2\mathrm{}_1d_0}}`$
$`+\kappa ^2[{\displaystyle \frac{5}{24}}T_1D^\mu D_\nu d_0{\displaystyle \frac{1}{6}}\delta _\nu ^\mu T_1D^2d_0`$
$`+{\displaystyle \frac{1}{2}}(T_{1\alpha }^\mu D^\alpha D_\nu d_0+T_{1\nu }^\alpha D_\alpha D^\mu d_0)`$
$`{\displaystyle \frac{5}{8}}\delta _\nu ^\mu T_{1\beta }^\alpha D^\beta D_\alpha d_0`$
$`+{\displaystyle \frac{1}{4}}T_{1\nu }^\alpha {\displaystyle \frac{D^\mu d_0D_\alpha d_0}{d_0}}+{\displaystyle \frac{1}{4}}T_{1\alpha }^\mu {\displaystyle \frac{D_\nu d_0D^\alpha d_0}{d_0}}`$
$`{\displaystyle \frac{1}{6}}T_1{\displaystyle \frac{D^\mu d_0D_\nu d_0}{d_0}}+{\displaystyle \frac{1}{6}}\delta _\nu ^\mu T_1{\displaystyle \frac{(Dd_0)^2}{d_0}}`$
$`{\displaystyle \frac{1}{4}}T_{1\nu }^\mu {\displaystyle \frac{(Dd_0)^2}{d_0}}{\displaystyle \frac{1}{2}}\delta _\nu ^\mu T_{1\beta }^\alpha {\displaystyle \frac{D^\beta d_0D_\alpha d_0}{d_0}}`$
$`{\displaystyle \frac{1}{4}}D^\alpha d_0(3D_\alpha T_{1\nu }^\mu 2D_\nu T_{1\alpha }^\mu 2D^\mu T_{1\nu \alpha })`$
$`+{\displaystyle \frac{1}{12}}(4\delta _\nu ^\mu D^\alpha d_0D_\alpha T_12D^\mu d_0D_\nu T_1`$
$`2D_\nu d_0D^\mu T_1)]`$
$`+\kappa ^4[{\displaystyle \frac{1}{12}}T_1T_{1\nu }^\mu +{\displaystyle \frac{1}{24}}\delta _\nu ^\mu T_1^2{\displaystyle \frac{1}{8}}\delta _\nu ^\mu T_{1\beta }^\alpha T_{1\alpha }^\beta `$
$`+{\displaystyle \frac{1}{8}}(T_{1\nu }^\alpha T_{2\alpha }^\mu T_{2\nu }^\alpha T_{1\alpha }^\mu )+{\displaystyle \frac{1}{16}}T_1(T_{1\nu }^\mu +T_{2\nu }^\mu )`$
$`+{\displaystyle \frac{3}{16}}\delta _\nu ^\mu T_{1\beta }^\alpha (T_{1\alpha }^\beta +T_{2\alpha }^\beta )].`$
Here note that we derived the above equation under the condition of $`d_0/\mathrm{},\kappa ^2T_1d_0,\kappa ^2T_2d_0<1`$.

The coupling of the matter with the radion is complicated. In the low energy limit such terms can be omitted and after taking of Randall-Sundrum tuning of tensions $`\mathrm{}=\mathrm{}_1=\mathrm{}_2`$, the gravitational equation is reduced to

$`{}_{}{}^{(4)}G_{\nu }^{\mu }`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{2d_0}}(T_{1\nu }^\mu +T_{2\nu }^\mu )+{\displaystyle \frac{\kappa ^2}{2\mathrm{}}}(T_{1\nu }^\mu 3T_{2\nu }^\mu )`$ (77)
$`+\left(1{\displaystyle \frac{d_0}{\mathrm{}}}\right){\displaystyle \frac{D^\mu D_\nu d_0\delta _\nu ^\mu D^2d_0}{d_0}}`$
$`+{\displaystyle \frac{2D^\mu d_0D_\nu d_0+\delta _\nu ^\mu (Dd_0)^2}{2\mathrm{}d_0}}.`$
Here reminded that the first term is much bigger than the second term in the first line of the right-hand side. It is easy to check that the above is identical to the close limit of the low energy effective theory obtained in Refs. .

## IV Applications

In this section we consider the applications of our formalism: the radion dynamics and the braneworld black holes.

### A Radion dynamics and bulk geometry

First, let us remember the work of Khoury and Zhang . They investigated the modified Friedman equation on the brane in empty bulk and observed that $`E_{00}`$ is proportional to the energy density on the brane. It corresponds to Eq (44) in our general analysis of the close limit because we can neglect the effect of bulk curvature for $`d_0\mathrm{}`$. At the leading order of $`d_0/\mathrm{},d_0\kappa ^2T`$, the effective equation on the positive tension brane becomes

$`G_\nu ^\mu `$ $`=`$ $`E_\nu ^\mu ,`$ (78)
$`E_\nu ^\mu `$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{2d_0}}(T_{1\nu }^\mu +T_{2\nu }^\mu )+{\displaystyle \frac{D^\mu D_\nu d_0\delta _\nu ^\mu D^2d_0}{d_0}},`$ (79)
$`D^2d_0`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{6}}(T_1+T_2).`$ (80)
From these effective equations, it is easy to derive the modified Friedmann equation on the brane

$$H^2=H\frac{\dot{d_0}}{d_0}+\frac{\kappa ^2(\rho _1+\rho _2)}{6d_0},$$
(81)
where $`H`$ is the Hubble parameter on the positive tension brane. $`\rho _1`$ and $`\rho _2`$ are energy densities on positive tension brane and negative tension brane, respectively. This agrees with Khoury and Zhang’s equation. From the evolution equation for radion and the conservation law for the energy-momentum tensor, the right-hand side of Eq. (81) can be rewritten as

$$H^2=Ca^4,$$
(82)
where $`a`$ is the scale factor on the positive tension brane and $`C`$ is the integration constant of the radion $`d_0`$. Thus we can understand the behavior of dark radiation in terms of the dynamics of the radion and that the bulk geometry is related to the initial condition for the radion. Hence the present formalism gives us a straightforward way to see what they observed.

### B Radion dynamics for deSitter-like branes

Next we discuss the radion dynamics for the deSitter-like branes. For simplicity we consider the vacuum energy cases

$`\mathrm{}=\mathrm{}_1=\mathrm{}_2`$ (83)
and

$`T_{1\nu }^\mu =\lambda _1\delta _\nu ^\mu T_{2\nu }^\mu =\lambda _2\delta _\nu ^\mu .`$ (84)
Then the branes will be deSitter spacetimes if the maximum symmetry is imposed. But, we do not restrict ourself to the consideration on the exact deSitter branes. In this case the equation for the radion is extremely simplified as

$`\left(1+{\displaystyle \frac{d_0}{L}}\right){\displaystyle \frac{D^2d_0}{d_0}}{\displaystyle \frac{(Dd_0)^2}{Ld_0}}{\displaystyle \frac{\kappa ^4}{9}}\lambda _1\lambda _2`$ (85)
$`+{\displaystyle \frac{2\kappa ^2}{3d_0}}\left[\left(1+{\displaystyle \frac{d_0}{\mathrm{}}}\right)\lambda _1+\left(1{\displaystyle \frac{d_0}{\mathrm{}}}\right)\lambda _2\right]=0,`$ (86)
where $`1/L=1/\mathrm{}+\kappa ^2\lambda _1/6`$. Defining $`\phi =d_0/\mathrm{}`$($`1`$) and the conformaly transformed metric $`g_{\mu \nu }=(1+\phi )^1q_{\mu \nu }`$, the equation becomes

$`D_g^2\phi V^{}(\phi )=0,`$ (87)
where

$`V(\phi )\alpha \phi +\beta \phi ^2+\mathrm{}.`$ (88)
The constants $`\alpha `$ and $`\beta `$ are defined by

$`\alpha ={\displaystyle \frac{2\kappa ^2}{3L}}(\lambda _1+\lambda _2)`$ (89)
and

$`\beta ={\displaystyle \frac{\kappa ^2}{3\mathrm{}}}(\lambda _1\lambda _2)+{\displaystyle \frac{\kappa ^4}{18}}\lambda _1\lambda _2.`$ (90)
The close limit corresponds to $`\phi 0`$ limit. At the very close limit, the potential is approximated by

$`V(\phi )\alpha \phi .`$ (91)
So the brane could collide if $`\alpha <0`$ due to the attractive force. If $`\alpha >0`$ the repulsive force acts between two branes. This feature is same as the result in Ref. . The term proportional to $`\beta `$ describes the correction in the next order. We should emphasize that high-energy terms such as $`\kappa ^2\lambda _1\lambda _2`$ can be included in our method. In the case of $`\alpha =0`$, the next order should be taken into account and the potential becomes

$`V(\phi ){\displaystyle \frac{2\kappa ^2}{3\mathrm{}}}\lambda _1\left(1+{\displaystyle \frac{1}{12}}\kappa ^2\lambda _1\mathrm{}\right)\phi ^2.`$ (92)
From the above we can see that the repulsive force will be induced between two branes when $`\lambda _1>0`$.

In general the potential at the close limit is given by

$`V(\phi ){\displaystyle \frac{\kappa ^2}{6L}}(T_1+T_2)\phi +\mathrm{}.`$ (93)
Thus the attractive force acts if $`T_1+T_2>0`$ and the repulsive force if $`T_1+T_2<0`$. If one thinks of the dust universe, the repulsive force does.

### C An implication into black holes

Let us finally consider brane world black holes in the case of $`\mathrm{}=\mathrm{}_1=\mathrm{}_2`$ and $`T_{1\nu }^\mu =T_{2\nu }^\mu =0`$. In this case the equation for the radion and the gravitational equation become

$`{\displaystyle \frac{1}{d_0}}\left(1+{\displaystyle \frac{d_0}{\mathrm{}}}\right)D^2d_0{\displaystyle \frac{(Dd_0)^2}{\mathrm{}d_0}}=0`$ (94)
and

$`{}_{}{}^{(4)}G_{\nu }^{\mu }`$ $`=`$ $`\left(1{\displaystyle \frac{d_0}{\mathrm{}}}\right){\displaystyle \frac{D^\mu D_\nu d_0\delta _\nu ^\mu D^2d_0}{d_0}}`$ (96)
$`+{\displaystyle \frac{2D^\mu d_0D_\nu d_0+\delta _\nu ^\mu (Dd_0)^2}{2\mathrm{}d_0}}.`$
Up to this order the higher order correction disappears because they are written by the energy-momentum tensor on the branes. The above system is equivalent to the scalar-tensor theory. In the present case it is easy to see that the static and spherical solution on the brane is only Schwarzshild spacetimes due to the Bekenstein theorem. At the close limit, thus, the black holes in the bulk are described by the adS black string solution. However, if the brane distance increase, it is natural to expect that there should be black holes solution confined on the each branes. That is, there might be a phase transition from black holes to the adS black string. Thus we hope we will be able to see this transition using our formalism. This expectation might be supported because the higher order corrections than that we considered here will contain the terms like $`D^\alpha {}_{}{}^{(4)}R_{\mu \nu \alpha \beta }^{}D^\beta \varphi `$ and $`{}_{}{}^{(4)}R_{\mu \alpha \nu \beta }^{}{}_{}{}^{(4)}R_{}^{\alpha \beta }`$ and then the modified solution from the Schwarzshild could be emerged such as the situation in the superstring theory.

## V Summary and discussion

In this paper we deriveed the effective theory at the close limit of two branes and discussed the radion dynamics/braneworld black holes. Our effective theory is complimentary with the low energy effective theory and can be regarded as the generalisation of the Khoury and Zhang’s recent work.

Finally we should stress on the limitation of our formulation. (i)To determine the bulk spacetime we used the Israel’s junction condition at the beginning. Obviously the derivation of the junction condition will be not correct at the just moment of collision of two branes. Our formalism is correct before the collision. (ii)In addition, we tacitly supposed that the coordinates to the extra dimensions works very well. (iii)In order to treat the collision we must think of the thick wall where the brane is described by a field theory and the brane-brane interactions is also given. Since the effective theory describing the thick wall might be reformulated in terms of the noncommutative geometry, the brane collision might be also treated in the effective theory with the noncommutative geometric aspects.

## Acknowledgments

We would like to thank Shoko Hayakawa, Daisuke Ida, Takashi Torii, and Toby Wiseman for fruitful discussions. This work is inspired by unpublished study on single brane models. The work of TS was supported by Grant-in-Aid for Scientific Research from Ministry of Education, Science, Sports and Culture of Japan(No. 13135208, No.14740155 and No.14102004). The work of KK were supported by the JSPS.