text:
# Lattice study on kaon nucleon scattering length in the 𝐼=1 channel

## 1 Introduction

Anisotropic lattices and improved actions have been extensively used in recent lattice QCD calculations. In our previous works, the following tadpole improved gluonic action on anisotropic lattices has been utilized:

$`S=`$ $``$ $`\beta {\displaystyle \underset{i>j}{}}\left[{\displaystyle \frac{5}{9}}{\displaystyle \frac{TrP_{ij}}{\xi u_s^4}}{\displaystyle \frac{1}{36}}{\displaystyle \frac{TrR_{ij}}{\xi u_s^6}}{\displaystyle \frac{1}{36}}{\displaystyle \frac{TrR_{ji}}{\xi u_s^6}}\right]`$ (1)
$``$ $`\beta {\displaystyle \underset{i}{}}\left[{\displaystyle \frac{4}{9}}{\displaystyle \frac{\xi TrP_{0i}}{u_s^2}}{\displaystyle \frac{1}{36}}{\displaystyle \frac{\xi TrR_{i0}}{u_s^4}}\right],`$
where $`P_{0i}`$ and $`P_{ij}`$ represents the usual temporal and spatial plaquette variable, respectively. $`R_{ij}`$ and $`R_{i0}`$ designates the $`2\times 1`$ spatial and temporal Wilson loops, where, in order to eliminate the spurious states, we have restricted the coupling of fields in the temporal direction to be within one lattice spacing. The parameter $`u_s`$, which is taken to be the fourth root of the average spatial plaquette value, implements the tadpole improvement. With the tadpole improvement factor, the renormalization of the anisotropy (or aspect ratio) $`\xi =a_s/a_t`$ would be small. In this work, we will not differentiate the bare aspect ratio and the renormalized one in this work. Using this action, glueball and hadron spectra have been studied within quenched approximation . In our previous works, configurations generated from this improved action have also been utilized to calculate the $`\pi \pi `$ scattering lengths in the $`I=2`$ channel . In this paper, we extend our study to the kaon nucleon scattering lengths within quenched approximation. Unlike the case in pion pion scattering where the chiral expansion is quite reliable, chiral predictions in $`KN`$ scattering suffer from large corrections . On the experimental side, Martin used dispersion relation fits to determine the $`KN`$ scattering parameters . Therefore, a lattice QCD calculation will offer an important and independent check on these results.

Lattice calculations of hadron scattering lengths have been performed by various authors using conventional Wilson fermions on symmetric lattices without the improvement . In Ref. , a systematic study of the hadron scattering was performed, including pion-pion, kaon-nucleon, pion-nucleon and nucleon-nucleon scattering. However, this study was done only at some specific parameters and the results were not extrapolated towards the continuum and chiral limit. In later works, chiral and continuum limit was studied systematically in the case of pion-pion scattering, especially in the $`I=2`$ channel  It is known that, using the symmetric lattices and Wilson action, large lattices have to be simulated which require substantial amount of computing resources. In Ref. , we have shown that, pion-pion scattering lengths can be obtained using relatively small lattices with the help of tadpole improved Wilson quarks on anisotropic lattices. It is also possible to perform the chiral, infinite volume and continuum limits extrapolation. The final extrapolated result on the pion-pion scattering length in the $`I=2`$ channel is in good agreement with Chiral perturbation result, the experimental result and other lattice results using conventional Wilson fermions . In this paper, similar techniques are applied and the kaon-nucleon scattering lengths in the $`I=1`$ channel is studied.

The fermion action used in this study is the tadpole improved clover Wilson action on anisotropic lattices :

$`_{xy}`$ $`=`$ $`\delta _{xy}\sigma +𝒜_{xy}`$
$`𝒜_{xy}`$ $`=`$ $`\delta _{xy}\left[{\displaystyle \frac{1}{2\kappa _{\mathrm{max}}}}+\rho _t{\displaystyle \underset{i=1}{\overset{3}{}}}\sigma _{0i}_{0i}+\rho _s(\sigma _{12}_{12}+\sigma _{23}_{23}+\sigma _{31}_{31})\right]`$ (2)
$``$ $`{\displaystyle \underset{\mu }{}}\eta _\mu \left[(1\gamma _\mu )U_\mu (x)\delta _{x+\mu ,y}+(1+\gamma _\mu )U_\mu ^{}(x\mu )\delta _{x\mu ,y}\right],`$
where various coefficients in the fermion matrix $``$ are given by:

$`\eta _i`$ $`=`$ $`\nu /(2u_s),\eta _0=\xi /2,\sigma =1/(2\kappa )1/(1\kappa _{max}),`$
$`\rho _t`$ $`=`$ $`c_{SW}(1+\xi )/(4u_s^2),\rho _s=c_{SW}/(2u_s^4).`$ (3)
Among the parameters which appear in the fermion matrix, $`c_{SW}`$ is the coefficient of the clover term and $`\nu `$ is the so-called bare velocity of light, which has to be tuned non-perturbatively using the single pion dispersion relations . Tuning the clover coefficients non-perturbatively is a difficult procedure. In this work, the tadpole improved tree-level value, namely $`c_{SW}=1`$, is used.

In the fermion matrix (1), the bare quark mass dependence is singled out into the parameter $`\sigma `$ and the matrix $`𝒜`$ remains unchanged if the bare quark mass is varied. This shifted structure of the matrix $``$ can be utilized to solve for quark propagators at various values of valance quark mass $`m_0`$ (or equivalently $`\kappa `$) at the cost of solving only the lightest valance quark mass value at $`\kappa =\kappa _{\mathrm{max}}`$, using the so-called Multi-mass Minimal Residual ($`M^3R`$ for short) algorithm . This is particularly advantageous in a quenched calculation since one needs the results at various quark mass values to perform the chiral extrapolation.

This paper is organized as follows. In Section 2, the method to calculate $`KN`$ scattering lengths is briefly reviewed. The $`KN`$ correlation functions used in our calculation are defined. In Section 3, some simulation details are described. In Section 4, our results of the $`KN`$ scattering lengths obtained on lattices of various sizes and lattice spacings are extrapolated towards the chiral limit. Physical strange quark mass is also determined from the $`K`$-input. Finite size effects are then studied and the infinite volume limit are taken. Finally, our lattice results are extrapolated towards the continuum limit. Comparisons with the experimental value and values from chiral perturbation theory at various orders are also discussed. In Section 5, we conclude with some general remarks.

## 2 Formulation to extract the scattering length

In order to calculate the elastic scattering lengths for hadron-hadron scattering on the lattice, or the scattering phase shifts in general, one uses Lüscher’s formula which relates the exact energy level of two hadron states in a finite box to the elastic scattering phase shift in the continuum. In the case of $`KN`$ scattering at zero relative three momentum, this formula amounts to a relation between the exact energy $`E_{KN}^{(I)}`$ of the $`KN`$ system in a finite box of size $`L`$ with isospin $`I`$, and the corresponding scattering length $`a_0^{(I)}`$ in the continuum. This formula reads :

$$E_{KN}^{(I)}(m_K+m_N)=\frac{2\pi a_0^{(I)}}{\mu _{KN}L^3}\left[1+c_1\frac{a_0^{(I)}}{L}+c_2(\frac{a_0^{(I)}}{L})^2\right]+O(L^6),$$
(4)
where $`c_1=2.837297`$, $`c_2=6.375183`$ are numerical constants and $`\mu _{KN}=m_Km_N/(m_K+m_N)`$ is the reduced mass of the $`KN`$ system. In this paper, the $`KN`$ scattering length $`a_0^{(1)}`$ in the $`I=1`$ channel will be studied.

To measure the hadron mass values $`m_K`$, $`m_N`$ and to extract the energy shift $`\delta E_{KN}^{(1)}`$, we construct the correlation functions from the corresponding operators in the $`I=1`$ channel. In this work, the hadron states that need to be measured include pion, kaon, rho, nucleon and $`KN`$ states. The local operators which create a single pion with appropriate isospin values are:

$`\pi ^+(𝐱,t)`$ $`=`$ $`\overline{d}(𝐱,t)\gamma _5u(𝐱,t),\pi ^{}(𝐱,t)=\overline{u}(𝐱,t)\gamma _5d(𝐱,t),`$
$`\pi ^0(𝐱,t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}[\overline{u}(𝐱,t)\gamma _5u(𝐱,t)\overline{d}(𝐱,t)\gamma _5d(𝐱,t)],`$ (5)
where $`u(𝐱,t)`$, $`d(𝐱,t)`$, $`\overline{u}(𝐱,t)`$ and $`\overline{d}(𝐱,t)`$ are basic local quark fields corresponding to $`u`$ and $`d`$ quarks, respectively. The operators which correspond to zero momentum single pions are:

$$\pi _0^a(t)=\frac{1}{L^{3/2}}\underset{𝐱}{}\pi ^a(𝐱,t),$$
(6)
where the flavor of pions $`a=+,,0`$ and $`L^3`$ is the three volume of the lattice. Zero momentum one pion correlation function can then be formed:

$$C_\pi (t)=<\pi _0^{\overline{a}}(t)\pi _0^a(0)>.$$
(7)
Similarly, one can construct the kaon correlation functions from kaon operators. In the $`I=1`$ channel, only $`K^+`$ and proton operators are needed:

$`K^+(𝐱,t)`$ $`=`$ $`\overline{s}(𝐱,t)\gamma _5u(𝐱,t),`$
$`p(𝐱,t)`$ $`=`$ $`ϵ_{abc}\left[u_a^T(𝐱,t)(C^1\gamma _5)d_b(𝐱,t)\right]u_c(𝐱,t).`$ (8)
The proton operator given above is itself a Dirac spinor. The zero momentum proton-proton correlation function is obtained by projecting out the positive parity part of the zero momentum spinor operator and taking the trace in Dirac space:

$$C_N(t)=\underset{𝐱,𝐲}{}p(𝐱,0)\left(\frac{1+\gamma _0}{2}\right)p(𝐲,t).$$
(9)
From the large $`t`$ behavior of the correlation functions, the pion mass $`m_\pi `$, the kaon mass $`m_K`$ and the nucleon mass $`m_N`$ are obtained.

Two hadron operators are used in the construction of two hadron correlation functions. For example, the $`KN`$ operator in the $`I=1`$ channel is given by:

$$O_{KN}^{I=1}(t)=K_0^+(t)p_0(t+1),\overline{O}_{KN}^{I=1}(t)=K_0^{}(t)\overline{p}_0(t+1),$$
(10)
where $`K_0^+(t)`$ and $`p_0(t)`$ are the zero momentum kaon and proton operators respectively. From these operators, two hadron correlation functions can be constructed. The $`KN`$ correlation function in the $`I=1`$ channel reads:

$$C_{KN}^{I=1}(t)=<\overline{O}_{KN}^{I=1}(0)\left(\frac{1+\gamma _0}{2}\right)O_{KN}^{I=1}(t)>.$$
(11)
Numerically, it is more advantageous to construct the ratio of the correlation functions defined above:

$$^{I=1}(t)=C_{KN}^{I=1}(t)/(C_K(t)C_N(t)).$$
(12)
This ratio thus exhibits the following asymptotic behavior for large $`t`$:

$$^{I=1}(t)\stackrel{t>>1}{}e^{\delta E_{KN}^{(1)}t},$$
(13)
with $`\delta E_{KN}^{(1)}=E_{KN}^{(1)}m_Km_N`$ is the energy shift in this channel which directly enters Lüscher’s formula (4). This formula is utilized for small $`\delta E_{KNi}^{(1)}t`$ values (large enough $`L`$) in the calculation where the signal to noise ratio is good. In this case, one could equally use the linear fitting function:

$$^{I=1}(t)\stackrel{T>>t>>1}{}1\delta E_{KN}^{(1)}t,$$
(14)
to determine the energy shift $`\delta E_{KN}^{(1)}`$. Similar procedure has also been used in lattice studies of pion-pion scattering lengths by various groups.

$`KN`$ correlation function, or equivalently, the ratio $`^{I=1}(t)`$ constructed above can be transformed into products of quark propagators using Wick’s theorem.

These contributions are represented by diagrams as shown in Fig. 1, where the solid, the dashed and the dashed-dotted lines represent the $`u`$, $`d`$ and $`s`$ quark propagators, respectively.111In this work, the $`u`$ and $`d`$ quark are assumed to be degenerate in mass.The $`KN`$ correlation function in the $`I=0`$ channel is, however, more complicated which involves vacuum diagrams that require to compute the quark propagators for wall sources placed at every time-slice, a procedure which is more time-consuming than in the $`I=1`$ channel.

## 3 Simulation details

Simulations are performed on a PC based cluster with $`30`$ nodes. Configurations are generated using the pure gauge action (1) for $`6^340`$, $`8^340`$ and $`10^350`$ lattices with the gauge coupling $`\beta =1.7`$, $`1.9`$, $`2.2`$ and $`2.4`$. The spatial lattice spacing $`a_s`$ is roughly between $`0.2`$fm and $`0.4`$fm while the physical size of the lattice ranges from $`1.2`$fm to $`4.0`$fm. For each set of parameters, several hundred decorrelated gauge field configurations are used to measure the fermionic quantities. Statistical errors are all analyzed using the usual jack-knife method. Basic information of these configurations is summarized in Table 1.

Quark propagators are measured using the Multi-mass Minimal Residue algorithm for several different values of bare quark mass. Periodic boundary condition is applied to all three spatial directions while in the temporal direction, Dirichlet boundary condition is utilized. In this calculation, it is advantageous to use the wall sources which enhance the signal . Values of the maximum hopping parameter $`\kappa _{\mathrm{max}}`$, which corresponds to the lowest valance quark mass, are also listed in Table 1. Typically, a few hundred Minimal Residual iterations are needed to obtain the solution vector for a given source. On small lattices, in particular those with a low value of $`\beta `$, the hopping parameter has to be kept relatively far away from its critical value in order to avoid the appearance of exceptional configurations. The parameter $`\nu `$, also known as the bare velocity of light, that enters the fermion matrix (1) is determined non-perturbatively using the single pion dispersion relations as described in Ref.  and Ref. . The optimal value of $`\nu `$ which corresponds to each parameter set of the simulation is also tabulated in Table 1

Single pion, kaon, rho and nucleon mass values are obtained from the plateau of their corresponding effective mass plots. The fitting interval is automatically chosen by minimal $`\chi ^2`$ per degree of freedom. Due to the usage of finer lattice spacing in the temporal direction, good plateau behavior was observed in these effective mass plots. Therefore, contaminations from excited states should be negligible. These mass values will be utilized in the chiral extrapolation.

$`KN`$ correlation functions and the ratio $`^{(1)}(t)`$ are constructed from products of suitable quark propagators according to Wick’s theorem. Both the Direct and the Cross contributions are included. For the ratio $`^{(1)}(t)`$, we obtain clear signal for all our data sets. Since a kaon consists of one up (down) and one strange quark whose mass is different from the up and down quarks, all physical quantities involving the kaon will depend on two quark masses. In terms of the hopping parameters, the $`KN`$ correlation functions depends on a pair of hopping parameters: $`(\kappa _u,\kappa _s)`$. In Fig. 2 we show the ratio $`^{(1)}(t)`$ for different pairs of $`\kappa `$ values as a function of the temporal separation $`t`$, together with the corresponding linear fit (14) for $`\beta =2.4`$ on $`8^340`$ lattices. The starting and ending positions of the straight lines indicate the appropriate fitting range which minimizes the $`\chi ^2`$ per degree of freedom.

After obtaining the energy shifts $`\delta E_{KN}^{(1)}`$, these values are substituted into Lüscher’s formula to solve for the scattering length $`a_0^{(1)}`$ for all possible hopping parameter pairs: $`(\kappa _u,\kappa _s)`$, which corresponds to the up (down) and the strange bare quark mass values respectively. This is done for lattices of all sizes being simulated and for all values of $`\beta `$. From these results, attempts are made to perform an extrapolation towards the chiral, infinite volume and continuum limit.

## 4 Extrapolation towards chiral, infinite volume and continuum limits

The chiral extrapolations of physical quantities involving both the up (down) and the strange quarks consists of two steps. In the first step, the bare strange quark mass, or equivalently the corresponding hopping parameter $`\kappa _s`$, is kept fixed while the hopping parameter of the up/down quark, $`\kappa _u`$, is brought to their physical value $`\kappa _u^{(phy)}`$. The precise value of $`\kappa _u^{(phy)}`$ can be obtained by inspecting the chiral behavior of the pseudo-scalar (pion) and the vector meson (rho) mass values. In the second step, one fixes the up/down quark hopping parameter at its physical value obtained in the first step and extrapolate/interpolate in the strange quark mass. The physical strange quark hopping parameter $`\kappa _s^{(phy)}`$ is determined by demanding the mass of the kaon being exactly at its physical value ($`493`$MeV). In practice, this could be accomplished by, for example, inspecting the kaon to rho mass ratio. This is also known as the $`K`$-input. After these two steps, we would be able to obtain physical quantities at the physical up/down and strange quark mass values. In Table 1, we have listed the physical hopping parameters $`\kappa _u^{(phy)}`$ and $`\kappa _s^{(phy)}`$ for the up/down quarks and the strange quark. Note that one only needs the pion, kaon and rho correlation functions to determine the physical hopping parameters of the quarks.

We now come to the chiral extrapolation of the scattering length. In the chiral limit, the $`KN`$ scattering length in the $`I=1`$ channel is given by the current algebra result:

$$a_0^{(1)}=\frac{1}{4\pi }\frac{\mu _{KN}}{f_K^2},$$
(15)
where $`a_0^{(1)}`$ is the $`KN`$ scattering length in the $`I=1`$ channel and $`f_K113`$MeV is the kaon decay constant. Chiral Perturbation Theory to one-loop order gives an expression which includes the next-to-leading order contributions . Unlike the chiral perturbation theory results for the mesons, the chiral expansion for baryons also receives corrections which are suppressed by the meson to baryon mass ratios. Since the mass of the kaon is not light, one typically obtains large corrections from the next chiral order compared with the previous one. In the case of the $`KN`$ scattering length $`a_0^{(1)}`$, the next order correction is $`114`$% with the opposite sign and the next-to-next order is about $`94`$% with the same sign of the tree-level. This is a signal that one really has no control over the expansion.

To perform the chiral extrapolation of the scattering length, it is more convenient to use the quantity $`F=a_0^{(1)}m_\rho ^2/\mu _{KN}`$, which in the chiral limit reads:

$$F\frac{a_0^{(1)}m_\rho ^2}{\mu _{KN}}=\frac{1}{4\pi }\frac{m_\rho ^2}{f_K^2}3.686,$$
(16)
where the final numerical value is obtained by substituting in the experimental values for $`m_\rho 770`$MeV and $`f_K113`$MeV. The factor $`F`$ can be calculated on the lattice with good precision without the lattice calculation of meson decay constants. The error of the factor $`F`$ obtained on the lattice will mainly come from the error of the scattering length $`a_0^{(1)}`$, or equivalently, the energy shift $`\delta E_{KN}^{(1)}`$ and nowhere else. Since we have calculated the factor $`F`$ for several different values of valance quark mass, we could make a chiral extrapolation and extract the corresponding results in the chiral limit.

First, we extrapolate the factor $`F`$ in the up/down quark mass values. In this step, we have adopted a simple linear extrapolation. the fitting range of the extrapolation is self-adjusted by the program to yield a minimal $`\chi ^2`$ per degree of freedom. Second, the extrapolated values for the factor $`F`$ are then extrapolate/interpolate in the strange quark mass. Our hopping parameters are chosen such that the physical strange quark mass is either within the range of our choice or very close to avoid a long extrapolation in the strange quark mass. In Fig. 3, we have shown the two step chiral extrapolation for one of our simulation points.

In the lower panel of Fig. 3, we have plotted the results for the factor $`F`$ for all possible pairs of quark mass values. They are arranged according to different up/down quark hopping parameter. The straight lines show the extrapolation in $`m_\pi ^2`$ (the first step extrapolation). The physical up/down quark hopping parameter are read from the pion and rho correlation function fits which are tabulated in Tab. 1. The resulting factor $`F`$ after the first extrapolation are plotted in the upper panel of Fig. 3 versus the strange quark mass parameter $`1/(2\kappa _s)`$. The straight line is the linear extrapolation/interpolation (second extrapolation) and the final result is also depicted as a solid square at the physical strange quark mass. Current algebra result and the result from chiral perturbation theory are also shown as a filled triangle (tree level), a filled square (next order) and a filled hexagon (next-to-next), respectively. The result from Ref.  is shown as a filled circle. To avoid crowdedness, these symbols are shifted to the left. The fitting quality for the factor $`F`$ are reasonable. The quality for the chiral extrapolation of other simulation points are similar. As is seen, the linear fit gives a reasonable modelling of the data. The divergent contributions and other higher order corrections from quenched chiral perturbation theory seem to be numerically small for the lattices being simulated.

After the chiral extrapolation, we now turn to study the finite volume effects of the simulation. According to formula (4), the quantity $`F`$ obtained from finite lattices differs from its infinite volume value by corrections of the form $`1/L^3`$. However, the situation might be different in a quenched calculation, which happened in the case of pion-pion scattering . In the case of $`KN`$ scattering length, a similar calculation using quenched chiral perturbation theory is still not available. Therefore, we will be using a linear extrapolation in $`1/L^3`$, as suggested by Lüscher’s formula. In Table. 1, we have listed the infinite volume extrapolation for our simulation points at $`\beta =2.4`$, $`2.2`$, $`1.9`$ and $`1.7`$.

Finally, we can make an extrapolation towards the continuum limit by eliminating the finite lattice spacing errors. Since we have used the tadpole improved clover Wilson action, all physical quantities differ from their continuum counterparts by terms that are proportional to $`a_s`$. The physical value of $`a_s`$ for each value of $`\beta `$ can be found from Ref. , which is also included in Table 1. This extrapolation is shown in Fig. 4 where the results from the chiral and infinite volume extrapolation discussed above are indicated as data points in the plot for all $`4`$ values of $`\beta `$ that have been simulated. The straight line shows the extrapolation towards the $`a_s=0`$ limit and the extrapolated results are also shown as solid squares together with the experimental result from Ref. which is shown as the filled circle. For comparison, the results from chiral perturbation theory ar various orders are also shown as the filled triangle, the filled square and the filled hexagon at $`a_s=0`$, respectively. It is seen from the figure that our extrapolated result is negative, with a large error bar. Therefore, we can only say that, from our preliminary study we obtain a result which is compatible with the experimental result and results from chiral perturbation theory. The error of our calculation can be reduced by further studies with more statistics and on lattices with smaller lattice spacing in the future. To summarize, we obtain from the linear extrapolation the following result for the quantity $`F=a_0^{(1)}m_\rho ^2/\mu _{KN}=(1.6\pm 5.4)`$. If we substitute in the physical values, we obtain the $`KN`$ scattering length in the $`I=1`$ channel: $`a_0^{(1)}=(0.17\pm 0.59)\text{fm}`$, which is to be compared with the experimental result of $`(0.33)\text{fm}`$ and the chiral result of $`(0.398+0.4570.389)\text{fm}`$ at the first, second and third chiral order, respectively.

## 5 Conclusions

In this paper, we have calculated kaon nucleon scattering lengths in isospin $`I=1`$ channel using quenched lattice QCD. It is shown that such a calculation is feasible using relatively small, coarse and anisotropic lattices. The calculation is done using the tadpole improved clover Wilson action on anisotropic lattices. Simulations are performed on lattices with various sizes, ranging from $`1.2`$fm to about $`4`$fm and with four different values of lattice spacing. Quark propagators are measured with different valence quark mass values. These enable us to explore the finite volume and the finite lattice spacing errors in a systematic fashion. The infinite volume extrapolation is made. The lattice result for the scattering length is extrapolated towards the chiral and continuum limit where a result consistent with the experiment and the Chiral Perturbation Theory is found. We believe that, using the method described in this exploratory study, more reliable results on $`KN`$ scattering lengths could be obtained with more statistics and with simulations on finer lattices.

## Acknowledgments

This work is supported by the National Natural Science Foundation (NFS) of China under grant No. 90103006, No. 10235040 and supported by the Trans-century fund from Chinese Ministry of Education. C. Liu would like to thank Prof. H. Q. Zheng and Prof. S. L. Zhu for helpful discussions.


format: text

text:
# High frequency features in the 1998 outburst of 4U 1630–47

## 1 Introduction

The soft X-ray transient 4U 1630–47 was discovered with *Uhuru*  but already detected with *Vela* 5B in 1969 . There is no optical counterpart known and the best estimate we have for the distance is about 10 kpc . The detection of X-ray absorption dips suggest the system has an inclination angle of $`60`$–75 . 4U 1630–47 is known for its regular outburst cycle: it was shown to have outbursts every $`600`$ days ( ) however, recently this period changed to about 690 days . Based on the X-ray spectral and timing characteristics 4U 1630–47 has been classified as a black hole candidate , and it was pointed out that 4U 1630–47 shows similarities to the galactic sources GRO J1655–40 and GRS 1915+105 (similar variability behaviour; ).

Currently, 8 black hole sources are known to show high frequency QPOs, including 4U 1630–47 . The other sources are: Cyg X-1 , GRO J1655-40 , XTE J1650-500 , XTE J1550-564 , GRS 1915+105 , XTE J1859+227 , and most recently H 1743–322 (IGR/XTE J1746–3213, ) was added to this list. In some of these cases (GRO J1655-40, XTE J1550-564 and GRS 1915+105) a doublet of high frequency features has been detected which seem to have 2:3 or 3:5 frequency ratios.

Here we report for the first time on the complete high frequency behaviour of 4U 1630–47 during its 1998 outburst. We show that 4U 1630–47 is capable of producing strong (2–6% rms), both broad and narrow (Q$``$0.3–5.0) high frequency features, that are comparable in frequency to those found in other black holes. However, we find that 4U 1630–47 shows about a factor $`3`$ more power at high frequencies compared to other sources. Moreover we find, for the first time in a black hole, simultaneous high frequency features that have a frequency ratio of 1:4. The complete timing analysis of this outburst covering the general outburst behaviour, the low- and high frequency behaviour will be presented by Klein-Wolt, Homan & van der Klis (2003, in prep.).

## 2 The 1998 outburst

### 2.1 Observations and data analysis

We analyse Proportional Counter Array (PCA) Rossi X-ray Timing Explorer (*RXTE*) observations of 4U 1630–47 made between February 11 1998 and June 8 1998 and for a total of 101 observations we have about $`300`$ ksec of data. For the timing study we use high time resolution data modes, and power spectra are created of data segments 128s using the standard Fast Fourier Transform techniques ( and references therein). Detector drop-outs are removed but no background or deadtime corrections are performed prior to the production of the power spectra. For each observation the power spectra are averaged.

We pay special attention to the presence of possible power spectral features at high frequencies. Note that in the following we refer to the interval above $`20`$ Hz as the high frequency, while the low frequency regime is the interval below $`20`$ Hz. As explained in Klein-Wolt, Homan & van der Klis (2003, in prep.), we estimate the Poisson level using the representation of Zhang  with an additional shift to fit the data in a frequency range where no source power is expected (above $`3000`$ Hz). After considerable experimentation we decided to produce the power spectra in the 7–60 keV range (channels 18–157, gain epoch 3111All energy ranges correspond to gain epoch 3) for observations belonging to proposals 30178 and 30188, while for observations belonging to proposal 30172 the 5.5–104 keV range (channels 14–254) was used. These ranges are chosen in such a way to eliminate an instrumental effect that appeared as a broad feature at high frequencies ($`1000`$ Hz) in our power spectra.

The power spectra are fitted with a multi-Lorentzian model (, ) using the approach of  to assign, to each component, a characteristic frequency which corresponds to the maximum in power for a Lorentzian feature in the $`\nu `$-$`P_\nu `$ representation. We only keep those Lorentzian components in the fits whose significance based on the error in the integrated power (from 0–$`\mathrm{}`$) is more than $`3\sigma `$, or whose inclusion gives a significant ($`>3\sigma `$) improvement of the $`\chi ^2`$ of the fit.

### 2.2 General outburst behaviour

In Fig. 1 we show the Hardness-Intensity Diagram (HID) of the outburst. The source starts at the top-middle of the diagram, moves to the upper-left as count rate increases and hardness decreases, then moves down as the outburst declines and finally ends up at the lower-right where the count rate is low and the spectrum is hard. During the movement in the HID the source occupies different locations in the diagram that are related to specific types of timing behaviour. We find that there are six different groups (A–F) in the HID (and also in the colour-colour diagram (CD)). The average power spectra of these groups are shown in Fig. 2; they are characteristic for the individual power spectra in each group. Note, that averaging power spectra together can result in extra features that are either due to small shifts in frequency of one component between different observations, or to the increase in significance of features that are too weak to be detected in individual observations. In group A for instance, we detect two or three features in the individual observations between 2–20 Hz, but the average spectrum in Fig. 2a clearly shows four features. These are the result of the aforementioned frequency shifts.

The power spectra are sometimes very complex, even in the individual observations, and show features from $`0.1`$$`300`$ Hz. Notice, for instance the so called ”dipping QPOs” around 0.1 Hz, already reported by , and the strong (up to $`7`$%) QPOs around 10 Hz (Fig. 2). The most remarkable features, however, are the high frequency features which are detected typically above $`20`$ Hz and are only present in groups A, B, C and D. We will discuss these in more detail in the next section.

During the outburst the source shows a number of flares, indicated with the roman numerals (I–IX) in the HID (Fig. 1). During these flares the spectrum hardens, the rms amplitude increases and features appear or disappear in the power spectrum. In cases when the source during a flare moves in the direction of another group in the HID (for instance during flare II and III from group D to B, Fig. 1), we find that the power spectral behaviour is very similar to that group. This, together with the fact that different groups have distinct different power spectral characteristics, suggest that the behaviour of 4U 1630–47 in the HID (and CD) is two-dimensional: the power spectral behaviour is determined by the position in the HID and not by the chronology by which the source moves through the various regions in the HID. Therefore, the power spectral behaviour depends on the two parameters that determine the position in the HID, the count rate and the hard colour (HC). This behaviour is very similar to that of XTE J1550–564  and it suggests that both sources share a common outburst mechanism.

## 3 The high frequency features

In Fig. 2 we show the high frequency features detected in groups A–D. In Table 1 we give the parameters of these components. For groups A and C we find one high frequency feature at 164 Hz and 262 Hz, respectively, that is relatively broad. In groups B and D we find two high frequency features each time, one below 100 Hz and one above. For both the twin high frequency features, we find that their ratio of frequencies is about 1:4.

In 22 of the 101 individual observations during the 1998 outburst of 4U 1630–47 we detect significant high frequency features. These features cover a range of widths with Q values of $`0`$$`6`$ and frequencies ranging from $`23`$$`300`$ Hz. For these features we find an increase in frequency during the rise of the outburst (first $`4`$ days) which is followed by a more or less stable region during which the frequencies hardly change (remain constant at $`180`$ Hz), and finally a decrease in frequencies during the decay of the outburst. This behaviour is reflected in a correlation with count rate and an anti-correlation with HC for all the high frequency features detected above 100 Hz. We find that the features detected below 100 Hz do not seem to follow that same relation. Note, that the same correlations between count rate and HC are also found for the features detected in the average power spectra. Finally, also in five of the individual observations in groups A, B, C and D we find a feature below 100 Hz. However, only on one of these occasions another feature above 100 Hz is detected: again we find two features (at 42 and 170 Hz) with a ratio in frequencies of 1:4.

From Table 1 and Fig. 2 it is clear that the high frequency features detected in 4U 1630–47 are usually very broad and strong: we find rms amplitudes from $`2`$$`7`$ % in both the average and individual observations. We compare the average power in the high frequency range (100–3000 Hz) with that of two black hole sources, Cyg X-1 and XTE J1550–564 and two neutron star sources, 4U 1608–522 and 4U 1728–34. We find that on average 4U 1630–47 shows $`3`$ times more power compared to the black hole sources and about 3 times less power compared to the neutron star sources. Also, when comparing the shape of the power spectra the difference between 4U 1630–47 and other black hole sources is evident. From Fig. 3 it is clear that for the two black hole sources (Cyg X-1 and XTE J1550–564) the power spectrum *declines* above $`100`$ Hz, while for both 4U 1630–47 and the neutron star 4U 1728–34 the power *increases*. This clearly shows that, like neutron stars this black hole candidate, too, can show relatively high power levels at high frequencies.

## 4 Discussion and conclusions

In Fig. 4 we show the Psaltis-Belloni-van der Klis relation (PBK relation) for a number of black hole sources and neutron star sources (taken from ). For the lower frequency QPO of 4U 1630–47 we always take the first Lorentzian of the complex of features around $`10`$ Hz, and we plot both the points for the individual and average power spectra. First of all, it is clear from this figure that 4U 1630–47 shows a correlation between the low frequency features and the high frequency features: 4U 1630–47 shows the same trend as the original PBK relation and falls on the lower PBK track. Notice that this is exactly at the same position as XTE J1550–564, indicating that both sources show features at similar frequencies. Secondly, we find that the high frequency features below 100 Hz (the big stars and diamonds in Fig. 4) do not fall on the PBK relation, nor do they show any consistent relation between the low and high frequency features. This, together with the fact that they also do not follow the same relations with HC and count rate as the features detected above 100 Hz and the fact that they are about a factor 2 lower in frequency, strongly suggest that the features below 100 Hz are different from the other high frequency features and might even have a different origin.

The agreement of 4U 1630–47 with the PBK relation raises the question if the high frequency features in this source are related to the kHz QPOs found in the neutron star sources? Although the relation between NS and BH features over a large range in frequencies shown by the PBK relation in itself suggests a similarity between these two types of sources, we would like to stress that the most important results that can be deduced from this relation are: *1) the correlation between the low and the high frequency features (PBK relation) for 4U 1630–47*; and *2) the overlap in frequencies between 4U 1630–47 and the black hole source XTE J1150-564*. Both these results clearly suggest that 4U 1630–47 can behave very similar to other black hole sources, showing similar frequencies and relations between them.

It is for the first time for a black hole X-ray binary system that twin high frequency features are detected in the same energy band, with a frequency ratio of 1:4, clearly different from the more common 2:3 or 3:5 ratios (see for instance ). The fact that the ratio is different in the case of 4U 1630–47 might be related to the fact that the lower of the two high frequency features seems to be different, as explained above. However, during the outburst of other black hole sources similar features are found. For instance during the 2001 outburst of XTE J1650-500 a feature around 50 Hz is detected . This feature is simultaneously present with features at higher frequencies (at twice or three times the frequency) and is broad but much stronger compared to the feature at higher frequencies . Also, similar features are reported for Cyg X-1, GX 339-4 and XTE J1118-480. It is not clear whether these features are all the same, however, they might be related to the Hecto Hz QPOs found in both neutron star and black hole sources (see  and references therein); the anomalous high frequency features found in 4U 1630–47 have similar frequencies, however, they are not as stable in frequency as the normal Hecto Hz QPOs.

Note that, although the 1:4 ratio is unconventional, it follows for instance from the discoseismic model ( and references therein) as a combination between a g- and c-mode oscillation or as a non-axisymmetric oscillation between two g-modes . Also, in the model by Abromowicz and Kluźniak it is perhaps possible to get a 1:4 ratio by putting the ratio between the relativistic epicyclic frequency and the relativistic Keplerian frequency equal to 1:4. These models do not exclusively predict the ratio of the frequencies, hence the ratio (2:3, 3:5 or 1:4) might be determined by the characteristics of the individual sources.

In conclusion, we find that during the 1998 outburst of 4U 1630–47 broad, strong high frequency features are present, which are comparable in frequency to the high frequency features detected in XTE J1550–564. We find that on average 4U 1630–47 shows about 3 times more power in the high frequency regime compared to other black hole sources, however, also about 3 times less power compared to neutron star sources. Also, the overall behaviour of 4U 1630–47 is similar to that of XTE J1550–564 and suggests a two-dimensional behaviour in the HID, which depends on both the count rate and the HC. Based on the behaviour described here we suggest that identifying black hole and neutron star sources based on the amount of power at high frequencies, as suggested by , is not straightforward: 4U 1630–47 is capable of showing an increase in power above $`100`$ Hz, like the neutron star sources. In fact, our findings seem to suggest that black holes and neutron stars can show very similar behaviour, indicating that the timing properties are mainly determined by processes that take place in the accretion disc and are less dependent on the presence of a surface, boundary layer, or strong magnetic field.


text:
# Inflating magnetically charged braneworlds

## I Introduction

In “braneworld” models, our universe is represented by a $`(3+1)`$-dimensional brane floating in a higher-dimensional bulk spacetime . The brane can be thought of as fundamental $`D`$-brane, or it can arise as a topological defect in a higher-dimensional field theory. It may well be that $`D`$-branes could also eventually find some dual field-theory description. In the simplest, codimension-one models, the brane can be pictured as a domain wall propagating in a 5D spacetime . Higher codimensions with both gauge and global defects have also been considered. In particular, models have been discussed where the field configuration in the directions orthogonal to the brane is that of a cosmic string , a monopole , or a texture . In these models, the bulk curvature produced by the defects plays a major role in localizing gravity on the brane and in solving the mass-hierarchy problem.

The physics of braneworld models crucially depends on the nature of the gravitational fields produced by defects in higher dimensions. These fields, however, are not yet fully understood. Exact analytic solutions have been obtained to the higher-dimensional Einstein equations, and to the combined system of Einstein, gauge and/or Goldstone field equations for branes carrying gauge or global charges . These solutions correspond to the limit of a vanishingly thin core and contain curvature singularities. One can hope that the singularities will be smoothed out in more realistic models including the defect core, and it has been verified in  that this is indeed the case for higher-dimensional global monopoles.

Apart from the spurious singularities resulting from the thin-core approximation, the defect spacetimes can also develop singularities in the exterior region. Well-known $`3`$-dimensional examples are static domain wall  and global string  solutions. In these cases, the singularities are removed by relaxing the requirement of staticity and allowing the defect worldsheets to inflate . Static 3-dimensional gauge strings are also known to be singular when their symmetry-breaking scale becomes greater than the Planck scale (see, e.g., Ref.  for discussion and references).

In Ref.  we showed that similar singularities arise in higher-dimensional monopole spacetimes, if the symmetry-breaking scale of the monopole exceeds certain critical value comparable to the Planck scale, and that the singularities can be removed if the monopole worldsheets are allowed to inflate. We argued that the inflation rate $`H`$ is uniquely fixed by the regularity condition.

Another motivation for looking into the super-Planckian regime is the interesting work by Dvali, Gabadadze and Shifman , suggesting that braneworld models of codimension $`n>2`$, with a very low bulk Planck scale, can provide a solution to the cosmological constant problem. It was conjectured in Ref.  that defect solutions in such models may exhibit de Sitter inflation of the worldsheet, with the expansion rate inversely proportional to the brane tension. The very low expansion rate observed today could then be obtained with a very large brane tension. Our analysis of global defects in  has shown that global higher-dimensional defects do not have this property.

The purpose of the present paper is to extend the analysis of Ref.  to the case of gauge defects. Specifically, we consider a codimension-3 gauge monopole. In the following Section, we introduce the model and review some asymptotic solutions found earlier by Gregory . The numerical solutions assuming a static (noninflating) worldsheet are presented in Section III, for both sub- and super-critical regimes. We find that sub-critical solutions are regular and asymptotically flat, while the super-critical ones have a curvature singularity at a finite distance from the core (as they do in the global monopole case). Inflating defect spacetimes are discussed in Sec. IV. We find analytic solutions describing the asymptotic behavior of these spacetimes and chart the rather complicated parameter space of these solutions. Our numerical results indicate that inflation of the worldsheet does remove the singularity and that the resulting spacetimes have a ‘cigar’ geometry. However, the numerical accuracy limitations do not allow us to conclude that the regularity requirement selects a unique value of the expansion rate $`H`$. Our conclusions are briefly summarized in Section V.

## II Field equations and some asymptotic solutions

The action for our model is

$$𝒮=d^7x\sqrt{g}\left(\frac{}{2\kappa ^2}+_m\right),$$
(1)
where $`\kappa ^2=1/M^{2+n}`$ and $`M`$ is the 7D Planck mass. We assume that the metric has de Sitter or Poincare symmetry along the brane and spherical symmetry in the extra three dimensions. This corresponds to the following ansatz:

$$ds^2=B^2(r)d\overline{s}_4^2+dr^2+C^2(r)r^2(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2),$$
(2)
where $`d\overline{s}_4^2`$ is the 4D world-volume metric.

The matter-field Lagrangian is given by

$$_m=\frac{1}{4}G_{MN}^aG^{aMN}\frac{1}{2}D_M\mathrm{\Phi }^aD^M\mathrm{\Phi }^a\frac{\lambda }{4}(\mathrm{\Phi }^a\mathrm{\Phi }^a\eta ^2)^2,$$
(3)
where

$`D_M\mathrm{\Phi }^a`$ $`=`$ $`_M\mathrm{\Phi }^a+eϵ^{abc}W_M^b\mathrm{\Phi }^c,`$ (4)
$`G_{MN}^a`$ $`=`$ $`_MW_N^a_NW_M^a+eϵ^{abc}W_M^bW_N^c,`$ (5)
$`D_M`$ is the covariant derivative, and $`e`$ is the gauge-coupling constant.Here, small Roman indices from the beginning of the alphabet run through the values $`a,b,c=1,2,3`$, and from the middle of the alphabet $`i,j=4,5,6`$. Capital Roman indices take values from 0 to 6, and Greek indices from 0 to 3. In the following sections, it will be convenient to use the notation

$$\beta \lambda /e^2$$
(6)
for the ratio of the scalar and gauge couplings.

We assume that the scalar triplet $`\mathrm{\Phi }^a`$ has a “hedgehog” configuration in extra dimensions,

$$\mathrm{\Phi }^a=\varphi (r)\frac{x^{3+a}}{r}.$$
(7)
The gauge field is given by

$`W_\mu ^a`$ $`=`$ $`0,`$ (8)
$`W_i^a`$ $`=`$ $`[1w(r)]ϵ^{a(i3)(j3)}{\displaystyle \frac{x^j}{er^2}},`$ (9)
or, in spherical coordinates (2),

$`W_r^a`$ $`=`$ $`0,`$ (10)
$`W_\theta ^a`$ $`=`$ $`{\displaystyle \frac{1w(r)}{er}}(\mathrm{sin}\phi ,\mathrm{cos}\phi ,0),`$ (11)
$`W_\phi ^a`$ $`=`$ $`{\displaystyle \frac{1w(r)}{er}}\mathrm{sin}\theta (\mathrm{cos}\theta \mathrm{cos}\phi ,\mathrm{cos}\theta \mathrm{sin}\phi ,\mathrm{sin}\theta ).`$ (12)
Einstein equations corresponding to the above ansätze are

$`G_\mu ^\mu `$ $`=`$ $`3{\displaystyle \frac{B^{\prime \prime }}{B}}3\left({\displaystyle \frac{B^{}}{B}}\right)^26{\displaystyle \frac{B^{}}{Br}}6{\displaystyle \frac{B^{}C^{}}{BC}}2{\displaystyle \frac{C^{\prime \prime }}{C}}\left({\displaystyle \frac{C^{}}{C}}\right)^26{\displaystyle \frac{C^{}}{Cr}}+{\displaystyle \frac{1}{C^2r^2}}{\displaystyle \frac{1}{r^2}}+{\displaystyle \frac{1}{4}}{\displaystyle \frac{\overline{R}^{(4)}}{B^2}}`$ (13)
$`=`$ $`\kappa ^2\left[{\displaystyle \frac{\varphi ^2}{2}}+{\displaystyle \frac{\varphi ^2w^2}{C^2r^2}}+{\displaystyle \frac{1}{e^2C^2r^2}}\left[w^2+{\displaystyle \frac{(1w^2)^2}{2C^2r^2}}\right]+{\displaystyle \frac{\lambda }{4}}(\varphi ^2\eta ^2)^2\right],`$ (14)
$`G_r^r`$ $`=`$ $`6\left({\displaystyle \frac{B^{}}{B}}\right)^28{\displaystyle \frac{B^{}}{Br}}8{\displaystyle \frac{B^{}C^{}}{BC}}\left({\displaystyle \frac{C^{}}{C}}\right)^22{\displaystyle \frac{C^{}}{Cr}}+{\displaystyle \frac{1}{C^2r^2}}{\displaystyle \frac{1}{r^2}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\overline{R}^{(4)}}{B^2}}`$ (15)
$`=`$ $`\kappa ^2\left[{\displaystyle \frac{\varphi ^2}{2}}+{\displaystyle \frac{\varphi ^2w^2}{C^2r^2}}+{\displaystyle \frac{1}{e^2C^2r^2}}\left[w^2+{\displaystyle \frac{(1w^2)^2}{2C^2r^2}}\right]+{\displaystyle \frac{\lambda }{4}}(\varphi ^2\eta ^2)^2\right],`$ (16)
$`G_{\theta _i}^{\theta _i}`$ $`=`$ $`4{\displaystyle \frac{B^{\prime \prime }}{B}}6\left({\displaystyle \frac{B^{}}{B}}\right)^24{\displaystyle \frac{B^{}}{Br}}4{\displaystyle \frac{B^{}C^{}}{BC}}{\displaystyle \frac{C^{\prime \prime }}{C}}2{\displaystyle \frac{C^{}}{Cr}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\overline{R}^{(4)}}{B^2}}`$ (17)
$`=`$ $`\kappa ^2\left[{\displaystyle \frac{\varphi ^2}{2}}{\displaystyle \frac{(1w^2)^2}{2e^2C^4r^4}}+{\displaystyle \frac{\lambda }{4}}(\varphi ^2\eta ^2)^2\right],`$ (18)
where $`\overline{R}^{(4)}`$ is the Ricci scalar of the 4D world volume.

The field equation for the scalar field is

$$\varphi ^{\prime \prime }+2\left(2\frac{B^{}}{B}+\frac{C^{}}{C}+\frac{1}{r}\right)\varphi ^{}2\frac{w^2\varphi }{C^2r^2}\lambda \varphi (\varphi ^2\eta ^2)=0.$$
(19)
The field equation for the gauge field is

$$w^{\prime \prime }+4\frac{B^{}}{B}w^{}+\frac{w(1w^2)}{C^2r^2}e^2\varphi ^2w=0.$$
(20)
At large distances from the monopole core, $`r(\sqrt{\lambda }\eta )^1`$, one can expect the scalar and gauge field to approach the asymptotic values

$$\varphi (r\mathrm{})=\eta ,w(r\mathrm{})=0.$$
(21)
In this regime, the monopole can be approximated as an infinitely thin charged 3-brane. The corresponding solution of Einstein equations has been found by Gregory  in the case of a flat 4D world-volume metric. We shall not reproduce this solution here and note only that it is asymptotically flat, $`B(r\mathrm{})=C(r\mathrm{})=1`$.

Gregory has also considered inflating branes with a de Sitter world-volume metric ,

$$d\overline{s}_4^2=dt^2+H^2\mathrm{cosh}^2(Ht)d\mathrm{\Omega }_3^2.$$
(22)
Using higher-dimensional vacuum Einstein’s equations, she found the asymptotic behavior of the metric at large $`r`$. In the case of a 3-brane in a 7D spacetime, her asymptotic solution has the form,

$$ds^2=\frac{3}{5}H^2r^2(dt^2+e^{2Ht}d\text{x}^2)+dr^2+\frac{1}{5}r^2d\mathrm{\Omega }_2^2.$$
(23)
Note that this is not globally flat at $`r\mathrm{}`$, but has a solid deficit angle in the extra dimensions. Although this asymptotic solution was obtained for an uncharged brane, one can expect it to give a good approximation, since the gauge field energy-momentum tensor vanishes at $`r\mathrm{}`$.

In the following sections, we shall see that Gregory’s solutions indeed describe the asymptotic behavior of the metric, but only if the symmetry-breaking scale $`\eta `$ is sufficiently small. For large values of $`\eta `$, the nature of the solutions is drastically different.

## III Static solutions

We numerically solved Einstein, scalar-field, and gauge-field equations (14)-(20) with boundary conditions, $`\varphi (0)=0`$, $`\varphi (\mathrm{})=\eta `$, $`w(0)=1`$, $`w(\mathrm{})=0`$, $`B(0)=C(0)=1`$, and $`B^{}(0)=C^{}(0)=0`$ . For static solutions, we assumed a flat worldsheet, $`d\overline{s}_4^2=\eta _{\mu \nu }dx^\mu dx^\nu `$, and $`\overline{R}^{(4)}=0`$.

For a sufficiently small symmetry-breaking scale $`\eta `$, we find regular, asymptotically flat solutions with $`B(\mathrm{})=\mathrm{const}`$ and $`C(\mathrm{})=1`$. The fields $`\varphi `$ and $`w`$ rapidly approach their asymptotic values outside the monopole core; see Fig. 1. As $`\eta `$ is increased, the character of the solution gradually changes: the asymptotic value of $`B(r)`$ grows and a ‘cylindrical’ region develops near the core, where $`C(r)r\mathrm{const}`$. This is illustrated in Figs. 2 and 3.

At some critical value $`\eta =\eta _c`$, the solution degenerates into a ‘cigar’, with $`Cr`$ approaching a constant at large $`r`$, and for $`\eta >\eta _c`$ all static solutions are singular. The singularity is at a finite value of $`r`$, where $`C(r)`$ vanishes and $`B(r)`$ diverges (see Figs. 4 and 5). This singular behavior is opposite to that of a global monopole, for which $`C`$ diverges and $`B`$ vanishes at the singular point .

Figure 6 shows the location of the singularity, $`r=r_s`$, as a function of $`\eta `$. As $`\eta `$ is decreased, the singularity moves away from the core, and $`r_s\mathrm{}`$ in the limit $`\eta \eta _c`$.

The asymptotic form of the critical ‘cigar’ solution can be found analytically. With the ansätze $`Cr=\mathrm{const}c_0`$, $`\varphi =\eta `$, $`w=0`$, Einstein equations take the form

$`G_\mu ^\mu `$ $`=`$ $`3{\displaystyle \frac{B^{\prime \prime }}{B}}3\left({\displaystyle \frac{B^{}}{B}}\right)^2+{\displaystyle \frac{1}{c_0^2}}={\displaystyle \frac{\kappa ^2}{2e^2c_0^4}},`$ (24)
$`G_r^r`$ $`=`$ $`6\left({\displaystyle \frac{B^{}}{B}}\right)^2+{\displaystyle \frac{1}{c_0^2}}={\displaystyle \frac{\kappa ^2}{2e^2c_0^4}},`$ (25)
$`G_{\theta _i}^{\theta _i}`$ $`=`$ $`4{\displaystyle \frac{B^{\prime \prime }}{B}}6\left({\displaystyle \frac{B^{}}{B}}\right)^2={\displaystyle \frac{\kappa ^2}{2e^2c_0^4}},`$ (26)
and the solution is easily found,

$$B=b_0e^{\frac{\sqrt{5}e}{8\kappa }r},c_0=\frac{2}{\sqrt{5}}\frac{\kappa }{e},$$
(27)
where $`b_0`$ is an integration constant. The value of $`\eta _c`$ depends on the parameter $`\beta =\lambda /e^2`$ and has to be determined numerically. The cigar solution in Figs. 4 and 5 is obtained for $`\beta =1`$, in which case we find $`\kappa \eta _c=1.28630998`$.

In Fig. 7, we show the critical value $`\kappa \eta _c`$ for several values of $`\beta `$. As $`\beta `$ increases, $`\kappa \eta _c`$ decreases and approaches $`\kappa \eta _c=1`$ in the global monopole limit $`\beta \mathrm{}`$.

## IV Inflating solutions

### A Sub-critical case: $`\eta <\eta _c`$

We shall now discuss monopoles with inflating 4D worldsheets. In this case, the worldsheet metric is a 4D de Sitter space (22), and the 4D curvature scalar is $`\overline{R}^{(4)}=12H^2`$.

Let us first consider sub-critical inflating monopoles with $`\eta `$ significantly below $`\eta _c`$. In this case, we found regular solutions for a range of expansion rates $`H`$ spanning almost two orders of magnitude. All these solutions approach Gregory’s asymptotic form (23) at large $`r`$. We note also that, as $`H`$ is increased, our solutions appear to approach Gregory’s metric (23) at all $`r`$ in the limit $`H\mathrm{}`$. This is illustrated in Figs. 8 and 9 for $`\kappa \eta =0.2`$.

A peculiar feature of these solutions is that the expansion rate $`H`$ is not determined by the energy density of the brane, as in the case of the “usual” 4D inflation, but is rather a free parameter. Also peculiar is the asymptotic geometry at large $`r`$, which exhibits a deficit solid angle $`\mathrm{\Delta }\mathrm{\Omega }=16\pi /5`$, also independent of the stress-energy of the brane. It thus appears that inflation is imposed on the brane by the choice of boundary conditions.

For larger values of $`\eta `$, $`1/\kappa \eta <\eta _c`$, we find a different pattern. As $`H`$ is increased, a cigar solution with $`Cr=\mathrm{const}`$ is obtained at some critical $`H_c(\eta )`$ . At still larger values of $`H`$, the solutions have a curvature singularity of the type similar to the $`H=0`$ case (see Figs. 10 and 11).

Based on this analysis, it appears reasonable to conjecture that the physical solutions for sub-critical branes ($`\eta <\eta _c`$) are the static regular solutions with $`H=0`$. The physical interpretation of inflating sub-critical solutions is not clear to us, and we are inclined to dismiss them as unphysical. We have not, therefore, attempted to determine the dependence $`H_c(\eta )`$ in any detail.

### B Super-critical case: $`\eta >\eta _c`$

In the super-critical case, $`\eta >\eta _c`$, regular static solutions no longer exist. It is well known that the singularity of some static defect solutions can be removed by relaxing the requirement of staticity and allowing the defect worldsheet to inflate. Examples are domain walls  and global strings  in $`(3+1)`$ dimensions. In both of these cases, the inflation rate $`H`$ is uniquely determined by the requirement of regularity, and the inflating solutions can be interpreted as true physical spacetimes for these sources. The same idea was applied to a 7D global monopole with a 3D core  (this corresponds to $`e=0`$, or $`\beta =\mathrm{}`$, in our model). In that paper, we found (i) that regular static solutions exist only for $`\eta <\eta _c=\kappa ^1`$, and (ii) that nonsingular super-critical solutions with $`\eta >\eta _c`$ exist only for inflating branes. In these solutions, the scalar field acquires an asymptotic value which is somewhat lower than the symmetry-breaking scale, $`\varphi =\varphi _0<\eta `$, the extra dimensions have a cigar geometry, and we have argued that the value of $`H`$ is uniquely fixed by the regularity requirement.

We have verified that regular solutions of the same type also exist for super-critical gauge monopoles. The asymptotic form of these solutions at large $`r`$ can be found analytically. With the ansätze $`\varphi =\varphi _0`$, $`w=w_0`$, $`\sqrt{\lambda }\eta Cr=C_0`$, where $`\varphi _0`$, $`w_0`$ and $`C_0`$ are constants, the combined field equations yield the following equation for the scalar field $`\varphi _0`$:

$$\frac{4\kappa ^2}{5\eta ^2}\left(\frac{1}{\beta }\frac{1}{2}\right)\varphi _0^4+\frac{1}{\eta ^2}\left(\frac{3\kappa ^2\eta ^2}{10}\frac{1}{\beta }+\frac{1}{2}\right)\varphi _0^2+\frac{\kappa ^2\eta ^2}{10}\frac{1}{2}=0.$$
(28)
The solution to this equation, as a function of $`\eta `$ and $`\beta `$, is given by

$$\left(\varphi _0^\pm \right)^2=\left[\left(\frac{3\kappa ^2\eta ^2}{10}\frac{1}{\beta }+\frac{1}{2}\right)\pm \sqrt{D}\right]\left[\frac{8\kappa ^2}{5}\left(\frac{1}{\beta }\frac{1}{2}\right)\right]^1,$$
(29)
where

$$D=\left(\frac{3\kappa ^2\eta ^2}{10}\frac{1}{\beta }+\frac{1}{2}\right)^2\frac{16\kappa ^2\eta ^2}{5}\left(\frac{1}{\beta }\frac{1}{2}\right)\left(\frac{\kappa ^2\eta ^2}{10}\frac{1}{2}\right).$$
(30)
With this expression for the scalar field, the other fields are

$`\sqrt{\lambda }\eta B`$ $`=`$ $`{\displaystyle \frac{H}{\sqrt{K}}}sinh(\sqrt{\lambda }\eta \sqrt{K}r)(\text{for }K>0),`$ (31)
$`\sqrt{\lambda }\eta B`$ $`=`$ $`{\displaystyle \frac{H}{\sqrt{K}}}\mathrm{sin}(\sqrt{\lambda }\eta \sqrt{K}r)(\text{for }K<0),`$ (32)
$`C_0^2`$ $`=`$ $`\left[{\displaystyle \frac{1}{2}}+\left({\displaystyle \frac{1}{\beta }}{\displaystyle \frac{1}{2}}\right){\displaystyle \frac{\varphi _0^2}{\eta ^2}}\right]^1,`$ (33)
$`w_0^2`$ $`=`$ $`{\displaystyle \frac{C_0^2}{2}}\left(1{\displaystyle \frac{\varphi _0^2}{\eta ^2}}\right),`$ (34)
where

$$K=\frac{\kappa ^2\varphi _0^21}{4C_0^2}.$$
(35)
Note that it follows from Eq. (34) that $`\varphi _0<\eta `$, and it follows from (35) that $`\varphi _0>\kappa ^1`$ for $`K>0`$ solutions and $`\varphi _0<\kappa ^1`$ for $`K<0`$ solutions.

In the asymptotic solutions (31) and (32), the expansion rate $`H`$ can be absorbed by rescaling the time coordinate $`t`$ in Eq. (22); then the solution takes the form

$`(\lambda \eta ^2)ds^2`$ $`=`$ $`K^1(sinh^2\chi d\overline{s}_+^2+d\chi ^2)+C_0^2d\mathrm{\Omega }_2^2(\text{for }K>0),`$ (36)
$`=`$ $`|K|^1(\mathrm{sin}^2\chi d\overline{s}_+^2+d\chi ^2)+C_0^2d\mathrm{\Omega }_2^2(\text{for }K<0),`$ (37)
where $`\chi =\sqrt{\lambda }\eta \sqrt{|K|}r`$ and $`d\overline{s}_+^2`$ is the 4D de Sitter metric with $`H=1`$. Note that, although $`H`$ drops out of the asymptotic solution, it is a meaningful parameter for the full spacetime. With our boundary conditions at $`r=0`$, it gives the inflation rate in the core of the monopole.

For $`K<0`$, the above metric is essentially the same as the one found in Ref.  as a solution for a global defect with $`\varphi =\eta `$ and a positive bulk cosmological constant. It was also obtained in Ref.  as a cigar solution for a global monopole with $`\varphi _0<\eta `$. In the latter case, the role of the cosmological constant is played by the scalar field potential, $`V(\varphi _0)>0`$. For a gauge monopole, there is in addition the gauge field energy density, which approaches a constant at large $`r`$.

The solution (37) has an apparent singularity at $`\chi =\pi `$, but as it was noted in Ref. , the first two terms in Eq. (37) describe a 5D de Sitter space. The 5D inflation rate is $`H_5=\sqrt{\lambda }\eta \sqrt{|K|}`$, and the surface $`\chi =\pi `$ corresponds to the de Sitter horizon. This shows that the space outside the monopole core inflates not only along the 4D worldsheet, but also in one of the transverse directions, while two of the extra dimensions remain compactified in a sphere. The effective 5D cosmological constant is

$$\mathrm{\Lambda }_5=6\lambda \eta ^2K,$$
(38)
so the geometry of the uncompactified dimensions is flat for $`K=0`$ and anti-de Sitter for $`K>0`$.

### C The parameter space of asymptotic solutions

We found in Ref.  that in the global monopole case, $`\eta _c=1/\kappa `$ and solutions of the form (37) exist only in the range $`1\kappa ^2\eta ^2<5`$. For a gauge monopole, we also find that asymptotic cigar solutions exist only in a restricted parameter space of $`\eta `$ and $`\beta `$. The conditions for a solution to exist are $`D0`$, $`\varphi _0^2>0`$, $`C_0^2>0`$, and $`w_0^2>0`$.

The $`\eta \beta `$ parameter space is partitioned by the lines $`D=0`$, $`\varphi _0=0`$, $`C_0=0`$, $`w_0=0`$, and $`K=0`$, as shown in Fig. 12. The vertical straight lines, $`\kappa ^2\eta ^2=5`$ and $`\kappa ^2\eta ^2=5/4`$ correspond to $`\varphi _0=0`$ and $`w_0=0`$, respectively. The curve $`D=0`$ is split into two sectors. The right sector tangents the $`\kappa ^2\eta ^2=5`$ line at $`\beta =1/2`$, while the left sector tangents the $`\kappa ^2\eta ^2=5/4`$ line at $`\beta =8`$ and asymptotes to $`\kappa ^2\eta ^2=1`$ as $`\beta \mathrm{}`$. Finally, the curve $`K=0`$ is given by

$$\beta =\frac{2}{(\kappa ^2\eta ^21)^2}$$
(39)
and also asymptotes $`\kappa ^2\eta ^2=1`$ at large $`\beta `$.

The parameter space is thus divided into several domains with different types of solutions. In the figure, the domains are marked (a)-(f) and the corresponding solutions are: (a) $`B_s^+`$, (b) $`B_s^+`$ and $`B_s^{}`$, (c) $`B_s^{}`$ and $`B_{sh}^+`$, (d) $`B_{sh}^+`$, (e) $`B_{sh}^+`$ and $`B_{sh}^{}`$, (f) $`B_{sh}^{}`$ and $`B_s^+`$. The superscripts $`\pm `$ indicate the choice of sign in the expression (29) for $`\varphi _0`$, and the subscripts indicate whether the corresponding solution for $`B`$ is a sine function ($`K<0`$) or a sinh function ($`K>0`$).

### D Nonsingular super-critical solutions

By analogy with the global monopole case, one can expect that, for a given value of $`\beta `$ and $`\eta >\eta _c(\beta )`$, there is a single value of $`H`$ that gives a nonsingular solution. We would then interpret that solution as the physical solution describing the spacetime of a super-critical gauge monopole. Of course, it is impossible to find the precise value of $`H`$ numerically, and all super-critical global monopole solutions we considered in Ref.  eventually developed a singularity. By adjusting the value of $`H`$, we were able to shift the onset of the singularity to larger and larger radii, so that the solution got closer and closer to the analytic asymptotic solution, presumably giving a better and better approximation to the nonsingular solution.

We tried to use the same strategy for gauge monopoles. We found, however, that numerical instabilities here are much more severe than in the global case. As we tried to extend the solutions from the origin towards larger values of $`r`$, they became unstable well before we got close to the singularity. We could not, therefore, get any information about the value of $`H`$ that yields the nonsingular solution.

Sample numerical solutions with the parameter values in the domains (a) and (d) are shown in Figs. 13-18 for several values of $`H`$. The solutions are extended in $`r`$ as far as the numerical accuracy allowed. All these solutions approach our analytical asymptotic solutions at large $`r`$, indicating that nonsingular super-critical inflating solutions do exist. However, our numerical solutions exhibit this behavior for a wide range of $`H`$. We cannot, therefore, conclude that the true nonsingular solution corresponds to a single value of $`H`$, and if it does, our results do not allow us to determine this value, even approximately. To make further progress, one will have to find ways to significantly improve the numerical accuracy. It may also be possible to prove the uniqueness of the nonsingular solution using the dynamical systems method employed in Refs. .

## V Conclusions

In this paper, we continued the investigation of the gravitational field of higher-dimensional defects that we started in Ref. . We have obtained numerical solutions of Einstein’s, scalar, and gauge field equations for a defect with a $`(3+1)`$-dimensional core, which has a monopole-like field configuration in the extra three dimensions. We have verified that, for symmetry-breaking scale $`\eta `$ below the critical value $`\eta _c`$, the spacetime of the defect worldsheet is flat, and the geometry of extra dimensions is asymptotically flat, in agreement with the asymptotic solution obtained in Ref. . For $`\eta =\eta _c`$, the extra dimensions degenerate into a ‘cigar’, and for $`\eta >\eta _c`$ all static solutions are singular. The critical value $`\eta _c`$ depends on the ratio of scalar and gauge couplings, $`\beta =\lambda /e^2`$, with $`\eta _c`$ approaching the Planck scale from above as $`\beta \mathrm{}`$.

In the super-critical regime, $`\eta >\eta _c`$, we found numerical solutions in which the defect core has the geometry of de Sitter space inflating at some rate $`H`$, while extra dimensions have a cigar geometry, with two of these dimensions compactified as a sphere of a fixed radius. At large distances from the core, the solutions have a very simple form, and we found exact analytic solutions of the field equations describing this asymptotic behavior.

The nature of super-critical monopole solutions is rather similar to that of super-critical global monopoles that we discussed in Ref. . There, our results indicated that, for a given $`\eta `$, there is a unique value of $`H`$ that gives a nonsingular solution, and that $`H`$ is a growing function of $`\eta `$. This behavior is opposite to that required in the Dvali, Gabadadze and Shifman (DGS) scenario for solving the cosmological constant problem. One of our goals in the present paper has been to find whether or not the situation is different for gauge monopoles. Unfortunately, we were unable to check this directly, since the numerical instabilities did not allow us to determine the dependence of $`H`$ on $`\eta `$. In fact, our numerical accuracy was not sufficient to conclude that the value of $`H`$ is fixed (or strongly constrained) by fixing $`\eta `$ and requiring regularity of the metric.

We do however have some indirect evidence indicating that higher-dimensional magnetic monopoles do not help to solve the cosmological constant problem. The DGS conjecture was originally based on the assumption that in singular static monopole solutions , the distance $`r_s`$ from the monopole center to the singularity is a growing function of $`\eta `$. Then, in the inflating solutions, where the singularity is replaced by a horizon, one expects that the expansion rate $`Hr_s^1`$ decreases with $`\eta `$. However, our analysis of static solutions shows that $`r_s`$ gets smaller as $`\eta `$ is increased for magnetic monopoles, just as it does for global monopoles. This suggests that the DGS scenario cannot be implemented using inflating defect solutions.

###### Acknowledgements.

We are grateful to Christos Charmousis, Ruth Gregory and Gia Dvali for useful discussions. The work of A.V. was supported in part by the National Science Foundation.


text:
# Crumpling of a stiff tethered membrane

## Abstract

A first-principles numerical simulation model for crumpling of a stiff tethered membrane is introduced. In our model membranes, wrinkles, ridge formation, ridge collapse, as well as the initiation of stiffness divergence, are observed. The ratio of the amplitude and wave length of the wrinkles, and the scaling exponent of the stiffness divergence, are consistent with both theory and experiment. We observe that close to the stiffness divergence there appears a crossover beyond which the elastic behavior of a tethered membrane becomes similar to that of dry granular media. This suggests that ridge formation in membranes and force-chain network formation in granular packings are different manifestations of a single phenomenon.

Crumpling of, e.g., sheets of paper, is an everyday phenomenon that one easily passes without much extra attention. For physicists, however, crumpling has for a long time been a topic of intense research. In addition to fundamental statistical mechanic properties Kantor:1986a , issues such as scaling of strength and energy t2 ; t3 , geometry and singularities Cerda:1999a ; t4 , and even acoustic emission Houle:1996a  and crumpling related phenomena in soft matter systems Helfer:2001a , have become topics of increasing interest. In this article, we present a first-principles numerical simulation model for crumpling of a stiff tethered membrane, and relate the problem to formation of force chain networks in granular media gn1 ; gn2 .

Crumpling of a thin elastic sheet or a stiff membrane demands deformation energy. At very small strains a membrane is uniformly compressed. Thin structures are, however, sensitive to buckling. At first buckling appears as wrinkles. The wave length and the amplitude of the wrinkles depend on the conditions of the external loading and on the dimensions of the membrane wrink . As crumpling proceeds the deformation energy begins to concentrate in narrow ridges and conical peaks t2 ; Cerda:1999a . Upon crumpling of a piece of paper, for example, the strain at the ridges and peaks become so large that irreversible plastic yielding takes place. After the paper is stretched out again, the ridge pattern can clearly be observed.

Another characteristic feature of membrane crumpling is that the effective stiffness of a crumpled membrane increases fast when the degree of crumpling becomes large. The ultimate limit of crumpling would be to press a membrane into a more or less cubic form with a volume corresponding to that of the volume of the membrane itself. In practice this limit can hardly be reached. Squeezing a piece of paper as hard as possible between the hands, results in a ball that typically contains about 75% air mylar . The behavior of the effective stiffness has been investigated by carefully monitored experiments, and for a crumpled membrane it displays a power-law divergence as the dense packing limit is approached mylar . The suggested explanation to the power-law divergence is that the total ridge length diverges in a scale invariant way as the porosity inside the crumpled membrane vanishes.

In order to study wrinkles, ridge formation, and the effective stiffness in crumpled membranes, we constructed a numerical model for stiff tethered membranes. The tethered membrane consists of frictionless spheres with stiffness $`Y_s`$ arranged in a triangular lattice. Each sphere has a mass $`m`$ and moment of inertia $`I`$, and is connected to its neighbors by massless elastic beams which have both bending and tensile stiffness. It is reasonable to use a higher Young’s modulus for the spheres ($`Y_s`$) than for the beams ($`Y_b`$). In the simulations here we chose $`Y_s=10^7`$ and $`Y_b=5\times 10^3`$. The membrane is placed horizontally (i.e., in the $`xy`$ plane) in the middle of a rigid cube of size $`X^3(t)`$, where $`X`$ is time dependent. Small random fluctuations of the $`z`$ component of the centers of mass of the spheres are introduced to avoid perfect symmetry. Newtonian dynamics is applied.

The environment inside the box is set to be strongly dissipative (the membrane can be thought of being immersed in a highly viscous fluid). This means that a force proportional to velocity, ($`d\dot{\stackrel{}{x}}`$), slows down the motion of each sphere. The membrane dimensions are considered to be large enough for temperature fluctuations to be negligible (e.g. a paper sheet). For computational reasons the membrane cannot contain more than about $`10^5`$ degrees of freedom, which sets a limit on its size. Wrinkles and the onset of ridge formation is nevertheless observed. Also the beginning of a power-law divergence of the effective stiffness with an exponent close to the one observed experimentally is seen. However, as porosity reaches a low enough value, the effective stiffness begins to deviate from the experimental power law. At low porosities the effective stiffness of a tethered membranes behaves in the same way as that of a packing of solid spheres makse .

Figure 1 shows a snapshot from a simulation displaying a few wrinkles at $`t=200`$. By minimizing the elastic energy it can be demonstrated that the relation between the wave-length $`(\lambda )`$ and the amplitude $`(A)`$ of the wrinkles is generally given by $`A/\lambda (\mathrm{\Delta }/W)^{1/2}`$, where $`W`$ is the width and $`\mathrm{\Delta }`$ the compression of the membrane perpendicular to the wrinkles wrink . This relation is tested in Fig. 2 for the model membranes. The square-root relation between $`A/\lambda `$ and $`\mathrm{\Delta }`$ holds up to $`\mathrm{\Delta }6`$ (here $`W23.4`$).

Deformation of ridges in controlled geometries has been thoroughly investigated both numerically and theoretically t2 ; t1 ; t3 ; t4 ; t5 ; t6 ; t7 ; t8 . The elastic deformation energy of a ridge of length $`L`$ scales as $`\kappa (L/\delta )^{1/3}`$, where $`\delta `$ is the thickness of the sheet and $`\kappa `$ its bending modulus. By taking advantage of this result, and following the derivation in Ref. mylar , the scaling of stiffness of a membrane can be estimated: A crumpled sheet is divided into facets by the ridges. Ridges surrounding a facet of size $`L^2`$ have a length proportional to $`L`$. Such facets should fill a space of size $`L^3`$. This leads to a bulk stiffness divergence of the form $`K=VdP/dVV^{11/3}`$, when the volume $`V`$ surrounding the sheet decreases under the influence of pressure $`P`$. This is not quite consistent with experimental findings. In Ref. mylar  $`K`$ was reported to scale as $`K(V_cV)^\alpha `$ with $`\alpha 2.85`$. Notice that this scaling relation can also be expressed in the forms $`F(\varphi _c\varphi )^{\alpha +1}`$ and $`F(XX_c)^{\alpha +1}`$, where $`F`$ is the external force applied on the membrane, $`\varphi `$ the solid volume fraction, and $`X(t)`$ the time dependent linear dimension of the crumpling membrane.

To investigate possible ridge formation in the simulations, we compare the local bending to the local energy. We define the sum of the dot products of nearest neighbor bonds along the three principal directions in a triangular lattice as a measure of the local bending of the membrane,

$$p(\stackrel{}{x},t)=\frac{1}{3}\underset{i=1}{\overset{3}{}}\frac{\stackrel{}{l}_i^n(t)\stackrel{}{l}_i^m(t)}{l_i^nl_i^m}.$$
(1)
This is then compared to the local elastic energies of the beams

$$E_b(\stackrel{}{x},t)=\underset{i=1}{\overset{6}{}}_0^t𝐅(Y_b,\stackrel{}{x}(t))𝑑\stackrel{}{x}(t).$$
(2)
In Eq. (1), $`l^m`$ and $`l^n`$ are the lengths of the beams on the opposite sides of a lattice site. In Eq. (2), the elastic force $`𝐅(Y_b,\stackrel{}{x}(t))`$ is integrated over the corresponding displacement and summed over the six degrees of freedom for each site.

In addition, the elastic deformation energy of spheres in contact can be used to identify ridges. It can be written as

$$E_s(\stackrel{}{x},t)=\underset{i=1}{\overset{n}{}}\frac{1}{2}Y_s\delta _s^2,$$
(3)
where $`Y_s`$ is the stiffness of the spheres and $`\delta _s`$ the virtual ’overlap’ of spheres in contact. The sum is taken over all contacts.

Figure 3 shows contour plots of the lowest values of $`p(\stackrel{}{x},t)`$ and the highest values of $`E_b(\stackrel{}{x},t)`$ and $`E_s(\stackrel{}{x},t)`$ at an early stage of the crumpling process ($`t=400`$). It is quite obvious that $`E_b`$ and $`E_s`$ reach their highest values at the locations of high local bending. The patterns are quite complex, but elongated structures, i.e. ridges, can be seen in the figure. As crumpling proceeds, ridges grow and collapse and the membrane folds. For later crumpling stages it is no longer possible to visibly detect ridges as the patterns become too complex.

The effective bulk stiffness of a membrane can be monitored, e.g., through its stress-strain curve. This is simply the total force $`F`$ on the walls of the box surrounding the membrane, divided by the linear size $`W`$ of the membrane ($`F/W`$) as a function of $`X(t)`$. Initially $`F(t)/W`$ increases fairly linearly until the membrane begins to buckle and fold. As the dense packing limit is approached, $`F/W`$ begins to increase rapidly. Finally, when the membrane almost fills the box, a saturation regime is entered.

It is important to notice that there are three types of dense packing limits: 1) The tightest packing is achieved when the membrane volume is equal to the volume of its container. In this limit, the solid volume fraction is $`\varphi _1=1.0`$. 2) Hard-core spheres. Here, the tightest packing is an FCC lattice with $`\varphi _2=\pi \sqrt{2}/60.74`$. 3) A random dense packing of spheres with $`\varphi _30.63`$. For membranes it seems that the stiffness divergence takes place at $`\varphi _c0.75`$ which is close to $`\varphi _2`$. The regime of saturating stiffness is, however, entered well before this limit is reached.

Figure 4A shows $`F/(Y_sW)`$ for membranes with $`10^4`$ and $`4\times 10^4`$ degrees of freedom (each simulation demands about $`10^610^7`$ time steps). The simulation results are compared to the experimental result reported by Matan et al. mylar , who found that $`F(X)(XX_c)^{\alpha +1}`$, with $`\alpha 2.85`$. Figure 4A shows good agreement between simulations and experiment until saturation sets in at $`(XX_c)/X_c0.03`$.

At saturation there is a crossover to a different type of stiffness behavior. Beyond saturation the beams no longer affect the stiffness of the membranes. Instead, the membranes begin to behave like granular packings makse . The stiffness of a 3D packing of solid spherical grains scales like $`(\varphi \varphi _c)^\beta `$, with $`\beta 1.62`$, and the critical volume is $`\varphi _c0.60.7`$ depending on the form of interaction between the grains. Figure 4B shows a comparison between simulations and the scaling behavior of $`F(\varphi )`$ in the two mentioned regimes.

By following the energy contents in a membrane from the initially flat membrane to a fully crumpled state, it is possible to determine the energy component responsible for the stiffness divergence seen in Fig. 4. The external energy $`_tF(t)\left[X(t)X(0)\right]𝑑t`$ is compared to the internal energy of the membrane in Fig. 5A. There is no difference between the external and internal energies as energy conservation also demands.

The internal energy can be split into five terms that are related to the force terms in the general equation of motion

$$m\ddot{\stackrel{}{x}}+d\dot{\stackrel{}{x}}+𝐊\stackrel{}{x}+Y_s\stackrel{}{\delta }_s=\stackrel{}{f}_{ext}.$$
(4)
Here, $`x`$ is the displacement, $`m`$ the mass, $`d`$ the damping coefficient, $`𝐊`$ the stiffness matrix, $`\stackrel{}{\delta }_s`$ the deformations (or virtual overlap) of spheres in contact, and $`\stackrel{}{f}_{ext}`$ are forces applied on the walls of the box surrounding the membrane. The related energy terms are the following: $`(\mathrm{𝟏})`$ The elastic energy of the beams, $`E_b`$. $`(\mathrm{𝟐})`$ The compression energy at the contacts between spheres, $`E_s`$. $`(\mathrm{𝟑})`$ The compression energy of the spheres against the external walls, $`E_{ext}=0.5Y_w\delta _w^2`$. $`(\mathrm{𝟒})`$ The energy dissipated through damping, $`E_d=_t\stackrel{}{F}_d𝑑\stackrel{}{x}`$. $`(\mathrm{𝟓})`$ The kinetic energy of the spheres, $`E_m=0.5m\dot{\stackrel{}{x}}^2`$. These energy components are shown separately as functions of time in Figs. 5B and 5C. The diverging terms are $`E_s`$ and $`E_{ext}`$. The divergence of $`E_{ext}`$ is rather trivial as it is just the deformation energy of the box surrounding the membrane. It is thus the energy $`E_s`$ that causes the stiffness divergence. This explains to some extent why there is a crossover to the stiffness behavior of granular packings. When the deformation energy at the contacts between spheres dominate over other energy terms, the beams no longer play a role in the membrane stiffness.

Since there are two different scaling behaviors of stiffness, and both regimes are dominated by $`E_s`$, one would expect that at the crossover point there appears a qualitative change in the function describing the distribution of contact energies of the spheres. There is, however, no indication of such a qualitative change as is demonstrated in Fig. 6, which means that, at the crossover, ridges are simultaneously force chains. When comparing the topology of the network of ridges in crumpled membranes and that of the force-chain network in dry granular materials gn1 ; gn2 , one notices that both are essentially one-dimensional lines of localized deformation energy, and form complicated networks in three-dimensional space. Based on these observations, it is plausible to draw the conclusion that the ridge formation and force-chain formation are not separate phenomena, but different manifestations of a single phenomenon. The crossover between the two then only marks the region where the membrane thickness begins to limit the density of network lines.

In conclusion, we have demonstrated for the first time that a numerical model of a tethered membrane can reproduce the theoretically predicted wrinkling and ridge formation, and experimentally observed stiffness divergence. The early-stage ridges can be traced both by the elastic deformations of the beams and by the contact deformations of spheres in a tethered membrane. The divergence of the effective stiffness is a result of the divergence of the contact energy of the spheres. This eventually leads to a crossover in the stiffness between two different scaling behaviors. When the solid volume fraction is far below that of the random dense limit of granular packings, the scaling of stiffness is ’membranic’, and when the random dense packing is approached, the stiffness scaling becomes ’granular’. There is no obvious difference in the distribution of deformation energy in the two regimes. It thus seems that deformation ridges in membranes and force chains in granular packings are different manifestations of the same phenomenon.


text:
# Formation of multi-solitons and vortex bright solitons in Bose-condensed alkali-metal atoms

## I Introduction

In the last years we have investigated various dynamical properties of pure Bose-Einstein condensates (BECs) made of alkali-metal atoms at ultra-low temperatures . The recent experimental observation of dark  and bright  solitons in Bose-Einstein condensates (BECs) has focused our interest on BEC bright solitons. In particular, by using the nonpolynomial Schrödinger equation (NPSE), an effective equation derived from the 3D Gross-Pitaevskii equation (GPE) , we have studied the stability properties of bright solitons, determining the critical strength beyond which bright solitons collapse . This collapse is a direct consequence of the non-constant size of the soliton in the transverse direction, as found by solving the 3D GPE and NPSE. Instead, the 1D GPE, for which the transverse size of the soliton is constant, admits only stable bright solitons .

In this paper we analyze the formation of bright solitons in a Bose-Einstein condensate of 7Li atoms induced by a sudden change in the sign of the scattering length from positive to negative, as reported in a recent experiment . In a previous paper  we have show that a number of bright solitons is produced and this can be interpreted in terms of the modulational instability of the time-dependent macroscopic wave function of the Bose condensate. From the 1D GPE we have derived a simple formula for the number of solitons that is in good agreement with the numerical results. Here we generalize such a formula taking into account the non-constant size of the BEC in the transverse direction. Then we investigate the properties of vortex bright solitons, which are bright solitons with a finite quantum number of circulation around the longitudinal axis.

## II Formation of BEC bright solitons

At zero temperature the macroscopic wave function of a BEC made of $`N`$ dilute 7Li atoms can be modeled by the 3D GPE. For a dilute gas the inter-atomic strength depend on the single parameter: the s-wave scattering length $`a_s`$. Note that only for $`a_s<0`$ there exist bright solitons, which are shape-invariant waves traveling without spreading due to the interplay between the repulsive kinetic energy and the attractive inter-atomic energy.

Following the experiment of Strecker et al. , the trapping potential of the BEC is modeled by an anisotropic harmonic potential with transverse frequency $`\omega _{}=2\pi \times 800`$ Hz and longitudinal frequency $`\omega _z=2\pi \times 4`$ Hz. In addition there is a box optical potential that initially confines the condensate in the longitudinal direction.

By using 3D GPE we first calculate the ground-state wave function of the condensate with a positive scattering length $`a_s=100a_B`$, where $`a_B=0.53\AA `$ is the Bohr radius. We choose a condensate with longitudinal width $`L=284.4`$ $`\mu `$m and $`N=10^4`$ atoms. Then the ground-state wave function is inserted as initial condition for the time-evolution of the 3D GPE with $`a_s=3a_B`$. We observe that a three-body dissipative term has been added to the 3D GPE. This term is important mainly during the collapse because it fixes the critical density at which the compression ceases (see also ).

In Figure 1 we plot the probability density in the longitudinal direction $`\rho (z)`$ of the evolving wave-function. The Bose-Einstein condensate shows the formation of 5 peaks. It is remarkable to stress that at the first instants the two external peaks become very high but the collapse is avoided by the dissipative term that quikly reduces their density. We have verified that without the dissipative term the Bose condensate collapses. The dissipative term does not alter the number of peaks but only their density. After the formation, the peaks start to separate each other showing a repulsive force between them. Eventually each peak evolves with small shape oscillations but without appreciable dispersion.

The formation of these solitonic peaks can be explained as due to MI of the time-dependent wave function of the BEC, driven by imaginary Bogoliubov excitations . Neglecting the dissipative term, in the 1D limit of very long cigar-shaped condensate, the Bogoliubov longitudinal elementary excitations $`ϵ_k`$ of the axially homogeneous BEC can be obtained from the 1D GPE. They are given by

$$ϵ_k=\sqrt{\frac{\mathrm{}^2k^2}{2m}\left(\frac{\mathrm{}^2k^2}{2m}+2g_{1d}n\right)},$$
(1)
where $`n=N/L`$ is the linear axial density and $`g_{1d}=2(\mathrm{}\omega _{})a_s`$ is the 1D interaction strength. By suddenly changing the scattering length $`a_s`$ to a negative value, the excitations frequencies corresponding to $`k<k_c=\left(16\pi |g_{1d}|n\right)^{1/2}`$ become imaginary and, as a result, small perturbations grow exponentially in time. It is easy to find that the maximum rate of growth is at $`k_0=k_c/\sqrt{2}`$. The wavelength of this mode is $`\lambda _0=2\pi /k_0`$ and the ratio $`L/\lambda _0`$ gives an estimate of the number $`N_s`$ of bright solitons which are generated:

$$N_s=\frac{\sqrt{N|a_s|L}}{\pi a_{}},$$
(2)
with $`a_{}=(\mathrm{}/(m\omega _{}))^{1/2}`$. The predicted number $`N_s`$ of solitons is in very good agreement with the numerical results of Figure 1 (see also ).

The Eq. (2) is valid in the 1D limit but can be extended to a generic 3D case by using NPSE instead of 1D GPE. In general, the Bogoliubov excitations can be written as

$$ϵ_k=\sqrt{\frac{\mathrm{}^2k^2}{2m}\left(\frac{\mathrm{}^2k^2}{2m}+2c^2\right)},$$
(3)
where $`c`$ is the first sound velocity in the longitudinal direction. The sound velocity $`c`$ satisfies the equation

$$c=\sqrt{n\frac{\mu }{n}},$$
(4)
where $`\mu `$ is the chemical potential of the Bose condensate. The chemical potential derived from 1D GPE is $`\mu =g_{1d}n+\mathrm{}\omega _{}`$, while using NPSE one finds

$$\mu =\frac{g_{1d}n}{\sqrt{1+2a_sn}}+\frac{\mathrm{}\omega _{}}{2}\left(\frac{1}{\sqrt{1+2a_sn}}+\sqrt{1+2a_sn}\right).$$
(5)
This chemical potential reduces to the 1D GPE one in the 1D limit, i.e. for $`2a_sn<<1`$. Moreover, if $`a_s<0`$ the chemical potential $`\mu `$ diverges at $`2|a_s|n=1`$, as a consequence of the collapse on the condensate . In conclusion, NPSE tells us that the number $`N_s`$ of bright solitons obtained via modulational instability is given by

$$N_s=\frac{g\sqrt{2}}{4\pi }\frac{\sqrt{43g\rho }}{\sqrt{g\rho }\left(1g\rho \right)^{3/4}},$$
(6)
where $`g=2|a_s|/a_{}`$ is the adimensional inter-atomic strength and $`\rho =na_{}`$ is the adimensional axial density of the condensate. In Figure 2 we plot $`N_s`$ as a function of $`g\rho `$ for a fixed value of $`\rho `$. For small values of $`g\rho `$ 1D GPE and NPSE give the same results but approaching $`g\rho =1`$ the results are strongly different. In fact, NPSE gives a divergent number of solitons at $`g\rho =1`$. As previously stated, only NPSE takes into account the non-constant width of the BEC in the transverse direction. Such a width shrinks to zero for a large attractive inter-atomic interaction.

## III Vortex bright solitons

The 3D GPE of a Bose-Einstein condensate of $`N`$ atoms confined in the transverse direction by a harmonic potential of frequency $`\omega _{}`$ and with an attractive inter-atomic interaction ($`a_s<0`$) can be derived by the following scaled GPE Lagrangian

$$L=\psi ^{}\left[i\frac{}{t}+\frac{1}{2}^2\frac{1}{2}(x^2+y^2)2\pi \gamma |\psi |^2\right]\psi $$
(7)
where $`\psi `$ is the macroscopic wave function of the condensate and $`\gamma =Na_s/a_{}`$, with $`a_{}=(\mathrm{}/(m\omega _{}))^{1/2}`$. Note that in the GPE Lagrangian lengths are in units $`a_{}`$, time in units $`\omega _{}^1`$, action in units $`\mathrm{}`$ and energy in units $`\mathrm{}\omega _{}`$.

Using the cylindrical coordinates $`r=(x^2+y^2)^{1/2}`$, $`\theta =arctg(y/x)`$ and $`z`$, we write the trial wave function of the vortex bright soliton as

$$\psi (r,z,\theta ,t)=r^ke^{ik\theta }\mathrm{exp}\left\{\frac{r^2}{2\sigma (t)^2}\frac{z^2}{2\eta (t)^2}+i\alpha (t)r^2+i\beta (t)z^2\right\},$$
(8)
where $`k`$ is the vortex quantum number and $`\sigma `$, $`\eta `$, $`\alpha `$ and $`\beta `$ are the time-dependent variational parameters. In  we have used this trial wave function to study vortex states with positive scattering length in a harmonic potential. Here we consider a negative scattering length. Moreover, in this case the (harmonic) confinement is only in the transverse direction. By inserting the trial wave function in the GPE Lagrangian, after spatial integration, one gets an effective Lagrangian depending on the variational parameters. By using the Euler-Lagrange equations one realizes that the time dependence of the phase parameters $`\alpha `$ and $`\beta `$ is fully determined by that of $`\sigma `$ and $`\eta `$ (see also ). Thus, one obtains a new effective Lagrangian depending on $`\sigma `$ and $`\eta `$ only. It is given by

$$L_{eff}=(k+1)\dot{\sigma }^2+\frac{1}{2}\dot{\eta }^2(k+1)\sigma ^2\frac{(k+1)}{\sigma ^2}\frac{1}{2\eta ^2}\sqrt{\frac{2}{\pi }}\frac{(2k)!}{2^{2k}(k!)^2}\frac{\gamma }{\sigma ^2\eta }.$$
(9)
In the case $`k=0`$ this effective Lagrangian reduces to the one we have successfully used to analyze the collective oscillations of the bright soliton in . In fact, the stationary solutions of the Euler-Lagrange equations of $`L_{eff}`$ give the widths of the bright soliton ($`k=0`$) and the widths of the vortex bright soliton ($`k0`$).

In Figure 3 we plot the transverse radial width $`\sigma `$ and the longitudinal axial width $`\eta `$ of the vortex bright soliton as a function of the inter-atomic strength $`\gamma `$. For values of $`\gamma `$ close to zero the axial width $`\eta `$ is very large with respect to the radial width $`\sigma `$: the soliton is a strongly elongated cigar. By decreasing $`\gamma `$ the soliton anisotropy is reduced up to the collapse, i.e. the critical point at which the width parameters are no more minima of the GPE potential energy. It is interesting to observe that the critical $`\gamma `$ changes with the vortex quantum number $`k`$: $`\gamma _c=0.78`$ for $`k=0`$, $`\gamma _c=2.20`$ for $`k=1`$, and $`\gamma _c=3.59`$ for $`k=2`$. This means that it is possible to avoid the collapse by creating vortex bright solitons with a large vortex quantum number $`k`$ (see also the recent numerical results  of Adhikari). Another remarkable effect shown in Figure 3 for $`k0`$ is the crossing of the widths by decreasing $`\gamma `$: the vortex soliton becomes slightly disc-shaped.

In the lower part of Figure 3 we plot the frequencies $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ of the collective oscillations of the vortex bright soliton along the axial and transverse directions. In cigar-shaped condensates $`\mathrm{\Omega }_1`$ is called monopole frequency and $`\mathrm{\Omega }_2`$ is called quadrupole frequency . These frequencies are obtained from the effective Lagrangian $`L_{eff}`$ linearizing around the minima of the potential energy (equilibrium solution). The transverse frequency $`\mathrm{\Omega }_2`$ remains practically constant while the axial frequency $`\mathrm{\Omega }_1`$ grows by decreasing $`\gamma `$ but goes to zero at the critical point, beyond which the frequencies are complex and the equilibrium solution becomes unstable. Note that for a fixed $`\gamma `$, the value of $`\mathrm{\Omega }_1`$ strongly depends on $`k`$.

## IV Conclusions

In this paper we have shown that a train of bright solitons can be created in a Bose-Einstein condensate by a sudden change of the scattering length from positive to negative. The number of bright solitons can be predicted by using an analytical formula we have derived from the nonpolynomial Schrödinger equation. We have also investigated vortex bright solitons, which are bright solitons with a finite angular momentum along the longitudinal axis. Our results suggest that vortex bright solitons can be created with a larger number of atoms than bright solitons without vorticity. Moreover, our calculations suggest that vortex bright solitons can be experimentally detected by analyzing their collective oscillations.

## Figure captions

Figure 1: Axial density profile $`\rho (z)`$ of the Bose-Einstein condensate made of $`10^4`$ 7Li atoms obtained by solving 3D GPE. For $`t<0`$ the scattering length is $`a_s=100a_B`$, while for $`t0`$ we set $`a_s=3a_B`$ with $`a_B`$ the Bohr radius. Length $`z`$ in units $`a_z=(\mathrm{}/m\omega _z)^{1/2}`$, density $`\rho `$ in units $`1/a_z`$ and time $`t`$ in units $`\omega _z^1`$. Adapted from .

Figure 2: Number $`N_s`$ of bright solitons created via modulational instability. $`\rho =Na_{}/L`$ is the initial (adimensional) axial density of the axially homogeneous Bose condensate. $`g=2|a_s|/a_{}`$ is the inter-atomic strength with $`a_{}=(\mathrm{}/m\omega _{})^{1/2}`$ and $`a_s<0`$ the final scattering length.

Figure 3: Top: transverse width $`\sigma `$ and axial width $`\eta `$ of the vortex bright soliton as a function of the interaction strength $`\gamma =Na_s/a_{}`$. Bottom: monopole frequency $`\mathrm{\Omega }_1`$ and quadrupole frequency $`\mathrm{\Omega }_2`$ of the vortex bright soliton as a function of the interaction strength $`\gamma `$. Lengths in units $`a_{}`$ and frequencies in units $`\omega _{}`$. $`k`$ is the vortex quantum number.