text:
# Current and Shot Noise in a Quantum Dot Coupled to Ferromagnetic Leads in the Large 𝑈 Limit

## Abstract

Using the Keldysh nonequilibrium Green function technique, we study the current and shot noise spectroscopy of a single interacting quantum dot coupled to two ferromagnetic leads with different polarizations. The polarizations of leads can be both parallel and antiparallel alignments. General formulas of current and shot noise are obtained, which can be applied in both the parallel and antiparallel alignment cases. We show that for large polarization value, the differential conductance and shot noise are completely diferent for spin up and spin down configurations in the parallel alignment case. However, the differential conductance and shot noise have the similar properties for parallel alignment case in the small polarization value and for antiparallel alignment case in any polarization value.

Spin-related phenomena in semiconductor quantum dot have attracted great interest recently as they are the crucial ingredient in the emerging field of spintronics Wolf01  and several quantum computation scheme Recher00 . In addition to their potential industrial applications, these devices also provide an ideal test ground for the study of basic physics including many-body effect, such as the Kondo effect Ng and Glazman88 . In this Letter, we use the Keldysh nonequilibrium Green function technique to study the current and shot noise through an interacting dot coupled to two ferromagnetic (FM) leads as a function of the applied bias voltage for parallel (P) and antiparallel (AP) lead-polarization alignments. Using the equation-of-motion approach, we obtain a general formula of current and shot noise for interacting quantum dot, which can be applied in studying the transport phenomena of dot coupled to FM leads with both the P and the AP alignments. Our results show that both the differential conductance and shot noise show the completely different behavior for spin up and spin down configurations in P alignment case with large polarization value. However, the differential conductance and shot noise show similar behavior for P alignment in small polarization value and for AP alignment in any polarization value.

The system Hamiltonian is written as

$$H=H_L+H_R+H_D+H_T.$$
(1)
The Hamiltonian for electrons in the left and right non-interacting metallic leads is

$$H_L+H_R=\underset{kL,R;\sigma }{}ϵ_{k\sigma }c_{k\sigma }^{}c_{k\sigma },$$
(2)
where the electron creation (annihilation) operators in the leads are denoted by $`c_{k\sigma }^{}`$ ($`c_{k\sigma }`$). The Hamiltonian of the dot is

$$H_D=\underset{\sigma }{}[ϵ_{d\sigma }d_\sigma ^{}d_\sigma +\frac{U}{2}n_{d,\sigma }n_{d,\overline{\sigma }}],$$
(3)
where $`d_\alpha ^{}`$ ($`d_\alpha `$) are the creation (annihilation) operators of dot electrons, and $`ϵ_{d\sigma }`$ is the resonance level of the dot which can be tuned by magnetic field. The coupling of the dot to the leads is

$$H_T=\underset{kL,R;\sigma ,\sigma ^{}}{}[V_{k\sigma ,\sigma ^{}}c_{k\sigma }^{}d_\sigma ^{}+H.c.],$$
(4)
where the tunneling matrix elements $`V_{k\sigma ,\alpha }`$ transfer electrons through an insulating barrier out of the dot.

Using the Keldysh nonequilibrium Green function formalism Caroli71 ; Keldysh65 , the terminal current is given by Meir92 ; Jauho94

$`I`$ $`=`$ $`{\displaystyle \frac{ie}{4\pi }}{\displaystyle }dϵ\{\text{Tr}[(f_L𝚪^Lf_R𝚪^R)(𝐆^r𝐆^a)]`$ (5)
$`+\text{Tr}[(𝚪^L𝚪^R)𝐆^<]\},`$
and the spectral density of shot noise in the zero-frequency is given by Zhu02

$`S(\omega 0)={\displaystyle \frac{e^2}{2\pi }}{\displaystyle }dϵ\{f_L(1f_L)(\text{Tr}[(𝚪^L𝐆^r)^2]`$
$`+\text{Tr}[(𝚪^L𝐆^a)^2])+if_L\text{Tr}[𝚪^L𝐆^>]i(1f_L)\text{Tr}[𝚪^L𝐆^<]`$
$`+f_L\text{Tr}[𝚪^L𝐆^>𝚪^L(𝐆^r𝐆^a)]+\text{Tr}[𝚪^L𝐆^>𝚪^L𝐆^<]`$
$`(1f_L)\text{Tr}[𝚪^L(𝐆^r𝐆^a)𝚪^L𝐆^<]\},`$ (6)
where $`f_{L(R)}`$ are the Fermi distribution function of the left and right leads, which has different chemical potential upon a voltage bias $`\mu _L\mu _R=eV`$. The coupling of the dot to the leads is characterized by the parameter

$$\mathrm{\Gamma }_{\sigma \sigma ^{}}^{L(R)}=2\pi \underset{\sigma ^{\prime \prime }}{}\rho _{L(R),\sigma ^{\prime \prime }}(ϵ)V_{\sigma ^{\prime \prime },\sigma }^{}(ϵ)V_{\sigma ^{\prime \prime },\sigma ^{}}(ϵ)$$
(7)
with $`\rho _{L(R),\sigma }`$ the spin-$`\sigma `$ band density of states in the two leads. $`G_{\sigma \sigma ^{}}^{r(a)}`$ and $`G_{\sigma \sigma ^{}}^{<(>)}`$ are the Fourier transform of the dot electron retarded (advanced) Green function $`G_{\sigma \sigma ^{}}^{r(a)}(t,t^{})=i\theta (\pm tt^{})\{d_\sigma (t),d_\sigma ^{}^{}(t^{})\}`$, the lesser Green function $`G_{\sigma \sigma ^{}}^<(t,t^{})=id_\sigma ^{}^{}(t^{})d_\sigma (t)`$, and the greater Green function $`G_{\alpha \alpha ^{}}^>(t,t^{})=id_\alpha (t)d_\alpha ^{}^{}(t^{})`$. It is noted that Eqs. (5) and (6) expresses the current and the fluctuations of current through the quantum dot, an interacting region, in terms of the distribution functions in the leads and local properties of the quantum dot, such as the occupation and density of states.

In order to compute the current and the shot noise, one has to compute the dot electron retarded Green function $`G^r`$ and Keldysh Green function $`G^<`$ in the presence of Coulomb interaction $`U`$. Without loss of physics we are considering here, we assume that the tunneling matrix elements are spin independent, $`V_{k\sigma ,\sigma ^{}}=V_k\delta _{\sigma \sigma ^{}}`$. Using the equation-of-motion approach, we obtain the retarded dot Green function in the large $`U`$ limit as:

$$G_\sigma ^r(\omega )=G_{\sigma \sigma }^r(\omega )=\frac{1n_{d,\overline{\sigma }}}{\omega ϵ_{d\sigma }\mathrm{\Sigma }_{0\sigma }^r\mathrm{\Sigma }_{1\sigma }^r+i0^+}.$$
(8)
Here

$`\mathrm{\Sigma }_{0\sigma }^r(\omega )`$ $`=`$ $`{\displaystyle \underset{kL,R}{}}{\displaystyle \frac{|V_k|^2}{\omega ϵ_{k\sigma }+i0^+}}`$ (9)
$`=`$ $`{\displaystyle \frac{dϵ}{2\pi }\frac{\mathrm{\Gamma }_\sigma ^L(ϵ)+\mathrm{\Gamma }_\sigma ^R(ϵ)}{\omega ϵ+i0^+}},`$
$`\mathrm{\Sigma }_{1\sigma }^r(\omega )`$ $`=`$ $`{\displaystyle \underset{kL,R}{}}{\displaystyle \frac{|V_k|^2f_{L/R}(ϵ_{k\overline{\sigma }})}{\omega ϵ_{d\sigma }+ϵ_{d\overline{\sigma }}ϵ_{k\overline{\sigma }}+i/2\tau _{\overline{\sigma }}}}`$ (10)
$`=`$ $`{\displaystyle \frac{dϵ}{2\pi }\frac{\mathrm{\Gamma }_{\overline{\sigma }}^L(ϵ)f_L(ϵ)+\mathrm{\Gamma }_{\overline{\sigma }}^R(ϵ)f_R(ϵ)}{\omega ϵ_{d\sigma }+ϵ_{d\overline{\sigma }}ϵ+i/2\tau _{\overline{\sigma }}}},`$
where $`\mathrm{\Gamma }_\sigma ^{L(R)}(ϵ)=\mathrm{\Gamma }_{\sigma \sigma }^{L(R)}(ϵ)`$, and the occupation number is subject to the self-consistency condition $`n_{d\sigma }=i\frac{d\omega }{2\pi }G_\sigma ^<(\omega )`$. The finite life-time in Eq. (10) for finite bias voltage and magnetic field can be obtained by using the second-order perturbation theory Meir93 .

The Green function $`G_{\sigma \sigma ^{}}^{<(>)}`$ cannot be obtained by directly using the above equation-of-motion approach without introducing additional assumptions. By instead applying the operational rules as given by Langreth Langreth76  to the Dyson equation for the contour-ordered Green function, one can show the following Keldysh equation for the lesser and greater functions Jauho94

$$𝐆^{<(>)}(\omega )=𝐆^r(\omega )𝚺_T^{<(>)}(\omega )𝐆^a(\omega ).$$
(11)
From Eq. (8), the retarded self-energy is given by

$$\mathrm{\Sigma }_{T,\sigma }^r(\omega )=\frac{n_{d\overline{\sigma }}(\omega ϵ_{d\sigma })+(\mathrm{\Sigma }_{0\sigma }^r+\mathrm{\Sigma }_{1\sigma }^r)}{1n_{d\overline{\sigma }}}.$$
(12)
We then arrive at

$`\mathrm{\Sigma }_{T,\sigma }^<(\omega )\mathrm{\Sigma }_{T,\sigma }^>(\omega )=\mathrm{\Sigma }_{T,\sigma }^a(\omega )\mathrm{\Sigma }_{T,\sigma }^r(\omega )`$ (13)
$`=i(1n_{d\overline{\sigma }})^1\{(\mathrm{\Gamma }_\sigma ^L(\omega )+\mathrm{\Gamma }_\sigma ^R(\omega ))`$
$`2Im\mathrm{\Sigma }_{1\sigma }^r\left(\omega \right)\}`$ $`.`$
To determine $`\mathrm{\Sigma }_{T,\sigma }^<`$ and $`\mathrm{\Sigma }_{T,\sigma }^>`$, we follow the ansatz proposed in Ref. Ng96  to assume further that these self-energies have the form

$$\mathrm{\Sigma }_{T,\sigma }^<=i[\mathrm{\Gamma }_\sigma ^L(\omega )f_L(\omega )+\mathrm{\Gamma }_\sigma ^R(\omega )f_R(\omega )]R_\sigma (\omega ),$$
(14a)
$$\mathrm{\Sigma }_{T,\sigma }^>=i[\mathrm{\Gamma }_\sigma ^L(\omega )(1f_L(\omega ))+\mathrm{\Gamma }_\sigma ^R(\omega )(1f_R(\omega ))]R_\sigma (\omega ),$$
(14b)
where $`R(\omega )`$ can be regarded as a renormalization factor due to strong electron-electron Coulomb repulsion in the dot. A little algebra yields

$$R_\sigma (\omega )=\frac{1}{1n_{d\overline{\sigma }}}\{1\frac{2}{\mathrm{\Gamma }_\sigma }Im\mathrm{\Sigma }_{1\sigma }^r\left(\omega \right)\}$$
(15)
where $`\mathrm{\Gamma }_\sigma =\mathrm{\Gamma }_\sigma ^L+\mathrm{\Gamma }_\sigma ^R`$.

As an extension, we can assume in general the “scattering in” and “scattering out” self-energies to have the form

$$𝚺_T^<=i[𝚪^L(\omega )f_L(\omega )+𝚪^R(\omega )f_R(\omega )]𝐑(\omega ),$$
(16a)
$$𝚺_T^>=i[𝚪^L(\omega )(1f_L(\omega ))+𝚪^R(\omega )(1f_R(\omega ))]𝐑(\omega ),$$
(16b)
Eqs. (16a) and (16b) leads to

$`𝐆^r𝐆^a`$ $`=`$ $`𝐆^>𝐆^<`$ (17)
$`=`$ $`i𝐆^r(𝚪^L+𝚪^R)\mathrm{𝐑𝐆}^a.`$
Substitution of Eq. (17) and Eq. (11) with Eq. (16b) into Eqs. (5) and (6) gives rise to

$`I`$ $`=`$ $`{\displaystyle \frac{e}{2\pi }}{\displaystyle 𝑑ϵ\text{Tr}\{𝚪^L[𝐆^<+f_L(𝐆^r𝐆^a)]\}}`$ (18)
$`=`$ $`{\displaystyle \frac{e}{2\pi }}{\displaystyle }dϵ\text{Tr}\{𝚪^L[i𝐆^r(f_L𝚪^L+f_R𝚪^R)\mathrm{𝐑𝐆}^a`$
$`if_L𝐆^r(𝚪^L+𝚪^R)\mathrm{𝐑𝐆}^a]\}`$
$`=`$ $`{\displaystyle \frac{e}{2\pi }}{\displaystyle 𝑑ϵ[f_Lf_R]\text{Tr}[𝐓]},`$
and

$`S`$ $`=`$ $`{\displaystyle \frac{e^2}{2\pi }}{\displaystyle }dϵ\{[f_L(1f_L)\text{Tr}\{𝐓[1+2(𝚪^R𝐑)^1(𝐑1)]\}`$ (19)
$`+f_R(1f_R)\text{Tr}[𝐓]+(f_Lf_R)^2\text{Tr}[(1𝐓)𝐓]\}.`$
where the transmission coefficient matrix is defined as $`𝐓=𝐆^a𝚪^L𝐆^r𝚪^R𝐑`$. We remark that Eqs. (18) and (19) are the general formulas for current and shot noise in the interacting dots with large $`U`$, which can be applied in studying the spin-transport of both the P and AP configurations. It is noted that the formula of proportional coupling (i. e., $`\mathrm{\Gamma }_\sigma ^L=\lambda \mathrm{\Gamma }_\sigma ^R`$) does not work in the AP alignment case at finite bias voltage and magnetic field. However, one can easily to compute the current and shot noise for AP alignment configuration by using our formalisms Eqs. (18) and (19). It is easy to find that for the noninteracting case, the above two expressions are just the Landauer-Büttiker formalisms developed based on the scattering matrix theory. The connection between the two formalisms was first established by Meir and Wingreen for the current Meir92 , and by Zhu and Balatsky for the shot noise Zhu02 .

In the wideband limit, we assume that $`\mathrm{\Gamma }_\sigma ^{L(R)}`$ to be energy independent, but polarization dependent with $`\mathrm{\Gamma }_{}^L=\mathrm{\Gamma }_{}^R=(1+P)\mathrm{\Gamma }_0`$, $`\mathrm{\Gamma }_{}^L=\mathrm{\Gamma }_{}^R=(1P)\mathrm{\Gamma }_0`$ for P alignment, while $`\mathrm{\Gamma }_{}^L=\mathrm{\Gamma }_{}^R=(1+P)\mathrm{\Gamma }_0`$, $`\mathrm{\Gamma }_{}^L=\mathrm{\Gamma }_{}^R=(1P)\mathrm{\Gamma }_0`$ for AP alignment, where $`\mathrm{\Gamma }_0`$ describes the tunneling coupling between the dot and the nonmagnetic leads, and $`0p<1`$ characterizes the polarization of leads. In Fig. 1 we plot the zero-temperature differential conductance of spin up (solid line) and spin down (dashed line) for FM leads with P alignment as a function of bias voltage. One interesting observation is that the behavior of differential conductance is quite different for small and large $`P`$ values. For small $`P=0.2`$, Fig. 1(a) shows that the zero-bias peak of differential conductance for nonmagnetic leads Meir93  splittes into two peaks even in the absence of magnetic field, and this splitting can be tuned away by applying an appropriate magnetic field \[Fig. 1(c)\] and can be recovered by increasing the magnetic field further \[Fig. 1(e)\]. This result agrees well with that in Ref. Martinek02 , which is due to the splitting of the dot levels renormalized by the spin-dependent interacting self-energy $`\stackrel{~}{ϵ}_{d\sigma }`$ satisfying the self-consistent equation $`\stackrel{~}{ϵ}_{d\sigma }=ϵ_{d\sigma }+Re\mathrm{\Sigma }_{1\sigma }^r(\stackrel{~}{ϵ}_{d\sigma },\stackrel{~}{ϵ}_{d\overline{\sigma }})`$. For large $`P=0.8`$ \[Figs. 2 (b,d,f)\], the behavior of spin up is similar to that of small $`P=0.2`$, which shows one peak at approriate value of magnetic field and two peaks when the field is away from the appropriate value. However for the spin down configuration, the differential conductance shows a small ”flat” near $`eV=0`$ when the magnetic field is smaller than the appropriate value, while it shows a ”dip” at $`eV=0`$ when the field is at the appropriate value, and finally it shows a large maximum when the field increases further. This interesting behavior can be understood with the help of DOS for spin down near the zero bias voltage (shown as insects of Figs. 1(b,d,f)). The DOS for spin down configuration shows a dip at $`\omega =\stackrel{~}{ϵ}_d\stackrel{~}{ϵ}_d`$ because $`Re\mathrm{\Sigma }_1^r\left(\omega \right)`$ grows logarithmically at this value, and the dip position can be tuned by magnetic field around the appropriate value corresponding to the single peak for spin up configuration. In Fig. 2 we plot the differential conductance for FM leads with AP alignment. In this case, the behavior of differential conductance looks similar for spin up and spin down configurations, which has one peak in the absence of magnetic field \[Fig. 2(a,b)\] and two shifted peaks in the presence of magnetic field \[Fig. 2(c,d)\]. The small difference between the up and down spin configurations is due to the small (compared to the P alignment case) splitting of renormalized dot levels for different spin configurations at finite voltage (as shown in the inset of Fig. 2(b)).

In Figs. 3 and 4, we plot the corresponding differential shot noise for FM leads with P and AP alignments as a function of the bias voltage. For P alignment case, the differential shot noise Fiq. 3 also shows the different behavior for spin up and down configurations in large $`P`$ values, while similar behavior for small $`P`$ value. However, for AP alignment case, the differential shot noise Fig. 4 shows the similar behavior for different spin configurations. Since the transmission probability corresponding to the conductance peaks (as shown in Figs. 1 and 2) is small, the differential noise is approximately proportional to the conductance and only shown one single peak rather than two peaks.

To summarize, using the Keldysh nonequilibrium Green function technique, we have studied the current and shot noise spectroscopy of a single dot with Coulomb interaction coupled to FM leads with P and AP polarization alignments. We have shown that the lead alignments affect both the current and current fluctuations. For large polarization value, the spin up and spin down configurations have the different behavior in the differential conductance and shot noise as a function of bias voltage in P alignment case. While the differential conductance and shot noise show similar behavior for different spin configurations in P alignment with small polarization value and in AP alignment case with any polarization value. The derived current and shot noise formulism can be applied to more complicated system. Work along this line is still in progress.

Acknowledgments: We thank Jian-Xin Zhu for many stimulating discussions which lead to this work. R. L. thanks Xiao-Bing Wang, M. Sindel, L. Borda, and J. von Delft for useful discussions.


format: text

text:
# A method for solving systems of non-linear differential equations with moving singularities

## 1 Introduction

When solving a physical problem, one usually encounters a set of coupled nonlinear differential equations, called dynamical equations or equations of motion, describing the dynamics of the system. This set of differential equations usually emanate from a general physical principle, and might include some subsidiary equations which can be categorized in two distinct classes, constraints or integrals of motion. The degrees of these subsidiary equations are usually at least one order lower than the dynamical equations. At the practical level these two classes of subsidiary conditions are treated very differently. The first class has to be solved simultaneously with the rest of the dynamical equations. However, the integrals of motion need only be used to put some constraints on the initial conditions. Moreover they can be used as a consistency check on the solutions.

For physically interesting cases, usually time evolution problems appear as initial valued ones, and static problems as boundary valued ones. While solving these problems, one often encounters various kinds of singularities. These singularities are usually indications of some profound physical laws or processes with significant implications. The most common type of singularities are those in which the coefficient functions of the differential equations have singularities at some fixed points, for example $`\delta `$-functions. These can be called the fixed type of singularities. In a previous work  we discussed several general methods to make an efficient numerical algorithm for boundary valued problems of this type.

In this paper we present a method for handling a class of initial valued, coupled, non-linear differential equations, whose solutions contain moving singularities. These singularities have the property that their positions and severities are apriori unknown and depend on the solutions yet to be obtained. Singularities of this type have been noted in such diverse fields as celestial mechanics, in particular the classical Kepler problem , and in the study of tensor fields defined on moving surfaces . We show that this type of singularity can be adequately treated numerically by establishing certain moving jump conditions across them. We establish the accuracy of our numerical solutions by showing that the equation representing the integral of motion is satisfied at all values of the independent variable including at the positions of the moving singularities.

For the integration algorithm we use the basic fourth-order Runge-Kutta method. It is worth mentioning that more accurate integration algorithms exist. For example, there are exponential and Bessel fitted variable step method of order 6 due to Raptis and Cash , and also a variable step P-stable method of order 6 and 8 with a phase lag of the same order due to Simos . For the embedded Runge-Kutta, formulae of order 6(5) and 8(7) have been developed by Prince and Dormand  and formulae of order 8(6) and 12(10) have been developed by Dormand et. al. . These integration algorithms should be more efficient for higher accuracy. Our choice of the integration algorithm is based on the following reasons. First, our main objective has been to find a solution to the problem of moving singularities and not the efficiency of the integration algorithm itself. Second, it turns out that in this problem the cumulative error of even the fourth order Runge-Kutta integration algorithm is negligible compared to the error introduced to the solutions at each jump across the singularities.

The set of equations that we discuss results from a classical model of gravitation in Robertson-Walker cosmology in which the signature of the metric undergoes a transition from a Euclidean to a Lorentzian domain. In section 2 we briefly discuss the physical origin of the problem and show its reduction to a set of ordinary differential equations. There, one sees an example in which this set automatically includes a subsidiary equation which is an integral of motion, along with the dynamical equations. It is worth mentioning that if an integral of motion is not directly included in the set of equations of motion, it can some times be derived directly from the dynamical equations. Also we employ a reparameterization transformation which allows one to seek continuous solutions across the hypersurface of signature change. Moreover, we employ a set of transformations which reduces the degree of severity of the moving singularities. The reader who is interested only in the numerical methods can skip to section 3, without loss of continuity.

## 2 Derivation of dynamical equations

Traditionally, one of the features of classical gravity is that the signature of the metric is usually considered as fixed. If one relaxes this condition, one may find solutions to the field equations which exhibit a signature transition . In the model that we study here a real scalar field is taken as the matter source interacting with gravity and itself in a Robertson-Walker geometry whose signature evolution is controlled by a preferred coordinate. In this model, we seek solutions to the dynamical equations which are smooth and continuous across the hypersurface of signature transition, where the metric is degenerate. The alternative would have been to find solutions by solving Einstein’s equations in disjoint regions next to the hypersurface, and then finding jump conditions to match them . For the spatially flat universes, the first approach yields exactly solvable Einstein’s equations . Here, we discuss the general case which includes the spatially flat as well as non-flat cases and solve the resulting dynamical equations numerically. For more details of the physical basis and significance of the problem, we refer the interested reader to the reference .

Consider gravity coupling to a scalar field through Einstein’s equation,

$`G_{\mu \nu }=\kappa T_{\mu \nu }[\varphi ],`$ (1)
where the scalar field $`\varphi `$ is a solution of the Klein-Gordon equation,

$`\mathrm{\Delta }\varphi {\displaystyle \frac{U}{\varphi }}=0.`$ (2)
Here, $`G_{\mu \nu }`$ is the Einstein tensor constructed from torsion-free connections compatible with the metric, and $`U(\varphi )`$ is the scalar potential for the real scalar field $`\varphi `$, which interacts with itself and gravity through the stress-energy tensor $`T[\varphi ]`$.

The above coupled equations are to be solved in a domain that would lead to Robertson-Walker cosmologies with Lorentzian signature. However, if the metric is suitably parametrized, one expects to see continuous transition to a Euclidean domain. As in , we adopt a chart with coordinate functions $`\{\beta ,x^1,x^2,x^3\}`$ where the hypersurface of signature change would be located at $`\beta =0`$. The metric can be parametrized to take the form

$`g=\beta d\beta d\beta +{\displaystyle \frac{R^2(\beta )}{[1+(k/4)r^2]^2}}{\displaystyle \underset{i}{}}dx^idx^i,`$ (3)
where $`r^2=_ix^ix^i`$. We seek solutions of the form $`R=R(\beta )`$ and $`\varphi =\varphi (\beta )`$. Now, it is apparent that the sign of $`\beta `$ determines the geometry, being Lorentzian if $`\beta >0`$ and Euclidean if $`\beta <0`$. For $`\beta >0`$, the traditional cosmic time can be recovered by the substitution $`t=(2/3)\beta ^{3/2}`$. Adopting the chart $`\{t,x^i\}`$ and using equations (1) through (3) with units in which $`\kappa =1`$, one finds

$`3{\displaystyle \frac{\dot{R}^2}{R^2}}3{\displaystyle \frac{k}{R^2}}`$ $`+`$ $`{\displaystyle \frac{\dot{\varphi }^2}{2}}+U(\varphi )=0,`$ (4)
$`2{\displaystyle \frac{\ddot{R}}{R}}+{\displaystyle \frac{\dot{R}^2}{R^2}}+{\displaystyle \frac{k}{R^2}}`$ $`=`$ $`{\displaystyle \frac{\dot{\varphi }^2}{2}}+U(\varphi ),`$ (5)
$`\ddot{\varphi }+3{\displaystyle \frac{\dot{R}}{R}}\dot{\varphi }+{\displaystyle \frac{U}{\varphi }}`$ $`=`$ $`0,`$ (6)
where a dot represents differentiation with respect to $`t`$ and $`\mathrm{}<\varphi <\mathrm{}`$, $`0R<\mathrm{}`$. Now a solution to the problem is furnished by finding $`R(t)`$ and $`\varphi (t)`$, for a given $`U(\varphi )`$. Note that these equations are not all independent. For example equation (5) can be obtained by combining equations (4) and (6). Upon a closer inspection we recognize that this is due to the fact that equation (4) is not a dynamical equation, rather it is actually an integral of motion representing a zero energy condition. That is, any solution of the dynamical equations (5) and (6) would yield a constant total energy (equation (4)). However, Einstein’s equations demand zero energy solutions only.

As is apparent from the dynamical equations, we have moving singularities at all times for which $`R=0`$. These moving singularities are potentially very severe and, as we shall see later, $`\varphi `$ actually diverges there. We can get an indication on the divergence of $`\varphi `$ from equations (4,5). These equations indicate that $`U(\varphi )`$ has to cancel the divergence of the $`k/R^2`$ terms, and for all physically relevant potentials, this implies that $`\varphi `$ has to diverge. We therefore need to use a set of transformations to reduce the severity of the divergence of solutions. We expect the following transformations to render the solutions more manageable, since it is formed of products of factors which go to zero and infinity at about the same strength,

$`X=R^{3/2}\mathrm{cosh}(\alpha \varphi ),`$ (7)
$`Y=R^{3/2}\mathrm{sinh}(\alpha \varphi ),`$ (8)
where $`\alpha ^2=\frac{3}{8}`$.

The above equations are considerably simplified if we take the potential to be

$`2\alpha ^2(X^2Y^2)U(\varphi (X,Y))=a_1X^2+a_2Y^2+2bXY,`$ (9)
where $`a_1`$, $`a_2`$ and $`b`$ are adjustable parameters. This choice for the potential stems from the fact that the left hand side of equation (9) directly appears in the Lagrangian, from which the dynamical eqations can be derived. The features of this potential and the physics involved in the choice of its parameters have been discussed in .

The dynamical equations (5) and (6) in terms of $`X`$ and $`Y`$ and the evolution variable $`\beta `$ now become,

$`Y^{\prime \prime }={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{\beta }}\right)Y^{}{\displaystyle \frac{3}{4}}\beta kY(X^2Y^2)^{2/3}\beta (a_2Y+bX),`$ (10)
$`X^{\prime \prime }={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{\beta }}\right)X^{}{\displaystyle \frac{3}{4}}\beta kX(X^2Y^2)^{2/3}+\beta (a_1X+bY),`$ (11)
subject to the subsidary ‘zero energy condition’, equation (4), which can be written in terms of the new variables as

$`\left({\displaystyle \frac{1}{\beta }}\right)(X^{^{\mathrm{\hspace{0.17em}2}}}+Y^{^{\mathrm{\hspace{0.17em}2}}}){\displaystyle \frac{9}{4}}k(X^2Y^2)^{1/3}+(a_1X^2+a_2Y^2+2bXY)=0.`$ (12)
Here, a prime represents differentiation with respect to $`\beta `$. Note that this equation does not contain any singularity, and equations (10) and (11) are actually less singular than equations (5) and (6). The coupled equations (10) and (11) must now be solved and, as explained before, equation (12) is merely a restriction on the initial conditions. However, it can also be used as a consistency check on the analytical or numerical solutions. These equations do not seem to have a closed form solution so a numerical treatment is necessary.

## 3 The numerical method

The dynamical equations that we have to solve are equations (10, 11). As mentioned in the last section, equation (12) is an integral of motion. That is, any true solution to equations (10, 11) automatically satisfies equation (12) for all values of the independent parameter, if it satisfies it at any one point. Therefore, we can use equation (12) to put a restriction on the initial conditions. More importantly, one can check the accuracy of the solutions by seeing how well equation (12) is satisfied, as the algorithm integrates the dynamical equations.

As a first step towards a numerical solution to the equations, we should study the restrictions imposed by the set of differential equations on the initial conditions. These restrictions are the result of the requirement of consistency of the initial conditions with the dynamical equations. However we can accomplish a more complete task by finding the general form of the analytic solutions close to the initial point. These solutions certainly include the complete information on the allowed set of the initial conditions111The question of the allowed set of the initial conditions, though interesting enough in its own right in all problems of this type, is of crucial importance for the problem at hand, as the determination of the correct initial conditions is an open problem in cosmology.. Moreover, the knowledge on the analytic solutions help with the first few steps of the integration algorithm.

### 3.1 Analytic solutions close to the initial point

In order to find analytic solutions which are valid near the initial point ($`\beta =0`$), we first study the restrictions imposed by equations (10)–(12) on the initial conditions. This is done by noting that in order to have well behaved solutions close to $`\beta =0`$, the first term of equation (12) shows that we must either have $`X^{}(\beta )\beta ^{n_x}`$ and $`Y^{}(\beta )\beta ^{n_y}`$, where $`n_x`$, $`n_y1/2`$, or $`|X^{}(0)|=|Y^{}(0)|`$. However, the first terms on the right hand sides of equations (10) and (11) impose a more severe restriction. These two equations admit solutions $`X^{}(\beta )\beta ^{1/2}`$ and $`Y^{}(\beta )\beta ^{1/2}`$ close to $`\beta =0`$, however, this class of solutions does not admit real or $`C^2`$ solutions across $`\beta =0`$. One can show that regular solutions close to $`\beta =0`$ are of the form

$`X(\beta )`$ $`=`$ $`A_x\beta ^3+X_0,\text{where}A_x={\displaystyle \frac{2}{9}}\left[{\displaystyle \frac{3}{4}}{\displaystyle \frac{kX_0}{(X_0^2Y_0^2)^{2/3}}}+a_1X_0+bY_0\right],`$ (13)
$`Y(\beta )`$ $`=`$ $`A_y\beta ^3+Y_0,\text{where}A_y={\displaystyle \frac{2}{9}}\left[{\displaystyle \frac{3}{4}}{\displaystyle \frac{kY_0}{(X_0^2Y_0^2)^{2/3}}}a_2Y_0bX_0\right],`$ (14)
with $`X_0X(0)`$, etc. Therefore, the initial conditions on the first and second derivatives must satisfy the relations

$`X^{}(0)=Y^{}(0)=0\text{and}X^{\prime \prime }(0)=Y^{\prime \prime }(0)=0.`$ (15)
Strictly speaking the conditions on the second derivatives are not initial conditions but rather consistency checks, since we have coupled second order equations. Therefore, the initial values for the functions $`X`$ and $`Y`$ must now satisfy, c.f. equation (12),

$`{\displaystyle \frac{9}{4}}k(X_0^2Y_0^2)^{1/3}+(a_1X_0^2+a_2Y_0^2+2bX_0Y_0)=0.`$ (16)
The contour plots of equation (16) for $`k=\pm 1`$ are given in figure 1. Along the contours, one finds the possible initial values for $`X`$ and $`Y`$. Although equation (16) is equivalent to a sixth order algebraic equation which cannot be directly solved analytically, we can solve it by going back to the original variables $`R`$ and $`\varphi `$. The solutions are either $`R(0)=0`$ giving $`\varphi (0)=\pm \mathrm{}`$, which we exclude because we have been seeking continuous solutions across $`\beta =0`$, or $`R(0)0`$ (it is a free parameter) with

$`\varphi (0)={\displaystyle \frac{1}{2\alpha }}\mathrm{cosh}^1\left[{\displaystyle \frac{DB\pm b\sqrt{D^2B^2+b^2}}{B^2b^2}}\right],`$ (17)
where

$`D={\displaystyle \frac{9k}{4R(0)^2}}{\displaystyle \frac{a_1a_2}{2}}\text{and}B={\displaystyle \frac{a_1+a_2}{2}}={\displaystyle \frac{m^2}{2}}.`$
Therefore, the acceptable values of $`X_0`$ and $`Y_0`$ ( $`X_0>|Y_0|`$) can also be obtained analytically from equation (17).

### 3.2 Integration algorithm

The important feature of equations (10) and (11) is that they are singular for all $`\beta `$ at which $`X=\pm Y`$. At these critical values of $`\beta `$ ($`\beta _c`$), the original variables take the values $`R(\beta _c)=0`$ and $`\varphi (\beta _c)=\pm \mathrm{}`$. We can directly infer from the differential equations that the solutions for the new variables and their first derivatives have to be continuous across the singularities. However, the second and higher derivatives will be singular at $`\beta _c`$. That is, the singularities of the new variables are considerably milder than those of the original variables.

Although the solutions and their first derivatives are continuous across $`\beta _c`$, they cause problems for the integration algorithm. Any attempt in solving these equations involves handling these moving singularities, as one encounters them when integrating the coupled equations. To proceed, we first establish jump conditions across these singular points as follows: close to $`\beta _c`$ we assume that the solutions have the following linear forms

$`X_\pm `$ $`=`$ $`a_\pm +b_\pm \beta ,`$ (18)
$`Y_\pm `$ $`=`$ $`c_\pm +d_\pm \beta ,`$ (19)
where $`\pm `$ refers to the right or left hand sides of the singularity, respectively. Substituting the above equations in (10) and (11) and dropping all non-singular terms, one can integrate these equations in the interval $`\beta _cϵ,\beta _c+ϵ`$, where $`2ϵ`$ is the distance across the jump. For the integration we have dropped all terms which would give rise to contributions $`𝒪(ϵ^{4/3})`$. One finds at $`Y_c=\pm X_c`$

$`b_+b_{}={\displaystyle \frac{9}{4}}k{\displaystyle \frac{(2X_cϵ)^{1/3}}{(b_{}d_{})^{2/3}}}\beta _c=(d_+d_{}),`$ (20)
where $`ϵ`$ can be taken as small a value as is desired for any required accuracy. Equation (20), together with the requirement of continuity of $`X`$ and $`Y`$, establish our jump condition for handling the singularities of the differential equations. It is apparent from equation (20) that the slopes $`X^{}(\beta )`$ and $`Y^{}(\beta )`$ are continuous at $`\beta _c`$.

Writing an actual algorithm for handling these singularities requires some care. Let us first write the original variable $`R(\beta )`$ in terms of the new variables,

$$R=(X^2Y^2)^{1/3}.$$

We recall that the differential equations become singular when $`R=0`$. As we integrate them, when $`|R|`$ becomes small (less than 1) we reduce the step size by one order of magnitude since the crossing of $`R(\beta )`$ through zero at $`\beta _c`$ is rather steep. Then at the first instant when $`|R|`$ becomes smaller than 0.1, henceforth called ‘the fixed point’, the algorithm records all the relevant values ($`X,Y,X^{},Y^{},R,\beta `$) and continues integrating towards the singular point with yet finer steps. Past the fixed point, the singular terms in the differential equations become too large and no integration algorithm can give reliable values for $`X`$ and $`Y`$. However, we can use the information obtained past this point to pinpoint $`\beta _c`$ as follows: We record the last two values of $`\beta `$ and $`R`$ right before the instant when the sign of $`R`$ changes. Then assuming $`R(\beta \beta _c)^{1/3}`$ (which is consistent with equations (18) and (19)) one obtains

$$\beta _c=\frac{R_2^3\beta _1R_1^3\beta _2}{R_2^3R_1^3}.$$
(21)
Having obtained $`\beta _c`$, we can calculate $`ϵ=\beta _{\text{fixed}}\beta _c`$ and use a linear extrapolation to obtain $`X_c`$ and $`Y_c`$ and see whether $`Y_c=\pm X_c`$ as a consistency check. We can then calculate the values of the slopes ($`X^{},Y^{}`$) on the other side of the singularity (at $`\beta _{\text{over}}=\beta _c+ϵ`$) using the jump conditions (equation (20)) and then using linear extrapolation on both sides of the singularity, the values of the functions can be calculated at $`\beta _{\text{over}}`$. The algorithm then continues integrating with fine steps until the value of $`|R|`$ increases beyond 1 and with regular steps until it approaches the next singularity.

We use a set of parameters ($`b=2`$, $`\lambda =0`$, $`m^2=4.5`$) in equations (10-12) which are physically relevant and choose our initial conditions consistent with equations (15) and (16). Recall that since equation (12) is a constant of motion, if it is satisfied at $`\beta =0`$, for a true solution it will be satisfied at all other values of $`\beta `$. Therefore if the initial values for the functions $`X`$ and $`Y`$ satisfy equation (16) at $`\beta =0`$, equation (12) should always be satisfied. For integrating equations (10) and (11), we have used the fourth order Runge-Kutta method. The resulting solutions for $`k=\pm 1`$ are shown in figure 2. It is apparent from figure 2 that at the singular points $`X=\pm Y`$, the solutions are continuous and the singularities are very mild. As a measure of the accuracy of the solutions, we have computed the ‘zero energy condition’, equation (12), as a function of $`\beta `$ for $`k=\pm 1`$ which are shown in figure 3. It is evident from figure 3 that the values of ‘total energy’ stay very close to zero, thus indicating the validity of the numerical solution. In figure 4 the variations of the original variables $`\varphi `$ and $`R`$ are shown. As can be seen, $`\varphi `$ actually diverges at the singular points. As a further check, we have numerically recovered the analytic solutions presented in  for $`k=0`$ in every detail.

## 4 Conclusions

We have shown how a particular class of initial valued coupled non-linear ordinary differential equations with moving singularities can be numerically solved. The main obstacle of having moving singularities can be overcome by establishing a set of jump conditions across them. These conditions are obtained by approximating the form of the solutions close to the singular points and directly integrating the differential equations in the neighbourhood of these points. We have found that a first order approximation close to these points is sufficiently accurate. Also since the main source of error in the solutions eminates from the jump conditions, we have found that the fourth order Runge-Kutta is sufficient for the integration algorithm.


text:
# 1 Coincidence Problem vs Equation of State

## 1 Coincidence Problem vs Equation of State

One of the most puzzling aspects of the current cosmological paradigm, is the so-called ‘cosmic coincidence problem’, i.e. the quasi equivalence between dark matter (DM) and dark energy (DE) densities today. Indeed, for a non-interacting fluid $`A`$, the continuity equation $`d(\rho _Aa^3)=p_Ada^3`$, relates the equation of state $`w_Ap_A/\rho _A`$ to the evolution of the energy density with the expansion, i.e., $`\rho _Aa^{3(w_A+1)}`$. Therefore, the ratio between dark matter ($`w_{DM}=0`$) and a dark-energy component with negative pressure ($`w_{DE}<0`$) scales as

$$\frac{\rho _{DM}}{\rho _{DE}}a^{3w_{DE}},$$
(1)
and could be only ‘coincidentally’ of order one today. This is the case for the cosmological constant ($`w_{DE}=1`$), for a scalar field with an inverse power law potential $`V(\phi )\phi ^n`$ ($`w_{DE}=2/(n+2)`$), or with any potential allowing solutions with negative pressure. On the other hand, exponential potentials exhibit attractor solutions in which there is no cosmic coincidence problem, since the scalar field energy density scales at a constant ratio with the background, be it matter or radiation. But, by virtue of eq. (1) this implies $`w_{DE}=0`$ after equivalence, thus making it unacceptable as a DE component. Thus, there seems to be an inextricable contradiction between a negative equation of state and the elimination of the cosmic coincidence.

The impasse may be escaped by relaxing the assumption leading to eq. (1), i.e. by allowing an interaction between DM and DE, besides the gravitational one.

## 2 Exponential Vamps

An interesting possibility is to consider dark matter particles whose masses depend on the quintessence field, or VArying Mass ParticleS (Vamps), which were considered in . Assuming exponential dependencies both in the potential and in the DM particle mass,

$$V(\phi )=V_0\mathrm{exp}(\beta \phi ),M=M_0\mathrm{exp}(\lambda \phi ),$$
(2)
an attractor solution exists on which both $`\rho _{\mathrm{V}amps}`$ and $`\rho _{DE}`$ scale as $`a^{3(W+1)}`$, with

$$W=\frac{\lambda }{\lambda +\beta },\mathrm{a}nd\mathrm{\Omega }_{DE}=1\mathrm{\Omega }_{\mathrm{V}amps}=\frac{3+\lambda (\lambda +\beta )}{(\beta +\lambda )^2},$$
(3)
provided $`\beta >\frac{\lambda +\sqrt{\lambda ^2+12}}{2}`$, thus solving the coincidence problem with a negative effective equation of state.

It should be noted that the crucial requirement in order to get an attractor with the above properties, is to have exponential dependencies both in the potential and in the mass, a case not considered in refs. . Moreover, the coupling $`\lambda `$ between Vamps and $`\phi `$ should be $`O(1)`$ in order to have a sensibly negative equation of state. If also baryons were coupled so strongly with the almost massless $`\phi `$, we would immediatly run into problem with the bounds coming from tests of the fifth force or the equivalence principle. In the remaining part of my talk, I will discuss how a model exhibiting these properties  can be formulated in scenarios with large extra dimensions .

Before doing that, let me mention that a similar model has been discussed recently in . Those authors implement the DE-DM interaction by assuming non-vanishing divergences for the energy momentum tensors

$$(T_{DM}^{\mu \nu })_{;\mu }=(T_{DE}^{\mu \nu })_{;\mu }=\lambda T_{DM}g^{\nu \rho }\phi ,_\rho ,(\mathrm{A}mendolaetal)$$
(4)
where $`T_{DM}`$ is the trace of the $`DM`$ energy-momentum tensor. The background behavior turns out to be the same as (3), but fluctuations evolve differently, since in the Vamps case we have

$$(T_{DM}^{\mu \nu })_{;\mu }=(T_{DE}^{\mu \nu })_{;\mu }=\lambda T_{DM}u^\nu u^\rho \phi ,_\rho (\mathrm{f}orVamps),$$
(5)
where $`u^\mu `$ is the four-velocity of Vamps particles. This will be studied in a forthcoming paper.

## 3 Vamps and Large Extra Dimensions

In principle, brane world scenarios offer a suggestive contact with the DE problem. Indeed, extra spatial dimensions compactified to a size as large as $`100\mu \mathrm{m}1\mathrm{m}\mathrm{m}`$ have been shown to be a viable possibility, provided that no Standard Model (SM) field propagates through them . It might then be tempting to attribute the observed value of the DE energy density to the Casimir energy associated to some field propagating in these large extra dimensions of size $`r`$, i.e.

$$\rho _{\mathrm{C}asimir}1/r^4.$$
(6)
If the radius $`r`$ is stabilized, then the energy density behaves exactly as a cosmological constant, and we have made no substantial progress with respect to the cosmic coincidence problem. It seems then necessary to make the radius $`r`$ dynamical, so that it may evolve on a cosmological time-scale. But this turns out to be quite dangerous. The ‘radion’ field, whose expectation value fixes $`r`$, is ultra-light ($`mH_010^{33}\mathrm{e}V`$) and couples to the trace of the SM energy momentum tensor with gravitational strength. It then behaves as a practically massless Brans-Dicke scalar with $`O(1)`$ couplings to matter, whereas present bounds are $`O(10^3)`$ .

However, besides large extra-dimensions in the $`100\mu \mathrm{m}`$ range there could also be smaller ones, and this may solve the problem. Considering for simplicity only two compact subspaces, each caracterized by a single radius, we note that the trace of the four-dimensional energy-momentum tensor couples with the total volume of the compact space, that is with the combination $`r_s^{n_s}r_l^{n_l}`$, where $`r_{s,l}`$ and $`n_{s,l}`$ are the radii and dimensionalities of the two subvolumes, respectively. If we assume a stabilization mechanism that fixes the total volume $`𝒱`$ of the compactified manifold but not its shape, the potentially dangerous combination of radion fields associated with volume fluctuations is made harmless, whereas the orthogonal one, associated with shape deformation, is decoupled from normal matter and may then be ultra-light. As a consequence, $`r_l`$ can grow on a cosmological time-scale (such that the associated Casimir energy $`1/r_l^4`$ decreases) and at the same time $`r_s`$ shrinks so that $`𝒱`$ keeps a fixed value.

If we now identify as DM candidate some stable state living on the $`4+n_s`$-dimensional brane, we may obtain a realization of the general Vamps model of eq. (2). Indeed, if the mass of the DM particle is inversely proportional to some power of $`r_s`$, as for a Kaluza-Klein (KK) state ($`m_{DM}=m_{KK}r_s^1`$), then its energy density scales as

$$\rho _{DM}1/(r_sa^3).$$
(7)
Assuming that the total volume of the compact subspace is fixed, the Casimir energy (6) and the DM energy density (7) give

$$V(\varphi )\mathrm{exp}\left(4\sqrt{\frac{n_s}{n_l𝒩}}\varphi \right),\mathrm{a}nd\rho _{DM}\mathrm{exp}\left(\sqrt{\frac{n_l}{n_s𝒩}}\varphi _2\right)a^3,$$
(8)
respectively, where $`\varphi `$ is the canonically normalized field associated to shape deformations of the compact sub-manifold . Moreover, ordinary matter such as baryons is only sensitive to the global volume of the compact manifold, not to its shape, and is thus decoupled from $`\varphi `$.

In conclusion, we have realized a Vamps model in which the DM-DE coupling is naturally $`O(1)`$, and the DE-baryon coupling is naturally suppressed. The parameters $`\lambda `$ and $`\beta `$ of eq. (2) are fixed by geometry, i.e. by $`n_l`$ and $`n_s`$, and then so are the values of $`W`$ and $`\mathrm{\Omega }_{DE}`$ on the attractor.

In Tab. 1 we list the possible values of $`\mathrm{\Omega }_{DE}`$, $`W`$ and of $`H_0t_0`$, $`t_0`$ being the present age of the Universe. We limited the dimensionality of the compact space according to the theoretical prejudice coming from string theory, i.e. $`n_l+n_s6`$.

| Table 1. Possible attractors for $`n_l+n_s6`$ | | | | |
| --- | --- | --- | --- | --- |
| $`n_s`$ | $`n_l`$ | $`\mathrm{\Omega }_{DE}`$ | $`W`$ | $`H_0t_0`$ |
| 1 - 5 | 1 | $`0.44`$ |  |  |
| 1 | 2 | 0.83 | $`1/3`$ | 1 |
| 2 | 2 | 0.68 | -0.20 | 0.83 |
| 3 | 2 | 0.60 | -0.14 | 0.78 |
| 4 | 2 | 0.56 | -0.11 | 0.75 |
| 1-3 | 3 | $`0.92`$ |  |  |
| 1-2 | 4 | no attractor |  |  |
| 1 | 5 | no attractor |  |  |

A noticeable fact about Tab. 1 is that the observed range for the dark energy density, 0.6

<

ΩDE

<

0.8

<

0.6subscriptΩ𝐷𝐸

<

0.80.6\mathrel{\vbox{\hbox{$<$}\nointerlineskip\hbox{$\sim$}}}\Omega_{DE}\mathrel{\vbox{\hbox{$<$}\nointerlineskip\hbox{$\sim$}}}0.8 uniquely selects the number of ‘large’ extra dimensions to be $`n_l=2`$, the same value that is required by the totally unrelated issue of solving the hierarchy problem with ‘millimeter’ size extra dimensions .

A model-independent prediction of the framework outlined in this paper is the presence of extra dimensions in the $`100\mu \mathrm{m}`$ range. Whereas for the solution of the hierarchy problem discussed in  this value was allowed, the link with the cosmological expansion discussed in this paper makes it a true prediction. Even more interestingly, the expected value for $`r_l`$ is quite close to the present bound of $`200\mu \mathrm{m}`$, obtained from measurements of the Newton’s law at small distances , and then will be quite soon tested experimentally.


text:
# Upper Bounds for Ropelength as a function of Crossing Number

## 1. Introduction

The ropelength of a space curve is defined to be the quotient of its length by its thickness, where thickness is the radius of the largest embedded tubular neighborhood around the curve. For a knot or link type $`L`$, we define the ropelength $`\mathrm{Rop}(L)`$ to be the minimum ropelength of all curves with the given link type. This minimum ropelength is a link invariant which measures the topological complexity of the link, much like crossing number, or bridge number, in classical knot theory.

It has been shown that every link type contains at least one $`C^{1,1}`$ tight representative which achieves this minimum ropelength . Much effort has been invested in the project of finding lower bounds for the ropelength of various link types in terms of classical topological invariants, such as the crossing number .

In this paper, we are interested in a converse problem: given a link type $`L`$ of crossing number $`\mathrm{c}(L)`$, can we guarantee the existence of a representative curve with ropelength less than some function of $`\mathrm{c}(L)`$? That is, can we find upper bounds on ropelength in terms of crossing number? Our main theorem states the following:

###### Theorem 1.

If $`L`$ is a non-split link, then

(1)
$$\mathrm{Rop}(L)1.64\mathrm{c}(L)^2+7.69\mathrm{c}(L)+6.74.$$

In particular, this bound holds for prime links.

Our Theorem 2 gives similar bounds for composite links.

Other groups () have attacked this problem by finding upper bounds on the number of edges required to embed a given link $`L`$ in the unit lattice (the lattice number $`k(L)`$ of the link), and then observing that $`\mathrm{Rop}(L)<2k(L)`$ . Both proofs rely on laying out a diagram of the knot as a graph in a planar grid and then adding bridges to form overcrossings. In this context, it has been observed that constructing a particular diagram of a link with crossing number $`\mathrm{c}(L)`$ may require ropelength $`O(\mathrm{c}(L)^2)`$ . These authors have obtained the weaker bounds $`\mathrm{Rop}(L)<24\mathrm{c}(L)^2`$ , and $`\mathrm{Rop}(L)<25\mathrm{c}(L)^2`$ . Johnston’s algorithm, like ours, produces an explicit realization of the knot in space, while the approach of  is less constructive.

By contrast, our methods are more three-dimensional and are not based on grid or lattice embeddings. Instead of using a planar diagram of a knot, we base our construction on Peter Cromwell’s idea of arc-presentations . It is curious that our methods, too, seem to be essentially of order $`\mathrm{c}(L)^2`$. While we believe that bounds with a slower order of growth must be attainable, it is becoming clear that the problem of constructing such bounds is likely to be challenging.

## 2. The Definition of Ropelength

The ropelength of a curve is defined to be the quotient of length by the radius of the largest embedded tubular neighborhood around the curve. This radius is called the thickness of the curve. For $`C^2`$ curves, this radius is locally controlled by curvature and globally controlled by distances of self-approach between various regions of the curve. Formally, we write

###### Definition 1.

The thickness of a $`C^2`$ curve $`c`$ is given by

(2)
$$\tau [c]:=\mathrm{min}\{\underset{s}{\mathrm{min}}\frac{1}{\kappa (s)},\frac{\mathrm{dcsd}(c)}{2}\},$$

where $`\kappa (s)`$ is the curvature of $`c`$ at $`s`$, and $`\mathrm{dcsd}(c)`$ is the shortest doubly-critical self-distance of $`c`$; that is, the length of the shortest chord of $`c`$ which is perpendicular to the tangent vector $`c^{}`$ at both endpoints.

We can extend this defintion to $`C^{1,1}`$ curves by adjusting our idea of the radius of curvature as follows (c.f. ):

###### Definition 2.

Let $`s`$ be a point on a $`C^{1,1}`$ curve. Consider a decreasing sequence of open neighborhoods $`U_n`$ of $`s`$. The infimal radius of curvature at $`s`$ is given by

(3)
$$\underset{U_n}{inf}\left\{\underset{tU_n}{inf}\frac{1}{\kappa (t)}\right\},$$

where the inner infimum is restricted to $`t`$ in $`U_n`$ such that $`\kappa (t)`$ exists.

Figure 1 shows examples of curves where thickness is controlled by curvature and by the doubly-critical self-distance.

Gonzalez and Maddocks have given another definition of thickness which looks somewhat less natural, but is often more useful. (See  for details). Another useful way to look at thickness comes from Federer’s notion of reach, which agrees with the thickness for curves .

###### Definition 3.

The reach of a set $`S`$ inside $`^n`$ is the greatest non-negative $`r`$ so that each point within distance $`r`$ of $`S`$ has a unique nearest neighbor in $`S`$.

## 3. Arc-presentations

We start with a definition:

###### Definition 4.

An arc-presentation of a link $`L`$ is an embedding of $`L`$ in a finite collection of $`\alpha `$ open half-planes arrayed around a common axis, or binding, so that the intersection of $`L`$ with each half-plane is a single simple arc. The number of half-planes $`\alpha `$ is called the arc-index of the arc-presentation. The minimal arc-index over all arc-presentations of a link $`L`$ is an invariant of the link type.

By isotopy, we can arrange that $`L`$ intersects the axis only at the points $`1,\mathrm{},\alpha `$. We call these the levels of the arc-presentation. Such an arc-presentation is then specified by combinatorial data: a collection of $`\alpha `$ triples in the form $`(x_i,y_i,\theta _i)`$, where each denotes an arc from level $`x_i`$ to level $`y_i`$ on the half-plane at angle $`\theta _i`$ around the axis.

Figure 2 shows an arc-presentation for the trefoil and the corresponding set of triples.

We will assemble our ropelength bounds from two ingredients. First, we define a notion of the total distance travelled by the arcs in an arc-presentation:

###### Definition 5.

The total skip of an arc-presentation $`A`$, denoted $`\mathrm{Skip}(A)`$, is

(4)
$$\mathrm{Skip}(A)=\underset{i=1}{\overset{\alpha }{}}|x_iy_i|.$$

For a given arc-presentation we can construct a realization of the knot in space with ropelength bounded in terms of $`\mathrm{Skip}(A)`$ and $`\alpha `$:

###### Proposition 1.

An arc-presentation $`A`$ composed of $`\alpha `$ half-planes can be realized with ropelength smaller than

(5)
$$\frac{2\alpha }{\mathrm{tan}(\pi /\alpha )}+(\pi 2)\alpha +2\mathrm{Skip}(A).$$

For the arc-presentation of the trefoil in Figure 2, we have $`\alpha =5`$ and $`\mathrm{Skip}(A)=12`$; so Proposition 5 yields an upper bound on the ropelength of the trefoil of about $`43.47`$. Numerical experiments estimate the ropelength of the tight trefoil to be about $`32.66`$ , so the slack in our estimate is about $`33\%`$ of the total value. Figure 3 shows the tubular neighborhoods of this trefoil knot and an arc-presentation of the knot $`7_1`$ as realized by the algorithm of Proposition 1.

Further, if we can bound $`\mathrm{Skip}(A)`$ in general, we will be able to draw conclusions about the ropelength of an arbitrary link. A combinatorial argument yields:

###### Proposition 2.

If an arc-presentation $`A`$ has arc-index $`\alpha `$, then

(6)
$$\mathrm{Skip}(A)\{\begin{array}{cc}\frac{\alpha ^21}{2}\hfill & \text{if }\alpha \text{ is odd,}\hfill \\ \frac{\alpha ^2}{2}\hfill & \text{if }\alpha \text{ is even.}\hfill \end{array}$$

This bound is sharp.

It is shown in  that any non-split link $`L`$ admits an arc-presentation with $`\alpha \mathrm{c}(L)+2`$. This result, when coupled with the previous two propositions, gives Theorem 1. We obtain an even stronger statement for composite links:

###### Theorem 2.

If $`L`$ is a non-split composite link with prime components $`L_1,L_2,\mathrm{},L_n`$, then

(7)
$$\mathrm{Rop}(L)1.64\underset{i=1}{\overset{n}{}}\mathrm{c}(L_i)^2+7.69\underset{i=1}{\overset{n}{}}\mathrm{c}(L_i)+6.74n.$$

## 4. Proofs of the key propositions and theorems

###### Proof of Proposition 1.

We would like to take an arc-presentation $`A`$ for $`L`$ as a template for constructing an embedding of $`L`$ with unit thickness. We will then bound the length of this embedding in terms of the arc-index and the total skip of $`A`$.

We begin by constructing a right regular polygonal prism $`P\times [0,2\alpha ]`$, where $`P`$ is a regular polygon with $`\alpha `$ sides of length $`2`$. This prism will serve as the binding of $`A`$; each vertical face of the prism will correspond to an open half-plane in the arc-presentation $`A`$. We divide the prism vertically into $`\alpha `$ floors, each a prism of height $`2`$, which will represent the $`\alpha `$ levels of the arc-presentation $`A`$.

We can now construct a link isotopic to $`L`$. First, represent the arcs of $`A`$ by $`\alpha `$ handles outside the prism which join different floors on the same vertical face. We will refer to these handles as fins. Next, add $`\alpha `$ circular sections inside the prism which join different vertical faces on the same floor. These sections represent the junctions between arcs on the binding of the open book described by $`A`$.

We must show that this construction can be accomplished with a unit thickness curve and then compute the length of that curve.

### 4.1. The Fins

Let us denote the fins $`_1,\mathrm{},_\alpha `$. Each fin consists of two quarter-circles of unit radius, joined by a straight vertical segment. Each fin joins two points on a vertical face of the prism and is contained in a $`2\alpha \times 2\times 2`$ rectangular box extending radially from a vertical face of the prism.

Since the $`_i`$’s stay outside the prism and each is contained in a different box, the tubes around the fins are pairwise disjoint, and disjoint from the tubes surrounding regions of the curve inside the prism. Given that each fin has curvature bounded above by $`1`$ and no doubly-critical chords, this means that the fins can be constructed with a unit-thickness curve.

###### Claim 1.

If $`\mathrm{Rop}(_i)`$ denotes the length of the segment of the curve $`_i`$, then

(8)
$$\underset{i=1}{\overset{\alpha }{}}\mathrm{Rop}(_i)=(\pi 2)\alpha +2\mathrm{Skip}(A).$$

###### Proof.

Suppose that $`_i`$ travels from floor $`x_i`$ to floor $`y_i`$ of the prism. The total vertical distance covered by the fin is $`2|x_iy_i|`$ (recall that each floor has height $`2`$). However, the quarter-circles on each end of the fin cover a vertical distance of $`2`$ units. Thus, the straight segment has length $`2|x_iy_i|2`$, and the total length of the fin is $`\pi 2+2|x_iy_i|`$. Summing over $`i=1,\mathrm{},\alpha `$ and using Definition 5 proves the claim. ∎

### 4.2. The Binding Prism

We denote the sections of the curve inside each floor of the binding prism by $`_1,\mathrm{},_\alpha `$. Each $`_i`$ is a circular arc joining the midpoints of two edges of the regular polygon which is the cross-section of the prism as shown in Figure 4.

Because the sides of the polygon have length $`2`$, each of these is an arc of a circle of radius at least one; so each arc has curvature bounded above by one. Further, since each floor has height $`2`$ and only one $`_i`$ lies in each floor, the tubes around each of the $`_i`$ are disjoint. Thus these $`_i`$ can be constructed with a tube of unit thickness.

###### Claim 2.

If $`\mathrm{Rop}(_i)`$ denotes the length of the segment of the curve $`_i`$, then

(9)
$$\underset{i=1}{\overset{\alpha }{}}\mathrm{Rop}(_i)\frac{2\alpha }{\mathrm{tan}(\pi /\alpha )}.$$

###### Proof.

Each of these circular arcs is contained in a sector of the circle inscribed within the polygonal cross-section of the prism as shown in Figure 5. Since each arc is convex, its length is bounded above by the diameter of the inscribed circle. This diameter is exactly $`2\mathrm{cot}(\pi /\alpha )`$. Summing over $`i=1,\mathrm{},\alpha `$ proves the claim. ∎

Combining Claims 1 and 2 yields the statement of Proposition 1. ∎

###### Proof of Proposition 2.

Our job is to find an upper bound for $`\mathrm{Skip}(A)=_{i=1}^\alpha |x_iy_i|.`$ We first observe that the difference $`|x_iy_i|`$ is one unit larger than the number of levels skipped over. For example, jumping from level $`3`$ to level $`6`$, a difference of $`3`$ levels, skips the fourth and fifth levels. Thus, we can rewrite the sum

(10)
$$\mathrm{Skip}(A)=\alpha +\underset{i=1}{\overset{\alpha }{}}\{\text{number of levels skipped by the arc }(x_i,y_i,\theta _i)\}.$$

Notice that any level $`j`$ contributes to the above sum exactly when it is skipped over. We can rewrite our sum in terms of $`j`$ as

(11) $`\mathrm{Skip}(A)`$ $`=\alpha +{\displaystyle \underset{j=1}{\overset{\alpha }{}}}\{\text{number of times level }j\text{ is skipped}\}`$
$`=\alpha +{\displaystyle \underset{j=1}{\overset{\alpha /2}{}}}\{\text{number of times level }j\text{ is skipped}\}`$
$`+{\displaystyle \underset{j=0}{\overset{\alpha \alpha /21}{}}}\{\text{number of times level }\alpha j\text{ is skipped}\},`$
where in the final equality we have split the second half of the sum off and let $`j\alpha j`$.

Now we bound the number of times level $`j`$ is skipped over. The only way to hop over $`j`$ from a higher level is to land on a lower level. There are $`j1`$ levels below the $`j`$th on which such a jump can land. Further, each of these levels can act as a launch pad for a jump back up which crosses the $`j`$th level again. This gives at most $`2(j1)`$ skips over level $`j`$. Similarly, the number of times we can skip over the $`\alpha j`$th level is twice the number of levels above it, or $`2j`$.

For even $`\alpha `$, these estimates are sharp (as we will see below). However, when level $`\alpha j`$ is the central level of an arc-presentation with $`2k+1`$ levels ($`j=k=\frac{\alpha 1}{2}`$), the situation is slightly different. Here all of the $`j`$ levels above the middle cannot be initial and terminal levels of arcs which skip level $`\alpha j`$. For if so, then no arcs land on level $`\alpha j`$, and we could have eliminated level $`\alpha j`$ from the original arc-presentation. Thus level $`\alpha j`$ is skipped at most $`2j1=\alpha 2`$ times.

Inserting these bounds into Equation 11, we apply the sum formulae for arithmatic progressions. When $`\alpha `$ is odd, we get

(12)
$$\mathrm{Skip}(A)\alpha +\underset{j=1}{\overset{\frac{\alpha 1}{2}}{}}2(j1)+\underset{j=0}{\overset{\frac{\alpha 3}{2}}{}}2j+(\alpha 2)=\frac{\alpha ^21}{2}.$$

If $`\alpha `$ is even, the proof is similar.

We now construct arc-presentations which show that these results are sharp. Consider the arc-presentation with even arc-index $`\alpha =2k`$ described by the data

$`(\alpha ,\alpha /2,\theta _1),(\alpha /2,\alpha `$ $`1,\theta _2),(\alpha 1,\alpha /21,\theta _3),(\alpha /21,\alpha 2,\theta _4),`$
$`\mathrm{},(\alpha /2+1,1,\theta _{2k1}),(1,\alpha ,\theta _{2k}).`$
If we add up the lengths of the jumps, we get

(13) $`\mathrm{Skip}(A)=\alpha ^2/2.`$
The same approach yields a realization of $`A`$ so that $`\mathrm{Skip}(A)=\frac{\alpha ^21}{2}`$ for odd $`\alpha `$.

###### Proof of Theorem 1.

Taylor’s theorem gives the approximation $`\frac{1}{\mathrm{tan}(x)}1/xx/3`$ for $`x>0`$. Via Propositions 1 and 2 we gather that

(14) $`\mathrm{Rop}(L)`$ $`{\displaystyle \frac{2\alpha }{\mathrm{tan}(\pi /\alpha )}}+(\pi 2)\alpha +\alpha ^2`$
$`(2/\pi +1)\alpha ^2+(\pi 2)\alpha 2\pi /3.`$
By Bae and Park , for any non-split link $`L`$ there exists an arc-presentation with $`\alpha \mathrm{c}(L)+2`$. Inserting this into the above bound for ropelength yields

(15)
$$\mathrm{Rop}(L)(2/\pi +1)\mathrm{c}(L)^2+(8/\pi +2+\pi )\mathrm{c}(L)+(8/\pi +4\pi /3),$$

and each of these constants evaluates to something smaller than the approximations given in the statement of the theorem. To gain the final remark in the theorem, we note that any prime link $`L`$ is non-split (otherwise it would consist of split compnents $`L_1`$ and $`L_2`$ and would admit the nontrivial factors $`L_1`$ and $`L_2`$ union a split unknot). ∎

###### Proof of Theorem 2.

The strategy for this proof is to arrange the prime components of our composite link so that we can make use of the bounds given by Theorem 1. So suppose that we have found arc-presentations with minimal arc-index for these components and embedded them as unit-thickness curves $`L_1,\mathrm{},L_n`$ according to the algorithm of Proposition 1.

We will now prove that for any links $`L_1`$ and $`L_2`$, constructed by the algorithm of Proposition 1, we can construct a curve $`L_1\mathrm{\#}L_2`$ with ropelength less than or equal to $`\mathrm{Rop}(L_1)+\mathrm{Rop}(L_2)`$. This is all that is required to complete the proof of our Theorem since the bound in the statement is just the sum of the bounds obtained for the $`L_i`$ by Theorem 1.

We begin by preparing $`L_1`$ and $`L_2`$. The top floor of $`L_1`$ contains only a single horizontal circular arc joining the centers of two sides of the binding prism. Since no fins jump over this level, we may rotate these quarter-circles to face one another and replace the horizontal circular arc with a horizontal line segment of shorter length without changing thickness or knot type. We do the same for the bottom floor of $`L_2`$. This procedure is shown in Figure 6.

We now arrange $`L_1`$ and $`L_2`$ in space so that the horizontal segments are colinear and share an endpoint. If we keep each oriented so that its floors are horizontal, the only overlap between the tubes surrounding each curve occurs on the shared floor. At the shared endpoint, we may delete two quarter-circles and replace them with a vertical line segment of length $`2`$. We could keep track of this savings and get a slightly better constant term in the statement of Theorem 2. For each prime component we add, we save $`\pi 2`$ in length.

Handling the other endpoints of the curve will prove to be a little more work. We may assume that both line segments lie along the x-axis with the shared endpoint at the origin. Suppose $`L_2`$’s segment has length $`\mathrm{}_2`$, while $`L_1`$’s segment has the smaller length $`\mathrm{}_1`$.

We now rotate the remaining vertical quarter-circle of $`L_1`$ to face the corresponding quarter-circle of $`L_2`$. If $`\mathrm{}_1\mathrm{}_22`$, we may replace both horizontal line segments with a single, shorter horizontal line segment joining the ends of these vertical quarter circles to obtain the desired curve. See Figure 7.

If $`\mathrm{}_1>\mathrm{}_22`$, we cannot simply connect the endpoints of the quarter-circles after rotating the lower quarter-circle to face right. The resulting curve would have cusps on both ends. We solve this problem by finding a line tangent to both circles and following the composite path shown in Figure 8.

It is less obvious that these changes reduce length. To see that they do, we consider the diagonal line tangent to both circles shown in Figure 8. Since both circles are also tangent to a horizontal line, by symmetry this horizontal line cuts the diagonal line in half. Consider Figure 9. We need only show that half of the diagonal line (labelled $`x`$ in the Figure) is shorter than the portion of the quarter-circle it replaces (twice the angle $`\theta `$).

Since the lower quarter-circle has unit radius, this amounts to proving that $`\mathrm{tan}\theta 2\theta `$ for $`0\theta \pi /4`$. This is shown by a simple computation.

Since the resulting curve remains $`C^{1,1}`$, is still of unit thickness, and has less length than the total length of the initial curves, this completes the proof. ∎

An example of this construction is shown in Figure 10.

## 5. Acknowledgements

We would like to thank many of our colleagues for helpful conversations on these topics over several years, including Yuanan Diao, Heather Johnston, Rob Kusner and John Sullivan. We are indebted to the VIGRE program at the University of Georgia, and to the other members of our VIGRE research group: Ted Ashton, Kenny Little, Heunggi Park, Darren Wolford, and Nancy Wrinkle. In addition, Cantarella would like to acknowledge the generous support of the National Science Foundation Postdoctoral Research Fellowship Program under grant number DMS-99-02397, and Faber and Mullikin would like to recognize the gracious support of the VIGRE fellowship program at UGA.