text:
# Analytical description of finite size effects for RNA secondary structures

## I Introduction

Folding of biopolymers is a fundamental process in molecular biology without which life as we know it would not be possible. In biopolymer folding, well-characterized interactions between individual monomers make a polymer fold into a specific structure believed to minimize the total interaction free energy. The apparent simplicity in the formulation of this biopolymer folding problem is in sharp contrast with the immense challenges faced in actually describing biopolymer folding quantitatively caused by the intricate interplay of monomer-monomer interactions and the constraint that the monomers are connected into a chain of a certain sequence. The biological importance of biopolymer folding paired with this immense intellectual challenge has sparked numerous computational and theoretical studies profold . These studies do not only attempt quantitative predictions of specific structures but also focus on more fundamental properties of the biopolymer folding problem such as its phase diagram.

While the bulk of the work concentrates on the folding of proteins due to its overwhelming importance in pharmaceutical applications, recently RNA folding has been identified as an ideal model system for biopolymer folding tinoco ; higgs00 . RNA is a biopolymer of four different bases G, C, A, and U. The most important interaction among these bases is the formation of Watson-Crick(WC) base pairs, i.e., A–U and G–C pairs. This comparatively simple interaction scheme makes the RNA folding problem very amenable to theoretical approaches without losing the overall flavor of the general biopolymer folding problem. Again, a lot of effort has been devoted to understanding fundamental properties of RNA folding such as the different thermodynamic phases an ensemble of RNA molecules can be in as a function of temperature, an external force acting on the molecules, and the sequence design muller ; higgs96 ; pagna ; hart ; bund99 ; degen .

All these theoretical approaches are concerned with the phase behavior of RNA molecules in the thermodynamic limit. In order to compare these theoretical predictions with numerical or actual biological experiments it is thus important to know which role finite size effects play, i.e., at which size of a molecule the universal predictions of the asymptotic theories are expected to hold. In this publication we precisely aim to answer this question. We study homogeneous RNA sequences which allows us to analytically solve for the universal asymptotic behavior as well as the cross-over length below which the universal theory is not applicable any more. We find that this cross-over length is very strongly dependent on the sequence of the molecule. For realistic energy parameters we find that the cross-over length can be as large as 10,000 bases. This is about the largest size of naturally occurring RNAs as well as the largest length of RNA molecules amenable to quantitative computational approaches. Thus, we conclude that finite size effects have to be seriously taken into account when describing RNA folding by asymptotic theories.

This article is organized as follows: In Sec. II, we briefly review the definition of RNA secondary structure. In Sec. III, we analytically derive the finite size effects of the simplest model of RNA folding namely a homogeneous sequence without loop entropies. While this model is mainly treated for pedagogical purposes, in Sec. IV, we sketch how the result can be generalized to more realistic models of RNA folding. In Sec. V, the behavior of the cross-over length $`N_0`$ is discussed. We find that $`N_0`$ depends mostly on the degree of branching of the RNA molecules and a simple approximate formula is derived. These results are shown to be consistent with the numerical values obtained using experimentally known energy parameters for specific sequences in Sec. VI. We point out how enormous finite size effects in the RNA secondary structure formation problem can be. The detailed derivations of the partition function and the cross-over length are relegated to various appendices.

## II Review of RNA secondary structures

### II.1 Definitions

RNA usually occurs as a single-stranded polymer with four types of monomers (bases) G, C, A, and U. The strand can bend back onto itself and form helices consisting of stacks of stable Watson-Crick pairs (A with U or G with C).

An RNA secondary structure describes which bases are bound and can be written as a set of binding pairs $`(i,j)`$, where $`i`$ and $`j`$ denote the $`i^{th}`$ and $`j^{th}`$ base of the RNA polymer respectively. For example, the secondary structure S of Fig. 1 is written as

$$S=\{(2,57),(3,56),\mathrm{},(i,j),\mathrm{},(k,l),(37,44)\}.$$

In this study, we apply the common approximation to exclude the so-called pseudo-knots tinoco , i.e., for two base pairs $`(i,j)`$ and $`(k,l)`$, the configurations $`i<k<j<l`$ and $`k<i<l<j`$ are not allowed. As a result, the analytical studies become more tractable.

This exclusion of pseudo-knots is reasonable. For long pseudo-knots, the double helix structure would require threading one end of the molecule through its loops many times so they are kinetically difficult to form. Thus, these pseudo-knots occur infrequently in natural RNA structures tinico; higgs00 . Short pseudo-knots, on the other hand, do not contribute much to the total free energy because of the relatively low binding energies for short sequences and the strong electrostatic repulsion of the backbone since the polymer backbone is highly charged. By excluding pseudo-knots, we will stay close to commonly used algorithms that compute the exact partition function which can be applied to test our model.

### II.2 Interaction energies

Since the tertiary interactions between structures are in general much weaker than the interactions among the secondary structures tinoco ; higgs00 , we will follow the common approaches and take into account only the energy contribution from the secondary structures.

If we assign a Gibbs free energy $`\mathrm{\Delta }G(𝐒)`$ to each secondary structure S, the partition function of the ensemble of all structures is given by

$$Z=\underset{𝐒}{}e^{\mathrm{\Delta }G(𝐒)}.$$
(1)
The Gibbs free energy is commonly used to describe the secondary structure since it contains entropic contributions from the formation of loops as well as enthalpic terms from the formation of base pairs. The total free energy $`G(𝐒)`$ is the sum of the energy contributed from each elementary piece such as the stacking of base pairs and the connecting loops. The largest contribution are the stacking energies between adjacent WC pairs, and these energies depend on the type of bases in the pairs. While the typical value of the stacking energy is on the order of $`k_BT`$ at room temperature, both the enthalpic and entropic terms are on the order of $`10k_BT`$. Thus, the stacking energy will become repulsive with a moderate increase of temperature to around $`80^{}C`$ and the RNA molecule denatures.

## III Molten Phase

In order to get a qualitative understanding of finite size effects, we first follow previous works higgs96 ; pagna ; bund99 ; higgs00 ; bund02  and assume that the Gibbs free energy is the sum of the binding free energies $`ϵ_{ij}`$ of each base pair in the structure,

$$\mathrm{\Delta }G(𝐒)=\underset{(i,j)𝐒}{}ϵ_{ij},$$
(2)
and neglect the entropic energies due to the formation of loops for the rest of this section.

The binding free energies $`ϵ_{ij}`$ in this model are differences between the gain in chemical binding energy and the loss in the configurational entropy associated with the formation of the base pairs. Since both contributions are large and comparable, the realistic values of the $`ϵ_{ij}`$ strongly depend on temperature. Since we do not describe spatial degrees of freedom in this model, it does not describe denaturation of the RNA molecule and we restrict ourselves from here on to a parameter regime where the majority of the bases is paired, i.e., where a significant fraction of the $`ϵ_{ij}`$ is negative.

In order to obtain analytical insights into the finite size effects, we additionally assume that the binding free energy $`ϵ_{ij}`$ is a constant $`ϵ_0`$, independent of the identities of the bases. Thus, in our simplified model $`\mathrm{\Delta }G(𝐒)=ϵ_0\times |𝐒|`$ where $`|𝐒|`$ stands for the number of pairs in S. This simplified energy model serves as the basis of our study for the more realistic energy model.

As it stands, this energy model and the more realistic energy model we will introduce later describe only homogeneous sequences. However, it has been argued that this energy model can be applied to random RNA sequences at high enough temperature when the disorder is sufficiently weak bund02 . Under this weak disorder, there exist many structures with nearly degenerate energies and the corresponding scaling laws match the predictions of the simplified energy model. Only as the temperature is lowered, a strong disorder phase arises. This low temperature phase is characterized by a small number of distinct low-energy structures, and is referred to as the “glass phase” in analogy with studies of other disorder system. However, this glass phase is not within the scope of this article.

The partition function of the molten phase model can be obtained through the recursive relation in Fig. 3. This figure shows how the possible ways of binding can be decomposed into two cases where the last base N is either unbound or bound to some base k. If we define $`G(N+1)`$ as the partition function of an RNA of length N, the relation reads

$$G(N+1)=G(N)+\underset{k=1}{\overset{N1}{}}G(k)qG(Nk),$$
(3)
where $`q=e^{\beta ϵ_0}`$. Together with the boundary condition $`G(1)=1`$, this equation allows calculation of the exact value of the partition function $`G(N)`$ in $`O(N^2)`$ time. The Vienna package hofac  is able to calculate this exact value with more complete sequence dependent energy parameters based on the similar scheme. This recursive equation (3) also leads us to the analytical expression for the partition function. By introducing the z-transform,

$$\widehat{G}(z)=\underset{N=1}{\overset{\mathrm{}}{}}G(N)z^N,$$
(4)
and applying it to Eq. (3), we get a quadratic equation for $`\widehat{G}(z)`$ as

$$z\widehat{G}(z)1=\widehat{G}(z)+q\widehat{G}^2(z),$$
(5)
from which $`\widehat{G}(z)`$ can be solved. For large sequence lengths, the partition function $`G(N)`$ is obtained by performing the inverse z-transform on $`\widehat{G}(z)`$, and can be approximated (see Appendix A) as

$`G(N)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi i}}{\displaystyle \widehat{G}(z)z^{N1}𝑑z}.`$ (6)
$``$ $`A(q)N^\theta z_c^N(q)\left[1{\displaystyle \frac{N_0(q)}{N}}+O({\displaystyle \frac{1}{N^2}})\right],`$ (7)
where $`z_c(q)=1+2\sqrt{q}`$ is the branch point of $`\widehat{G}(z)`$, $`\theta =3/2`$, $`A(q)=[(1+2\sqrt{q})/4\pi q^{3/2}]^{1/2}`$ and $`N_0(q)=3(1+4\sqrt{q})/16\sqrt{q}`$. This asymptotic analytical formula is only determined by the behavior of $`\widehat{G}(z)`$ near the branch point $`z_c`$. The exponent $`\theta =3/2`$ indicates the characteristic universal behavior of this partition function for long sequences. The non-universal cross-over length $`N_0(q)`$ characterizes how long a sequence has to be for the universal laws to hold. Here, we find an explicit analytical formula for $`N_0`$ as a function of parameter q.

From the formula of $`N_0(q)`$, we can see that the cross-over length in this simple model is on the order of 1 for all values of q. However, in the following sections we will show that the cross-over length may vary over several order of magnitudes from several bases to 1,000,000 bases. Thus, the loop entropies of the more realistic energy model would greatly modify the behavior of cross-over length.

## IV Including loop entropies

To get a more quantitative understanding of the cross-over length, we now take into account the loop entropies and introduce Boltzmann factors, $`s`$, $`b`$, $`i`$, $`h`$, and $`m`$, for the different types of loops (see Fig. 1). The values of these free energy parameters have been carefully measured freier  such that our model can be applied quantitatively to realistic RNA molecules. Typically, the free energy of a stacking loop is large and negative ($`s1`$) while the free energies for all the other loops tend to be large and positive, leading to Boltzmann factors much less than one. The binding energy $`ϵ_0`$ of the simple model introduced above is now absorbed into these loop free energies. As mentioned in Sec. III, we still restrict ourselves to a temperature regime below denaturation which is determined by the true energy model.

Again, we want to calculate the partition function of the structure ensemble and derive the cross-over length as a function of the loop parameters. This calculation in principle follows along the lines of Sec. III, but is technically much more elaborate because of the more complicated energy model. A reader more interested in the final results than in the technical details is advised to directly skip to the Sec. V.

In order to calculate the partition function, we separate the secondary structure into two categories as shown in Fig. 4. One is the bubble structure which contains only hairpins and multiloops. The other is the stem structure, which connects the bubbles, containing only stacking loops, bulges, and interior loops. We will study each of them individually and later combine them together.

### IV.1 Stem structure

In principle, the loop free energy depends on the length of a loop. Thus, unbound bases also contribute to the total free energy. This contribution has been experimentally measured for small loops and behaves logarithmically with length when the loop is large. However, in the following we show that the free base energy of unbound bases provides only a negligible effect on the behavior of a stem and thus on the cross-over length if the stacking energy is relatively large.

To explore all possible ways of pair bindings, again a graphical recursion relation is helpful. Such a recursion relation is shown in Fig. 5. Starting from a closed pair on the left, the following loop can be either a stacking loop, a bulge or an interior loop which correspond to the terms on the right hand side. To study the influence of a free energy for unbound bases, we assign the Boltzmann factors, $`\stackrel{~}{b}`$ and $`\stackrel{~}{i}`$, to each unbound base in a bulge and an interior loop respectively. If we define the partition function of stem structures with N bases and the first and the last bases of which is paired as $`S(N1)`$, which corresponds to the left hand term, this relation in Fig. 5 is formulated as

$`S(N1)`$ $`=`$ $`sS(N3)+2{\displaystyle \underset{k=3}{\overset{N2}{}}}b\stackrel{~}{b}^{Nk1}S(k2)`$ (8)
$`+`$ $`{\displaystyle \underset{a=3}{\overset{N3}{}}}{\displaystyle \underset{b=a+1}{\overset{N2}{}}}i\stackrel{~}{i}^{N(ba)3}S(ba).`$
To perform the z-transform, we have to consider the initial conditions, $`S(1)=1`$, $`S(2)=0`$, $`S(3)=s`$, and $`S(4)=2b\stackrel{~}{b}`$. These initial conditions arise because certain structures are not allowed, e.g., since a base can not be shared in two base pairings, a stem with 3 bases does not exist and this leads to $`S(2)=0`$. Also when the length of a stem is small, certain loops which require many bases are not allowed, e.g., the only available structure for a stem with 4 bases is a stacking loop. Including these conditions, we apply the z-transform on Eq. (8) with the definition

$$\widehat{S}(z)=\underset{N=1}{\overset{\mathrm{}}{}}S(N)z^N,$$
(9)
to resolve the convolution. Solving for $`\widehat{S}(z)`$ gives us

$$\widehat{S}(z)=\frac{1}{z}\left[1\frac{s}{z^2}\frac{2b\stackrel{~}{b}}{z^2(z\stackrel{~}{b})}\frac{i\stackrel{~}{i}^2}{z^2(z\stackrel{~}{i})^2}\right]^1.$$
(10)
Again the partition function $`S(N)`$ can be obtained by applying the inverse z-transform on $`\widehat{S}(z)`$.

To illustrate the effect from the unbound bases, we show how their Boltzmann factors affect the average behavior of the long stem structure. The average quantity of a certain type of loop or unbound base can be calculated as $`\tau _\tau \mathrm{ln}S(N)`$ where $`\tau `$ is the corresponding Boltzmann factor. Here, we specifically calculate the average number of unbound bases per interior loop, which is defined by $`\stackrel{~}{i}_{\stackrel{~}{i}}\mathrm{ln}S(N)/i_i\mathrm{ln}S(N)`$ in the large N limit, as a function of $`\stackrel{~}{i}`$ (see Appendix B for detailed calculations).

From Fig. 6, we can see that this average number is barely affected by $`\stackrel{~}{i}`$ when the Boltzmann factor $`s`$ for a stacking loops is large. This can be easily understood as follows: For an interior loop, the unbound bases always introduce a strong energy penalty since less bases are available for stacking. This penalty is much larger than the penalty due to loss of the degree of freedom by one more free base. Thus when the binding energy is large, the free base entropic penalty can be neglected. From Fig. 6, we can see that interior loops tend to stay at the smallest length (2 unbound bases) when $`s`$ is large, independent of the value of $`\stackrel{~}{i}`$. Since the same argument applies to bulges as well, we will set $`\stackrel{~}{i}`$ and $`\stackrel{~}{b}`$ to 1 for the rest of this publication.

### IV.2 Bubble structure

In a similar fashion, a recursive relation for the bubble structure is found graphically as shown in Fig. 7. In the first relation for a closed bubble structure, we can have either a hairpin loop or a multiloop following from the closed pair at the end. In the second relation, the multiloop structure can be decomposed into two cases where the last base is either unbound or bound. Since a multiloop has to have at least 3 branches, we have a term with two more bubble structures; the last recursive term produces more branches.

By defining the partition function for the closed bubble structure and multiloop structure with N bases as $`B(N1)`$ and $`M(N+1)`$ respectively, the recursive relations read

$`B(N1)`$ $`=`$ $`t\left(h+mM(N1)\right)`$ (11)
$`M(N+1)`$ $`=`$ $`M(N)+{\displaystyle \underset{k=1}{\overset{N1}{}}}M(k)B(Nk)`$
$`+`$ $`{\displaystyle \underset{a=1}{\overset{N3}{}}}{\displaystyle \underset{b=a+1}{\overset{N2}{}}}{\displaystyle \underset{k=b+1}{\overset{N1}{}}}B(ba)B(Nk).`$
Here an additional Boltzmann factor $`t`$ is introduced in the first relation at the position where two bubbles are connected. Later we will insert stems into the bubble structure by replacing $`t`$ with the partition function of the stem structure. We also neglect the free base energy for unbound bases in hairpins and multiloops for similar arguments as above.

In this recursive relation, the smallest multioop should have at least 4 bases such that two branchings can be connected. Thus, we set the initial conditions as $`M(1)=M(2)=M(3)=M(4)=0`$, which forbids a multiloop with length less than 4. With the initial conditions, the recursive relations result in the following quadratic equation for the z-transformed $`\widehat{B}(z)`$ :

$$\left(\frac{1}{t}+\frac{m}{(z1)^2}\right)\widehat{B}^2\left(\frac{z1}{t}+\frac{h}{z1}\right)\widehat{B}+h=0.$$
(13)
### IV.3 Complete structure

To combine the stem and bubble structures, we insert a stem structure at each position represented by $`t`$ which is a placeholder for the connections between multiloops and hairpin loops in the bubble structure. In this case the first relation in Fig. (7) is modified as indicated in Fig. 8.

By defining the partition function of the closed structure on the left hand side as $`C(N1)`$, the Fig. 8 reads

$$C(N1)=\underset{k=1}{\overset{N1}{}}S(k)B(Nk).$$
(14)
After z-transform, this relation results in $`\widehat{C}=(z\widehat{S})\widehat{B}`$. Thus, the following replacement

$$tz\widehat{S},\widehat{B}\widehat{C}$$
(15)
combines stem and bubble structures together. In order to complete all possible structures, the single strands outside the closed end pair also have to be included. This can be done by going back to the first recursion relation, Eq. (3), and replacing $`G(Nk)`$ by the closed structures $`C(Nk)`$, which relates $`\widehat{G}(z)`$ to $`\widehat{C}(z)`$ as

$$\widehat{G}=\frac{1}{(z1)\widehat{C}}.$$
(16)
Putting everything together, we obtain

$`\widehat{G}`$ $`=`$ $`{\displaystyle \frac{1}{2m}}[({\displaystyle \frac{z1}{z\widehat{S}}}+{\displaystyle \frac{2mh}{z1}})`$
$``$ $`\sqrt{\left({\displaystyle \frac{z1}{z\widehat{S}}}+{\displaystyle \frac{h}{z1}}\right)^2\left({\displaystyle \frac{4h}{z\widehat{S}}}+{\displaystyle \frac{4hm}{(z1)^2}}\right)}].`$
Notice that the leading singularity in $`\widehat{G}`$ is again from the branch cut induced by the square root. Thus, as expected, the inverse z-transform leads to the same universal behavior (7) with an exponent $`\theta =3/2`$ as the simple model we studied first. However, non-universal quantities such as the cross-over length $`N_0`$ will depend on the parameters of the extended model.

### IV.4 Minimum hairpin length constraint

For natural RNA molecules, a hairpin loop needs to have at least 3 unbound bases due to the width of the double helix (which implies $`ji>3`$ in the secondary structure). This minimum hairpin length constraint is easily taken into account in our calculations. Under the constraint, a bubble structure which contains at least one hairpin must have at least 5 bases. Thus, we adopt the initial condition, $`B(1)=B(2)=B(3)=0`$, when we perform the z-transform on Eq. (11). The summation range in Eq. (11) is then changed and it simply leads to a substitution of $`h`$ by $`h/z^3`$ in all subsequent equations.

This substitution is reasonable since $`z`$ represents the Boltzmann factor for the free energy of one single base. The minimum hairpin length constraint reduces the number of available bases for binding by 3, so it introduces an energy penalty $`3\mathrm{ln}z`$ which causes the Boltzmann factor $`h`$ to be divided by $`z^3`$. In this way, we can easily introduce any kind of minimum length constraint via a similar substitution. For example, if we require a bulge loop to have at least two unbound bases instead of one, the replacement of $`b`$ by $`b/z`$ will include this constraint.

Note, this principle also helps us to understand equation (10) for the stem structure. The terms with different powers of z, namely $`s/z^2`$, $`b/z^3`$, and $`i/z^4`$, arise because a stacking loop reduces the number of available bases for pairing by 2, a bulge loop has at least one unbound base, and an interior loop has a minimum of 2 unbound bases, respectively.

## V Behavior of the cross-over length

In this section, we will use the general results of Sec. IV in order to calculate the cross-over length in our model for sequences with loop entropies.

### V.1 Large stacking energy approximation

To derive the cross-over length, we need to solve for the branch point $`z_c`$ that is defined by the vanishing of the term under the square root in Eq. (17). In principle, we can always obtain the numerical value of the branch point and expand $`\widehat{G}(z)`$ around this point to obtain a numerical value for the cross-over length (see Appendix A). We will refer to this numerical value as the homogeneous cross-over length since it is obtained under the homogeneous molten phase model.

However, since this calculation involves finding the root of a fourth order polynomial, no meaningful analytical expression can be found in general. Thus, we resort to a large $`s`$ approximation in order to obtain an analytical expression. This approximation is justified since the Boltzmann factor of a stacking loop, $`s`$, is usually much larger than 1 while the loop Boltzmann factors $`b,i,h,`$ and $`m`$ are less than one. In this approximation, from Eq. (17) we find the branch point $`z_c\sqrt{s}`$, i.e., the free energy per base is $`f=kT\mathrm{ln}(z_c)\frac{1}{2}\mathrm{\Delta }G(s)`$. This can be easily interpreted since we expect most bases to form pairs due to the favorable stacking loops such that the free energy per base is half of the free energy of a stacking loop.

### V.2 Cross-over length

Including the minimum hairpin length constraint introduced in Sec. IV.D, we expand the branch point $`z_c`$ near $`\sqrt{s}`$. Then, the approximated analytical formula for the cross-over length becomes

$`N_0`$ $`=`$ $`{\displaystyle \frac{3s^{3/4}}{8\sqrt{hm}}}[(\sqrt{s}1)^2+h`$ (18)
$`+`$ $`{\displaystyle \frac{9b}{2\sqrt{s}}}+{\displaystyle \frac{11i6b}{4s}}+O(s^{3/2})].`$
It has a straight forward interpretation: The simplest possible structure is a long stem with one hairpin at the end. Every additional branching of the structure requires formation of one hairpin and one multiloop. Since upon formation of the hairpin and mulitloop, at least 3 bases become unbound, the Boltzmann factor for a branching is $`hm/s^{3/2}`$. The prefactor in Eq. (18) is up to the numerical factor of $`3/8`$ the inverse square root of this expression. Thus, we conclude that the cross-over length which becomes larger as the Boltzmann factor for a single branching becomes smaller can be interpreted as the minimum length that allows a certain degree of branchings.

Since the Boltzmann factors $`b`$ and $`i`$ appear only in the higher order terms, they barely modify the cross-over length and can be neglected altogether. This is consistent with the fact that $`b`$ and $`i`$ only play roles in the stem structure, but not in the bubble structure. The leading term of the approximated analytical formula (18) will be referred to as the large s cross-over length.

### V.3 Reliability of the large s approximation

For the analytical large s cross-over length to be useful, we have to know how good they agree with the numerical value of the homogeneous cross-over length. Here, we compare these two values for many different choices of energy parameters covering the whole range of realistic values. Fig. 9 shows how the large s cross-over length approaches the homogeneous one as $`s`$ gets large. Typical values for the Boltzmann factor of a stacking loop $`s`$ involving GC pairs are $`s30`$ freier , so the approximation is very good in this region. For stacking loops involving AU pairs, $`s`$ is around 5 so a deviation from the approximated formula in the large $`s`$ limit can be seen. However, since $`N_0`$ only sets the order of magnitude of the length beyond which the asymptotic theory is applicable, the large s cross-over length with a deviation by a factor of 2 at $`s=5`$ is still a good estimation.

## VI Numerical verification

While for a generic RNA sequence, the molten phase model is believed to only apply at sufficiently high temperature, it can be applied to repeated sequences at all temperatures below denaturation since each repeated unit can be viewed as the equivalent of a base in the molten phase (Fig. 10). To illustrate the correctness of the calculations shown in Secs. IV and V, and to get a feeling for typical cross-over lengths, we now compare our large s cross-over length of repeated sequences with the full numerical results. The full numerical result is obtained by using the Vienna package hofac  which can calculate the exact value of the partition function and other observables for arbitrary sequences using a realistic energy model.

As an observable, we choose the average size $`l`$ of a structure. This quantity is defined as

$$l=\underset{k=1}{\overset{N/2}{}}\underset{k^{}=N/2+1}{\overset{N}{}}P_{k,k^{}},$$
(19)
where $`P_{k,k^{}}`$ is the probability that bases $`k`$ and $`k^{}`$ are paired. The latter probability is also calculated exactly by the Vienna package. This size measures the average number of base pairs to be crossed when connecting the $`N/2^{th}`$ base to base 1 (Fig. 10), which captures the size of the secondary structure. We expect $`l`$ to obey

$$lN^{1/2}\left(1(\frac{N_0}{N})^{1/2}+O(\frac{1}{N})\right),$$
(20)
where the leading term is the asymptotic behavior bund02  and the next term reflects the first expected correction which is a constant independent of N. We determine $`l`$ for sequences of different lengths and extract the full numerical cross-over length $`N_0`$ by fitting data obtained via the Vienna package to Eq. (20). This is then compared to our large s cross-over length.

### VI.1 $`(GCA)_n`$ sequence

We apply this scheme to a repeated $`(GCA)_n`$ sequence. Such sequences naturally occur in the gene for Huntington’s disease and their secondary structures are believed to play a role in this disease. Since the GC/CG stack in the secondary structure of the $`(GCA)_n`$ sequence is much more favorable than any other combination, we can exclude the possibility to have binding pairs other than GC/CG. Thus, GCA is viewed as one unit base in the molten phase. With this equivalence we can use the experimentally determined parameters mathews  and calculate the equivalent energy parameters $`s`$, $`b`$, $`i`$, $`h`$, and $`m`$ for the molten phase model. For example, the stacking energy of the molten phase is the sum of GC/CG stacking energy and the free energy for the interior loop of length 2 in the $`(GCA)_n`$ sequence.

Fig. 11 shows the full numerical results for the average size of the secondary structure for a $`(GCA)_n`$ sequence as a function of the sequence length $`N`$. By fitting the result to the Eq. (22), we get a full numerical value of the cross-over length of 6.9 bases.

To compare this full numerical value with our large s cross-over length, we plug the equivalent energy parameters of the molten phase model into the approximated formula $`3s^{1/4}(\sqrt{s}1)^2/8\sqrt{hm}`$. This formula is different from equation (18) because the minimum hairpin length of the corresponding molten phase model only 1 instead of 3. The resulting large s cross-over length is about 2.3 repeat units which corresponds to 6.9 bases in the $`(GCA)_n`$ sequence. This agrees very well with the full numerical value.

### VI.2 $`(AU)_n`$ sequences

Repeated AU sequences have already been suggested as models for the molten phase by de Gennes in 1968 degen . For such sequences, we exclude the possibility of AA or UU binding pairs since they are not favorable at all. In a similar fashion as for the $`(GCA)_n`$ sequence in Fig.10, the smallest bulge has two free bases. However, since the minimum hairpin length is 4 instead of 3, in order to match AU or UA closing pairs, the large s cross-over length is $`3s(\sqrt{s}1)^2/8\sqrt{hm}`$.

Plugging in the correct values for the parameters at body temperature results in a large s cross-over length of $`N_07700`$. A verification of this value is beyond the reach of the numerical procedure using the Vienna package. Since the cross-over length is expected to decrease as the denaturation transition is reached, the full numerical verification could be performed at a higher temperature. Thus, we first study the the temperature dependence of the cross-over length $`N_0`$. To this end we choose to study the homogeneous cross-over length instead of large s cross-over length because it is not clear whether the large s approximation is appropriate when approaching the denaturation temperature. When studying the temperature dependence, it is important to note that the free energies themselves have a very strong temperature dependence which changes Boltzmann factors much more than the explicit $`\beta `$ in the exponent. By looking up the table of enthalpies corresponding to different loops, we determine this temperature dependence of the Boltzmann factors $`s`$, $`b`$, $`i`$, $`h`$, and $`m`$. The numerical values of the homogeneous cross-over lengths thus obtained are shown in Fig. 12. From this figure, we choose to numerically verify the estimate of the cross-over length at $`57^{}C`$.

The full numerical verification is done in a similar way as Sec. VI.A. Fitting numerical data at $`57^{}C`$ gives a full cross-over length $`N_0900`$. Our homogeneous cross-over length predicts a value $`N_01500`$ which is of the same order of magnitude as the full numerical result. From the above numerical studies, we conclude that our theory for the large s cross-over length is a good estimate. Thus, the cross-over length of $`(AU)_n`$ at body temperature is truly around 7,700 bases, which is almost the largest size of naturally occurred RNA sequences. For $`(GC)_n`$ sequences at body temperature, the situation is even more dramatic with a predicted cross-over length about 105,000 bases. This is beyond the limit of most natural RNAs and suggests that the structure ensemble of $`(AU)_n`$ and $`(GC)_n`$ molecules can never be described by the asymptotic theory for any naturally available molecules.

### VI.3 Distribution of cross-over length for short repeated sequences

Lastly, we want to explore the sequence dependence of the cross-over length for short repeated sequences. We study all possible partially self-complementary repeated sequences of unit length two and three and evaluate their cross-over lengths using the formula for the large s cross-over length.

For non self-complementary repeated sequences, e.g., the $`(GCC)_n`$ sequence, it is not a priori clear that we can map the repeated unit GCC onto one single base in the molten phase (Fig. 13) since base G can pair with either one of the two bases C. However, we verified numerically that more than 99% of the G’s are paired in the stacked configuration indicated by the dashed lines in Fig. 13. This is due to the large loss in free energy as one stacking loop is lost and a bulge loop is created. We verified this behavior for all the partially self-complementary sequences in our discussion. Thus, we can safely use the same approach as used for the $`(GCA)_n`$ sequence and calculate the large s cross-over length. The results are shown in the following table.

We can see that the cross-over lengths spread over the whole range from only several bases to $`10^6`$ bases. This non universal behavior has to be treated individually for different sequences. For longer repeated units, we expect similar behavior but we did not verify this due to the complicated entropic effects.

## VII Conclusion

For our asymptotic description to be appropriate, an RNA sequence must be long enough to allow a certain degree of branching. The cross-over length at which this happens for homogeneous sequences is strongly sequence dependent. Specifically for cases when stacking is very favorable and hairpins and multiloops are unfavorable, the cross-over length $`N_0`$ can be very large. In this region, finite size effects are very important and have to be taken into account in any numerical simulation or interpretation of experimental data in terms of asymptotic formulas. Our analytical results provide an easy way to estimate this cross-over length for realistic energy models.

It remains as an interesting and relevant question how large the cross-over lengths for generic, non-homogeneous sequences in the glassy phase are. Unfortunately, this question can at this stage neither be addressed numerically nor analytically since even the expected asymptotic behavior is elusive. Nevertheless, the existence of large cross-over length in the molten phase as established here suggests that also in the glass phase large cross-over length have to be considered as a serious possibility.

## VIII Appendix

### VIII.1 z-transform

For a series of numbers {G(1), G(2), …, G(N), …}, the z-transform of G(N) is defined as

$$\widehat{G}(z)=\underset{N=1}{\overset{\mathrm{}}{}}G(N)z^N.$$

To recover $`G(N)`$, we apply the inverse z-transform:

$$G(N)=\frac{1}{2\pi i}_CG(z)z^{N1}𝑑z.$$

As an example, here we solve for the partition function for the molten phase model. From Eq. (5), $`\widehat{G}(z)`$ is calculated as

$$\widehat{G}(z)=\frac{1}{2q}(z1\sqrt{(z1)^24q}),$$

so the partition function $`G(N)`$ can be obtained by plugging $`\widehat{G}(z)`$ into the inverse z-transform.

The analytical part $`(z1)`$ can be dropped since it does not contribute to the integral. The square root introduces a branch cut with two branch points at the ends. By taking the contour C as the smallest loop around the branch cut, the integral becomes

$`G(N)`$ $``$ $`{\displaystyle \frac{1}{4\pi qi}}{\displaystyle _{z_0}^{z_c}}(\sqrt{(z1)^24q+iϵ}`$
$`\sqrt{(z1)^24qiϵ})z^{N1}dz`$
$`=`$ $`{\displaystyle \frac{1}{2\pi q}}{\displaystyle _{z_0}^{z_c}}\sqrt{4q(z1)^2}z^{N1}𝑑z,`$
where $`z_c=1+2\sqrt{q}`$ is the branch point with greatest real part and $`z_0`$ is the other branch point. In the large N limit, due to the exponent factor $`z^N`$, we expect that only the behavior near $`z_c`$ is important, so we can expand the integrand around $`z_c`$ and perform the approximation.

To do this integral, we first consider a more general case where $`\widehat{G}(z)=\sqrt{f(z)}`$. The behavior near the greatest value of $`z`$ is obtained correctly if we replace $`z`$ by $`e^\mu `$ and expand $`\sqrt{f(\mu )}`$ in the power terms of $`(\mu _c\mu )`$ as

$`\sqrt{f(\mu )}`$ $`=`$ $`\sqrt{f^{}(\mu _c)}(\mu _c\mu )^{1/2}`$
$`\left[1{\displaystyle \frac{f^{\prime \prime }(\mu _c)}{4f^{}(\mu _c)}}(\mu _c\mu )+O((\mu _c\mu )^2)\right].`$
Together with the following approximation valid in the large N limit,

$$_{\mu _0}^{\mu _c}(\mu _c\mu )^\alpha e^{\mu N}𝑑\mu \mathrm{\Gamma }(1+\alpha )N^{(1+\alpha )}e^{\mu _cN},$$

the inverse z-transform of $`\widehat{G}(z)`$ can be expressed as

$`G(N)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _{\mu _0}^{\mu _c}}\sqrt{f(\mu )}e^{\mu N}𝑑\mu `$
$`=`$ $`{\displaystyle \frac{\sqrt{f^{}(\mu _c)}}{\pi }}\mathrm{\Gamma }({\displaystyle \frac{3}{2}})N^{3/2}e^{\mu _cN}`$
$`\left[1{\displaystyle \frac{f^{\prime \prime }(\mu _c)}{4f^{}(\mu _c)}}{\displaystyle \frac{\mathrm{\Gamma }(5/2)}{\mathrm{\Gamma }(3/2)}}{\displaystyle \frac{1}{N}}+O({\displaystyle \frac{1}{N^2}})\right]`$
$`=`$ $`A_0N^{3/2}z_c^N\left[1{\displaystyle \frac{N_0}{N}}+O({\displaystyle \frac{1}{N^2}})\right].`$
Here, the partition function of the molten phase model (7) is obtained by setting $`f(\mu )=4q(e^\mu 1)^2`$. For the general case, we can put the result in terms of $`z`$ with the substitutions

$`{\displaystyle \frac{df(\mu )}{d\mu }}|_{\mu _c}`$ $`=`$ $`z_c{\displaystyle \frac{df(z)}{dz}}|_{z_c}`$
$`{\displaystyle \frac{d^2f(\mu )}{d\mu ^2}}|_{\mu _c}`$ $`=`$ $`z_c{\displaystyle \frac{df(z)}{dz}}|_{z_c}+z_c^2{\displaystyle \frac{d^2f(z)}{dz^2}}|_{z_c}.`$
After little algebra, the cross-over length is obtained as

$$N_0=\frac{3}{8}\left[1+\frac{z_c\frac{d^2f(z)}{dz^2}|_{z_c}}{\frac{df(z)}{dz}|_{z_c}}\right]$$

for arbitrary $`f(z)`$.

### VIII.2 average length of interior loops

The calculation of partition function for the stem structure $`S(N)`$ is different from that for the molten phase model because there is no branch cut, so we do not expect the universal behavior $`N^{3/2}`$. However, the fact that only the pole with the greatest real part is important in the large N limit is preserved. From Eq. (8), the z-transform of the partition function for the stem structure $`\widehat{S}(z)`$ has the form $`f(z)/g(z)`$ where

$`g(z)`$ $`=`$ $`z^2(z\stackrel{~}{b})(z\stackrel{~}{i})^2s(z\stackrel{~}{b})(z\stackrel{~}{i})^2`$
$`2b\stackrel{~}{b}(z\stackrel{~}{i})^2i\stackrel{~}{i}^2(z\stackrel{~}{b}).`$
In the large N limit, the inverse z-transform can be approximated by

$`S(N)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi i}}{\displaystyle \widehat{S}(z)z^{N1}𝑑z}`$
$``$ $`{\displaystyle \frac{1}{2\pi i}}{\displaystyle \frac{f(z)}{g^{}(z_c)(zz_c)}z^{N1}𝑑z}`$
$``$ $`{\displaystyle \frac{f(z_c)}{z_cg^{}(z_c)}}z_c^N,`$
where $`z_c=z_c(s,b,\stackrel{~}{b},i,\stackrel{~}{i})`$ is the pole with greatest real part.

Since $`z_c^N`$ is the dominating term for large N, the average quantity can be approximated as

$`\tau `$ $`=`$ $`\tau _\tau \mathrm{ln}(S(N))`$
$``$ $`N\tau _\tau \mathrm{ln}(z_c)=Nz_c^1\tau _\tau z_c.`$
So the average number of unbound bases per interior loop can be obtained by $`(\stackrel{~}{i}_{\stackrel{~}{i}}z_c)/(i_iz_c)`$. Here, the pole $`z_c`$ can only be solved numerically.


format: text

text:
# Loss cone: past, present and future

## I Introduction

The characteristic frequency to which LISA  is sensitive is comparable to the orbital frequency at the innermost stable orbit of Schwarzschild black holes of mass $`M_{BH}10^6\text{M}_{}`$ and associated Schwarzschild radius, $`r_S`$. One of the more promising “guaranteed” sources for LISA is the gravitational radiation from the final stages of coalescence of low mass ($`1`$$`1000\text{M}_{}`$) compact objects with low mass massive black holes ($`M_{BH}10^{6\pm 1}\text{M}_{}`$). White dwarfs, with masses of $`0.5`$-$`1.3\text{M}_{}`$; neutron stars, with masses of $`1.4\text{M}_{}`$; and, stellar mass black holes, with masses from order solar masses upwards of $`100\text{M}_{}`$, may all coalesce with central supermassive black holes. The stellar population on the most tightly bound, low angular momentum orbits in the cusp around the central supermassive black hole is depleted by prompt swallowing of the stars “initially” on such orbits. The net flux of stars into the black hole is then set by the rate at which stars (stellar remnants) can enter the “loss cone”  on orbits such that gravitational radiation will shrink the semi-major axis more rapidly than other processes can increase the energy of the object on that orbit. Stellar dynamical processes, and possibly star formation processes, determine the instantaneous and mean rate at which (compact) stars can enter the gravitational radiation loss cone and evolve dynamically to coalescence with the central black hole.

Previous papers  considered the likely rate for detectable signals from degenerate compact objects in the cusps of normal galaxies, coalescing with central black holes such as the one inferred to be present in the Milky Way. Estimates for detectable signal rates are of the order one per year and higher, from these sources, but the estimates are sensitive to systematic uncertainties in the population contributing to the signal, with very large (many orders of magnitude) cumulative formal uncertainties in the expected signal rate. Since high rates of coalescence are self–limiting (through depletion of low mass compact objects at high coalescence rates; or the growth of the primary through coalescences, to the point where the frequency of the innermost stable orbit is too low for LISA to be sensitive to further coalescences), the formal uncertainty in the detectable rate is skewed to lower rates of coalescence than the “canonical” estimates in the literature. The true rate of detectable coalescences is sensitive to the assumptions about the mass function of the primary black hole, the number density of stars in regions where the appropriate mass black holes are found, the mass and number of compact remnants available for coalescence over a Hubble time in such regions, and the large scale structure of the stellar system in which the central supermassive black hole is embedded.

### A Assumptions

A number of similar assumptions are made by authors in estimating coalescence rates:

* Galactic cusps: Galaxies tend to have “cuspy”, non-isothermal centres. Observations and theory suggest that real galaxies have rising density profiles at small radii with:

$$\underset{r0}{lim}\rho r^\gamma $$

and $`\gamma =3/2+p`$, $`p\{1,+1\}`$ . Cosmological simulations suggest that a moderate cusp forms from initial condition cold collapse in hierarchical models , with the cusp possibly modified by subsequent black hole growth and interaction . Alternatively, one may assume that the cusp grows with the black hole from an initially dense, relaxed cluster of stars, with constant density profile at small radii.
* Relaxation: An implicit assumption is that central SMBH are near ubiquitous, with $`M_{BH}M_{GalSph}\sigma _{gal}^n`$ , where $`\sigma `$ is the dispersion of the stellar spheroid, defined at some suitable radius far from the black hole. Inside the black hole effective radius, $`r_hGM_{BH}/\sigma ^2`$, the velocity dispersion is rising like $`lim_{r0}\sigma (r)r^{1/2}`$, due simply to the Keplerian potential of the central black hole. Additional assumptions may be made about the stellar velocity anisotropy, $`\beta (r)`$, and the flattening or possible triaxiality of the stellar density profile at small and large radii .

Note that the census of low mass supermassive black holes $`M_{BH}10^510^6\text{M}_{}`$ of particular interest for LISA is not available, the existence and properties of this populations is extrapolated from the local population (Milky Way and M32), from the higher mass SMBH observed, and from circumstantial evidence from low luminosity active galactic nuclei (AGN), consistent with accretion onto lower mass SMBH.
* Stellar population: The relaxation time of the stellar population is $`t_r\sigma ^3(r)/\rho (r)`$, constant for $`p=0`$, which thus defines the critical cusp.

Most relevant galaxies are not relaxed in the centre, from observations. However, low mass galaxies are more likely to have high stellar densities at $`r_h`$ and steep cusps, and therefore may have short relaxation times inside $`r_h`$. The characteristic galaxy has luminosity $`L^{}`$. The galaxies most likely to have relaxed cusps with central supermassive black holes in the right mass range for LISA are sub-$`L^{}`$ galaxies, less luminous than the characteristic luminosity. It is possible that the stellar population was initially relaxed as densities were higher and high mass stars were present in galactic centers at early times. The current relaxation time may be longer than the initial relaxation time, and relaxation processes may have imprinted on the initial conditions.
* Mass function: A general assumption that is made is that the stellar population is normal, with a Salpeter or Scalo initial mass function. This may of course be wrong, and it has been suggested that the initial mass function (IMF) in galactic centres or nuclear star forming regions is biased to high masses (see Kroupa 2001 for review), which would substantially increase the estimated LISA rates. A further assumption is made about the white dwarf (WD) mass function, as a function of the progenitor zero-age mains sequence (ZAMS) mass (usually assumed to be independent of metallicity); the cut-off ZAMS mass for neutron star formation (also usually assumed to be independent of metallicity); and the cut-off ZAMS mass and mass function of stellar mass black holes formed (which ought not to be independent of metallicity, see Figer, this Proceedings). Binarity may in principle be different in nuclear star clusters than is observed in the field, either due to initial conditions, or due to dynamical modifications of the binary population. Since mass transfer processes in close binaries affect the formation channels for compact objects, this will affect the inferred rate of coalescence. In particular, binary interactions may reopen formation channels for classes of compact objects which are not available for single stars, reinjecting a population that might naively be assumed to be depleted. Assuming a Salpeter IMF and some upper mass truncation, the number fractions of different compact object populations are then of the order: $`f_{WD}0.10.3`$, $`f_{NS}10^3f_{kick}`$, $`f_{BH}10^4`$. $`f_{kick}`$ is the fraction of neutron stars receiving small enough a kick as not to be displaced from the centre of galaxies containing low mass SMBH, estimating from the globular cluster population of neutron stars, we might infer $`f_{kick}0.1`$.
* Compact Remnants: Compact remnants can then coalesce with SMBH through gravitational radiation or direct capture. Capture into eccentric orbits with pericentres $``$ few $`r_S`$ which can “grind down” (whereby the orbital semi-major axis shrinks gradually) through episodic emissions of gravitational radiation, before they are upscattered by a subsequent stellar encounter, provides the sought after periodic signal in the final stages of coalescence (see Sigurdsson & Rees 1997 for discussion).
* The loss cone: The loss cone \- $`\theta _{lc}(r)`$ \- is, by definition, the region where the time scale for coalescence through gravitational radiation emission is shorter than time scale to diffuse out of loss cone through random walks or large scatterings in angular momentum $`J`$ (or, rarely, energy, $`E`$).

The rate is sensitive to number density profile of stellar population in the center of galaxies, $`\rho r^{3/2p}`$, to the mass function and to any anisotropies in the stellar distribution function.

### B Past calculations

A number of authors have estimated the rate of coalescences of compact objects into central supermassive black holes, using a number of different formalisms and assumptions about the underlying “initial conditions”:

* Formalisms: The primary formalisms used follow Frank & Rees (1976), Lightman & Shapiro (1977, see also ) and Hénon (1973):

+ Hils & Bender (1995) used the Lightman & Shapiro formalism calculating the coalescence rate in M32 like galaxies (nucleated dwarf ellipticals with high density cores, short core relaxation times and an evolved stellar population). The rate calculated by Hils & Bender provides both a useful benchmark for other people’s calculations, and a good estimate for the range of event rates LISA may observe.

    The Hils & Bender rate for white dwarf-SMBH coalescence in M32 is $`1.8\times 10^8`$ events per year.
    + Sigurdsson & Rees 1997, and Sigurdsson 1997, 1998 used the Frank & Rees formalism, with the refinement of estimating also the large angle scattering rate, which may contribute significantly to the total rate in the “pinhole” regime, as opposed to the diffusion regime implicitly assumed by other formalisms. The pinhole regime assumes that change in angular momentum of individual stars is large per unit dynamical time, compared to the width of the loss cone in angular momentum space; the diffusion regime applies to the opposite limiting case of small changes in angular momentum per star per unit crossing time. S & R ’97 found a rate for M32 like galaxies consistent with Hils and Bender, with a somewhat higher net coalescence rate if large angle scattering is included. S & R found that the cosmological rate is dominated by low mass nucleated spiral galaxies, like the Milky Way, extending to lower (bulge) mass galaxies. Such galaxies are more numerous than M32 analogues , may have higher central densities, lower SMBH masses (and surprisingly higher coalescence rates) and ongoing star formation that replenishes the compact object population.

    Including large angle scattering, S & R found an event rate of $`3\times 10^8`$ per year for white dwarf-SMBH coalescence in M32. Excluding the large angle coalescence their rate agrees with Hils and Bender.
    + Miralda-Escude & Gould (2000) used essentially the Frank & Rees formalism with the assumptions of Bahcall & Wolf. They focused on low mass black hole coalescence, with the assumption of a relaxed central populations, finding a steady coalescence rate for Milky Way–like galaxies of $`10^6`$ per year. In comparison, S & R found an “initial” rate for LMBH-SMBH coalescence of $`10^4`$ per year in cuspy spirals, declining to a sustained merger rate of $`10^610^8`$ per year, depending on the assumptions about the cusp structure and total black hole population. With high initial LMBH coalescence rates, the initial black hole population is rapidly depleted (predicting a burst of coalescence at high–moderate redshift) and then settles to an asymptotic lower rate as replenishment of the inner population matches the steady state depletion rate.
    + Ivanov (2002) used the Lightman & Shapiro formalism, and correcting for the different WD mass function, gets a rate lower than but consistent with Hils & Bender for M32 like galaxies. Note that the statement in that paper that S & R get a much lower rate than Hils & Bender is wrong and is directly contradicted by the S & R paper. Some care must be made in comparing the rates since the models are scaled to different structural parameters, which are not independent. Also, differential and integrated rates are of course not directly comparable.
    + The current state–of–the–art models are by Freitag (2000,2002), using a time explicit Monte Carlo realisation of the stellar cluster surrounding the central black hole, doing a Fokker-Planck approximation evolution of the stellar dynamics, including SMBH coalescence and growth, following the Hénon formalism. The assumptions in Freitag’s work are typically for a very massive and compact isothermal cluster, with substantial structural evolution due to black hole growth and stellar evolution. As a consequence, the rates in Freitag’s models tend to be strongly time dependent, but they are consistent with those of other authors, scaling Freitag’s initial conditions to the appropriate late–time quasi–steady state values typically approximated by other work.

All the calculations assume some $`\gamma `$ (or $`p`$), scale to $`M_{BH}`$ and one of $`\sigma `$, or the number density of stars at the “radius of influence” of the black hole $`n_{}(r_h)`$. Inspection of the published papers reveals that the dominant source of the difference in the rates is the choice of “initial conditions”; principally the choice of evolved or unevolved stellar population (IMF, BH mass function), initial density and whether the population is initially relaxed or not. All assume isotropy and sphericity.
* Normal galaxies:Unrelaxed, nucleated sub-$`L^{}`$ spirals dominate the rate.
* Sustainability:Very high coalescence rates are self-limiting - they deplete the compact object population or grow the SMBH to larger mass, moving out of the LISA sensitivity range. The coalescence rate declines with higher mass SMBH, primarily due to the anti-correlation with stellar density and correlation with stellar dispersion.
* Aside - runaway growth: we can check for runaway growth under these assumptions, to see whether the growth of the loss cone with SMBH mass can trigger rapid mass growth, since adding some mass $`\delta MM_{BH}`$ causes the loss cone to expand. Runaway growth occurs if the time averaged occupancy number of additional stars in the additional slice of loss cone, $`\delta \theta (\delta M)`$, is greater than one. Trivially we find $`\delta \theta _{lc}/\delta M\theta _{lc}`$, implying exponential growth, but on inspection, with very low growth constant, so effectively there is linear growth of the SMBH due to the expansion of the loss cone, with a low growth rate.

The exception is for high $`m_{}/M_{BH}`$ where $`\delta \theta (\delta M)`$ can be large and therefore approach runaway conditions. So, in this approximation, we expect runaway growth for the lowest mass seed black hole only, which saturates when $`M_{BH}m_{}`$, or, if there is a substantial population of very massive compact objects in the central regions ($`m_{}>100\text{M}_{})`$.

### C Predictions

WD merger rates are uncertain by about a factor of three. The primary sources of uncertainty are the relative number density of WDs given the observed stellar luminosity density; and, whether the inner cusps are typically relaxed or not.

NS merger rates are bounded above to be no more than approximately few % of the WD rate, unless nuclear IMFs are very skewed to high masses, since the number fraction of NS from a normal evolved stellar population is small, and the NS are not massive enough to mass segregate efficiently in the cusp. Since they are more massive than typical WDs, they may still contribute significantly to the observed rate. However, if NS natal kicks are ubiquitous, then the NS population is displaced from the inner cusp and the number density in the inner cusp is small and the rate negligible at all times.

BH merger rates are uncertain by at least two orders of magnitude and in reality are likely highly time dependent. Critical issues is whether the cusp is relaxed, allowing an infusion of low mass BH from the outer cusp, what the number fraction of low mass BH formed is, and whether it varies strongly with the metallicity of the star forming gas, and whether star formation persists in the inner nuclear region at late times (as seems likely given observations of the Milky Way).

* LISA detection: All predictions in the literature are consistent with LISA detections. That is, despite the range of coalescence rates, even the most conservative rate estimates predict LISA will see multiple events during its planned operating lifetime.
* Source confusion: The highest rate predictions imply possible source confusion.
* The uncertainties in the rate estimates are primarily assumption driven, dominated by the initial conditions assumed, not the calculation method. In principle, the true initial conditions are constrainable by observations.
* Other processes: Some significant physics are still not included in the models. There is some reason to believe that the missing physical processes will not substantially affect the rate, and will if anything tend to lead to enhanced coalescence rates, with one possible exception. Significant tangential anisotropy of the inner stellar cusp could suppress the coalescence rate.

Other issues of concern include:

Direct capture vs spiral–in.

Gravitational radiation detection requires orbits that spiral in, gradually shrinking the semi-major axis over $`O(10^5)`$ orbits. The capture is into orbits with peribothrons of few $`r_S`$ ($`e0.9990.9999999`$), but not directly into event horizon, or rapidly perturbed into event horizon. Some fraction of those orbits are scattered down, across the event horizon, by stellar encounters near apobothron.

Estimates of direct loss fraction varies substantially ($`1/3`$ or more). There is some concern that this fraction could be very large, eliminating a large fraction of the estimated sources. Calculations suggest that the loss fraction is in fact modest .

Tidal disruption estimates. similar coalescence processes apply to main-sequence stars, but those stars generally undergo tidal disruption far outside the Schwarzschild radius (except for the most massive SMBH) and we should observe stellar tidal disruption events . The tidal disruption rate scales approximately with the compact object coalescence rate (but see Miralda-Escude & Gould 2000), typically $`r_{coll}`$ for stars is much larger than for compact stellar remnants, so the density distribution of stars “breaks” at larger radii in the cusp.

Observationally, the inferred tidal disruption rate is consistent with maximum predicted rate Donley et al (2002). Naively this would imply that the loss cone is generally full. Donley et al also find that the inferred tidal disruption rate is significantly higher for AGNs, suggesting that dynamical processes refilling the loss cone, or nuclear star formation, may substantially enhance the rate, contributing strongly to the integrated cosmological rate.

## II Future issues

A number of open issues remain to be explored in more detail.

### A Brownian wandering

The formalisms used generally assume the SMBH is fixed centrally.

In reality there can be substantial “Brownian motion” of the central SMBH, particularly for the lower mass SMBH of particular interest to LISA. In general we have some black hole velocity, $`v_{BH}`$; mean black hole displacement from the center of mass of the stellar population, $`r_{BH}`$; and some time scale on which the Brownian motion occurs, $`t_{BH}`$.

If the stellar population is isothermal and relaxed, then the problem is trivial, the SMBH is in equipartition and moving as a thermally excited massive particle in a harmonic oscillator potential with, $`v_{isoBH}=\sqrt{m_{}/M_{BH}}\sigma `$. Clearly there are two limits:

* If $`M_{BH}m_{}`$ there is no loss cone.
* If $`M_{BH}\mathrm{}`$ then loss cone is fixed.

For non-isothermal distributions, this is a hard problem. Even for flat density, non-isothermal distributions the problem becomes non-trivial .

Clearly for a real cusp composed of a distinct population of roughly stellar mass objects, the black hole moves. We want to estimate the mean free path, time scale for wandering and whether the black holes “carry the cusp with them” allowing an adiabatic response of the stars to the black hole displacement. We know that non-isothermal cusps in general have $`v_{BH}=\eta v_{isoBH}`$, $`\eta >1`$ (cf Chatterjee et al 2002 ).

We can assume ballistic loss cone exit. ie consider the problem in the limit where the black hole moves on some ballistic trajectory relative to the centre of mass of the stars, changing directions and velocity amplitude only at discrete intervals, with the change occurring rapidly compared to the duration of ballistic motion. A compact stellar remnant is generally captured if it comes within $`(10r_S)`$ of the black hole. If the the black hole has been displaced by a comparable amount, then a “new loss cone” has been entered. The time scale for this to occur is:

$$t_{lc}10r_S/\eta v_{isoBH}M_{BH}^{3/2}$$

For a $`10^6\text{M}_{}`$ BH embedded in a $`\sigma 100`$ km/sec cusp, we find $`t_{lc}1/\eta 10^8`$ seconds, or about 1 year. The time scale is longer if the mean-free path is shorter, since in that case the black hole has to random walk out of loss cone, rather than being displaced ballistically.

The motion should be ballistic on time scales comparable to some fraction of the cusp dynamical time, very roughly, $`10^{23}`$ years.

So, orbits with time scale less than 1 year are carried with the BH as it wanders. Which are orbits inside $`100AU`$ for our canonical $`10^6\text{M}_{}`$SMBH.

This is $`10^4r_S`$, which is precisely where the dominant influx of compact objects is coming from.

Therefore BH wandering is significant, but possibly not dominant, with a critical mass close to $`10^6\text{M}_{}`$ \- more massive SMBH wander slowly and carry cusp with them, less massive SMBH leave loss cone ballistically.

It is critical to calculate $`\eta `$ for non-isothermal cusps in order to estimate a robust coalescence rate. This may require numerical simulations.

We also want the RMS displacement $`r_{BH}`$, as well as the mean free path $`r_{BH}`$. The black hole is effectively fixed to the centre when $`r_{BH}r_S`$.

### B Future broken symmetries

* Triaxiality.
If we can’t move the SMBH, can we move the stars?
If the nucleus is triaxial, $`J`$ is not conserved for individual stars on orbits - box orbits and boxlets walk through centre. Close to the SMBH, orbits are perturbed Keplerian orbits, precessing ellipses, the transition is inside $`r_h`$, at what radius the perturbed Keplerian limit is appropriate is important.
* Chaos:
Scattering of stars by SMBH can lead to dynamical chaos where the only integral of motion conserved for stars moving in the mean potential of the galaxy is the energy, ignoring long time scale relaxation processes.
Ergodic orbits $``$ spherical potential (cf Gerhard & Binney 1985).
However, triaxial systems can be constructed with a central SMBH in which chaos is suppressed and triaxiality persist to small radii, well inside $`r_h`$ (see Holley–Bockelmann et al 2002).

Regular–regular scattering, allows continued centrophilic orbits, and a persistent flux of stars to very small radii on dynamical rather than diffusive time scales. This is a hard problem, and in particular we don’t know what real galaxies do yet, the existence of theoretical models with suppressed chaos is not sufficient to show that real galaxies can achieve this limit, and chaos may dominate the evolution of the centres of real galaxies with supermassive black holes.

We need to determine flux to low $`J`$ from high $`E`$ (large $`r`$) and the time scales for evolution of potential shape (cf Zhao et al 2002 ).
* Mergers:
A significant fraction of galaxies has undergone mergers with other galaxies, including low mass satellite galaxies. Merger processes provide for strongly non-axisymmetric time dependent central potentials which leads to refilling of the loss cone on dynamical time scales, for many cusp dynamical times. Most of the integrated rate may be from post-merger galaxies rather than quasi-static passively evolving galaxies. After the stellar population involved in a galaxy merger has settled into a stationary distribution, the SMBH brought together by the collision may evolve towards coalescence through dynamical processes. At late stages of SMBH–SMBH coalescence is driven by gravitational radiation, and therefore the inner stellar cusp around each of the SMBH coming together will see a strongly time dependent potential. This may drive significant numbers of stars into the gravitational radiation loss-cone leading to substantial enhancement of the rate of coalescence of low mass objects during the process of binary merger.

AGN activity may correlate with galaxy mergers, and coalescence may thus be enhanced during AGN phases. This may lead to significant fraction of coalescences occurring when the SMBH has a substantial gaseous disk, and the torquing of the inner accretion stream during late stages of in-spiral, producing a low amplitude periodic flickering of the AGN, which would provide a valuable electro–magnetic signature for the counterpart to the gravitational radiation emission.
* Spin:
Wilson & Colbert  noted that major mergers may both lead to near maximal spin black holes, and cause abrupt changes in the spin parameter of the black holes. Black hole spin will change the shape and extent of the loss-cone, but at a level small compared to other uncertainties. However, the waveform of gravitational radiation emitted during inspiral is sensitive to the spin of the SMBH .

Two issues arise: one is whether we can predict the fraction of near maximal Kerr black holes from observed statistics of accreting systems and models of major mergers; and, whether the observations of spin parameters from inspiraling compact objects can test models of SMBH evolution. Both issues are largely open and need to be explored further.
* Star Formation:
The largest outstanding uncertainty are the IMF, the low mass BH mass function, and whether star formation is ongoing in galactic nuclei, and if so, with what duty cycle.,
* Pop III:
If the Pop III IMF is peaked strongly to high mass , then at high z, intermediate mass black holes are strong contributors to the coalescence rate and observable LISA event rate.
* Black Hole Mass Function:
We need the zero mass cut-off masses, frequency and distribution of masses for low mass black holes as function of the IMF and metallicity (see Fryer this proceedings and references therein). Knowing the binarity and initial spin of massive stars and their compact remnant descendants would be nice too.
* Ongoing star formation:
If there is ongoing “normal” nuclear star formation in inner 1-2 pc of normal spirals, with as little as 1% duty cycle, then compact objects are replenished and high coalescence rate sustained. That is, if there is star formation in galactic nuclei, and if the mass function of the stars formed is similar to that of stars formed in the field, or possibly biased towards high mass stars, then the nuclear star formation is important for generating a new population of compact stars in the nuclear regions. This is a particularly important process if low mass black holes continue to form in the metal rich gas found in galactic centres at late times.

There is some observational evidence for ongoing central star formation in the Milky Way . An open question is whether there is in situ formation, or off centre clusters that then sink to the centre . We need to understand the formation channel and IMF of stars formed from metal rich gas deep in galactic potentials, in possibly intense radiation environments, near the central black holes.

## III Conclusions

* Coalescence of low mass compact objects with SMBH is certain.
* Rates are uncertain but consistent and optimistically high.
* The outstanding model uncertainties are resolvable in principle.

The coalescence of low mass compact stellar remnants with central supermassive black holes in galaxies is one of the major hoped for LISA sources, and detection of such events with reasonable S/N can provide strong tests of both astrophysical processes and fundamental physics.

LISA will hopefully test some of our models, some model parameters will probably remain free until LISA constrains them.

###### Acknowledgements.

This research was supported in part by NSF grant PHY-0203046 and the Center for Gravitational Wave Physics, an NSF funded Physics Frontier Center at Penn State.


text:
# On Nagata’s conjecture

## 1 Introduction

Let $`X`$ be a normal projective variety of dimension $`n`$ over an algebraically closed field $`k`$, and $`L`$ an ample divisor. Given $`r`$ points $`p_1,\mathrm{},p_rX`$ and an integer $`1dn`$, the $`d`$-dimensional Seshadri constant of $`L`$ at the points $`p_1,\mathrm{},p_r`$ is the real number

$$\epsilon _d(L,p_1,\mathrm{},p_r)=\sqrt[d]{\underset{Z}{inf}\left\{\frac{L^dZ}{\mathrm{mult}_{p_i}Z}\right\}},$$

where $`Z`$ runs over all positive $`d`$-dimensional cycles (see \[21, 1.1\], or  for the original definition). As $`L`$ is ample, we have $`L^dZ>0`$ for all $`Z`$, so $`\frac{L^dZ}{{\scriptscriptstyle \mathrm{mult}_{p_i}Z}}_+\mathrm{}`$. Moreover, there exist $`Z`$ which contain some point $`p_i`$, and therefore the Seshadri constant is indeed a finite real number.

Most work on Seshadri constants deals with the $`d=1`$ case, and usually one writes $`\epsilon (L,p_1,\mathrm{},p_r)=\epsilon _1(L,p_1,\mathrm{},p_r)`$; we shall be concerned here with the codimension 1 case $`(d=n1)`$. Also, we use the shorthand notation $`\epsilon _d(L,r)=\epsilon _d(L,p_1,\mathrm{},p_r)`$ for very general points $`p_1,\mathrm{},p_r`$ (i.e., in the intersection of countably many Zariski open subsets of $`X^r`$) which is the case we are mostly interested in.

In connection with his solution to the fourteenth problem of Hilbert, Nagata posed in  (in different terminology) the following conjecture concerning Seshadri constants of the plane:

###### Conjecture 1 (Nagata).

If $`r9`$, then

$$\epsilon (𝒪_{P^2}(1),r)=1/\sqrt{r}.$$

If $`r=s^2`$ is a square, then it is not hard to prove that the conjecture is true, and in fact Nagata proved a slightly stronger result in this case. In a variety of dimension $`n`$ it is also not difficult to prove that

$$\epsilon (L,p_1,\mathrm{},p_r)\sqrt[n]{\frac{L^n}{r}}$$

for every set of $`r`$ points (see \[20, Remark 1\], for example), so Nagata’s conjecture claims that the Seshadri constant of a very general set of $`r9`$ points in the plane is maximal. All available information on Seshadri constants (see , , , , , ,  for the case of surfaces, , ,  for dimension $`n>2`$) suggests that, in fact, in an arbitrary variety, for $`r`$ large enough, the Seshadri constant of $`r`$ very general points is maximal. This led Szemberg to propose in  the following generalization:

###### Conjecture 2 (Nagata-Szemberg).

Given a smooth variety of dimension $`n`$ and $`L`$ an ample divisor on $`X`$, there exists a number $`r_0=r_0(X,L)`$ such that for every $`rr_0`$

$$\epsilon (L,r)=\sqrt[n]{\frac{L^n}{r}}.$$

This note is devoted to the following result, that gives a lower bound for $`(n1)`$-dimensional Seshadri constants of $`r`$ very general points in a variety, relating them to the analogous constants of one point in the same variety and of $`r`$ points in projective space:

###### Theorem 3.

Let $`X`$ be a normal projective variety of dimension $`n2`$, $`L`$ an ample divisor. Then for every smooth point $`pX`$ and $`r1`$,

$$\epsilon _{n1}(L,r)\epsilon _{n1}(L,p)\epsilon _{n1}(𝒪_^n(1),r).$$

Combining it with known results on the value of the Seshadri constants in projective space, Theorem 3 implies more explicit relations between $`r`$-point and 1-point Seshadri constants. The consequences support the Nagata-Szemberg conjecture, especially in the case of surfaces (note that for surfaces, $`(n1)`$-dimensional Seshadri constants are the usual Seshadri constants).

###### Corollary 4.

Suppose $`r=s^n`$ for some integer $`s`$. Then for every smooth point $`pX`$ and $`r1`$,

$$\epsilon _{n1}(L,r)\epsilon _{n1}(L,p)/s.$$

###### Proof.

Use G. V. Choodnovsky’s result  that $`\epsilon _{n1}(𝒪_^n(1),s^n)=1/s`$. ∎

If $`r`$ is not the $`n`$-th power of an integer, then we do not know the exact value of $`\epsilon _{n1}(𝒪_^n(1),r)`$ (there is a conjecture similar to Nagata’s, posed by Choodnovsky in the same paper , and by A. Iarrobino in ). However, B. Harbourne pointed out that, using results of J. Alexander and A. Hirschowitz and of M. Hochster and C. Huneke, one can prove an asymptotically optimal bound for $`\epsilon _{n1}(𝒪_^n(1),r)`$. Combining it with Theorem 3, in section 2 we prove the following asymptotic bound for $`\epsilon _{n1}(L,r)`$ that depends only on the $`(n1)`$-dimensional Seshadri constant of $`L`$ at a smooth point $`pX`$:

###### Corollary 5.

For every $`ϵ>0`$ there exists an integer $`s`$, depending only on $`ϵ`$ and $`n`$, such that for every smooth point $`pX`$ and $`rs`$,

$$\epsilon _{n1}(L,r)\sqrt[n]{r}\epsilon _{n1}(L,p)ϵ.$$

###### Corollary 6.

If $`X`$ is a normal projective surface, $`L`$ an ample divisor on $`X`$, then for every smooth point $`pX`$ and $`r9`$, Nagata’s conjecture implies that

$$\epsilon (L,r)\frac{\epsilon (L,p)}{\sqrt{r}}.$$

Observe that this tells us that Nagata’s conjecture (on the plane) implies the Nagata-Szemberg conjecture on a large family of surfaces. Indeed, the obtained bound is equal to $`\epsilon (L,p)/\sqrt{L^2}`$ times the conjectured value of $`\epsilon (L,r)`$, so we have the following:

###### Corollary 7.

If $`X`$ is a normal projective surface, $`L`$ an ample divisor on $`X`$ and $`pX`$ is a smooth point such that $`\epsilon (L,p)=\sqrt{L^2}`$, then Nagata’s conjecture implies the Nagata-Szemberg conjecture on $`(X,L)`$ with $`r_0(X,L)9`$.

In particular, we can apply this to complex surfaces with Picard number 1, using A. Steffens’ result \[20, Proposition 1\], which says that if $`L`$ is an ample generator of $`NS(X)`$ then $`\epsilon (L,p)\sqrt{L^2}`$ for very general points:

###### Corollary 8.

Let $`X`$ be a smooth projective surface defined over $``$, $`L`$ an ample generator of $`NS(X)`$ and assume $`L^2=d^2`$ is a square. Then Nagata’s conjecture implies the Nagata-Szemberg conjecture on $`X`$ with $`r_0(X,L)9`$.

In particular, the Seshadri constant of $`r9`$ very general points on a complex surface $`X`$ with Picard number equal to 1, is maximal if both $`L^2`$ and the number of points $`r`$ are squares. This can be compared to B. Harbourne’s result \[7, I.2.b(i)\] (over an arbitrary base field and with no assumption on the Picard number) that the Seshadri constant is maximal whenever $`L`$ is very ample, $`rL^2`$ is a square, and $`rL^2`$.

Also, known bounds approximating Nagata’s conjecture give new bounds on surfaces; for instance, H. Tutaj-Gasińska’s bound in  showing that $`\epsilon (𝒪__{}^2(1),r)\frac{1}{\sqrt{r+\frac{1}{12}}}`$ gives the following:

###### Corollary 9.

If $`X`$ is a normal projective surface defined over $``$, then for every smooth point $`pX`$ and $`r>9`$,

$$\epsilon (L,r)\frac{\epsilon (L,p)}{\sqrt{r+\frac{1}{12}}}.$$

In a similar vein, Harbourne’s bounds on Seshadri constants of $`^2`$ in  and  imply that

###### Corollary 10.

Let $`X`$ be a normal projective surface, $`pX`$ a smooth point and $`r1`$. Then for every pair of integers $`1sr`$, $`1d`$, it holds

$$\epsilon (L,r)\{\begin{array}{cc}\frac{s}{rd}\epsilon (L,p)\hfill & \text{ if }s^2rd^2,\hfill \\ \frac{d}{s}\epsilon (L,p)\hfill & \text{ if }s^2rd^2.\hfill \end{array}$$

It should be mentioned, however, that both Corollaries 9 and 10 are usually weaker (but stronger for some surfaces and numbers of points) than Harbourne’s results on algebraic surfaces of , where he gets the bounds

$$\epsilon (L,r)\{\begin{array}{cc}\frac{s}{rd}\hfill & \text{ if }s^2rd^2L^2,\hfill \\ \frac{dL^2}{s}\hfill & \text{ if }s^2rd^2L^2\hfill \end{array}$$

for very ample $`L`$, assuming moreover that $`rL^2`$.

The proof of Theorem 3 is based on the idea, due to L. Évain (see ) that $`r`$ point Seshadri constants of the plane can be computed by means of homotetic collisions of fat points. For convenience of the exposition we express this, generalized to $`n`$-dimensional projective space, in terms of the order of a nonreduced curve singularity, rather than collisions. Then, we observe that it is enough to know the formal germ of such a singularity, and the fact that the completion of the local ring at a smooth point of a variety is a ring of formal power series that only depends on the dimension of the variety, to obtain the bound.

It might be interesting to note that the Viro method developed by E. Shustin in  can also be used to relate the existence of singular curves on $`^2`$ with the existence of singular curves in algebraic surfaces, see for instance \[12, §5\] or \[3, 3.A\].

We thank B. Harbourne for many valuable comments, which largely improved the paper.

## 2 Singularities of arrangements of multiple lines

To lighten notations for rings of polynomials and of power series, we write $`𝐱=(x_0,`$ …, $`x_n)`$ for a collection of variables, so $`k[𝐱]`$ and $`k[[𝐱]]`$ denote $`k[x_0,\mathrm{},x_n]`$ and $`k[[x_0,\mathrm{},x_n]]`$, respectively. Also, $`𝔪`$ and $`\widehat{𝔪}`$ will be the maximal ideals generated by $`x_0,\mathrm{},x_n`$ in $`k[𝐱]`$ and $`k[[𝐱]]`$ respectively, and for every point $`p=[\xi _0:\mathrm{}:\xi _n]`$ in projective $`n`$-space, we use the notation $`I_p`$ (respectively $`\widehat{I}_p`$) for the ideal generated by the $`2\times 2`$ minors of the matrix

$$\left(\begin{array}{ccc}\xi _0& \mathrm{}& \xi _n\\ x_0& \mathrm{}& x_n\end{array}\right)$$

in $`k[𝐱]`$ (respectively in $`k[[𝐱]]`$).

Given distinct points $`p_1,\mathrm{},p_r^n`$ and $`𝐦=(m_1,\mathrm{},m_r)_0^r`$, we define $`\alpha _𝐦(p_1,\mathrm{},p_r)`$ to be the minimal degree of a homogeneous polynomial vanishing to order $`m_i`$ at $`p_i`$. As $`I_p`$ is homogeneous for all $`p`$, this number coincides with the maximal integer $`\alpha `$ such that

$$I=\underset{i=1}{\overset{r}{}}I_{p_i}^{m_i}𝔪^\alpha ,$$

or equivalently, such that $`\widehat{I}=\widehat{I}_{p_i}^{m_i}\widehat{𝔪}^\alpha `$. In other words, $`\alpha _𝐦(p_1,\mathrm{},p_r)`$ is the order at the origin of $`𝔸^{n+1}`$ of the arrangement of multiple lines defined by $`I`$ (which is the affine cone over the fat point scheme consisting of the points $`p_i`$ with multiplicities $`m_i`$).

###### Remark 11.

The definition of $`\alpha `$ and the $`(n1)`$-dimensional Seshadri constants immediately give that $`𝐦`$,

$$\alpha _𝐦(p_1,\mathrm{},p_r)\left(\epsilon _{n1}(𝒪_^n(1),p_1,\mathrm{},p_r)\right)^{n1}\underset{i=1}{\overset{r}{}}m_i.$$

Let $`X`$ be a variety of dimension $`n2`$, $`qX`$ a smooth point, and fix uniformizing parameters $`x_1,\mathrm{},x_n`$ in some neighborhood $`V`$ centered at $`q`$. So (see \[17, §III.6\], for example) the germs of $`x_1,\mathrm{},x_n`$ at $`q`$ generate the maximal ideal in the local ring $`𝒪_{X,q}`$ and the morphism of $`k`$-algebras $`k[𝐱]𝒪_X(V)`$ determines an étale morphism $`\phi :V𝔸^nT_qX`$ (where $`VX`$ is open) and an isomorphism $`k[[𝐱]]\widehat{𝒪}_{X,q}`$ .

To every $`p=[\xi _0:\mathrm{}:\xi _n]^n`$ such that $`\xi _00`$ and $`p[1:0:\mathrm{}:0]`$ (so that $`p^{}=[\xi _1:\mathrm{}:\xi _n]^{n1}`$), we shall assign an irreducible curve $`C_p`$, smooth at $`q`$, and a regular parameter $`\overline{x}_pO_{C_p,q}`$. Define $`C_p`$ as the closure of the component through $`q`$ of $`\phi ^1(L_p^{})`$, where $`L_p^{}𝔸^n`$ is the affine cone over $`p^{}`$, i.e., the line through the origin in the direction determined by $`p^{}`$. As $`\phi `$ is étale, $`C_p`$ is smooth at $`q`$, and the ideal $`I_{C_p}𝒪_{X,q}`$ of its germ at $`q`$ is generated by the $`2\times 2`$ minors of the matrix

$$\left(\begin{array}{ccc}\xi _1& \mathrm{}& \xi _n\\ x_1& \mathrm{}& x_n\end{array}\right).$$

Now consider $`x_p=x_i\xi _i/\xi _0`$ for some $`\xi _i0`$, $`i1`$. It is easy to see that the restriction $`\overline{x}_p`$ of $`x_p`$ to $`C_p`$ does not depend on the choice of $`i`$, and that it is a uniformizing parameter for the curve. Also, abusing notation, in $`𝒪_{C_p\times X,(q,q)}=𝒪_{C_p,q}𝒪_{X,q}`$ we write $`x_p=\overline{x}_p1`$, $`x_1=1x_1`$, …, and $`x_n=1x_n`$; then the ideal of the germ of the diagonal $`\mathrm{\Delta }(C_p)C_p\times X`$ is generated by the $`2\times 2`$ minors of the matrix

$$\left(\begin{array}{cccc}\xi _0& \xi _1& \mathrm{}& \xi _n\\ x_p& x_1& \mathrm{}& x_n\end{array}\right),$$

in other words, $`\mathrm{\Delta }(C_p)`$ is the closure of the component through the origin of $`\psi ^1(L_p)`$, where $`L_p`$ is the affine cone over $`p^n`$, i.e., a line $`L_p𝔸^{n+1}=T_{(q,q)}(C_p\times X)`$, and $`\psi `$ is the étale morphism given by the parameters $`x_p`$, $`x_1`$, …, $`x_n`$.

With these notations, Theorem 3 follows from the more precise proposition:

###### Proposition 12.

Let $`X`$ be a variety of dimension $`n2`$, $`L`$ an ample divisor, $`qX`$ a smooth point, $`p_1,\mathrm{},p_r^n[1:0:\mathrm{}:0]`$ distinct points not on the hyperplane $`\xi _0=0`$. Then for very general points $`q_kC_{p_k}`$ it holds

$$\epsilon _{n1}(L,q_1,\mathrm{},q_r)\epsilon _{n1}(L,q)\epsilon _{n1}(𝒪_^n(1),p_1,\mathrm{},p_r).$$

###### Proof.

First note that due to the semicontinuity of multiplicity (see \[15, §8\] or \[16, §3\]), for each component $``$ of the Hilbert scheme of hypersurfaces in $`X`$, and each system of multiplicities $`𝐦`$, the sets of points $`(q_1,\mathrm{},q_r)`$ such that there is $`Y`$ with multiplicity $`m_i`$ at $`q_i`$ form a Zariski-closed subset of $`X^r`$. Thus, it will be enough to prove that, given $``$ and $`𝐦`$, the existence of $`Y`$ with multiplicity $`m_k`$ at $`q_k`$ for general $`q_kC_{p_k}`$ implies

$$YL^{n1}\left(\epsilon _{n1}(L,q)\epsilon _{n1}(𝒪_^n(1),p_1,\mathrm{},p_r)\right)^{n1}\underset{k=1}{\overset{r}{}}m_k.$$

So, fix $``$ and $`𝐦`$, and assume that for general points $`q_kC_{p_k}`$ there is a hypersurface $`Y`$ going through $`q_k`$ with multiplicity at least $`m_k`$.

In the local ring of $`C_{p_k}`$ at $`\mathrm{\Delta }(q)=(q,\mathrm{},q)`$, the $`1\mathrm{}\overline{x}_{p_k}\mathrm{}1`$, $`k=1,\mathrm{},r`$, form a regular system of parameters; abusing notation we call them simply $`x_{p_k}`$. Let $`\mathrm{\Gamma }C_{p_k}`$ be the irreducible curve defined locally by the equations $`x_{p_1}=\mathrm{}=x_{p_r}`$, which is obviously smooth at $`\mathrm{\Delta }(q)`$, and admits $`x_0:=\overline{x}_{p_1}𝒪_{\mathrm{\Gamma },\mathrm{\Delta }(q)}`$ as a local parameter. For every $`p_k=[\xi _{k,0}:\mathrm{}:\xi _{k,n}]`$, the ideal $`J_{p_k}𝒪_{\mathrm{\Gamma }\times X,\mathrm{\Delta }(q)}`$ generated by the $`2\times 2`$ minors of the matrix

$$\left(\begin{array}{ccc}\xi _{k,0}& \mathrm{}& \xi _{k,n}\\ x_0& \mathrm{}& x_n\end{array}\right),$$

defines the germ of a curve $`C_{p_k}^{}`$ (the preimage of a line in $`𝔸^{n+1}T_{\mathrm{\Delta }(q)}(\mathrm{\Gamma }\times X)`$) whose projection to $`X`$ is exactly $`C_{p_k}`$; more precisely, the fiber of $`C_{p_k}^{}`$ over $`\gamma =(q_1,\mathrm{},q_n)\mathrm{\Gamma }`$ is $`q_kC_{p_k}`$.

Consider now the diagonals $`\mathrm{\Delta }_{i,j}=\{(q_1,\mathrm{},q_r)C_{p_k}|q_i=q_j\}`$, $`\mathrm{\Delta }=\mathrm{\Delta }_{i,j}`$, and $`U=\mathrm{\Gamma }\mathrm{\Delta }`$. As the points $`p_1`$, …, $`p_r`$ are distinct, $`U`$ is a nonempty open subset of $`\mathrm{\Gamma }`$. The assumption on the existence of $`Y`$ tells us that there is an effective relative Weil divisor $`𝐘_UU\times X`$ (flat over $`U`$) whose fiber over $`(q_1,\mathrm{},q_r)U`$ belongs to $``$ and has multiplicity at least $`m_i`$ at $`q_i`$. By the smoothness of $`\mathrm{\Gamma }`$ at $`\mathrm{\Delta }(q)`$, $`𝐘_U`$ can be extended to a flat family $`𝐘U\{\mathrm{\Delta }(q)\}\times X`$, and then the condition on the multiplicity of the fibers of $`𝐘`$ means that $`𝐘`$ contains the arrangement of multiple curves whose germ at $`\mathrm{\Delta }(q)`$ is defined by the ideal

$$J=\underset{k=1}{\overset{r}{}}J_{p_k}^{m_k},$$

and this implies that the fiber $`Y_q`$ of $`𝐘`$ over $`\mathrm{\Delta }(q)`$ has multiplicity at least equal to the order of this arrangement at $`\mathrm{\Delta }(q)`$, i.e., at least equal to the maximal integer $`\alpha `$ such that $`J(x_0,\mathrm{},x_n)^\alpha `$. This can be computed equivalently as the order of the completion $`\widehat{J}`$ in $`\widehat{𝒪}_{𝔸^1\times X,(0,q)}k[[𝐱]]`$ but, by construction of the curves $`C_{p_k}^{}`$, one has $`\widehat{J}=\widehat{I}`$ as defined above, so this order is exactly $`\alpha _𝐦(p_1,\mathrm{},p_r)`$.

We remark also that, by the smoothness of $`X`$ in a neighborhood of $`q`$, $`𝐘`$ is defined by a principal ideal at $`q`$, so we get a Weil divisor $`Y_q^{}Y_q`$ with

$`L^{n1}Y_q^{}`$ $`=L^{n1}Y,`$
$`\mathrm{mult}_qY_q^{}`$ $`\alpha _𝐦(p_1,\mathrm{},p_r),`$
which together with the definition of the 1-point Seshadri constant, gives $`L^{n1}Y\epsilon _{n1}(X,q)\alpha _𝐦(p_1,\mathrm{},p_r)`$, and then it is enough to apply the bound of Remark 11. ∎

###### Proof of Corollary 5.

By Theorem 3 it is enough to see that

$$\underset{r\mathrm{}}{lim}\epsilon _{n1}(𝒪_^n(1),r)\sqrt[n]{r}1$$

(which in fact means that one has an equality, the converse inequality being well-known). More precisely, we shall prove that given $`k>0`$ there exists $`s_k=s_k(n)`$ such that if $`rs_k`$ then for all $`𝐦=(m_1,\mathrm{},m_r)`$ and general points $`p_1`$, …, $`p_r`$,

$$\alpha _𝐦(p_1,\mathrm{},p_r)\frac{m_i}{\sqrt[n]{r^{n1}}}\frac{k+1}{k+n}.$$

So let $`F`$ be a homogeneous polynomial defining a hypersurface of degree $`d`$ in $`^n`$ which has multiplicity $`m_i`$ at $`p_i`$ for general points $`p_1`$, …, $`p_r`$. Then, by the genericity of the points, for every permutation $`\sigma 𝒮_r`$ there is a polynomial $`F_\sigma `$ which has multiplicity $`m_i`$ at the point $`p_{\sigma (i)}`$. Therefore $`G=_{\sigma 𝒮_r}F_\sigma `$ is a polynomial of degree $`D=r!d`$ which has (the same) multiplicity $`M=(r1)!m_i`$ at $`p_1`$, …, $`p_r`$, and $`G^{k+n}`$ has degree $`(k+n)D`$ and multiplicity $`(k+n)M`$ at each point. By \[9, Theorem 1.1(a)\], applied to the ideal $`I`$ of the (reduced) scheme $`\{p_1,\mathrm{},p_r\}`$, this implies the existence of hypersurfaces of degree $`t(k+n)D/M`$ with multiplicity at least $`k+1`$ at $`p_1`$, …, $`p_r`$. Now write $`𝐦^{}=(k+1,\mathrm{},k+1)`$; by \[1, Corollary 1.2\], there is $`s_k(n)`$ such that if $`rs_k(n)`$ then $`\alpha _𝐦^{}(p_1,\mathrm{},p_r)(k+1)\sqrt[n]{r}`$ (again, because the points are general). Therefore we get $`(k+n)D/M(k+1)\sqrt[n]{r}`$ and

$$\frac{d}{m_i}=\frac{D}{rM}\frac{1}{\sqrt[n]{r^{n1}}}\frac{k+1}{k+n},$$

as wanted. ∎