text:
# A SHORTER PROOF FOR THE TRANSITIVITY OF TRANSFINITE CONNECTEDNESS

Abstract — A criterion is established for the transitivity of connectedness in a transfinite graph. Its proof is much shorter than a prior argument published previously for that criterion.

Key Words: Transfinite graphs, transfinite connectedness, nondisconnectable tips.

The objective of this work is provide a proof of a criterion for the transitivity of transfinite connectedness that is much shorter than a previously published proof \[1, pages 57-63\] for that criterion.

We shall use the terminology defined in that book  and will add a few more definitions here. $`𝒢^\nu `$ will denote a transfinite graph of rank $`\nu `$ with $`0\nu \omega `$. Two branches in $`𝒢^\nu `$ are said to be $`\rho `$-connected if there exists a $`\gamma `$-path $`(0\gamma \rho \nu )`$ that terminates at a 0-node of each branch. Then, $`𝒢^\nu `$ is called transfinitely connected if every two of its branches are $`\rho `$-connected for some $`\rho `$ depending on the branches. Transfinite connectedness is reflexive and symmetric binary relationship between the branches of $`𝒢^\nu `$, but it need not be transitive. Examples of such nontransitivity are given in \[1, Examples 3.1-5 and 3.1-6\] and \[2, Example 3.1-1\].

However, the transitivity of $`\rho `$-connectedness can be assured by restricting the graph $`𝒢^\nu `$ appropriately. Two criteria, each of which suffices for such transitivity, are \[1, Condition 3.2-1 and Condition 3.5-1\]. The first criterion is rather lengthy and of limited use. The second one is quite simple to state and not overly restrictive, but the proof of it given in  is quite long and complicated.111Please note the correction for \[1, pages 63-69\] stated in the Errata for that book given in the URL: www.ee.sunysb.edu/~zeman. Namely, every elementary set of a sequence of any rank is required to be node-distinct. This was implied by the fact that the argument was based on paths, but it should have been explicitly stated. This latter case can be simplified by assuming that every node of $`𝒢^\nu `$ is pristine, that is, each node does not embrace a node of lower rank—and consequently is not embraced by a node of higher rank either.222See \[2, Sec. 3.1\] in this regard. Under this condition, the second criterion \[1, Condition 3.5-1\] works again and with a much shorter proof. It turns out that that latter proof can be modified to apply to nonpristine graphs as well, yielding thereby a more elegant argument for the transitivity of $`\rho `$-connectedness. This is the subject of the present work.

We need some definitions. Two nonelementary tips (not necessarily of the the same rank) are said to be disconnectable if one can find two representatives, one for each tip, that are totally disjoint. (Any two tips, at least one of which is elementary, are simply taken to be disconnectable.) As the negation of “disconnectable,” we say that two tips are nondisconnectable if every representative of one of them meets every representative of the other tip infinitely often, that is, at infinitely many nodes. We can state this somewhat differently, as follows. Two tips, $`t_1^\gamma `$ and $`t_2^\delta `$, are called nondisconnectable if $`P_1^\gamma `$ and $`P_2^\delta `$ meet at least once whenever $`P_1^\gamma `$ is a representative of $`t_1^\gamma `$ and $`P_2^\delta `$ is a representative of $`t_2^\delta `$.

Here is the sufficient condition that will ensure the transitivity of connectedness.

Condition 1. If two tips of ranks less than $`\nu `$ (and not necessarily of the same rank) are nondisconnectable, then either they are shorted together or at least one of them is open.

Before turning to our principle result, we need to introduce some more specificity into our terminology. The maximal nodes partition all the nodes of $`𝒢^\nu `$; indeed, two nodes are in the same set of the partition if they are embraced by the same maximal node. Thus, corresponding to any path, there is a unique set of maximal nodes that together embrace all the nodes of the path, and that set of maximal nodes becomes totally ordered in accordance with an orientation assigned to the path. Moreover, two maximal nodes are perforce totally disjoint. In order to avoid explicating repeatedly which nodes embrace which nodes, it is convenient to deal only with the maximal nodes of our graph $`𝒢^\nu `$.

Let us expand upon this point because it is the key idea that allows us to replace the long proof in \[1, pages 57-73\] by a shorter, modified version of the proof in \[2, pages 31-35\]. Given any path $`P^\rho `$ of any rank $`\rho `$ in $`𝒢^\nu `$, let $`𝒳`$ be the set of (possibly nonmaximal) nodes encountered in the recursive construction of $`P^\rho `$; that is, $`𝒳`$ consists of the nodes in the sequential representation of $`P^\rho `$ along with the nodes of lower ranks in the sequential representations of the paths between those nodes, and also along with the nodes of still lower ranks in the sequential representations of the paths of lower ranks in those latter paths, and so on down to the 0-nodes of the 0-paths embraced by $`P^\rho `$. Each node of $`𝒳`$ need not be maximal, but it has a unique maximal node embracing it. Let $`𝒳_m`$ denote the set of those maximal nodes. So, given $`P^\rho `$ and thereby $`𝒳`$ and $`𝒳_m`$, each node of $`𝒳_m`$ corresponds to a unique set of nodes in $`𝒳`$ such that, for every two nodes in that latter set, one embraces the other. The union of those sets is $`𝒳`$. In this way, there is a bijection between $`𝒳_m`$ and the collection of the said sets, which in fact comprise a partition of $`𝒳`$. Moreover, when $`P^\rho `$ has an orientation, that orientation induces a total ordering of $`𝒳_m`$. Thus, when dealing with the nodes of $`P^\rho `$, we can fix our attention on the nodes of $`𝒳_m`$. This we shall do. Any node $`x^\beta `$ of $`𝒳_m`$ will be called a maximal node for $`P^\rho `$ to distinguish it from the corresponding subset of nodes in $`𝒳`$. In general, we may have $`\beta >\rho `$, but not necessarily always. Also, we can transfer our terminology for the nodes of a path, such as “incident to a node,” “meets a node,” or “terminates at a node,” to those maximal nodes. For example, a path is incident to a node $`x𝒳_m`$ if one of the path’s terminal tips is embraced by $`x`$; this will be so if it has a terminal node embraced by $`x`$. In the latter case, we say that the path meets x and also terminates at $`x`$. Thus, two paths meet at $`x𝒳_m`$ if $`x`$ embraces a node $`y`$ of one path and a node $`z`$ of the second path such that $`y`$ embraces or is embraced by $`z`$.

The proof of our main result (Theorem 3) requires another idea, namely, “path cuts.” Let $`P^\rho `$ be a $`\rho `$-path with an orientation. Let $`𝒴`$ be the set of all branches and all (not necessarily maximal) nodes of all ranks in $`P^\rho `$. The orientation of $`P^\rho `$ totally orders $`𝒴`$. With $`y_1`$ and $`y_2`$ being two members of $`𝒴`$, we say that $`y_1`$ is before $`y_2`$ and that $`y_2`$ is after $`y_1`$ if in a tracing of $`P^\rho `$ in the direction of its orientation $`y_1`$ is met before $`y_2`$ is met. A path cut $`\{_1,_2\}`$ for $`P^\rho `$ is a partitioning of the set of branches of $`P^\rho `$ into two nonempty subsets, $`_1`$ and $`_2`$, such that every branch of $`_1`$ is before every branch of $`_2`$. Another way of stating this is as follows. The partition $`\{_1,_2\}`$ of the branch set of $`P^\rho `$ comprises a path cut for $`P^\rho `$ if and only if, for each branch $`b_1`$, every branch of $`P^\rho `$ before $`b`$ is also a member of $`_1`$.

Here is a result that is easily proven through induction on ranks.

Lemma 2. For each path cut $`\{_1,_2\}`$ for $`P^\rho `$, there is a unique maximal node $`x^\gamma `$ $`(\gamma \rho )`$ for $`P^\rho `$ such that every branch $`b_1_1`$ is before $`x^\gamma `$ and every branch $`b_2_2`$ is after $`x^\gamma `$.

We will say that the path cut occurs at the maximal node $`x^\gamma `$. As was stated before, a node is called maximal if it is not embraced by a node of higher rank. On the other hand, a maximal node for $`P^\rho `$ may embrace many other nodes. It follows that all the nodes of $`P^\rho `$ other than the nodes of $`P^\rho `$ embraced by $`x^\gamma `$ are also partitioned into two sets, the nodes of one set being before $`x^\gamma `$ and the nodes of the other set being after $`x^\gamma `$.

Here now is our main result.

Theorem 3. Let $`𝒢^\nu `$ $`(0\nu \omega )`$ be a $`\nu `$-graph for which Condition 1 is satisfied. Let $`x_a`$, $`x_b`$, and $`x_c`$ be three different maximal nodes in $`𝒢^\nu `$ such that, if any one of them is a singleton, its sole tip is disconnectable from every tip in the other two nodes. If $`x_a`$ and $`x_b`$ are $`\rho `$-connected and if $`x_b`$ and $`x_c`$ are $`\rho `$-connected $`(0\rho \nu )`$, then $`x_a`$ and $`x_c`$ are $`\rho `$-connected.

Proof. Note that, if any one of $`x_a`$, $`x_b`$, and $`x_c`$ is a nonsingleton, then each of its embraced tips (perforce, nonopen) must be disconnectable from every tip embraced by the other two nodes; indeed, otherwise, the tips of that node would be shorted to the tips of another one of those three nodes, according to Condition 1 and our hypothesis, and thus $`x_a`$, $`x_b`$, and $`x_c`$ could not be three different maximal nodes.

Let $`P_{ab}^\alpha `$ $`(\alpha \rho )`$ be a two-ended $`\alpha `$-path that terminates at the maximal nodes $`x_a`$ and $`x_b`$ and is oriented from $`x_a`$ to $`x_b`$, and let $`P_{bc}^\beta `$ $`(\beta \rho )`$ be a two-ended $`\beta `$-path that terminates at the maximal nodes $`x_b`$ and $`x_c`$ and is oriented from $`x_b`$ to $`x_c`$. Let $`P_{ba}^\alpha `$ be $`P_{ab}^\alpha `$ but with the reverse orientation. $`P_{ba}^\alpha `$ cannot have infinitely many $`\alpha `$-nodes because it is two-ended.

Let $`\{x_i\}_{iI}`$ be the set of maximal nodes at which $`P_{ba}^\alpha `$ and $`P_{bc}^\beta `$ meet, and let $`𝒳_1`$ be that set of nodes with the order induced by the orientation of $`P_{ba}^\alpha `$. If $`𝒳_1`$ has a last node $`x_l`$, then a tracing along $`P_{ab}^\alpha `$ from $`x_a`$ to $`x_l`$ followed by a tracing along $`P_{bc}^\beta `$ from $`x_l`$ to $`x_c`$ yields a path of rank no larger than $`\rho `$ that connects $`x_a`$ and $`x_c`$. Thus, $`x_a`$ and $`x_c`$ are $`\rho `$-connected in this case. This will certainly be so when $`\{x_i\}_{iI}`$ is a finite set.

So, assume $`𝒳_1`$ is an infinite, ordered set (ordered as stated). We shall show that $`𝒳_1`$ has a last node. Let $`Q_1`$ be the path induced by those branches of $`P_{ba}^\alpha `$ that lie between nodes of $`𝒳_1`$ (i.e., as $`P_{ba}^\alpha `$ is traced from $`x_b`$ onward, such a branch is traced after some node of $`𝒳_1`$ and before another node of $`𝒳_1`$, those nodes depending upon the choice of the branch.) Let $`_1`$ be the set of those branches. We can take it that $`P_{ba}^\alpha `$ extends beyond the nodes of $`𝒳_1`$, for otherwise $`𝒳_1`$ would have $`x_a`$ as its last node, and $`x_a`$ and $`x_c`$ would be $`\rho `$-connected. Therefore, we also have a nonempty set $`_2`$ consisting of those branches in $`P_{ba}^\alpha `$ that are not in $`_1`$. $`\{_1,_2\}`$ is a path cut for $`P_{ba}^\alpha `$. Consequently, by Lemma 2., there is a unique maximal node $`x_1^{\gamma _1}`$ at which that path cut occurs. Thus, $`Q_1`$ terminates at $`x_1^{\gamma _1}`$. Let $`t_1^{\rho _1}`$ be the $`\rho _1`$-tip through which $`Q_1`$ reaches $`x_1^{\gamma _1}`$. Every representative of $`t_1^{\rho _1}`$ contains infinitely many nodes of $`𝒳_1`$ (otherwise, $`𝒳_1`$ would have a last node).

Now, consider $`P_{bc}^\beta `$. We can take it that there is a maximal node $`x_d`$ in $`P_{bc}^\beta `$ different from $`x_c`$ such that the subpath of $`P_{bc}^\beta `$ between $`x_d`$ and $`x_c`$ is totally disjoint from $`P_{ba}^\alpha `$, for otherwise $`x_a`$ and $`x_c`$ would have to be the same node according to Condition 1 and our hypothesis again. We can partition the branches of $`P_{bc}^\beta `$ into two sets, $`_3`$ and $`_4`$, as follows. Each branch of the first set $`_3`$ is such that it lies before (according to the orientation of $`P_{bc}^\beta `$) at least one node in $`𝒳_1`$ of each representative of $`t_1^{\rho _1}`$, this being so for all such representatives.333Let us note here a correction for \[2, page 34, line 19 up\]. Replace the sentence on that line by the following sentence: “Let $`𝒩_2`$ be the set of those nodes in $`\{n_i\}_{iI}`$ that lie before (according to the orientation of $`P_{bc}^\beta `$) at least one node of $`\{n_i\}_{iI}`$ in each representative of $`t_1^{\rho _1}`$, this being so for all such representatives.” The second set $`_4`$ consists of all the branches of $`P_{bc}^\beta `$ that are not in $`_3`$. No branch of $`_4`$ can precede a branch of $`_3`$. Thus, we have a path cut $`\{_3,_4\}`$ for $`P_{bc}^\beta `$ and thereby (according to Lemma 2.) a unique maximal node $`x_2^{\gamma _2}`$ lying after the branches of $`_3`$ and before the branches of $`_4`$. Let $`Q_2`$ be the path induced by $`_3`$. It reaches $`x_2^{\gamma _2}`$ through some tip $`t_2^{\rho _2}`$. Furthermore, each representative of $`t_2^{\rho _2}`$ must meet each representative of $`t_1^{\rho _1}`$ at least once because they meet at at least one node of $`𝒳_1`$. Thus, $`t_1^{\rho _1}`$ and $`t_2^{\rho _2}`$ are nondisconnectable. Moreover, neither of those tips can be open (i.e., be in a singleton node) because the paths $`P_{ba}^\alpha `$ and $`P_{bc}^\beta `$ pass through and beyond their respective nodes $`x_1^{\gamma _1}`$ and $`x_2^{\gamma _2}`$. So, by Condition 1, $`x_1^{\gamma _1}`$ and $`x_2^{\gamma _2}`$ must be the same node because the tips $`t_1^{\rho _1}`$ and $`t_2^{\rho _2}`$ are shorted together.

This means that $`𝒳_1`$ has a last node, namely, $`x_1^{\gamma _1}=x_2^{\gamma _2}`$. It follows now that $`x_a`$ and $`x_c`$ are $`\rho `$-connected. $`\mathrm{}`$

The last proof has established the following two result.

Corollary 4. Under the hypothesis of Theorem 3, let $`P_{ab}^\alpha `$ be a two-ended $`\alpha `$-path connecting nodes $`x_a`$ and $`x_b`$, and let $`P_{bc}^\beta `$ be a two-ended $`\beta `$-path connecting nodes $`x_b`$ and $`x_c`$. Then, there is a two-ended $`\gamma `$-path ($`\gamma \mathrm{max}\{\alpha ,\beta \}`$) connecting $`x_a`$ and $`x_c`$ that lies in $`P_{ab}^\alpha P_{bc}^\beta `$.

Corollary 5. Under the hypothesis of Theorem 3, let $`\{x_i\}_{iI}`$ be the set of maximal nodes at which two two-ended paths meet. Assume that, if any terminal node of either path is a singleton, its sole tip is disconnectable from every tip in the two terminal nodes of the other path. Assign to $`\{x_i\}_{iI}`$ the total ordering induced by an orientation of one of these two paths. Then, $`\{x_i\}_{iI}`$ has both a first node and a last node.


format: text

text:
# Untitled Document

Supermassive black holes may be limited by the Holographic Bound

Paulo Sergio Custodio and J.E. Horvath

Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo

R. do Matão 1226, Cidade Universitária, 05508-900 São Paulo, Brazil

Supermassive Black Holes are the most entropic objects found in the universe. The Holographic Bound (HB) to the entropy is used to constrain their formation time with initial masses $`10^{68}M_{}`$, as inferred from observations. We find that the entropy considerations are more limiting than causality for this ”direct” formation. Later we analyze the possibility of SMBHs growing from seed black holes. The growth of the initial mass is studied in the case of accretion of pure radiation and quintessence fields, and we find that there is a class of models that may allow this metamorphosis. Our analysis generalizes recent work for some models of quintessence capable of producing a substantial growth in a short time, while simultaneously obeying the causal and Holographic Bound limits.

Supermassive Black Holes in the Universe

The continued observations of galaxies has revealed a hidden population of huge massive objects in compact nuclear regions of size $``$ few pc. Dynamical measurements taken along the last decade point out that the compact objects have masses in the range $`10^610^8M_{}`$, and perhaps more importantly, that every galaxy seems to host a central massive object . Even though some exotic alternatives have been proposed for their nature (e.g. neutrinoballs, see ), the simplest explanation is that the central parsecs of the galaxies are sites of residence of supermassive black holes (SMBH). Among the possible formation scenarios a hierarchical merging of smaller black holes has been suggested , although it is not guaranteed that the efficiency of the merging process is high enough to provide large masses. A likely alternative is that the SMBHs are primordial, i.e. preexist the galaxies, and perhaps are important for their very formation .

While the mass budget of the universe is not likely to be affected by the presence of the nuclear SMBHs, the total entropy will certainly be, since the entropy content of the black holes is huge. This suggests a connection between the formation of the SMBHs and the total entropy, possibly limited by the Holographic Bound, which has been proposed to limit the entropy enclosed in a given volume and may be deeply related to the fundamental theories .

We discuss in this work the issue of entropic limitations to the formation of SMBHs with masses $`10^6M_{}`$. After a brief presentation of the Holographic Bound and related concepts in Section 2, direct formation of SMBHs is addressed in Section 3. We analyze the possibility of growing the ”seed” black holes to those large values is addressed in Section 4. Conditions for fast growth due to accretion of a quintessence scalar field are addressed in Section 5. Section 6 discusses the role of causality in the process of accretion. Some general conclusions are given in Section 7.

The Holographic Bound

The Holographic Bound may be formulated by asserting that for a given volume $`V`$, the state of maximal entropy is the one containing the largest black hole that fits inside $`V`$, and this maximum is given by the finite area that encloses this volume. This idea generalizes a conjecture made by Bekenstein  in which this maximum is fixed by the non-gravitational energy within a sphere of size $`R`$, i.e. $`S<\frac{2\pi ER}{hc}`$ (now being properly called the Bekenstein limit).

Several analysis made in recent years reformulated this conjecture and proposed slightly different forms for the HB, but rather than discussing which one is correct we will base our argument on the very existence of some entropy bound, yet to be definitively identified.

To be concrete we shall assume the entropy $`S`$ to be bounded by the Bekenstein-Hawking value

$$S\frac{A}{4}$$
(1)
where $`A`$ is the area of the enclosed system under consideration. Unless explicitly indicated, we shall use natural units throughout this paper, then the Planck length $`L_{planck}^2=1`$ in eq.(1) above and so on.

Verlinde  observed some time ago that this bound must be modified in a cosmology with an arbitrary number of dimensions. Considering the Einstein space-time with the metric

$$ds^2=dt^2+R^2d\mathrm{\Omega }_{n}^{}{}_{}{}^{2}$$
(2)
where $`d\mathrm{\Omega }_{n}^{}{}_{}{}^{2}`$ is the line element of a unit n-dimensional sphere, the entropy of the conformal field in this space-time can be expressed in terms of its total energy $`E`$ and the Casimir $`E_C`$ by a generalized form of the Cardy-Verlinde formula as

$$S=\frac{2\pi }{n}R\sqrt{E_C(2EE_C)}$$
(3)
For a $`(n+1)`$ dimensional closed universe, the FRW equations are

$$H^2=\frac{16\pi G_n}{n(n1)}\frac{E}{V}\frac{1}{R^2}$$
(4)
$$\dot{H}=\frac{8\pi G_n}{(n1)}\left(\frac{E}{V}+P\right)+\frac{1}{R^2}$$
(5)
where $`H(t)=\frac{\dot{a}}{a}`$ is the Hubble parameter (describing the expansion/contraction of the universe), the dot stands for differentiation with respect to the proper time, $`E`$ is the total energy of matter filling the universe, and $`G_n`$ is the Newton constant in $`(n+1)`$ dimensions. $`a(t)`$ describes the scale factor of the Universe and $`R(t)a(t)`$ its physical size.

The FRW equation can then be related to three cosmological entropy bounds; the Bekenstein-Verlinde bound $`S_{BV}=\frac{2\pi }{n}ER`$, the Bekenstein-Hawking bound $`S_{BH}=(n1)\frac{V}{4G_nR}`$ (expressing that the black hole entropy is bounded by the area of the cosmological model), and the Hubble bound $`S_H=(n1)\frac{HV}{4G_n}`$ (a reflection of the fact that the maximal entropy is produced by black holes of the size of Hubble horizon). At a critical point defined by $`HR=1`$, all these three entropy bounds coincide with each other. Let us define $`E_{BH}`$ such that $`S_{BH}=\frac{(n1)V}{4G_nR}=\frac{2\pi }{n}E_{BH}R`$. Then, the first FRW equation takes the form

$$S_H=\frac{2\pi R}{n}\sqrt{E_{BH}(2EE_{BH})}$$
(6)
which is precisely of the same form as the Cardy-Verlinde formula. Its maximum reproduces the Hubble bound

$$S_H\frac{2\pi R}{n}E$$
(7)
Therefore, in some sense it may be said that the FRW dynamics ”knows” the value of the maximum entropy filling the universe. This connection between geometry and dynamics is a consequence of the Holographic Principle. However, and suggestive as these arguments are, we do not intend to further analyze them. It is just enough to keep in mind that independently of its specific form, an Holographic Bound is likely to hold for the actual 3+1 universe.

Another important concept needed for the study of SMBHs is the generalized second law of thermodynamics, formulated by J.Bekenstein  using a series of gedanken experiments. The generalized second law attempted to cure serious problems with the matter + radiation entropy as the latter was absorbed onto black holes (thus causing a growth of the black hole mass). Given that black holes need a few macroscopic parameters (mass, angular momentum and charge) for their description, the absorption of matter+radiation seemed to lead to a decrease of the entropy of the universe, since the matter+radiation entropy ended hidden behind the horizon. This was very problematic, since that this kind of Geroch process seems to go against the second law of thermodynamics $`\mathrm{\Delta }S>0`$. Bekenstein conjectured that the total entropy of the universe plus $`N`$ black holes is given by the sum of the matter+radiation entropy, plus the black hole entropy (which is proportional to the horizon area) in what is now known as the generalized second law (GSL). The GSL takes the form

$$S_{total}=S_{m+r}+\frac{1}{4}\underset{i}{\overset{N}{}}A_i$$
(8)
where the first contribution is the entropy associated to usual matter and radiation, and the second term describes the black hole contribution to the total entropy. Note that the entropy of just one black hole is numerically huge, $`S_{bh}10^{77}(M/M_{})^2`$, and this fact will be very important to set astrophysical constraints. According to the GSL, as long as we deal with classical process involving black holes and matter, the total variation of entropy must be positive

$$\mathrm{\Delta }S_{total}>0$$
(9)
In the next Sections we will evaluate some constraints to the mass and formation time of supermassive black holes using the concepts of HB and the GSL.

SMBH and Primordial Black Holes: direct formation

While the exact origin of the SMBHs is not known, it is possible that either a primordial process contributed to form them as they are, and that they have grown explosively from a seed population.

Let us discuss direct formation first. It is well-known that big black holes have a huge entropy, and if the bounds to the cosmic entropy apply, restrictions would arise for a ”direct formation”mechanism. Actually, if we impose

$$S_{smbh}(t_i)2\times 10^{77}N(M/M_{})^2<S_H(t_i)$$
(10)
where $`t_i`$ stands for the formation time hereafter, and it is further assumed that all black holes form more or less simultaneously, we may find an upper bound to the SMBH directly formed allowed by the HB. We start by evaluating $`S_H(t)`$ from eq.(7) for a 3+1 FRW universe. We identify $`R`$ with the particle horizon $`R=R_{hp}(t)t`$, and the total energy contained within this radius $`E=\frac{4\pi }{3}\varrho (t)R_{hp}^3(t)`$. We further restrict the analysis to the epochs in which $`a(t)t^n`$. Inserting all the quantities, we obtain $`E_{hp}(t)a^4(t)t^3t`$, for $`n=1/2`$. We do not need consider other forms of entropy in eq.(10), since black holes actually dominate the entropy budget by a large factor.

Then, multiplying by the particle horizon, we express the HB in the radiation-dominated era as $`S_{hp}(t)S_{hp}(t_D)(t/t_D)^2`$. For the matter-dominated era $`n=2/3`$ an analogous procedure yields $`E_{hp}(t)\varrho _m(t)R_{hp}^3(t)t`$. Therefore, the entropy contents become $`S_{hp}(t)S_{hp}(t_0)(t/t_0)^2`$, with $`t_0H_{0}^{}{}_{}{}^{1}`$ the present age of the universe ($`10^{17}h_{0}^{}{}_{}{}^{1}s`$) and $`t_D`$ ($`10^{13}s`$) is the radiation-matter decoupling time. For the sake of generality the dimensionless Hubble constant $`h_0=(H_0/100kms^1Mpc^1)`$ has not been fixed, although recent measurements suggest $`h_00.65`$. We also know that $`S_{hp}(t_0)8\times 10^{121}`$, therefore $`S_{hp}(t_D)\frac{8\times 10^{121}}{(t_0/t_D)^2}`$. Then, for the radiation-dominated era the condition $`S_{smbh}(N,M)<S_{hp}(t)`$ yields the maximum SMBH mass allowed by the entropy bound

$$M_{smbh}(t_i)<5.6\times 10^3h_0(t_i/s)(10^{11}/N)^{1/2}M_{}$$
(11)
if their number $`N`$ is equal to the number of galaxies and all them have been assumed to be of the same mass. The SMBHs can form directly only after

$$t_{min}(M)>5\times 10^7(M/10^8M_{})(N/10^{11})^{1/2}s$$
(12)
because before $`t_{min}`$ the entropy of these $`N`$ black holes would be larger than the entropy allowed by the Holographic Bound.

The same reasoning as above can be applied to the formation in the matter-dominated era with the result

$$M_{smbh}(t_i)<0.6h_0^2(t_i/s)(10^{11}/N)^{1/2}M_{}$$
(13)
and an earliest formation time

$$t_{min}(M)>1.8\times 10^9h_0^1(M/10^8M_{})(N/10^{11})^{1/2}s$$
(14)
Since the mass inside the horizon in the radiation-dominated era is just $`M_{hor}(t)7.6\times 10^{37}(t/1s)g`$, and achieves $`10^{15}M_{}`$ at its very end, we conclude that the availability of entropy is more restrictive than the demand of a causal formation. In other words, it is not sufficient to have a large horizon in which the SMBH can fit, to be allowed by the entropy of the HB seems to be even more important than that primary requirement.

A realistic and complete model would take into account an Initial Mass Function (IMF) for these black holes, with a general form given by $`(\frac{dN}{dM})`$. In this case we would need to replace the formula above by an integral of the form $`S_{smbh}(N,M)2\times 10^{77}_{M_1}^{M_2}𝑑M(\frac{dN}{dM})(M/M_{})^2`$, with $`M_1`$, $`M_2`$ the lower and upper limits to the masses determined by the specific physical conditions at formation. We shall not address this complex model here and leave it for future work. We have just considered a delta-type $`(\frac{dN}{dM})\delta (MM_{})`$, leading to the simple model given by $`S_{smbh}(N,M)N(M/M_{})^2`$ in this paper.

SMBH and Primordial Black Holes: grow of seeds to the SMBH scale

If the direct formation of SMBHs is difficult, one may wonder if there is still a possibility of starting with black holes of masses $`M_i10^6M_{}`$ which grow subsequently by accretion. Considering the absorption-evaporation processes of PBHs, we can identify epochs in which these objects grow or evaporate. The complete evolution of PBH mass is given semiclassically by

$$\left(\frac{dM}{dt}\right)=\frac{D}{M^2}+BM^2c[\underset{i}{}\varrho _i(t)]$$
(15)
where the term $`c[_i\varrho _i(t)]`$ describes the flux of whatever component flows through the black hole horizon, $`D10^{26}g^3s^1`$ is the evaporation constant (which, strictly speaking, depends somewhat on the number of degrees of freedom of the incoming material) and $`B=\frac{27\pi G^2}{c^4}4.6\times 10^{55}cm^2g^2`$ is the absorption constant (related to the cross-section of the black hole).

When writing down the semiclassical eq.(15) we have not ruled out any ”fuel” contributing to the growth of the black hole (quintessence has been proposed by Bean and Magueijo in ), provided their flux is large enough to contribute to the mass balance. If the black holes formed in the radiation-dominated or matter-dominated era, the main contribution to the second term is the flux of background radiation (quintessence will be explicitly addressed later). If we consider the radiation only, the balance of the r.h.s. terms define the critical mass, (see Refs.\[11- 13\]) for the instantaneous equilibrium between black holes and radiation, and its value is

$$M_c(t)\frac{10^{26}g}{T_{rad}(t)/T_0}$$
(16)
where the radiation temperature $`T_{rad}(t)a(t)^1`$ falls along the cosmological expansion and $`T_0`$ is the present temperature of the CMBR. If, say, $`t_i1s`$, the critical mass is then very small, and all the PBHs candidates to grow to SMBH were well above this instantaneous equilibrium mass value. Therefore, the Hawking radiation was negligible for them . If the PBHs are feed, they may grow with time, and the question is whether they can gain mass until the supermassive regime $`M10^6M_{}`$ is reached.

It is generally agreed that the gas and dust accretion are not likely to be important at early times (and do not have simple behavior with the cosmological time either). However, the radiation is always absorbed and might be important if the accretion rate is high enough . Neglecting the gas and dust fuels, and keeping only the radiation we may obtain analytic lower limits to the growth.

A first specific question we want to address is the following : if one black hole formed initially at $`t_i`$ satisfies the HB, i.e. $`S(M(t_i))10^{77}M(t_i)^2<S_{hp}(t_i)8\times 10^{121}(t_i/t_0)^2`$, is it automatically guaranteed that this object will satisfy the HB all the time?

To answer this question let us consider the most general flux $`F(\varrho (t))`$ absorbed by the event horizon by this black hole. Then, the mass accretion rate is given by

$$\left(\frac{dM}{dt}\right)=\frac{27\pi }{4}r_{g}^{}{}_{}{}^{2}F(\varrho (t))$$
(17)
If we choose $`\varrho =\varrho _{rad}`$, then $`F=c\varrho _{rad}(t)`$; for quintessence accretion we use $`F(\phi )=\frac{\dot{\phi }^2}{2}`$, see .

Solving formally the eq.(17) above, yields

$$M(t)=\frac{M_i}{[1\frac{27\pi M_i}{M_{pl}^4}_{t_i}^t𝑑t^{}F(\varrho (t^{}))]}$$
(18)
On the other hand, these black holes evolve obeying the HB if the local flux satisfies

$$_{t_i}^t𝑑t^{}F(\varrho (t^{}))<\frac{M_{pl}^{}{}_{}{}^{4}}{27\pi M_i}\left[1\frac{(M_i/M_{})(t_0/t)}{2\times 10^{22}}\right]$$
(19)
For the radiation flux, $`F(t)=c\varrho _{rad}(t)`$ and therefore eq.(19) requires

$$\varrho _{rad}(t_i)<6.3\times 10^{28}gcm^3(10^{15}g/M_i)(t_i/s)$$
(20)
For $`t_i1s`$ and $`M_i<1M_{}`$, typical of the radiation-dominated era, this condition is indeed satisfied.

The above general results suggest that the global constraint $`S_{total}<S_{hp}t^2`$, (i.e, the HB) implies some kind of restriction to the cross-section for the local accretion onto the black hole, as if the total flux through the horizon event would have to be modified. A detailed study of the flow into the black hole is needed to address this issue.

Quintessence models and the SMBH growth

Recent work by R.Bean and J.Magueijo  suggested an important growth of seed PBHs when a quintessence scalar field $`\phi `$ dominates the accretion. The key new ingredient is the role played by the kinetic term $`\frac{\dot{\phi }^2}{2}`$ in the flux onto the PBHs, which is absent in the case of absorption of pure radiation. We shall describe some growth solutions of the Bean-Magueijo model of quintessence around SMBHs. The evolution of the mass of these objects is given by the following formula

$$\left(\frac{dM}{dt}\right)=\eta CM^2\dot{\phi }^2$$
(21)
where $`C=\frac{27\pi }{2M_{pl}^{}{}_{}{}^{4}}`$ and $`\eta 1`$ is a parameter that measures the efficiency of the accretion process. Because of the uncertainties on the value of $`\eta `$ (related to the details of the depletion of material nearby the black hole, see Ref. for a recent discussion), we have left it free in our calculations, so that the adoption of a different value can be easily done. The connection between the potential $`V(\phi )`$ and $`\dot{\phi }`$ comes because the latter regulates the expansion rate and hence the behavior of the flux $`\dot{\phi }^2`$. The complete system to solve is given by eq.(21) above together with the dynamical equations

$$\ddot{\phi }+3H\dot{\phi }+V^{}=0$$
(22)
and

$$H^2=\frac{8\pi }{3M_{pl}^2}[V(\phi )+\varrho _{pbh}]$$
(23)
This set of differential equations is very difficult to solve, even in the approximation $`\varrho _{pbh}<<V(\phi )`$, which is relevant for this work. Bean and Magueijo analyzed one particular model of quintessence in which $`V(\phi )=\lambda \mathrm{exp}[\lambda \phi ]`$, implying $`\dot{\phi }^2t^2`$. Within this model, in which the quintessence flux around the SMBH decreases with time, they claim that the black holes grow in time. However, the decreasing quintessence flux around these objects casts doubts on this result, since the mass gain term decreases accordingly. We pointed out elsewhere that in the radiation-dominated era the growth of black holes is actually quenched when the background flux decreases with time as $`t^2`$ or faster (see Ref.).

The question is whether there are PBHs that gain substantial mass at asymptotic times for a given potential $`V(\phi )`$ which determines the quintessence flux. Let us show a class of solutions involving quintessence accretion only which can make the PBHs grow, as an example of this general behavior. The quintessence models satisfying

$$\dot{\phi }^2=M_{pl}^4(t/t_{})^n$$
(24)
constitute a class of interesting growing models (we used natural units and $`t_{}=E_{}^{}{}_{}{}^{1}`$ is a time constant). For this choice, the kinetic energy $`\frac{\dot{\phi }^2}{2}=Kt^n`$ with $`K=\frac{M_{pl}^{}{}_{}{}^{4}}{2}E_{}^n`$ and $`E_{}`$ measures directly the kinetic contribution of the field.

Inserting into eq.(21) above, and solving for $`M(t)`$ we obtain, assuming $`n>0`$

$$M(t)=\frac{M_i}{1F(t,t_i)}$$
(25)
where

$$F(t,t_i)=\frac{27\pi \eta t_{}M_i}{2(n+1)}\left[(t/t_{})^{n+1}(t_i/t_{})^{n+1}\right]$$
(26)
It is easy to show that a set of solutions parametrized by the initial masses, time constants $`t_{}`$ and $`n`$ exist where a huge growing of the seed black holes is possible, provided the constraint given by eq.(19) holds all the time. SMBHs will arise (from initially small PBHs) at a final time $`t_f`$ if

$$0<1F(t_f,t_i)\mathrm{\hspace{0.17em}1}$$
(27)
Thus, all those PBHs with initial masses of order

$$M_i(t_fE_{})^{(n+1)}\left[\frac{2(n+1)E_{}}{27\pi \eta }\right]$$
(28)
would end with large masses ($`M10^6M_{}`$ or bigger) at final times $`t_ft_i`$.

Numerically, the mass is

$$M_i\frac{(n+1)}{\eta }\left[10^{(24+\frac{2}{n+1})}\right]^{(n+1)}(t_f/s)^{(n+1)}(E_{}/GeV)^nGeV$$
(29)
where we have absorbed a coefficient $`O(1)`$ into the efficiency $`\eta `$. Initial masses masses may be large only if the scale $`E_{}`$ is extremely small when measured in $`GeV`$ if $`t_f`$ is inside the radiation-dominated era, according to eq.(29).

We may invert the reasoning above and assert that if we had some initial black hole formed with $`M_i`$ at $`t_i`$, then, the constant $`E_{}`$ need to be larger than

$$E_{}10^{\frac{(26+24n)}{n}}\left(\frac{t_f}{s}\right)^{(1+1/n)}\left(\frac{(n+1)}{(\omega _\eta M_i/GeV)}\right)^{1/n}GeV\mathrm{\Theta }_1$$
(30)
(with $`\omega _\eta (27\pi \eta /2)`$) for the black hole to grow to the SMBH regime. Note that when $`E_{}`$ is larger, we need smaller initial masses in order to obtain larger SMBHs at the final time $`t>>t_f`$, as expected.

Eq.(28) also says that our approximations to the actual physical accretion are valid if and only if

$$M_i<\frac{2(n+1)}{27\pi \eta }E_{}(tE_{})^{(n+1)}$$
(31)
and the formulae above stay valid only if the parameter $`E_{}`$ does not change with time. The bottomline of eq.(29) is that if HB+GSL hold for all times, then seed PBHs can not have arbitrary initial masses (independently of the details of their formation) if they had to grow by accreting a quintessential field within the proposed class. Mutatis mutandis the same conclusions hold for other fuels for accretion, where we must use eq.(17) instead), in order to obtain supermassive black holes at the final time $`t_f`$. Note that this constraint gets weaker with time because $`S_{hp}t^2`$, and this constraint will become at some time weaker than the geometric causal condition $`r_g<R_{hp}`$. The HB is quite restrictive for large masses at black hole formation, as discussed in Section 3. Considerations on the accretion before formation must be added to this picture (see next Section).

Finally, according to eq.(19), $`E_{}`$ must also satisfy

$$E_{}6.6\times 10^{2/n}10^{(25+24/n)}\left(\frac{t_f}{s}\right)^{(1+1/n)}\left(\frac{(n+1)}{(M_i/GeV)}\right)^{1/n}GeV\mathrm{\Theta }_2$$
(32)
Therefore, only the PBHs contained within the range defined by $`\mathrm{\Theta }_1(M_i)<E_{}<\mathrm{\Theta }_2(M_i)`$ will satisfy the HB and become SMBHs at late times $`t`$ simultaneously. This leads to the with the following constraint on a positive $`n`$

$$n<0.65\mathrm{log}(27\pi \eta /2)$$
(33)
A careful examination of the $`n<0`$ case leads us to conclude that the lower limit thus obtained is irrelevant when compared to the $`n=1`$ case already discussed, that is, the index is actually limited by $`1`$ from below. The case of a constant flux $`n=0`$ can be also worked out without complications. We conclude that a window of indexes $`n`$ exist for quintessence to cause the growth of seed PBHs to the SMBH regime. Such a window is independent of $`M_i`$. Other physical effects may be important, for example, generally speaking, the depletion of the quintessence flux around the black hole can not be ignored for large masses, an effect that has to affect the parameter $`\eta `$.

Causality and Holographic requirements

As shown in the previous Section, when eqs.(28-30) are satisfied, the energy input by the accretion of quintessence would be enough to drive a black hole with initial mass $`M_i`$ to values $`M_f>10^6M_{}`$.

This energy input strongly depends on the initial mass $`M_i`$ and needs to be very large if the black hole was initially very small. In addition, if the accretion rate is very high, the black hole that was initially below the HB will blow that bound at some point. Then, to keep these black holes below the HB they must also obey the constraints $`\dot{S}_{smbh}(t)<\dot{S}_{hp}(t)`$. However, according to the eq.(10), this inequality is equivalent to

$$\dot{M}(t)<\frac{3.4h_0^2M_{pl}^{}{}_{}{}^{4}}{M(t)}t$$
(34)
We solve the evolution from $`t_i=0`$ until $`t_f>>t_i`$ as above and then

$$M(t_f)<M_i\left(1+(t_f/\tau )^2\right)^{1/2}$$
(35)
For $`t_f\tau =\frac{M_ih_0^1}{2.6M_{pl}^2}2.5\times 10^6s(\frac{M_i}{M_{}})h_{0}^{}{}_{}{}^{1}`$, eq.(35) becomes $`M(t_f)<1.84h_0M_{pl}^2t_f`$.

In order to enforce a strictly causal growth $`\dot{r_g}<c`$, the PBHs must obey also the condition $`\dot{M}<0.5M_{pl}^{}{}_{}{}^{2}`$ at all times. Then, solving for $`M(t_f)`$ we obtain

$$M(t_f)<M_i\left(1+8\times 10^4(t_f/s)(M_{}/M_i)\right)$$
(36)
In other words, any black hole with initial mass satisfying the Holographic Bound at $`t_i`$ will eventually be superholographic at $`t=t_f`$ (that is, $`S_{bh}(t_f)>S_{total}(t_f)`$) unless causality holds.

Since that the causality requirement is very strong, we rule out superholographic black holes in normal circumstances of physical accretion. If we impose that the solution given by eqs.(25-26) must satisfy the Holographic Bound for the maximal rate of energy gain, then combining with eq.(35)

$$\frac{M_i}{1F}<M_i\left(1+(t_f/\tau )^2\right)^{1/2}$$
(37)
Using eq.(37), we may obtain an upper bound for the cosmological time at which our approximations must break down. Defining the dimensionless quantities $`\theta _1=(t_f/\tau )`$, $`\theta _2=(t_f/\zeta t_{})`$ and $`\zeta ^{n+1}=2(n+1)/27\pi \eta `$, eq.(37) becomes

$$\theta _{1}^{}{}_{}{}^{2}(2\theta _{1}^{}{}_{}{}^{2}+1)\theta _{2}^{}{}_{}{}^{n+1}+(\theta _{1}^{}{}_{}{}^{2}+1)\theta _{2}^{}{}_{}{}^{2(n+1)}>0$$
(38)
to be solved for a given $`M_i`$, $`\eta `$ and $`n`$ set of values.

We must acknowledge that, in principle, seed black holes may still reach the SMBH regime without blowing the HB even if the condition $`\dot{S}_{smbh}(t)<\dot{S}_{hp}(t)`$ is not satisfied, provided its growth effectively stopped while still below the HB value. These models, however, must be analyzed in a one-by-one basis to check their viability.

We close this Section with the observation that both the HB requirement and the relativistic bound $`\dot{r_g}<c`$ on $`\dot{M}`$ lead to essentially the same value $`\dot{M}<2\times 10^{38}gs^1`$ within a numerical factor of the order of one.

Conclusions

We have discussed a possible form to limit the formation times of primordial black holes formed directly or grown by accretion which may be residing at the center of most galaxies as recently identified by a variety of observations. The huge entropy contained in these SMBH allows to limit their formation quite efficiently, since the total content of entropy of the universe is likely to be bounded by the HB. Even if preliminar,our analysis of the quintessence models for the growth of seed PBHs has been found to leave room for their formation and further growth, although not for arbitrary fluxes. The general argument developed in Ref.  against fast-growing solutions for PBH growth with radiation flux $`\varrho _{rad}(t)t^2`$ can be directly applied to the the particular model involving quintessence flux $`\dot{\phi }^2t^2`$. Generally speaking, the quintessence model must allow the flux to decrease slower than $`t^2`$ for PBHs to grow at all, and to stay constant or increase for substantial accretion to occur, as needed for achieving the SMBH condition in a short time. General conditions on $`\dot{M}`$ have been obtained by a combination of causal and holographic arguments and are amenable of specific applications.

Other models can be constructed to produce a population of SMBH starting from seed PBHs. For example, accretion in a brane-world high-energy phase has been recently studied  and shown to allow a substantial growth in which $`\dot{M}M/t`$. It may be possible to arrive to the end of the high-energy phase with very massive black holes, although the full consequences of this scenarios are yet to be explored.

Acknowledgments

Both authors wish to thank the São Paulo State Agency FAPESP for financial support through grants and fellowships. J.E.H. has been partially supported by CNPq (Brazil). An anonymous referee is acknowledge for criticisms that helped to improved the original draft.

References

 B.M.Peterson, An introduction to Active Galactic Nuclei, Cambridge, 1997.

 F. De Paolis, G. Ingrosso, A.A. Nucita, D. Orlando, S. Capozziello and G. Iovane, Astron. Astrophys.376, 853 (2002).

 T.Kephart and Y.Jack Ng, gr-qc/0204081 (2002).

 L. Ferrarese, Ap.J.578, 90 (2002) ; J. Silk, Space Sci. Rev.100, 41 (2002).

 D.Bigatti and L.Susskind, hep-th/0002044 (2000) and N.Ohta, gr-qc/0205036 (2002).

 J.Bekenstein, Phys.Rev.D23, 287 (1981).

 E.Verlinde, hep-th/0008140 (2000).

 J.Bekenstein, Phys.Rev.D9, 3292 (1974).

 R. Bean and J.Magueijo, astro-ph/0204486 (2002).

 R. Guedens, D. Clancy and A.R. Liddle, astro-ph/0208299

 J.D.Barrow, E.J. Copeland and A.R. Liddle, Phys. Rev. D46, 645 (1992); see also MNRAS 253, 675 (1991).

 P.Custodio and J.E. Horvath, gr-qc/0203031 (2002).

 P.Custodio and J.E. Horvath, Phys.Rev.D58, 023504 (1998).

Figure captions

Fig. 1. The window of indexes $`n`$ that allow a growth of small PBHs to the SMBH regime. As explained in the text, the upper bound is set by the HB requirement on $`E_{}`$ (eq.32); while the lower bound $`n=1`$ is imposed by the requirement of having enough flux of fuel to complete the process (Ref.12). The possible modes of growth include the constant quintessence flux case $`n=0`$, detailed conditions to be satisfied by these particular power-law models are given in the text.


text:
# The Roton in a Bose-Einstein Condensate

## Abstract

The roton in a Bose-Einstein condensate is computed, both near and far from a Feshbach resonance. A low-density approximation is made, allowing for an analytic result. A Monte Carlo calculation shows that the roton is larger than predicted by the low-density approximation, for the upper range of densities considered here. The low-density approximation is applied to superfluid 4He, roughly reproducing the results of previous Monte Carlo calculations.

A roton, first predicted by Landau landau47  for superfluid 4He, is an excitation in a Bose-Einstein condensed fluid, characterized by a minimum in the excitation spectrum $`\omega (k)`$. By Feynman’s relation feynman54 , this minimum also corresponds to a maximum in the static structure factor $`S(k)`$. This peak in $`S(k)`$ occurs for excitations whose wavelength $`2\pi /k`$ is equal to the characteristic wavelength of density fluctuations in the ground-state wave function of the quantum fluid. The peak in $`S(k)`$ exceeds unity, the value for an uncorrelated gas.

Feynman feynman54  described several possible microscopic pictures of a roton in superfluid 4He. He found that the most likely description is that the roton is analogous to a single atom moving through the condensate, with wave number $`k`$ close to $`2\pi /n^{1/3}`$, where $`n^{1/3}`$ is the mean atomic spacing.

The roton in superfluid 4He was calculated by a Monte Carlo technique mcmillan65 ; schiff67  using a many-body Jastrow wave function jastrow55 , which is of the form

$$\psi =\underset{j>i=1}{\overset{N}{}}f(|\text{r}_i\text{r}_j|),$$
(1)
where the pair function $`f(r)`$ should be determined for the quantum fluid under consideration. The wave function (1) deviates significantly from unity if any two atoms become very close to one another. For superfluid 4He, the two-particle correlation function $`g(r)`$ for this wave function has fluctuations with a preferential length scale of $`n^{1/3}`$, which results in a peak greater than unity in $`S(k)`$ near $`k2\pi /n^{1/3}`$, in agreement with Feynman’s picture.

Both Feynman’s result and the results of the Jastrow wave function roughly agree with measurements henshaw60  of $`S(k)`$ and $`\omega (k)`$.

The Jastrow wave function (1) has also been considered for a low-density Bose-Einstein condensate (BEC) mcmillan65 , and has been used to compute various properties of a high-density BEC cowell02 . We employ a Jastrow wave function to compute the roton for a low-density BEC. By an analytic calculation, we find that the hard-sphere size $`r_o`$ of the atoms determines the location of the roton, rather than $`n^{1/3}`$. We see that Feynman’s view that $`n^{1/3}`$ is the relevant length scale is the special case of high-density, for which $`r_on^{1/3}`$.

In a low-density BEC, Feynman’s view of the roton as a single atom moving through the condensate seems natural. In general, the excitation spectrum of a low-density BEC is of the Bogoliubov form bogo47 ; steinhauer02 , which consists of phonons and single-particle excitations. We find that the roton occurs on the single-particle part of the spectrum, at $`k8/\pi a`$, where $`a`$ is the $`s`$-wave scattering length.

Although the general form of the wave function (1) can be applied to both superfluid 4He and BEC, the wave function is fundamentally different for these two quantum fluids. For both of these fluids, below some temperature $`T_{swave}`$, the thermal de-Broglie wavelength $`\lambda _{dB}`$ becomes longer than the characteristic length scale $`R`$ of the interparticle potential. Below $`T_{swave}`$, all scattering processes except for $`s`$-wave scattering become negligible. Since a BEC is a dilute gas, $`n^{1/3}`$, which is typically about $`1500\AA `$, is much larger than $`R`$. Therefore, the critical temperature $`T_c`$ for quantum degeneracy, at which $`\lambda _{dB}`$ becomes comparable to $`n^{1/3}`$, is much lower than $`T_{swave}`$. Therefore, only $`s`$-wave scattering plays a role here for a BEC. In contrast to BEC, superfluid 4He is a relatively dense liquid, for which $`n^{1/3}`$ is comparable to $`R`$. Therefore, the temperature $`T_{swave}`$ also marks the transition $`T_c`$ to quantum degeneracy. Even below $`T_c`$ for superfluid 4He, partial waves other than the $`s`$-wave contribute to the wave function.

For many properties of a BEC, $`a`$ is sufficient to describe the interparticle potential, and the form of the potential does not play a role huang ; dalfovo99 . To quantitatively describe the roton however, the $`s`$-wave scattering wave function must be known to interparticle separations somewhat smaller than $`R`$, where the details of the potential are relevant.

Rotons have been predicted in a BEC in the presence of optical fields stamperkurn02 ; odell03 , as well as in a dipolar BEC shlyapnikov03 . In this work, we find that a roton occurs in an unperturbed BEC without any dipole interaction. We consider the enhancement of the roton near a Feshbach resonance feshbach ; inouye , but the physics is qualitatively the same as in the unperturbed case.

When we refer to a roton, we are referring to a peak in $`S(k)`$, rather than a minimum in $`\omega (k)`$. By the Feynman relation, the peak in $`S(k)`$ computed here for a BEC is not steep enough to produce a minimum in $`\omega (k)`$, for the range of densities considered here.

We use a low-density approximation jastrow55  to compute $`g(r)`$ and the roton for a BEC with a positive scattering length. The low-density approximation is compared to a Monte Carlo calculation.

For a BEC of alkali atoms, the potential can be taken as the van der Waals potential $`C_6/r^6`$ for atomic separations $`r`$ greater than a few Angstroms flambaum93 ; leggett01 . For this potential, $`R`$ is given by $`(mC_6/\mathrm{}^2)^{1/4}`$, which ranges from about $`50`$ to $`100\AA `$ leggett01 ; marinescu94 . The wave function for $`s`$-wave scattering between two particles in the limit of zero energy is given by $`\psi _6(r)`$, which is the solution of the radial Schrödinger equation with the van der Walls potential. $`\psi _6(r)`$ can be written in terms of elliptic integrals flambaum93 . For $`r>R`$, this wave function is very close in form to $`1a/r`$ huang ; comment1 . For the case of $`aR`$, $`\psi _6(r)`$ therefore obtains large negative values for $`ra`$ dalibard .

We cannot use $`\psi _6(r)`$ as $`f(r)`$ in (1), because even for $`rn^{1/3}`$, $`\psi _6(r)`$ is significantly less than unity. This is non-physical, in the sense that no matter how large the gas for fixed density, the value of the many-body wave function (1) depends on the size of the gas. To account for many-body effects, we use the following pair function which goes to unity for $`r>n^{1/3}`$ leggett01 ; cowell02 .

$$f(r)=\{\begin{array}{cc}\psi _6(r)/\psi _6(n^{1/3})\hfill & (rn^{1/3})\hfill \\ 1\hfill & (r>n^{1/3}).\hfill \end{array}$$
(2)
To compute (2), $`R`$ is needed. Throughout this work, we use $`R=0.05`$ $`n^{1/3}`$. This is a typical experimental value, and the results here are rather insensitive to $`R`$.

Equation (1) with (2) is shown schematically in Fig. 1. To aid in visualization, Fig. 1 shows the wave function squared in one dimension, as a function of the position $`x_1`$ of atom number $`1`$. The positions of all of the other atoms, such as atoms $`b`$ through $`e`$, are fixed. As long as atom number $`1`$ is far from the other atoms, the wave function has the constant value $`\psi _o`$. This value is determined by the positions of the atoms other than $`x_1`$. The rapid oscillations in $`\psi _6(r)`$ for $`r<R`$ appear in Fig. 1 as dark vertical bands, when atom number 1 is very close to another atom.

The correlation function $`g(r)`$ gives the unconditional probability of two atoms being at a distance $`r`$. $`g(r)`$ is related to the pair function $`f(r)`$ by jastrow55

$$g(|\text{r}_1\text{r}_2|)=V^2\frac{𝑑\text{r}_3\mathrm{}𝑑\text{r}_N\underset{j>i=1}{\overset{N}{}}f^2(|\text{r}_i\text{r}_j|)}{𝑑\text{r}_1\mathrm{}𝑑\text{r}_N\underset{j>i=1}{\overset{N}{}}f^2(|\text{r}_i\text{r}_j|)},$$
(3)
where the integrals are over the volume $`V`$. (3) can be evaluated starting with the first integral in the denominator, $`𝑑\text{r}_1\underset{j>1}{\overset{N}{}}f^2(|\text{r}_1\text{r}_j|)`$. This can be visualized as the integral over the function shown in Fig. 1 for one dimension. Neglecting three-body interactions greatly simplifies this integral. Three-body interactions are rare for small values of the gas parameter $`na^3`$, assuming that the range of three-body interactions is of the same order of magnitude as the range of two-body interactions. A three-body interaction is represented in Fig. 1 by the points $`x_c`$ and $`x_d`$, where atoms $`1`$, $`c`$, and $`d`$ interact. Neglecting such interactions, the integral is a function of the volume $`v`$ indicated by the shaded region in Fig. 1. The integral is then given by $`V(N1)v`$, where $`v=V𝑑\text{r}f^2(r)`$. This result is independent of the positions of the $`\text{r}_j`$.

Evaluating all of the integrals in (3) similarly to the first yields $`g(|\text{r}_1\text{r}_2|)V^2f^2(|\text{r}_1\text{r}_2|)[(Vv)V]^1`$. Since $`Vv`$, we obtain the result of the low-density approximation jastrow55

$$g(r)f^2(r).$$
(4)
The result (4) with (2) is indicated by the solid curve in Fig. 2 for $`na^3=2.2\times 10^4`$. This value of $`na^3`$ is an order of magnitude greater than typical experimental values without a Feshbach resonance. The result for $`na^3=0.011`$ is also shown in the figure.

While (4) is useful for obtaining analytic results, a more accurate computation can be made by the Monte Carlo technique described in Ref. mcmillan65 . This technique effectively evaluates (3) by using a Metropolis algorithm to randomly choose configurations of the $`N`$ atoms ($`N=100`$ here), according to the probability distribution given by (1) with (2), and computing the distribution of distances between the atoms, with periodic boundary conditions. This distribution, averaged over many likely configurations, is proportional to $`r^2g(r)`$. The result of the Monte Carlo computation is shown in Fig. 2 for $`na^3=2.2\times 10^4`$ and 0.011. These results were obtained with $`9\times 10^8`$ and $`2\times 10^7`$ iterations, respectively.

The small oscillations for small $`r`$ shown in Fig. 2 are negligible in the computation of the static structure factor. These are the oscillations for r values less than that of the first large peak below $`r/a=1`$, in the solid and dashed curves. To save computing time these oscillations are not included in $`f^2(r)`$ in the Monte Carlo computation.

For $`a1.3R`$ (corresponding to $`na^3>2.7\times 10^4`$ in this work), The peak in $`f^2(r)`$ located at $`r<R`$ is greater than unity. The dashed line of the inset of Fig. 2 shows an example of such a peak, for $`a=4.4R`$. Because this peak is greater than unity, a cluster of atoms is the most likely configuration. Since we are interested in the un-clustered phase, we do not let the Monte Carlo computation proceed long enough for the clear transition to clusters to occur.

The static structure factor $`S(k)`$ is given by $`1+n[g(r)1]e^{i\text{k}\text{r}}𝑑\text{r}`$, which in general can be written

$$S(k)=1+4\pi n\underset{0}{\overset{\mathrm{}}{}}𝑑rr^2[g(r)1]\frac{sin(kr)}{kr}.$$
(5)
We use (5) to compute $`S(k)`$ for the low-density approximation (4) with (2), as indicated in Fig. 3. The height $`S(k_r)`$ and location $`k_r`$ of the roton are indicated by the solid curves of Fig. 4. $`S(k)`$ in Fig. 3 has several maxima, the tallest of which is taken as the roton. As $`na^3`$ is varied, the peak which is taken as the roton varies, resulting in the jagged appearance of the solid curves of Fig. 4.

For the Monte Carlo computation, $`S(k)`$ is found by inserting $`g(r)`$ such as is shown in Fig. 2 into (5). The results are indicated in Fig. 3 by the dash dotted and solid curves. For small $`na^3`$, the height and location of the roton are seen in Fig. 4 to be the same for the low-density approximation and the Monte Carlo result. The low-density approximation is therefore valid for small $`na^3`$. For the larger values of $`na^3`$, the height of the roton in the Monte Carlo calculation is larger than that of the low-density approximation, as seen in Fig. 4a.

By applying the low-density approximation to the limit of $`aR`$, we can obtain an analytic expression for $`S(k)`$ for $`Ran^{1/3}`$. In this range, $`f(r)`$ can be taken as $`1a/r`$. By (4) and (5),

$$S(k)=1+4\pi na^3[\pi (ka)^1/22(ka)^2].$$
(6)
This small-R limit is indicated by the dotted curve in Fig. 3. The height of the roton in this limit is $`S(k_r)=1+\pi ^3na^3/8`$, and the location is given by $`k_r=8/\pi a`$. These values are indicated by the dashed curves of Fig. 4, which as expected, agree with the low-density approximation (solid curves) for $`aR`$, where $`nR^3`$ is indicated by the dotted lines. As seen in the above expression for $`k_r`$, the location of the roton is determined by $`a`$, rather than by $`n^{1/3}`$.

For $`na^3=0.011`$, the height of the roton given by the Monte Carlo calculation is $`S(k_r)=1.08`$, as shown in Figs. 3 and 4a. Measuring this $`8\%`$ effect could be experimentally feasible. $`na^3=0.011`$ could be attained by a Feshbach resonance comment1 . It should be noted though, that for this relatively large value of $`na^3`$, the form of (2) is only approximate, so the expressions for the height and location of the roton should be considered as estimates only.

Due to phonons, the true $`S(k)`$ is proportional to $`k`$ for $`k\xi ^1`$, where $`\xi ^1=a^1\sqrt{8\pi na^3}`$ is the inverse healing length. The curves shown in Fig. 3 do not show this linear behavior for small $`k`$ because the wave function (1) with (2) does not have the long-range correlations of a phonon mcmillan65 ; schiff67 .

For superfluid 4He the three inverse length scales $`2\pi /n^{1/3}`$, $`k_r`$, and $`\xi ^1`$, are roughly equal. This is also true for a strong roton in a BEC. As $`na^3`$ increases, both $`k_r`$ and $`\xi ^1`$ approach $`2\pi /n^{1/3}`$. More precisely, the ratios of $`\xi ^1`$ and $`k_r`$ to $`2\pi /n^{1/3}`$ are $`\sqrt{2/\pi }(na^3)^{1/6}`$ and approximately $`4/\pi ^2(na^3)^{1/3}`$, respectively. The latter ratio implies that an appropriate measurement system for measuring a roton should be able to probe wavelengths somewhat shorter than $`n^{1/3}`$.

As a demonstration of the validity of the low-density approximation (4), we use (4) to compute various properties of superfluid 4He found in mcmillan65  and schiff67  by means of Monte Carlo calculations. For superfluid 4He, $`nb^30.4`$ where $`b`$ is the hard-sphere size. Therefore, this is a stringent test of the low-density approximation. We find that the location of the roton for superfluid 4He is given by $`k_r=5.3/b`$, which differs from the Monte Carlo result by 6% mcmillan65 ; schiff67 . We find that the height of the roton for superfluid 4He is $`S(k_r)1+0.3nb^3`$, which gives a roton of $`S(k_r)=1.1`$, compared to $`S(k_r)=1.2`$ from the Monte Carlo calculation mcmillan65 ; schiff67 , and $`S(k_r)=1.5`$ from experiment henshaw60 . The order-of-magnitude agreement between the low-density approximation and the Monte Carlo result for superfluid 4He suggests that the low-density approximation preserves the essence of (3).

In conclusion, we find the height $`S(k_r)`$ and location $`k_r`$ of a roton in a BEC, for a range of densities. A low-density approximation is compared to a Monte Carlo calculation. The values of $`S(k_r)`$ and $`k_r`$ given by the two methods agree for the lowest densities. For higher densities, the Monte Carlo calculation predicts an enhancement in $`S(k_r)`$ of almost a factor of 2 over the low-density approximation.

In contrast to the Monte Carlo calculation for superfluid 4He, the small-R limit gives explicit expressions for the height and location of the roton.

We thank Servaas Kokkelmans, Ady Stern, Yoseph Imry, Daniel Kandel, Johnny Vogels, and Ananth Chikkatur for helpful discussions. This work was supported by the Israel Science Foundation. J. S. is a Landau Fellow, supported by the Taub and Shalom Foundations.


text:
# Time and energy-resolved two photon-photoemission of the Cu(100) and Cu(111) metal surfaces

## I Introduction

The presence of a metal surface creates electron states that do not exist in the bulk metal. For example, an electron situated at a distance $`z`$ from the metal surface experiments an attractive force, $`F(z)=e^2/(2z)^2`$, equivalent to that produced by its image charge situated at distance $`z`$ inside the metal. For large $`z`$, the potential generated by this surface induced charge approaches the classical image potential, $`V(z)=e^2/4z`$111This expression corresponds to assuming an infinite repulsive barrier at the surface. In the more general case we have $`V(z)=\frac{e^2}{4(zz_{\mathrm{im}})}\frac{ϵ1}{ϵ+1}`$, where $`z_{\mathrm{im}}`$ is the position of the image plane (that depends on the electronic and structural properties of the surface), and $`ϵ`$ is the static dielectric constant of the surface. . When the metal has a surface band gap near the vacuum level, the electrons below the vacuum level are trapped between the well of the image potential and the surface barrier. These quantized states correspond to image potential states, and are spatially localized in the region in front of the surface (see Fig. 1).

Metal image states typically have a large spatial extension, but barely penetrate the bulk metal. They form a Rydberg-like series with energies, $`E_n`$, approximately given byPedro :

$$E_n=\frac{0.85}{(n+a)^2}[\mathrm{eV}],$$
(1)
where $`n`$ is a positive integer, and $`0a0.5`$ represents the quantum defect (which depends on the position and width of the gap). The energy of an electron trapped in an image state and moving parallel to the surface with momentum $`𝐤_{}`$ is therefore:

$$E=\frac{\mathrm{}^2k_{}^2}{2m_e}+E_n$$
(2)
In contrast to bulk states, the small overlap between surface states and bulk states reduces considerably the inelastic electron scattering, leading to a long lifetime, $`\tau _n`$, of the image states, that scales asymptotically with the quantum number $`n`$ as $`\tau _nn^3`$Chulkov2 . Due to this fact, these states are very interesting for the study of electron correlations. Furthermore, image states play an important role in the laser induced chemical control of reactivity at metal surfaces. It is therefore not surprising the considerable number of studies focusing on this subject (see Ref. Chemphys  and references within).

Image potential states of different surfaces have been observed experimentallyexperiment . A powerful technique to measure their energy and lifetime is two-photon photoemission spectroscopy (2PPE)Giesen . Recently, and with the advent of ultrafast laser technology it became possible to perform time-resolved 2PPE spectroscopy (TR2PPE), which allowed, with the technique of quantum beat spectroscopyHofer ; Klamroth , the direct measurement of the lifetime of image states, even for states with large quantum number $`n`$. The 2PPE technique is depicted schematically in Fig. 1. One photon (pump), of energy $`\mathrm{}\omega _a`$, excites an electron from an occupied state below the Fermi energy ($`E_\mathrm{F}`$) to an image-potential state of quantum number $`n`$; Then a second photon (probe), with energy $`\mathrm{}\omega _b`$, ejects the electron out of the surface, above the vacuum energy. This electron has a kinetic energy $`E_k=\mathrm{}\omega _bE_n`$, that is measured in a detector far away from the surface. By varying the delay between the pump and probe pulses, the intensity of the 2PPE signal as a function of the delay will reflect the evolution of the population of the state $`n`$. From this information it is then possible to extract the lifetime of the state $`n`$.

Lifetimes of image states for different metal surfaces have been calculated in the framework of the self-energy formalism in the GW approximationOsma ; Chemphys . The results obtained happen to be in rather good agreement with experiments. For example, the calculated lifetime of the first three image potential states of the Cu(100) surface are $`\tau _1^{\mathrm{theor}}=30`$ fs, $`\tau _2^{\mathrm{theor}}=132`$ fs, and $`\tau _3^{\mathrm{theor}}=367`$ fsChemphys , while the experimental values are $`\tau _1^{\mathrm{exp}}=40\pm 6`$ fs, $`\tau _2^{exp}=110\pm 10`$ fs, and $`\tau _3^{exp}=300\pm 15`$ fsHofer , while for the first image state of Cu(111) surface the theoretical lifetime is $`\tau _1^{\mathrm{theor}}=17.5`$ fsChemphys  and the experimental is $`\tau _1^{\mathrm{exp}}=15\pm 5`$ fsMcNeil .

In the present paper we present calculations on energy-resolved 2PPE spectra of Cu(100) and Cu(111) surfaces. Our theoretical approach is based on the propagation of an electron wave-packet in an one-dimensional model potentialChulkov2  simulating the Copper surface. The finite electronic lifetimes of the image potential states are taken into account by an empirical self-energy term. The photoemission spectra are obtained through a geometrical scheme described in Ref. 2bepub . Note that in this approach the potential is fixed during the whole simulations, i.e., the change of the electronic screening during the excitation is not taken into account. This could be important for short laser excitations as done in TR2PPE.

To assert the relevance of this approximation we also performed time-dependent simulations of the Copper surface, but including the change of the Hartree potential due to the electronic excitations.

## II Method

To model the Copper surface we used a one-dimensional slab of 48 Copper layers, surrounded by 290 a.u. of vacuum on each side. The large portion of vacuum is necessary in order to describe the four first image potential wave functions (see Fig. 1). The bulk and image potentials are modeled by the 1D potential model of Ref. Chulkov2 . This model reproduces the position and the width of the energy gap, as well as the energy of the surface state and of the first image potential state. Furthermore, it provides a good description of the electronic structure of simple and noble metal surfacesChulkov2 . In this model the number of bulk states is discrete for each value of $`q_{}`$, and is related to the number of layers in the simulation. We do not believe that this is a major limitation, essentially because all the relevant physics occurs close to $`E_\mathrm{F}`$, while the remaining states mainly act as a polarizable background. The electrons are allowed to move and interact only in the z-direction, perpendicular to the surface, and we assume a parabolic dispersion in the parallel direction. Within this approximation, the Hamiltonian describing such system can be written as

$$\widehat{H}_0=\frac{d^2}{dz^2}+\widehat{v}_{\mathrm{model}}.$$
(3)
This operator is discretized in real-space using a grid-spacing of 0.2 a.u, which is sufficiently small to allow for a proper description of the relevant electronic states.

At $`t=0`$ we assume the system to be in the ground-state of the Hamiltonian (3). We then propagate this state with the new Hamiltonian

$$\widehat{H}(t)=\widehat{H}_0+\widehat{v}_{\mathrm{laser}}(z,t)+\widehat{\mathrm{\Sigma }}(z),$$
(4)
where $`\widehat{v}_{\mathrm{laser}}(z,t)`$ describes a laser polarized in the $`z`$ direction that reads, in the length gauge, $`\widehat{v}_{\mathrm{laser}}(z,t)=zE(t)`$. The electric field is composed of a pump and a probe pulses,

$$E(t)=E_{\mathrm{pump}}(t)\mathrm{cos}(\omega _{\mathrm{pump}}t)+E_{\mathrm{probe}}(t)\mathrm{cos}(\omega _{\mathrm{probe}}t).$$
(5)
the functions $`E_{\mathrm{pump}/\mathrm{probe}}(t)`$ are envelope functions of the form

$$E_i(t)=E_i\mathrm{cos}^2\left[\frac{\pi }{\sigma _i}(tt_i)\right]\mathrm{\Theta }(tt_i+\sigma _i/2)\mathrm{\Theta }(t+t_i+\sigma _i/2),$$
(6)
where $`\sigma _i`$ is the width of the pulse and $`t_i`$ the center of the pulse. For the pump pulse, $`t_{\mathrm{pump}}`$ is simply $`\sigma _{\mathrm{pump}}/2`$, while for the probe pulse $`t_{\mathrm{probe}}=t_\mathrm{D}+\sigma _{\mathrm{pump}}/2`$, with $`t_\mathrm{D}`$ the delay time between the pump and probe lasers. In order to take into account the finite lifetime of the image potential states we add to the Hamiltonian (3) a time-independent self-energy operator of the form:

$$\mathrm{\Sigma }(z)=\mathrm{i}\underset{k=1}{\overset{N}{}}\mathrm{\Gamma }_k|n_kn_k|$$
(7)
where $`|n_i`$ are the image potential states and $`\mathrm{\Gamma }_k`$ are the inverse experimental lifetimes of each state taken from Ref.Chemphys . Note that $`\widehat{H}(0)=\widehat{H}_0`$. The time-propagation is performed in real-time following the method of Ref.octopus . This propagation scheme preserves unitarity for a hermitian Hamiltonian, and has been proven very robust and stable in diverse applications.

## III Energy-resolved spectra

In order to obtain the energy-resolved photo-electron spectra we follow the technique of Ref.2bepub  and divide the simulation box in two regions: one region, A, containing the slab and the region where the first image potential states lie, and another, B, defined as the complement of A (see Fig. 2). In region B the electrons are considered free outgoing particles, and are treated in momentum space. The separation between the two regions is achieved through a smooth masking function, $`M`$, defined as one in the interior of region A and zero outside. The method consists in evolving the wave-packet in the first region with the Hamiltonian (4), and mask the orbitals at each time-step. The electrons that “leave” region A are then treated as free-particles and accumulated in momentum space. In mathematical terms, the evolution is performed in the following way: We start by propagating the wave-functions in regions A and B using

$`\psi _k^\mathrm{A}(z,t+\mathrm{\Delta }t)`$ $`=`$ $`M(z)\mathrm{exp}\left(\mathrm{i}\widehat{H}\mathrm{\Delta }t\right)\psi _k^\mathrm{A}(z,t)`$ (8)
$`\psi _k^\mathrm{B}(p,t+\mathrm{\Delta }t)`$ $`=`$ $`\mathrm{exp}\left\{\mathrm{i}{\displaystyle \frac{[pA(t)]^2}{2}}\mathrm{\Delta }t\right\}\psi _k^\mathrm{B}(p,t)+\stackrel{~}{\psi }_k^\mathrm{A}(p,t+\mathrm{\Delta }t),`$
where $`\stackrel{~}{\psi }_k^A(p,t+\mathrm{\Delta }t)`$ is the Fourier transform of the part of the wave function $`\psi _k^A`$ that left region $`A`$ during the time step $`\mathrm{\Delta }t`$, i.e,

$$\stackrel{~}{\psi }_k^\mathrm{A}(p,t+\mathrm{\Delta }t)=𝑑z\mathrm{exp}(\mathrm{i}pz)[1M(z)]\mathrm{exp}(\mathrm{i}\widehat{H}\mathrm{\Delta }t)\psi _k^\mathrm{A}(z,t),$$
(9)
where $`p`$ denotes the momentum. If we wait long enough after that the laser has been turned off, region B will contain those electrons that were ionized. The photo-electron spectra is then identified with:

$$P(\sqrt{2mE})=\underset{k=1}{\overset{N}{}}|\psi _k^\mathrm{B}(p,t\mathrm{})|^2$$
(10)
where $`N`$ is the total number of electrons. This method can be derived from the interpretation of the Wigner transform of the one-body density matrix as a probability density2bepub .

In Fig. 3 we depict the energy resolved spectra of a Cu(100) and a Cu(111) surfaces for both zero and finite pump-probe delay. For the (100) surface, the binding energies of the first image potential states are -0.57, -0.18, and -0.08 eV (see Fig. 1), while for the (111) surface the first states are at -0.82, -0.22 and -0.009 eV. Note that the second and third image states of the Cu(111) surface are resonances. Furthermore, the laser parameters used in these calculations were, for the (100) surface, $`\mathrm{}\omega _{\mathrm{pump}}=4.7`$ eV, $`\mathrm{}\omega _{\mathrm{probe}}=1.8`$ eV, $`\sigma _{\mathrm{pump}}=95`$ fs, and $`\sigma _{\mathrm{probe}}=54`$ fs, and for the (111) surface, $`\mathrm{}\omega _{\mathrm{pump}}=5.1`$ eV, $`\mathrm{}\omega _{\mathrm{probe}}=1.8`$ eV, $`\sigma _{\mathrm{pump}}=87`$ fs, and $`\sigma _{\mathrm{probe}}=50`$ fs. From these results we observe that the peaks corresponding to the first image potential state are obtained at the correct energy position. The signal due to the image states of higher quantum number cannot be distinguished due to the small difference in energy that cannot be resolved with our grid resolution. We also note that, for zero delay, the broadening of the peak relative to the first image potential state is larger in the case of the Cu(100) surface than in the Cu(111). This is due to the fact that in the former case, the probe pulse can ’pump’ low lying energy electrons from the bulk into energy levels near the gap which can be then ejected from the surface by the pump pulse. As the Cu(111) surface has a different electronic structure, such situation does not occurbands . Note that Cu(100) shows an intrinsic surface state resonance, while Cu(111) has a surface state below the Fermi level. In the last case, the dominant transition happens between the occupied surface state and the first image potential state. This can be rationalized in terms of the larger spacial overlap between surface and image states as compared to bulk states.

## IV Time-resolved spectra

The typical energy width of the laser pulses used in TR-2PPE varies between about 10 and 30 meV. With this energy resolution it is possible to excite separately each one of the first three image potential states of Cu(100). By following the TR-2PPE intensity in function of the pump-probe delay it is then possible to follow the evolution of the population of these states. In Fig.4 we present calculations of the populations of the first three image states of Cu(100). The curves were obtained by propagating the system with the Hamiltonian (4). For comparison we also show the experimental data from Ref. Hofer . Clearly, there is a good qualitative agreement between the theoretical and experimental curves.

In these calculations, the electron-electron interaction was taken into account by the model potential of Ref. Chulkov2 . However, the variation of this term with time was mostly neglected. Part of the time-dependence is taken indirectly into account through the self-energy operator, but one can question how reasonable this approximation really is. In order to answer to this concern, we performed simulations in which we allowed the classical part of the electron-electron interaction (i.e. the Hartree term) to change with time. This was achieved by adding to the Hamiltonian (4) a new term, $`\delta v_{\mathrm{Hartree}}(z,t)`$, that corresponds to the variation of the Hartree potential due to the change of electronic density in the direction perpendicular to the surface,

$$\frac{^2\delta v_{\mathrm{Hartree}}(z,t)}{z^2}=4\pi \delta \rho (z,t).$$
(11)
where $`\delta \rho (z,t)=\rho (z,t)\rho (z,0)`$. Exchange and correlation effects were treated as before, i.e. through the self-energy term (7) used to simulate the lifetime of the states.

The evolution of the population of the first image state with and without this term is shown in Fig 5. We observe only a minor shift of the peak to lower times from 24fs to 21.3fs and the lifetime is reduced by around 2% when the term of Eq.(11) is included. This result validates previous studies that have neglected this effectOsma .

## V Conclusions

We have performed dynamical simulations in order to calculate energy-resolved photoemission spectra of the Cu(100) and Cu(111) metal surfaces. Using the same technique we also obtained time-resolved photoemission spectra for individual excitations of the first three image potential states of Cu(100). This technique is quite general, and can be used to gain insight on the dynamics of image states with large quantum number and when more than one eigenstate is excited coherently (the situation in quantum beat spectroscopy). Furthermore, we studied the influence of the variation of the Hartree potential during the excitation. This effects turns out to be small (around 2%).

## Acknowledgments

This work was supported by the Basque Country University, DGES and European Community under the research training network NANOPHASE (HPRN-CT-2000-00167). We wish to thank H. Appel, E. K. U. Gross and P. M. Echenique, for enlightening discussions and to the high-performance computing facilities of CEPBA where most of the calculations were done.


text:
# Complex behavior from a simple physical system: A numerical investigation

## I Introduction

For many centuries now, students of introductory physics have studied Newton’s Laws of Motion. The two systems done ad nausea are a particle in free-fall and the mass-spring system (inclines being a special case of the former). Additionally, it is well know in the chaos community, over the past few decades, that coupled oscillators can exhibit chaotic behavior i.e. motion that is very sensitive to initial conditions.tl75  Examples of this kind include the driven damped oscillator, the double pendulum, and the magnetic pendulum. A few of these are illustrated in Marion and Thornton’s text on classical dynamics.mt95

The purpose of this paper is to answer the question, “What is the behavior of a mass-spring system which is subject to free-fall motion above a rigid floor without damping?”. This system is a variant of a class of problems know as impact oscillators. One of the earliest was put forward by Fermi to explain the origin of cosmic rays.ef49  His model was an elastic ball which bounced off of a vertically oscillating floor. The interested reader is referred to the bibliography for descriptions of similar systems.ws96 ; ld98

The free-fall spring subject to impact (FSI) constitutes a coupled oscillator problem and the resulting motion is chaotic.

## II The Physical System: Coupled Oscillators

Consider the motion of a point particle released from rest at an initial height $`y`$ subject to the acceleration due to gravity $`g`$. If the particle suffers a perfectly elastic collision with the floor (as is done in billiard theory) it will return to height $`y`$ after a time $`T_{freefall}`$ given by

$$T_{freefall}=\sqrt{\frac{8y}{g}}$$
(1)
If the system is allowed to continually bounce, then it behaves like an oscillator with frequency given by

$$\nu _{freefall}=\frac{1}{T_{freefall}}$$
(2)
Next, consider two equal masses at either end of a perfect spring with an equilibrium position of $`y_0`$ and Hooke’s constant $`k`$. The motion of the difference coordinate i.e. $`Y=y_2y_1`$ will be SHM with frequency

$$\nu _{SHM}=\frac{1}{2\pi }\sqrt{\frac{2k}{m}}$$
(3)
The motion of the two systems can be coupled by releasing the spring system vertically and having the lower mass reflect off of a perfectly elastic surface. The system is illustrated in figure 1 and represents a coupled oscillator. Examples in nature of such a system would be a diatomic molecule of gas which suffers collisions with the floor of its container or a child on a pogo stick.

## III Simulation: The Route from Physics to Computer

In this section, the steps to get from physical model to a computer simulation are outlined. Dimensionless coordinates are also introduced.

### III.1 Equations of Motion:Model to Theory

Newton’s second law is used to write down the equations of motion. Taking the rigid floor to be at $`y=0`$, $`m_1=m_2=m`$, and the rest length of the spring to be $`y_0`$ the accelerations of the respective particles are

$`\ddot{y}_1`$ $`=`$ $`g+\omega _{0}^{}{}_{}{}^{2}(y_2y_1y_0)`$ (4a)
$`\ddot{y}_2`$ $`=`$ $`g\omega _{0}^{}{}_{}{}^{2}(y_2y_1y_0)`$ (4b)
The conventional dot notation has been used to represent the temporal derivatives and $`\omega _{0}^{}{}_{}{}^{2}=k/m`$ . To simplify the above equations 4, dimensionless coordinates are used

$`{\displaystyle \frac{y_i}{y_0}}y_i`$ (5a)
$`t\omega _0t`$ (5b)
$`\beta ={\displaystyle \frac{g}{y_0\omega _{0}^{}{}_{}{}^{2}}}`$ (5c)
The dimensionless quantity $`\beta `$ contains all of the physical parameters of the problem.

The two second order equations can be rewritten into a system of four first order equations by noting that velocity is the time derivative of position. In this way, and using system 5, the dimensionless equations of motion for the FSI are

$`\dot{y}_1`$ $`=`$ $`v_1`$ (6a)
$`\dot{v}_1`$ $`=`$ $`\beta +(y_2y_11)`$ (6b)
$`\dot{y}_2`$ $`=`$ $`v_2`$ (6c)
$`\dot{v}_2`$ $`=`$ $`\beta (y_2y_11)`$ (6d)
The above system of equations ignores the interaction of the lower mass $`y_1`$ with the rigid floor. The lower mass receives a delta function impulse which switches its velocity $`v_1v_1`$ upon impact. Care must be taken on handling this impact and will be discussed in the numerics section.

### III.2 Energy Bounds

There are energy bounds that must be taken into consideration for the free-fall spring subject to impact(FSI).First, the energy is written down in dimensionless form

$$E=\frac{1}{2\beta }(v_{1}^{}{}_{}{}^{2}+v_{2}^{}{}_{}{}^{2})+\frac{1}{2\beta }(y_2y_11)^2+y_2+y_1$$
(7)
The terms are in the order of kinetic, spring potential, and gravitational potential energy. The minimum of energy is chosen to be the system at rest with the lower mass on the floor at $`y_1=0`$ and the upper mass to be at $`y_2=1\beta `$ which corresponds to

$$E_{min}=1\frac{\beta }{2}$$
(8)
Since negative energy is unphysical in this problem, this relation requires that $`\beta <2`$. Moreover, $`\beta `$ must be chosen so that $`y_2>0`$ i.e. $`\beta <1`$. By choosing the spring to start off at its equilibrium length, the minimum of energy is $`E_{min}=1`$.

The upper bound on energy is found by requiring that the spring not be compressed to zero length i.e. $`y_1=y_2=0`$. This is to ensure that particle two does not suffer a collision with the floor and that the two particles don’t pass through each other. Substituting these initial conditions into 7 gives

$$E_{max}=\frac{1}{2\beta }$$
(9)
The energy and parameter requirements for the FSI are

$`1<E<{\displaystyle \frac{1}{2\beta }}`$ (10a)
$`\beta <{\displaystyle \frac{1}{2}}`$ (10b)
### III.3 Theory to Numerics

The system was integrated in time with a fourth order Runge-Kutta method. The time step chosen was $`h=0.001`$ which results in a local truncation error of approximately $`𝒪(h^5)=10^{15}`$. The energy of the system was monitored and noted to be conserved to eight decimal places in double precision arithmetic. For all the simulations $`\beta `$ was chosen to be $`\beta =0.25`$ making the energy range $`1<E<2`$.

The subtle part of the simulation is the actual impact of the lower mass $`y_1`$ with the floor at $`y=0`$. The time of impact can not be predetermined with a static discrete time step. As an example, let $`y_1(t_i)`$ equal the location of the lower mass right before impact at time $`t_i`$. At the next time step $`y_1(t_i+h)`$ , the lower mass has moved through the floor and $`y_1(t_i+h)<0`$ which is unphysical and a consequence of using discrete time steps. The situation can be remedied by switching the independent variable from time to $`y_1`$ and integrating the system back to the surface. Applying the chain rule from calculus, the system is transformed to

$`{\displaystyle \frac{dy_1}{dy_1}}=1`$ (11a)
$`{\displaystyle \frac{dv_1}{dy_1}}={\displaystyle \frac{\dot{v_1}}{v_1}}`$ (11b)
$`{\displaystyle \frac{dy_2}{dy_1}}={\displaystyle \frac{\dot{y_2}}{v_1}}`$ (11c)
$`{\displaystyle \frac{dv_2}{dy_1}}={\displaystyle \frac{\dot{v_2}}{v_1}}`$ (11d)
and integrated with h replaced by $`hy_1`$. The system , in this way, is moved back by the amount of the overshoot. At this point, the velocity of the lower particle is switched $`v_1v_1`$ and the system is integrated forward again with the original time step. A similar example can be found in the references as well as a description of Runge-Kutta routines.km90

## IV The Indicators of Chaos

In this section, the system is investigated for the existence of chaos and two popular chaos indicators are described.

### IV.1 Sensitivity to Initial Conditions: Lyapunov Exponents

Chaos has been described as a system which exhibits strong sensitivity to initial conditions. The famous “butterfly effect” of Edward Lorenz is an example.el63  In figure 2.a two time series of $`y_2`$ are plotted on the same graph. Trajectory 1 is at an energy of $`E=1.1`$ and trajectory two at $`E=1.1+\delta E`$ where $`\delta E=0.001`$. The small separation in the two trajectories is slightly noticeable and the two paths follow each other closely. Also, the trajectory looks periodic. When the energy is increased to $`E=1.75`$ and $`E=1.75+\delta E`$ something very different happens. The two paths in figure 2.b start out following each other but around $`t=10`$ they start to diverge apart. Also, the motion of neither looks periodic implying that higher harmonics have been introduced. The proper way to quantify this phenomena is with the Lyapunov characteristic exponent.

If a system is chaotic, then two trajectories separated slightly in phase space will, on average, diverge away from each other with an exponential behavior given by

$$d(t)d_0e^{\lambda t}$$
(12)
Where $`d_0`$ is the initial phase space separation at $`t=0`$ and $`\lambda `$ is known as the Lyapunov characteristic exponent or the largest Lyapunov exponent. Equation 12 should not be taken too literally but only represents an average behavior of dynamical systems. The calculation of $`\lambda `$ is actually an iterative process.

To start off, two trajectories initially separated in phase space by an amount $`d_0`$ are integrate forward for some time $`\tau `$. The system is now separated in phase space by an amount $`d_1`$. At this point the two trajectories are renormalized so that the separation is again $`d_0`$ and then integrated forward again for $`\tau `$. Continuing in this manner, a set of N $`d_i`$ ’s is used to calculate the Nth approximate of the Lyapunov exponent, given by

$$\lambda _N=\frac{1}{N\tau }\underset{i=1}{\overset{N}{}}\mathrm{ln}\frac{d_i}{d_0}$$
(13)
The largest Lyapunov exponent is then given by taking the limit as $`N\mathrm{}`$ of equation 13. The interested reader should follow the references for a more detailed description. mt95 ; gb76 ; gb80 ; sd90 ; je93

In the presence of chaos, $`\lambda `$ tends toward some positive finite value. In the absence of chaos, $`\lambda 0`$. In figure 3, the results of $`\lambda `$ for the above two values of energy separated in phase space by $`d_0=10^6`$ are shown. After $`10^6`$ iterations the results were $`\lambda _{E=1.1}0.004`$ and $`\lambda _{E=1.75}0.187`$. The results imply that for the smaller energy the motion of the FSI is regular and for the larger energy the motion is chaotic.

### IV.2 Poincaré maps

The system under consideration is four dimensional in phase space given by the set $`(y_1,y_2,v_1,v_2)`$. Since energy is a conserved quantity the degrees of freedom can be reduced to three. A trajectory through three dimensional phase space of the system would be impossible to visualize and too complicated to follow. Information about the system is still possible by cutting phase space with a plane and looking at the crossing through this plane by trajectories at the same energy. The pattern left behind on the plane is known as a Poincaré map or a Poincaré surface of section.km90

In figure four, four maps for various values of energy are shown. The surface of section plane was chosen by plotting points $`(v_2,y_2)`$ whenever $`y_1`$ made an impact with the floor. For the lowest value of energy, $`E=1.1`$, the motion is regular and periodic. As the energy is increased to $`E=1.75`$ the motion becomes more complex and the map becomes dense with points. This phase space filling behavior is indicative of chaos.

## V Conclusion

The free-fall spring subject to impact is a simple model which exhibits rich behavior. The reader is invited to carry out other investigations such as higher energies or setting $`m_1m_2`$. The interested reader should see the references for non-technicaljg87 ; fb89  and technical introductions to chaos.eo93

###### Acknowledgements.

I would like to thank the Imperatore School of Science and Arts and the Department of Physics and Engineering Physics at Stevens Institute of Technology. In particular, I would like to thank Prof. George Schmidt for helpful suggestions and Petra Sauer for proof reading this manuscript.

## The Physical System