text:
# Small Cycles in Small Worlds

## Abstract

We characterize the distributions of short cycles in a large metabolic network previously shown to have small world characteristics and a power law degree distribution. Compared with three classes of random networks, including Erdős-Rényi random graphs and synthetic small world networks of the same connectivity, the metabolic network has a particularly large number of triangles and a deficit in large cycles. Short cycles reduce the length of detours when a connection is clipped, so we propose that long cycles in metabolism may have been selected against in order to shorten transition times and reduce the likelihood of oscillations in response to external perturbations.

Systems as diverse as the Western US power grid, metabolic networks of a cell, or the World Wide Web are well described as graphs with characteristic topology. Small world networks have received considerable attention since the seminal paper by Watts and Strogatz Watts:98 .

Most of the existing literature discusses small world networks in terms of the average path length between two vertices Newman:00a  or of the network’s clustering coefficient Herzel:98 ; Barrat:00a  which measures how close the neighborhood of a each vertex comes on average to being a complete subgraph (clique) Watts:98 . Barabási et al. Barabasi:99a ; Barabasi:99b  focussed on the degree distributions, finding a power law in a suite of real world examples including the world wide web or the US power-grid. Recent work on the spread of epidemics on a small world network Pandit:99  emphasizes the importance of “far-reaching” edges. The idea is that clipping a far edge will force a (relatively) long detour in the network. Hence it is these edges that are responsible for the small diameter of the graph $`G`$.

Let us look at detours in graphs in more systematic way. Throughout this paper we will represent a network as a simple (unweighted, undirected) graph $`G(V,E)`$ with vertex set $`V`$ and edge set $`E`$. A cycle in $`G`$ is a closed path which meets each of its vertices and edges exactly once. The length of a cycle $`C`$, i.e., the number of its vertices or edges, is denoted by $`|C|`$. With each edge $`eE`$ we can associate the set $`𝒮(e)`$ containing the shortest cycles in $`G`$ that go through $`e`$. It is easily verified that a far edge in the sense of Pandit:99  is an edge that is not contained in a triangle. In other words, $`e`$ is a far edge if and only if $`𝒮(e)`$ does not contain a triangle. The cycles $`C𝒮(e)`$ determine the shortest detours (which have length $`|C|1`$) when $`e`$ is removed from the graph.

It seems natural to consider the set $`𝒮(G)=_{eE}𝒮(e)`$ of shortest cycles of all edges in $`G`$ and to study e.g. their length distribution. However, as the example in Fig.1 shows, the shortest cycles $`𝒮(G)`$ do not convey the complete information about the graph. Additional cycles appear to be relevant, such as the hexagon in Figure 1. In order to extend $`𝒮(G)`$ to a more complete collection of cycles we need some more information on the cycle structure of graphs. Recall that the set of all subsets of $`E`$ forms an $`|E|`$-dimensional vector space over $`\{0,1\}`$ (with addition and multiplication modulo $`2`$). Vector addition in this *edge space* is given by symmetric difference $`XY=(XY)(XY)`$. The *cycle space* $`𝒞`$ consisting of all cycles and edge-disjoint unions of cycles in $`G`$ is a particularly important subspace of the edge space Chen:71 . Its dimension is the *cyclomatic number* $`\nu (G)=|E||V|+c(G)`$, where $`c(G)`$ is the number of connected components of $`G`$. The length $`\mathrm{}()`$ of a basis $``$ of the cycle space (*cycle basis* for short) is the sum of the lengths of its cycles: $`\mathrm{}()=_C|C|`$. A *minimum cycle basis* (MCB) is a cycle basis with minimum length. MCBs have the property that their longest cycle is at most as long as the longest cycle of any basis of $`𝒞`$ Chickering:95 . A MCB therefore contains the salient information about the cycle structure of a graph in its most compressed form. Most graphs, however, do not have a unique MCB. On the other hand, the distribution of cycle lengths is the same in all MCBs of a given graph Stepanec:64 . The way to avoid ambiguities is to consider the union of all minimum cycles bases, also known as the set $`(G)`$ of  *relevant cycles*. The term “relevant” is justified by two important properties of $`(G)`$: (i) a cycle is relevant if and only if it cannot be written as an $``$-sum of shorter cycles Vismara:97 , and (ii) the shortest cycles through an edge are relevant, i.e., $`𝒮(G)(G)`$ Stepanec:64 ; Zykov:69 . Consequently, the composition of $`(G)`$ in terms of number and length distribution of cycles is an important characteristic of a graph. The numerical studies below make use of Vismara’s Vismara:97  algorithm for computing $`(G)`$, which is based on Horton’s MCB algorithm Horton:87 .

The most common model of graph evolution, introduced by Erdős and Rényi Erdoes:60 , assumes a fixed number $`n=|V|`$ of vertices and assigns edges independently with a certain probability $`p`$ Bollobas:85 . In many cases ER random graphs turn out the be quite different from a network of interest. The Watts-Strogatz Watts:98  model of small world networks starts with a deterministic graph, usually a circular arrangement of vertices in which each vertex is connected to $`k`$ nearest neighbors on each side. Then edges are “rewired” (in the original version) or added Newman:99a ; Newman:00a  with probability $`p`$. We shall consider the latter model for $`k=1`$, denoted SW1 below, which corresponds to adding random edges to a Hamiltonian cycle. Both ER and SW1 graphs exhibit an approximately Gaussian degree distribution.

In many real networks, however, the degree distribution follows a power law. Barabási et al. Barabasi:99a ; Barabasi:99b  show that the scale invariant behavior of the degree distributions can be explained in terms of simple graph evolution model (AB model): Starting from a small core graph, at each time step a vertex is added together with $`m`$ edges that are connected to each previously present vertex $`k`$ with probability $`\mathrm{\Pi }(k)=d(k)/_jd(j)`$, where $`d(j)`$ is the degree of vertex $`j`$. In this contribution we will focus mostly on the AB model instead of Watt’s original construction, because we will apply the analysis of the cycle structure to an empirical network for which a power-law like degree distribution has been established. This network is the system of all chemical reactions required for the synthesis of small-molecule building blocks and energy in the bacterium Escherichia coli. Its structure described in ref. Wagner:00a . Such chemical reaction networks are often referred to as metabolic networks.

It is clear that all triangles in a graph are relevant, since a triangle is for sure a shortest cycle through each edge. Hence $`|(G)|\mathrm{\Delta }`$, where $`\mathrm{\Delta }`$ denotes the number of triangles in $`G`$. We expect $`\mathrm{\Delta }_{\mathrm{ER}}=\left(\genfrac{}{}{0pt}{}{n}{3}\right)p^3`$ triangles in an ER random graph with edge-drawing probability $`p`$. For the SW1 graphs we obtain a similar expression:

$$\mathrm{\Delta }_{\mathrm{SW1}}=np+n(n4)p^2+\frac{1}{6}n(n^29n+20)p^3.$$
(1)
The MCB will therefore consist almost exclusively of triangles if $`\mathrm{\Delta }\nu (G)`$. The average vertex degree is $`d=2|E|/n=p(n1)`$ for ER and $`d=2+p(n3)`$ for SW1, resp. Assuming that $`n`$ is large we expect to find only triangles in $`(G)`$ for $`d\sqrt{3n}`$. Numerical simulations show that this is indeed the case, Fig 2. In this regime, we have $`|(G)|d^3/6`$, and the graph contains no far edges. Not surprisingly, there is little difference between SW1 and ER random graphs for large $`n`$.

Since the AB model is based on a fixed vertex degree $`d`$, it should be compared to random graph models with given vertex degree $`d`$, not with given edge drawing probabilities $`p`$. We have an asymptotically constant number of triangles for both ER and SW1: $`\mathrm{\Delta }_{\mathrm{ER}}d^3/6`$ and $`\mathrm{\Delta }_{\mathrm{SW1}}d^3/6d+2/3`$, resp. Note that as a consequence the clustering coefficient vanishes asymptotically. In SW networks with a priori connectivity $`k>1`$ we find of course a number of triangles that grows at least linearly with $`n`$, since the initial ($`p=0`$) networks already contains $`(k1)n`$ triangles. The clustering coefficient stays finite for large $`n`$ Watts:99 .

The large vertex degree of the “early” vertices in the AB model suggests that there should be many more triangles than in ER or SW1 models. The expected degree of vertex $`s`$ at “time” $`t`$ is known Dorogovtsev:00 : $`d(s|t)=m[\sqrt{t/s}1]`$. The probability of an edge between $`s`$ and $`t`$, $`t>s`$, is therefore $`p_{st}=md(s|t1)/2(t1)m`$, where $`2(t1)m`$ is the sum of the vertex degrees at “time” $`t1`$. Thus $`\mathrm{\Delta }=_{r<s<t}p_{rs}p_{st}p_{rt}`$. This can be approximated by

$`\mathrm{\Delta }`$ $``$ $`{\displaystyle \frac{m^3}{8}}{\displaystyle _{1<r<s<t}^n}(1/st^2)(\sqrt{{\displaystyle \frac{s}{r}}}1)(\sqrt{{\displaystyle \frac{t}{r}}}1)\left(\sqrt{{\displaystyle \frac{t}{s}}}1\right)`$ (2)
$``$ $`Cm^3\mathrm{ln}^3n+𝒪(\mathrm{ln}^2n)`$
Fig. 3 shows $`\mathrm{\Delta }`$ for typical AB-random graphs with $`m=2,\mathrm{},8`$ as a function of “time’. The behavior of $`\mathrm{\Delta }`$ in a individual growing network is well represented by equ.(2).

The number $`||\mathrm{\Delta }`$ of non-trivial relevant cycles has its maximum around $`|E|0.74n^{3/2}`$ independent of the model. The scaling of their number is consistent with $`||\mathrm{\Delta }Cn^{5/2}`$, where the constant $`C0.036`$ is the same for ER and SW1 random graphs and $`C0.016`$ for the AB models. For small vertex degrees, $`d|V|^{1/2}`$ we find $`(G)\nu (G)`$, i.e., the MCB is (almost) unique.

The cyclomatic number of a AB random graph is $`\nu (G)(m/21)n`$; Hence eq.(2) implies that almost all relevant cycles must be long. Fig. 4 shows that the average length of a relevant cycle grows logarithmically with $`n`$. Not surprisingly, the slopes decrease with $`m`$.

Let us now turn to an example of metabolic networks. Because it is germane to their functional analysis, we first point out a nexus between graph representations of metabolic network, and metabolic flux analysis (MFA), the most generic framework to analyze the biological function of metabolic networks. A graph representation of metabolic networks was introduced as a substrate graph $`\mathrm{\Sigma }`$ in Wagner:00a . Its vertices are the molecular compounds (substrates); two substrates $`k`$ and $`l`$ are adjacent in $`\mathrm{\Sigma }`$ if they participate in the same reaction $`r`$. Substrate graphs are undirected because directed graphs would not properly represent the propagation of perturbations: even for irreversible reactions the product concentration may affect the the reaction rate, for instance by product occupancy of the enzyme’s active site; this in turn affects the substrate concentration. Thus, perturbations may travel backwards even from irreversible reactions. A similar argument for considering undirected graphs can be derived from metabolic control theory Sen:91 .

The key ingredient of MFA is the stoichiometric matrix $`𝐒`$. Its entries are the stoichiometric coefficients $`s_{kr}`$, i.e., the number of molecules of species $`k`$ produced ($`s_{kr}>0`$) or consumed ($`s_{kr}<0`$) in each reaction $`r`$. Reversible reactions are entered as two separate reactions in most references. In general, additional “pseudo-reactions” are added to describe the interface of the metabolic reaction network with its environment. Stationary flux vectors $`\stackrel{}{f}`$ in the network satisfy $`𝐒f=\stackrel{}{o}`$ and $`f_r0`$ for each reaction $`r`$, see e.g. Clarke:88 ; Heinrich:96 ; Fell:97 ; Schilling:00 ; Edwards:00a ; SchusterS:00a . It is not hard to see that if all reactions are mono-molecular, then $`𝐒`$ is the incidence matrix of a directed graph; The stationary flux vectors span the cycle space of this graph. The close connection between the cycle space of a directed graph and its underlying undirected graph Bollobas:98  allows us to use the relevant cycles of the substrate graph $`\mathrm{\Sigma }`$ to describe the structure of the metabolic network in a way complementary to that provided by MFA.

For our analysis of metabolic graphs, we use the substrate graph of the Ecoli1 core metabolism, a set of chemical reactions representing the central routes of energy metabolism and small-molecule building block synthesis. Similar to Wagner:00a , we omit the following substrates from the graph: $`\mathrm{𝖢𝖮}_\mathrm{𝟤}`$, $`\mathrm{𝖭𝖧}_\mathrm{𝟥}`$, $`\mathrm{𝖲𝖮}_\mathrm{𝟦}`$, $`\mathrm{𝖠𝖬𝖯}`$,$`\mathrm{𝖠𝖣𝖯}`$, and $`\mathrm{𝖠𝖳𝖯}`$, their deoxy-derivatives, both the oxidized and reduced form of thioredoxine, organic phosphate and pyrophosphate. The resulting graph has $`n=272`$ vertices and $`|E|=652`$ edges. It is analyzed below.

Table 1 shows that the three random models AB, SW1, and ER agree at least qualitatively with each other. The AB random graphs exhibit a much broader distribution of cycle sizes (not shown) than the ER and SW1 models. As a consequence, the average cycle numbers for ER and SW1 have statistical uncertainty of about 2%, while the uncertainty of the AB values is 5 to 10 times higher. Note that ER and SW1 have a similar number of relevant cycles, but the cycles are slightly longer in SW1. Two features distinguish the metabolic network Ecoli1 from all random networks: (1) The number $`\mathrm{\Delta }`$ of triangles is almost 10 times larger than expected. This can be explained at least in part as a consequence of the substrate graph representation: multi-molecular reactions translate to cliques and hence a large number of triangles. The ratio $`282/3790.744`$ indicates that in fact almost all triangles are contained in 4-cliques, since in each 4-clique we have three triangles that belong to a particular MCB, while the fourth face of the tetrahedron is their $``$-sum Gleiss:00a . (2) There is a much smaller number of relevant pentagons and hexagons, which results in an overall somewhat reduced number of relevant cycles: $`723`$ compared to about $`1060`$ (AB), $`904`$ (ER), and $`805`$ (SW1).

Strictly speaking, we do not know the biological significance of this relative paucity of longer cycles. However, we would like to venture a speculation. Organisms are constantly exposed to environmental fluctuations requiring transitions in metabolic states. That is, a metabolic network needs to produce different outputs depending on the environment. Environments may vary rapidly, requiring rapid transition between metabolic states. Possibly, networks with long cycles have longer transition times, because environmental perturbations may lead to prolonged oscillations in such networks. The dynamical system representation of metabolic networks required to test this idea rigorously lies beyond the scope of this article.

The longest relevant cycles in a metabolic networks are of particular interest since they reflect parts of the network that cannot easily be replaced by alternative routes. In Fig. 5 we show the largest such cycle in Ecoli1. We emphasize that the cycles in our analysis represent routes for transmission of perturbations, but not necessarily of mass, as it is commonly considered in MFA. This is apparent from Fig.5 , which does not correspond to a pathway from a biochemical chart, but links serval pathways together.


format: text

text:
# Scales of the Extra Dimensions and their Gravitational Wave Backgrounds

## I Relic Backgrounds from Brane Formation in “the Desert”

If quantum gravity lies fundamentally in ten spatial dimensions, then the seven “extra” dimensions must somehow be hidden. The traditional approach, going back to Kaluza and Klein, is to make extra dimensions very small, close the Planck scale. Recently there has been a resurgence of interest in another approach where many of the extra dimensions are of much larger size. The Standard Model fields are confined to a three-dimensional wall or “brane”, a 3-D defect in the higher-dimensional space; only gravity can propagate in the other dimensions or “bulk”.

These ideas are observationally constrained from the particle-physics side and from the gravity side. Roughly speaking, the success of the Standard Model limits the unification energy scale to energies well above 1 TeV, although some of the most interesting brane+bulk scenarios are those where the unification scale is not many orders of magnitude more than a TeV, providing a solution to to the “hierarchy problem.” Laboratory measurements of gravity limit the size (or radius of curvature) of the extra dimensions to be smaller than about 0.3 mm. Constraints of a model-dependent nature, including astrophysical ones, extend these limits, since the models generally add new particles and interactions to the Standard Model (e.g.,). A broad range of unification scales and brane configurations is consistent with current constraints; the unification scale and extra dimension size can be anywhere between the TeV/mm limits and the Planck scale.

This broad range of possibilities of course just reflects the relative lack of data in the physics/cosmology “desert”. Standard unification ideas often show running couplings extrapolated from actual experiments at 100 GeV by 14 orders of magnitude to the unification scale at $`10^{16}`$ GeV, but there are few direct particle data beyond the TeV scale. A similar lacuna appears in the cosmological data; we can cite direct information from cosmic abundances about microphysical processes as early as 1 MeV, which can be affected by inhomogeneities originating as early as 100 GeV (e.g. ) and about earlier inflationary effects on very large scales (such the microwave background anisotropy, e.g. ), but we have almost no direct data about early “mesoscopic” structure on comoving scales smaller than a few centimeters at 100 GeV (a quantity of energy about equal to the mass of the Earth) at times earlier than $`10^9`$ second (except possibly for the mean baryon-to-photon ratio and the mean density of some forms of dark matter.) Scalar metric perturbations are very effectively erased on these small scales by neutrino diffusion long before they can have any effect on observables such as element abundances. Small-scale perturbations in baryon/photon ratio $`\eta `$ on small scales, even of large amplitude, are erased by nucleon diffusion before nucleosynthesis. From most of cosmic history (in log space) since the Planck time, and on most scales of structure, only the tensor modes— gravitational waves— survive to the present.

New gravitational wave observatories now under development will soon achieve a critical level of sensitivity at frequencies where they can plausibly detect the relic stochastic background of gravitational waves from cosmological events far into the early, previously unobservable “mesoscopic” era. One of the most interesting possible sources of radiation, from the point of view of both physics and cosmology, is the formation of our 3-brane. This paper surveys broadly the circumstances and the values of parameters— the number and the size of large extra dimensions, and the scale of unification— for which the background might be detected with current technology.

To create a detectable background, the metric during the unification era must be far from equilibrium— it must be close (within a few orders of magnitude) to chaotic on the mesoscopic scale of the horizon at unification. Such circumstances can indeed be more generic than the usual assumption of a uniform equlibrium system. The brane+bulk configurations being contemplated are classical setups: the vacua of the fundamental theory have effective potentials that tend to drive them to forming classical three-dimensional defects. The details of the stabilization of the setup are not known but the cosmological formation of such defects is in general not an equilibrium achieved by a microscopic quantum process— it involves macroscopic (or at least mesoscopic) collective dynamics, in the same way that spontaneous symmetry breaking occurs in the standard Higgs mechanism. In the cosmological context, symmetry breaking during the setting up of the brane is likely to create gravitational waves via mechanisms analogous to such intense classical sources as cosmic strings or first-order phase transitions. The properties of the waves can be informative about both the physics of the brane and its early cosmological history.

## II Characteristic Frequency and Amplitude

The formation and stabilization of the brane are associated with two new geometrical degrees of freedom, radion modes governing the size or curvature of the extra dimensions and Nambu-Goldstone modes corresponding to spatial variations in brane displacement. The characteristic scales of the resulting broadband background spectra can be estimated from general scaling considerations. The characteristic amplitude of the metric perturbations induced by brane condensation is of the the order of unity and the characteristic frequency is the horizon scale. We will use the “maximal amplitude” and “Hubble frequency” to estimate the extra-dimensional scales accessible to gravitational wave detectors.

The characteristic scale for setting up our 3+1 world is the gravitational timescale $`H^1`$ in standard 3+1 General Relativity. This is cleanly determined by General Relativity and thermodynamics except for a weak dependence on the particle-physics uncertainties encapsulated in the number of effective relativistic degrees of freedom $`g_{}`$ (e.g.). The characteristic frequency of observed gravitational waves is thus the Hubbble frequency redshifted to the present day, $`f_{H0}(T)H(T)a(T)/a_0`$:

$$f_{H0}=7.65\times 10^5\mathrm{Hz}(T/TeV)g_{}^{1/6}(g_{}/g_S)^{1/3}T_{2.728}=9.37\times 10^5\mathrm{Hz}(H\times 1\mathrm{m}\mathrm{m})^{1/2}g_{}^{1/12}(g_{}/g_S)^{1/3}T_{2.728}.$$
(1)
The estimate is valid back to the threshold of the extradimensional dynamics: $`H^1b_0`$, (where $`b_0`$ is the size or curvature radius of the largest extra dimension), or $`TM_{}`$, whichever comes first.

The fiducial “maximal” energy density is set by the mean energy density in all relativistic species (photons and three massless neutrinos), $`\mathrm{\Omega }_{rel}=8.51\times 10^5h_{70}^2T_{2.728}^4`$, where $`h_{70}`$ refers to the Hubble constant. The energy density of relativistic matter redshifts in the same way as gravitational waves so they have the same ratio of energy densities today as when the waves were generated. (Constraints from nucleosynthesis limit the gravitational waves to about 10% or less of the relativistic matter density, or rms strain about three times smaller than the maximal value.)

At the low frequencies observed by spacecraft interferometers (such as LISA) there are many astrophysical foregrounds, including known sources such as galactic binaries. A stochastic background can however be distinguished from other sources of noise and astrophysical wave sources by resolving the background in frequency. For the projected sensitivity of LISA, a maximal background is detectable above the instrument noise or the other likely astrophysical foregrounds over a frequency range from about $`10^1`$ to $`10^4`$Hz. In the center of this range, $`10^4`$Hz, backgrounds are detectable with

$$\left(\frac{\mathrm{\Omega }_{GW}(\mathrm{\Delta }f=f)}{\mathrm{\Omega }_{rel}}\right)10^6.$$
(2)
At the higher frequencies observed by ground-based observatories (e.g., LIGO, VIRGO, TAMA300, GEO600) it could well be that other stochastic astrophysical backgrounds are relatively weak, and that the observable sources are so limited in time duration and frequency that the main contribution to the rms noise comes from the instrument itself. In this case we can use the entire bandwidth to measure the background and reach a level which places meaningful constraints on backgrounds. The most promising technique is to correlate the signals from two interferometers within a wavelength of each other. Combining two early LIGO or VIRGO systems will likely achieve a sensitivity

$$\left(\frac{\mathrm{\Omega }_{GW}(\mathrm{\Delta }f=f)}{\mathrm{\Omega }_{rel}}\right)10^2,$$
(3)
and LIGO II may reach another two orders of magnitude below this— close to the level where other stochastic astrophysical backgrounds (such as neutron star emission) may be expected. Sufficient sensitivity to constrain primordial stochastic backgrounds may be expected for early ground-based interferometers roughly from 50 to 500 Hz.

## III Excitation of Classical Displacement and Radion Perturbations

There are many possible “brane worlds” with various configurations of large extra dimensions. To simplify the discussion we will assume that there are $`n`$ equal-size largest extra dimensions of size $`b_0`$ and all the others are much smaller, of size $`M_{}^1`$ where $`M_{}`$ is the true fundamental unification scale. The apparent (usual) Planck mass in 3+1-D, $`M_{Planck}`$, is approximately given by $`M_{Planck}^2M_{}^2(M_{}^nV_n)`$ where $`V_n`$ is the volume of the extra dimensions. In our situation,

$$M_{Planck}^2M_{}^2(M_{}^nb_0^n)$$
(4)
where $`b_0^n`$ is the volume of $`n`$ extra large dimensions. Thus each choice of $`n`$ defines a relation between $`M_{}`$ and $`b_0`$, shown in Figure 1. The various relations intersect at $`M_{}=M_{Planck}`$:

$$b_0(M_{Planck}/M_{})^{2/n}M_{}^1.$$
(5)
Gravity is “normal”***Gauss’ law shows that gravity obeys a $`r^2`$ law for $`r>b_0`$, and $`r^{2n}`$ for $`r<b_0`$. on scales larger than $`b_0`$; standard field theory holds on the brane at energies below $`M_{}`$.

The cosmological formation of the brane setup is regulated on the gravitational timescale. In the 3+1-D era, this is just given by the Hubble scale discussed above,

$$H^1/1mm(T/1TeV)^2$$
(6)
which defines the usual gravitational relationship between length and energy. If we follow 3+1-D cosmology back using this equation, we come first to either (1) the time when $`T=M_{}`$ and extra-dimensional effects come into play on a small scale, possibly associated with a first-order phase transition; or (2) the time when $`H=b_0^1`$, and the geometrical degrees of freedom of the extra dimensions are important even on large scales. In general the excitation of the extra-dimensional modes has a different character depending on whether $`b_0(M_{})`$ lies above or below the Hubble relation $`H(T=M_{})`$. The former case produces displacement modes at with a scale determined by $`b_0`$; the latter case produces radion modes with a scale determined by $`M_{}`$.

The Hubble relation $`H(T)`$ is degenerate with the $`n=2`$ relation between $`b_0`$ and $`M_{}`$, meaning that in this case macroscopic and microscopic departures from the standard picture happen at about the same time. The transition to new microscopic unification physics and to 5-dimensional cosmological expansion would both happen at the same point on these lines, also corresponding to the condensation of the brane and the creation of the gravitational waves.

In the $`n=1`$ case where $`b_0>H(T=M_{})^1`$, the condensation of the brane (at $`T=M_{}`$) takes place while the Hubble length is still much smaller than $`b_0`$. In this case a long evolution takes place after the brane forms (and non-gravitational fields are confined to 3D as today), during which the cosmological evolution involves one classically large extra dimension and the usual relativistic cosmology does not apply. In particular there is not time for signals to propagate as far as $`b_0`$ in an expansion time. If the formation of the brane spontaneously breaks Poincaré invariance of the whole spacetime, the position of the brane when it condenses is therefore uncorrelated on scales larger than $`H(T=M_{})^1`$ and less than $`b_0`$. The Nambu-Goldstone modes corresponding to brane displacement are therefore substantially excited by the Kibble mechanism up to wavelength $`b_0`$. These large-amplitude scalar perturbations dynamically couple to the tensor modes to produce a gravitational wave background, in much the same way that gravitational waves are generated by defects confined to a 3+1-D space. The spectrum is peaked at the corresponding redshifted Hubble frequency $`f_{H0}(H^1=b_0)`$ and falls off as a power law at higher and lower frequencies. The spectrum at high frequencies depends on details of the $`4+1`$-D era of cosmic evolution, but may be as strong as $`\mathrm{\Omega }_{GW}(\mathrm{\Delta }f=f)=`$ constant; at lower frequencies, it is estimated to fall off as $`\mathrm{\Omega }_{GW}(\delta f=f)f^7`$.

In the $`3n7`$ case where $`b_0<H(T=M_{})^1`$, the universe can find itself in a macroscopic space of conventional 3-D dimensionality (on the Hubble scale) at the temperature where new extra-dimensional unification physics comes into play on a microscopic scale. Conventional cosmology can therefore be used as a framework. The stabilization of the $`3n7`$ extra dimensions appears as a change of physics, nongravitational and gravitational, which can be described macroscopically as a change in the energy-momentum tensor of cosmic matter, and may be (3D-) spatially inhomogeneous as a result of symmetry breaking. The excitation of these “radion modes” can be described as a change in the order parameter of a vacuum state in 3D, as is familiar in cosmological phase transitions. If this transition described by the effective radion potential is first order, free energy is released in macroscopic flows created by bubble nucleation and collisions, and can have a substantial gravitational-wave component, although for most plausible situations these backgrounds are several orders of magnitude weaker than $`\mathrm{\Omega }_{rel}`$. Figure 1 illustrates that the phase-transition bubbles are smaller than the horizon and the peak of the spectrum lies significantly higher than the redshifted Hubble frequency.Account is not taken of the narrower instrumental bandpasses for the less-than-maximal intensity of the background in these cases.

## IV Range of accessible parameters

Although based on a simplified picture, Figure 1 illustrates accurately the wide range of $`n,b,M_{}`$ potentially accessible via this technique. The possibilities include both the radion and Nambu-Goldstone displacement modes, and include, for some $`M_{}`$ and $`b_0`$, all likely values of $`n`$. The observatories constrain theories with values of these parameters including the limits accessible by other means, near $`M_{}=1`$ TeV and $`b_0=1`$ mm, and extending beyond this range by many orders of magnitude. For example, direct laboratory gravitational experiments become extremely challenging below the mm scale, whereas primordial gravitational wave backgrounds at the highest detectable frequencies (about 1 kHz) probe new classical spacetime dimensions directly (via the displacement modes) on scales 14 orders of magnitude smaller— aided by the cosmological redshift which stretches the waves to detectable wavelengths. Even a nondetection of backgrounds will constrain the behavior of theories with these numbers and sizes of extra dimensions.

Of particular interest are the models which potentially solve the hierarchy problem, where $`M_{}`$ is not many orders of magnitude more than 1 TeV. For these models, backgrounds in the LISA band can arise from all the viable numbers of extra dimensions ($`n=2`$ to $`n=7`$), and a large range of dimension sizes— $`10^{12}\mathrm{mm}b_01`$mm, depending on $`n`$. (It is already known that the $`n=1`$ case requires a unification scale of at least $`10^6`$ TeV to be consistent with the EotWash limits).

The LIGO/VIRGO frequencies probe much higher energy unification scales, from $`M_{}10^4`$ to $`10^{10}`$ TeV, and $`b_0`$ less than $`10^{10}`$mm. The intermediate range of $`b_0`$ and $`M_{}`$, corresponding to the frequency gap between these two sets of experiments, would be obervable with a space interferometer with a smaller baseline than LISA.

The smallest dimension accessible is determined by the highest frequency, which is set by the photon noise limit of LIGO or its successors. It is reasonable to expect high sensitivity with this technology up to about 1 kHz, which probes displacement modes down to $`b_0H^110^{14}`$mm. The largest $`M_{}`$ probed (corresponding to $`n=1`$) can be as high as about $`10^{10}`$ TeV; the smallest $`b_0`$ (corresponding to $`n=7`$) can be as short as $`10^{18}`$ mm. There are no known strong astrophysical sources of gravitational radiation above $`10^4`$Hz, so new types of experiments at higher frequency, with sufficient sensitivity (i.e. $`\mathrm{\Omega }_{GW}\mathrm{\Omega }_{rel}`$), would probe even more extreme physics with little astrophysical contamination.

As an aside, it is interesting to consider the traditional case with all the extra dimensions at the classical Planck scale. This point corresponds to $`M_{Planck}^1b_0M_{}^1`$ for any value of $`n`$. (The relation $`b_0M_{}^1`$ shown in Figure 1 for other values of $`M_{}`$ corresponds to the $`n=\mathrm{}`$ limit of equation 5.) The gravitational wave background in this model has a peak frequency of about $`10^{12}`$ Hz, or a peak wavelength of about 0.3 mm! Although the maximal classical background typically has $`10^{10}`$ times more energy density than tensor perturbations from inflationary quantum effects,Most discussions of primordial tensor modes concentrate on waves excited by quantum effects at inflation. These much weaker waves are probably unobservable at LIGO and LISA frequencies, but may already be observable on much larger scales through their effect on microwave background anisotropy. A gravitational wave background is produced by quantum fluctuations of the graviton during inflation, with rms tensor-mode amplitude (e.g. ) $`h_t(H_{inflation}/M_{Planck})^2`$. This is most easily detected on comoving scales close to the present-day Hubble length, $`fH_02\times 10^{18}`$ Hz, where background radiation anisotropy implies an amplitude $`h_t10^5`$. Of course the amplitude might be very different on the scales of the $`10^4`$ to 1000 Hz waves considered here, but since $`H_{inflation}`$ generally decreases slowly during inflation the “tilt” from a scale-free spectrum generally makes $`h_t`$ smaller at higher frequencies. The energy density of this background is $`h_t^2`$ times smaller than the maximal value ($`h_t=1`$) and is likely to be unobservable. there is no plausible technique for detecting a gravitational wave background at such high frequencies.

In models enclosed by the boxes in Figure 1, properties of these backgrounds can in some circumstances provide measurements of the main parameters of the largest extra dimensions: $`n`$, $`b_0`$ and $`M_{}`$. For example, suppose that a stochastic background is detected with various signatures (such as anisotropy statistics) leading us to suspect that it comes from high redshift. If the background is within one or two orders of magnitude of $`\mathrm{\Omega }_{rel}`$ it is likely from a displacement mode. This in turn indicates that $`n=1`$ or $`n=2`$— one or two extra dimensions are much larger than the rest. Since only $`n=1`$ undergoes a substantial period of classical expansion with $`b_0>H^1`$, the background spectra in these two situations are likely to be qualitatively different: the $`n=2`$ spectrum is more sharply peaked at $`f_{H0}`$ and the $`n=1`$ spectrum includes a more intense high-frequency tail from modes excited before $`H_0^1=b_0`$. Knowing $`n`$, the frequency of the spectral peak then provides an estimate of $`b_0`$ and $`M_{}`$. A nondection of a displacement-mode background conveys a significant constraint on parameters and on cosmology since the conditions for generating a background from symmetry breaking in this regime are relatively generic.

In the case of the $`n3`$ radion modes, it is difficult to break the degeneracies using just the information available from the background measurements. For example, the peak frequency and intensity depend on the details of bubble nucleation dynamics as well as on $`f_{H0}`$; these in turn depend on the details of the effective radion potential. That is, a weaker transition leads to a smaller nucleation scale, a higher frequency peak and a less intense background. (The interpretation of the background might in the case of LISA-band backgrounds be supplemented by other data, since the energy scale $`M_{}`$ is not too far above energies accessible by direct experiment.) A nondetection in this regime will only constrain the corresponding parameters in particular models; lack of a background could just mean that the relevant transition is not strongly first order, but occurs gently and does not generate intense gravitational waves.

###### Acknowledgements.

I am grateful for useful discussions with E. Adelberger, B. Allen, B. Barish, P. Bender, D. Kaplan, A. Nelson, and K. Thorne. I thank the Albert Einstein Institute, Golm, the LIGO Hanford Observatory, and JPL for hospitality. This work was supported at the University of Washington by NSF.


text:
# Shot-noise suppression in Schottky barrier diodes.

## Abstract

We give a theoretical interpretation of the noise properties of Schottky barrier diodes based on the role played by the long range Coulomb interaction. We show that at low bias Schottky diodes display shot noise because the presence of the depletion layer makes negligible the effects of the Coulomb interaction on the current fluctuations. When the device passes from barrier to flat band conditions, the Coulomb interaction becomes active, thus introducing correlation between different current fluctuations. Therefore, the cross-over between shot and thermal noise represents the suppression due to long range Coulomb interaction of the otherwise full shot-noise. Similar ideas can be used to interpret the noise properties of others semiconductor devices.

It is well known that, in the absence of $`1/f`$ contributions, excess noise in Schottky barrier diodes (SBD) exhibits shot noise followed by thermal noise (neglecting hot-carrier effects). Shot noise is attributed to the presence of a barrier which controls the exponential increase of current with applied voltage. By contrast, the thermal behavior is attributed to the series resistance which controls the more or less linear increase of the current at the highest voltages. A phenomenological approach to describe the cross-over from shot to thermal behaviors is to express the current spectral density, $`S_I`$, as the partition between two noise generators as:

$$S_I=\frac{2q\overline{I}R_j^2}{(R_j+R_s)^2}+\frac{4k_BTR_s}{(R_j+R_s)^2}\text{,}$$
(1)
where $`\overline{I}`$ is the average electric current, $`q`$ the electron charge, $`R_j`$ the junction resistance, $`R_s`$ the series resistance, $`k_B`$ Boltzman’s constants and $`T`$ the temperature.

Microscopic approaches at hydrodynamic and kinetic levels confirm the above expression without resorting to two independent noise sources, but providing a self-consistent solution of the dynamics together with the Poisson equation. However, a physical interpretation of the noise properties of SBD, and in particular of the cross-over between shot and thermal noise, is still lacking in the current literature and constitutes an intriguing question.

The aim of this paper is to address this issue by showing that the cross-over between shot and thermal noise behaviors represents the suppression due to long range Coulomb interaction of the otherwise full shot-noise. This suppression occurs in SBD when the device passes from barrier to flat band conditions because of the following reason. By increasing the free carrier concentration the originally negligible screening effects due to the Coulomb interaction become relevant, thus acting as a shot noise suppressor. Accordingly, when controlled by screening effects carrier number fluctuations give a negligible contribution to the total noise power, which is now practically given by the velocity fluctuations only, thus providing thermal noise.

The starting point of our approach is the non-degenerate fluctuating drift-diffusion equation which, for the case of one dimensional geometries and low frequencies, reads,

$$\frac{I(t)}{A}=qn(x,t)\mu E(x,t)+qD\frac{dn(x,t)}{dx}+\frac{\delta I_x(t)}{A}\text{,}$$
(2)
where $`I`$ is the instantaneous electric current, $`A`$ the cross-sectional area, $`n`$ the number density, $`\mu `$ the mobility, $`E`$ the electric field, and $`D`$ the diffusion coefficient. Moreover, $`\delta I_x(t)`$ is a Langevin diffusion noise source that accounts for current fluctuations, and has zero mean and correlation function

$$\overline{\delta I_x(t)\delta I_x^{}(t^{})}=2Aq^2D\overline{n}(x)\delta (xx^{})\delta (tt^{})\text{.}$$
(3)
Henceforth, a bar denotes an averaged quantity. In previous equations the time variable just specifies that the quantities are changing with time due to the existence of fluctuations. Equation (2) has to be completed with Poisson equation and with appropriate boundary conditions to take into account the long range Coulomb interaction and the thermionic emission processes, respectively. For the rectifying contact the boundary condition reads

$$n_0(t)n^{eq}=\frac{I(t)\delta I_{TE}(t)}{qAv_r}\text{,}$$
(4)
where $`v_r`$ is the recombination velocity of the contact, $`n^{eq}=N_Ce^{q\beta \varphi _{bn}}`$ the contact equilibrium density, with $`N_C`$ being the effective density of states in the conduction band, $`\varphi _{bn}`$ the barrier height, and $`\beta =1/k_BT`$, and $`\delta I_{TE}(t)`$ is an additional Langevin noise source that accounts for the fluctuations due to the random emission of carriers. It has zero mean and spectral density

$$S_I^{TE}=2q\left[\overline{I}+2I_C\right],$$
(5)
with $`I_C=qn^{eq}v_rA`$. For the ohmic contact, if the sample is long enough we can neglect their effects and take

$$n_L(t)=N_D\text{.}$$
(6)
By dividing Eq.(2) by the wronskian

$$W(x,t)=e^{q\beta _0^xE(x^{},t)𝑑x^{}}\text{,}$$
(7)
and integrating the result between $`x=0`$ and $`x=L`$, where $`L`$ is the sample length, one obtains (see also Refs.  )

$$I(t)=I_L(t)I_0(t)+\delta I_D(t)\text{,}$$
(8)
where

$$I_0(t)=qAv_D(t)n_0(t)\text{}I_L(t)=qAv_D(t)\frac{n_L(t)}{W_L(t)}\text{,}$$
(9)
and

$$\delta I_D(t)=\overline{v}_D_0^L\frac{\delta I_x(t)}{D\overline{W}(x)}𝑑x\text{.}$$
(10)
To arrive at the previous results we have made use of Einstein’s relation $`D=\mu k_BT/q`$, and have introduced the diffusion velocity as

$$\frac{1}{v_D(t)}=_0^L\frac{dx}{DW(x,t)}\text{.}$$
(11)
From Eq.(8), the average current is given by

$$\overline{I}=\overline{I}_L\overline{I}_0=I_s\left(e^{q\beta \overline{V}}1\right)\text{,}$$
(12)
where $`I_s=qn^{eq}\frac{v_r\overline{v}_D}{v_r+\overline{v}_D}A`$ and $`\overline{V}`$ is the applied bias. In deriving Eq.(12) we have used the averaged version of Eqs. (4), (6) and (9), and the fact that $`\overline{W}_L=\mathrm{exp}\left(q\beta (\overline{V}V_{bi})\right)`$, with $`V_{bi}=k_BT/q\mathrm{ln}(N_D/n^{eq})`$ being the built-in voltage.

Notice that in the definition of $`\overline{v}_D`$ we have included all the sample (and not only the depletion layer as usual). As a result, the formal expression for the $`IV`$ characteristics holds in the whole range of applied bias. For low applied bias $`I_s`$ is nearly constant, and the $`IV`$ characteristics is exponential. At larger applied bias (above flat band conditions), $`I_s`$ depends on $`\overline{V}`$ in such a way as to cancel the exponential behavior and give rise to a linear characteristics, as should be.

The current fluctuations can be obtained by linearizing Eq.(8) around the stationary state. After some algebra one arrives at

$$\delta I(t)=\delta I_L^{}(t)\delta I_0^{}(t)+\delta I_E(t)+\delta I_D(t)\text{,}$$
(13)
where

$$\delta I_0^{}(t)=qA\overline{v}_D\delta n_0(t)\text{}\delta I_L^{}(t)=qA\frac{\overline{v}_D}{\overline{W}_L}\delta n_L(t)\text{,}$$
(14)
and

$$\delta I_E(t)=_0^LH(x)\delta E_x(t)\text{,}$$
(15)
with

$$H(x)=Aq^2\beta \overline{v}_D\frac{\overline{n}(x)}{\overline{W}(x)}=q^2\beta A\overline{v}_DN_Ce^{\beta (F_n(x)E_C(0))}\text{.}$$
(16)
Here $`E_C(0)`$ is the edge of the conduction band and $`F_n(x)`$ the quasi-Fermi level. By substituting the linearized version of Eqs. (4) and (6) into Eq. (14), and substituting the result into Eq.(13), one finally obtains

$`\delta I(t)`$ $`=`$ $`{\displaystyle \frac{\overline{v}_D/v_r}{1+\overline{v}_D/v_r}}\delta I_{TE}(t)+{\displaystyle \frac{1}{1+\overline{v}_D/v_r}}\delta I_D(t)+`$ (18)
$`{\displaystyle \frac{1}{1+\overline{v}_D/v_r}}\delta I_E(t)\text{,}`$
The meaning of the different contributions to the current fluctuations in Eq.(18) is the following:

(i) The term involving $`\delta I_{TE}(t)`$ represents the contribution to the current fluctuations of the thermionic emission processes.

(ii) The term involving $`\delta I_E(t)`$ describes the current fluctuations resulting from a response to an electric field fluctuation. This term accounts for the effects of the long range Coulomb interaction on the current fluctuations. To obtain their statistical properties one needs to compute explicitly the electric field fluctuations. For the present paper it is enough to notice that $`\delta I_E(t)`$ is negligible when $`H(x)`$ is constant (or nearly constant) , since then $`\delta I_E(t)H_0^L𝑑x\delta E_x(t)=H\delta V(t)=0`$, where we have used that current fluctuations are calculated under fixed bias conditions.

(iii) The term involving $`\delta I_D(t)`$ represents the current fluctuations due only to diffusion processes since it is directly related to the fundamental diffusion noise source, $`\delta I_x(t)`$. From Eqs.(3) and (10) it can be shown that the current fluctuations associated to $`\delta I_D(t)`$ always display a low frequency spectral density, $`S_I^D=2_{\mathrm{}}^+\mathrm{}𝑑t\overline{\delta I_D(t)\delta I_D(0)}`$ , of the shot noise type (see Refs.  )

$`S_I^D`$ $`=`$ $`2q\left[\overline{I}_0+\overline{I}_L\right]=2q\overline{I}\left({\displaystyle \frac{\overline{I}_L+\overline{I}_0}{\overline{I}_L\overline{I}_0}}\right)=`$ (20)
$`2q\overline{I}\mathrm{coth}\left({\displaystyle \frac{\mathrm{\Delta }F_n}{2k_BT}}\right)\text{,}`$
where $`\mathrm{\Delta }F_n=F_n(L)F_n(0)`$ is the quasi-Fermi level drop across the sample.

Since thermionic emission and diffusion are independent processes, one can compute their joint contribution to the current fluctuations. We obtain

$`S_I^{TE,D}`$ $`=`$ $`\left({\displaystyle \frac{\overline{v}_D/v_r}{1+\overline{v}_D/v_r}}\right)^2S_I^{TE}+\left({\displaystyle \frac{1}{1+\overline{v}_D/v_r}}\right)^2S_I^D`$ (21)
$`=`$ $`2q\left[\overline{I}+2I_s\right]=2q\overline{I}\mathrm{coth}\left({\displaystyle \frac{q\overline{V}}{2k_BT}}\right)\text{,}`$ (22)
which is of the full shot noise type. This result is just a manifestation of the fact that both diffusion and thermionic emission processes are in origin uncorrelated processes. Notice, that from the present derivation, this result is valid for the whole range of applied bias.

Now we are in a position to develope a physical interpretation of the noise properties of SBD, in particular of the cross-over between shot and thermal noise behaviors. Let us start considering the case of low applied bias. As is well known, in this range of bias the device is controlled by the depletion layer present near the rectifying contact, where almost of the applied potential drops. This implies that the quasi-Fermi level is nearly constant all over the sample, and hence, according to Eq.(16), $`H(x)`$ is also (nearly) constant. Therefore, we may neglect the contribution of the long range Coulomb interaction term $`\delta I_E(t)`$. Accordingly, the current fluctuations are determined by the diffusion and thermionic emission terms, thus giving rise to a spectrum of the shot noise type (see Eq.(22)). On the contrary, for higher applied bias the depletion layer disappears and the applied bias drops all over the sample. In this situation the quasi-Fermi level is no longer constant, and the self-consistent term $`\delta I_E(t)`$ plays a relevant role. This term introduces correlations between the otherwise uncorrelated thermionic emission and diffusion processes, thus suppressing the shot noise. Therefore, the thermal noise associated to series resistance may be interpreted as just suppressed shot noise due to long-range Coulomb interaction.

The previous results can be understood physically by considering the density profile corresponding to these devices. At low voltages these devices display a density distribution with a deep depletion layer near the interface followed by a quasi-neutral region that extends up to the back ohmic contact. Since the device is controlled by the depletion layer, one only needs to take into account this part in the discussion. But, in the depletion layer the number density becomes very small. As a consequence, screening effects are very small and long range Coulomb interaction can be neglected, in agreement with the mathematical prove presented above. On the other hand at higher voltages, the density distribution is nearly uniform and equal to the doping density, thus providing appreciable screening effects, which are the responsible of suppressing the shot noise.

In summary, we have developed a theoretical approach to interpret physically the noise properties of Schottky diodes. We have shown that the occurrence of shot noise at low bias has to be attributed to the fact that the density distribution is such as to make negligible the effects of the long range Coulomb interaction. As a result the originally uncorrelated diffusion and thermionic emission processes remain uncorrelated, thus giving rise to a shot noise spectrum. On the other hand, at higher applied bias the nearly uniform carrier density profile induces considerable screening effects, which introduce correlation between the diffusion and thermionic emission process. The resulting thermal noise spectrum is thus the result of the suppression of the intrinsic shot noise of the device. The main ideas developed in the present paper can be extended to other devices, like $`pn`$ diodes, tunneling diodes, ballistic devices or mesoscopic devices, thus offering new perspectives on what concerns the understanding of the noise properties of semiconductor devices.

Acknowledgments

Drs. T. González and O. M. Bulashenko are thanked for stimulating discussions on the subject. Partial support from the Italian MURST ”Physics of Nanostructures”, from the Spanish SEUID and from the NATO linkage grant HTECH.LG 974610, is gratefully acknowledged.


text:
# 1 Introduction

## 1 Introduction

In 1990, a complete global analysis of an SIRS model of disease transmission in a population with varying size was given by Busenberg and van den Driessch . In an SIRS epidemiological model, we divide the population into three groups, Susceptibles, Infectives and Removeds, and the problem is to examine the behavior of the size of each group when the time goes to infinity. They considered a disease with horizonal transmission, that is a disease which is transmitted by contact between an infective and a susceptible individual. We intend to generalize their results for a vertically transmitted disease, that is a disease which is also transmitted from infective parents to their newborns. The assumption of vertical transmission has two consequences. First it causes some newborns to die and forces us to assume that the birth rate differs from one group to another. (See the demographic assumption in .) The second fact is that some newborns are infected, hence a group of newborns enter to the infective class . We also assume that a part of these infected newborns are known and removed after their birth. We shall show that the latter parameter can play an important role in the epidemic process.

We consider an SIRS epidemiological model for a vertically transmitted disease. In our model, the incidence function is of proportionate mixing type introduced by Nold . Natural births and deaths are assumed to be proportional to the class numbers with different rates. We also assume that a proportion of the infected newborns are known and removed. We follow  to examine our model equation which is homogeneous of degree one. We consider the proportions system and show that this system has no periodic orbit in its feasibility region. This reduces our analysis to the discussion of existence and stability of rest points of a palnar system. The technique used here to show the uniquness of endemic equilibria is based on the Poincaré index. This technique has no hard analysis and can be easily applied to other similar systems . The reader can verify that our results hold for similar SIRI systems as well .

We first in the next section state the model and a result concerning the non-existence of certain types of solutions . We consider the proportions system and prove that every solution in the feasibility region tends to a rest point of this system. In Section 3., we introduce a quadratic planar system with the same dynamics as the proportions system and then we discuss the existence and stability of rest points of this quadratic planar system. This gives a complete global analysis of the proportions system which is used to provide a global analysis of the original system is Section 4.

## 2 The model

In order to derive our model equations, we divide the population into three classes, the susceptible, the infective and the removed individuals with total numbers $`S`$, $`I`$ and $`R`$ respectively,. We set $`N=S+I+R`$ which is the total size of the population. The following parameters are used in our model equations:

* per capita birth rate of susceptible individuals,
* per capita birth rate of infective individuals born uninfected,
* per capita birth rate of infective individuals born infected,
* per capita birth rate of removed individuals,
* per capita disease free death rate,
* excess per capita death rate of infected individuals,
* excess per capita death rate of removed individuals,
* per capita removal rate of infective individuals,
* per capita recovery rate of removed individuals,
* effective per capita contact rate of infective individuals.

As mentioned before, we assume that the infected newborns enter the classes $`I`$ and $`R`$ of proportions $`\beta _1`$ and $`\beta _2`$ respectively, hence $`\beta =\beta _1+\beta _2`$. In this paper, all the above parameters, are positive, however some of them can also be zero. These hypotheses yield the following system of differential equation in $`_3^+`$, where “  ” denotes the derivatives with respect to $`t`$, the time.

$`\{\begin{array}{c}S^{}=(b_0d)S+b_1I+(b_2+\gamma )R\lambda {\displaystyle \frac{IS}{N}}(21)\hfill \\ I^{}=(\beta _1d\epsilon _1\alpha )I+\lambda {\displaystyle \frac{IS}{N}}(22)\hfill \\ R^{}=(\beta _2+\alpha )I(d+\epsilon _2+\gamma )R(23)\hfill \end{array}`$

where $`\lambda {\displaystyle \frac{IS}{N}}`$ is of the proportionate or random mixing type . By adding the above three equations, the total population equation is

$$N=b_0S+(b_1+\beta \epsilon _1)I+(b_2\epsilon _2)RdN.$$

If we consider the proportions $`s={\displaystyle \frac{S}{N}}`$, $`i={\displaystyle \frac{I}{N}}`$, $`r={\displaystyle \frac{R}{N}}`$, we get the following system of equations.

$`\{\begin{array}{c}s^{}=b_0s+b_1i+(b_2+\gamma )rb_0s^2(b_1+\beta +\lambda \epsilon _1)is(b_2\epsilon _2)sr(21)^{}\hfill \\ i^{}=(\beta _1\epsilon _1\alpha )i+(\lambda b_0)is(b_1+\beta \epsilon _1)i^2(b_2\epsilon _2)ir(22)^{}\hfill \\ r^{}=(\beta _2+\alpha )i(\epsilon _2+\gamma )rb_0sr(b_1+\beta \epsilon _1)ir(b_2\epsilon _2)r^2(23)^{}\hfill \end{array}`$

This system is called proportions system and the feasibility region of this system is the triangle

$$D=\{(s,i,r)|s0,i0,r0,s+i+r=1\}.$$

If we set $`\mathrm{\Sigma }=s+i+r`$ then $`\mathrm{\Sigma }^{}=(1\mathrm{\Sigma })(b_0s+(b_1+\beta \epsilon _1)i+(b_2\epsilon _2)r)`$. Thus the plane $`\mathrm{\Sigma }=1`$ is invariant. Moreover, on the sides of $`D`$, we have:

$`\{\begin{array}{c}s=0s^{}=b_1i+(b_2+\gamma )r0,\hfill \\ i=0i^{}=0\text{hence the line}\{i=0\}\text{is invariant},\hfill \\ r=0r^{}=(\beta _2+\alpha )i.\hfill \end{array}`$

Therefore, $`D`$ is positively invariant. On the invariant line $`\{\mathrm{\Sigma }=1\}\{i=0\}`$, we have

$$r^{}=(\epsilon _2+\gamma )rb_0r(1r)(b_2\epsilon _2)r^2=(b_0+\epsilon _2+\gamma )r+(b_0b_2+\epsilon _2)r^2.$$

It follows that this invariant line contains two rest points, the Disease-Free Equilibrium $`(1,0,0)`$ and possibly another one which is outside of D. It is easy to see that the DFE attracts the side $`D\{i=0\}`$. Furthermore, our vector field is strictly inward on the other sides of $`D`$. Thus $`\stackrel{}{D}`$, the interior of $`D`$, is positively invariant too. The following theorem reduces our problem to the discussion of existence and stability of rest points in $`D`$.

### Theorem 2.1.

The $`\omega `$-limit set of any solution for the system $`(21)^{}(23)^{}`$ with initial point in $`D`$ is a rest point in $`D`$.

### Proof.

Since $`D`$ is compact and positively invariant, the $`\omega `$-limit set of any solution with initial point in $`D`$ is a compact nonempty invariant subset of $`D`$. Here we use the Poincaré-Bendixon theorem and the terminology used to prove it . In the next section, we will see that this system has at most three rest points in $`D`$. Thus it satisfies the assumptions of Poincaré-Bendixon theorem. We follow  and define the vector field $`g=(g_1,g_2,g_3)`$ on $`D`$ by

$`g_1(i,r)=[0,{\displaystyle \frac{f_3(i,r)}{ir}},{\displaystyle \frac{f_2(i,r)}{ir}}],`$
$`g_2(s,r)=[{\displaystyle \frac{f_3(s,r)}{sr}},0,{\displaystyle \frac{f_1(s,r)}{sr}}],`$
$`g_3(s,i)=[{\displaystyle \frac{f_2(s,i)}{si}},{\displaystyle \frac{f_1(s,i)}{si}},0],`$
where $`f_1,f_2`$ and $`f_3`$ are the right hand side of $`(21)^{}`$, $`(22)^{}`$ and $`(23)^{}`$ reduced to functions of two variables by using $`=1`$ respectively. Clearly $`g.f=0`$ in $`\stackrel{}{D}`$ and after some computations , we get

$$(\text{curl}g).(1,1,1)=\left(\frac{b_1}{s^2r}+\frac{b_2+\gamma }{s^2i}+\frac{\beta _2+\alpha }{sr^2}\right).$$

Since the DFE is the only invariant subset of $`D`$ (i.e. the boundary of $`D`$), the $`\omega `$-limit set must have some regular point in $`\stackrel{}{D}`$ if it is not a rest point. Let $`x`$ be such a regular point and $`h`$ be the first return map (Poincaré map) defined on a tranversal at $`x`$. For a point $`y`$ near $`x`$ on the transversal, Let $`V`$ be the region surrounded by the orbit $`\mathrm{\Gamma }`$ from $`y`$ to $`h(y)`$ and the segment between them. (This region is known as Bendixon Sack, Sec Fig 2.1) Now by Stokes’ theorem

$$_V(\text{curl}g).(1,1,1)d\sigma =_\mathrm{\Gamma }g.fdt+_0^1g(ty+(1t)h(y)).(yh(y))dt$$

Since $`g.f=0`$ and $`h(x)=x`$, the right hand side of the above equality tends to zero when $`y`$ tends to $`x`$, but the left hand side tends to the integral over the region bounded by the $`\omega `$-limit set, this is a contradiction since $`(\text{curl}g).(1,1,1)<0`$ in $`\stackrel{}{D}`$. $`\mathrm{}`$

## 3 The planar system

using the relation $`s+i+r=1`$, we see that our system is essentially two dimensional. Thus we can eliminate one of the variable to arrive at the following quadratic planar system:

$$\{\begin{array}{c}i^{}=(\lambda +\beta _1b_0\epsilon _1\alpha )i+(b_0+\epsilon _1\lambda b_1\beta )i^2+(b_0+\epsilon _2\lambda b_2)ir(31)\hfill \\ r^{}=(\beta _2+\alpha )i(b_0+\epsilon _2+\gamma )r+(b_0+\epsilon _1b_1\beta )ir+(b_0+\epsilon _2b_2)r^2(32)\hfill \end{array}$$

The dynamics of the system $`(21)^{}(23^{})`$ in $`D`$ is equivalent to the dynamics of this planar system in the positively invariant region $`D_1=\{(i,r):i0,r0,i+r1\}`$. This quadratic system has at most four rest points and since $`D_1`$ misses a rest point on $`i=0`$, there are at most three rest points in $`D_1`$. One of these rest points is the origin which comes from the DFE. The matrix of the linearization of the system $`(31),(32)`$ at the origin is:

$$\left[\begin{array}{cc}\lambda +\beta _1b_0\epsilon _1\alpha \hfill & 0\\ \beta _2+\alpha \hfill & (b_0+\epsilon _2+\gamma )\end{array}\right]$$

with the eigenvalues $`\lambda +\beta _1b_0\epsilon _1\alpha `$ and $`(b_0+\epsilon _2+\gamma )`$. Now we define the first threshold parameter $`R_0={\displaystyle \frac{\lambda +\beta _1}{b_0+\epsilon _1+\alpha }}`$ which governs the stability of the origin.

### Theorem 3.1.

The origin is globally asymptotically stable in the feasibility region $`D_1`$ when $`R_01`$ and it is a saddle point when $`R_0>1`$.

### Proof.

By Theorem 2.1. it is enough to prove that if $`R_01`$, then the origin is the only rest point in $`D_1`$. If there exists a rest point in $`\stackrel{}{D_1}`$, we have $`i^{}=0`$ and $`i0`$ at this point. Thus it belongs to the line

$$(\lambda +\beta _1b_0\epsilon _1\alpha )+(b_0+\epsilon _1\lambda b_1\beta )i+(b_0+\epsilon _2\lambda b_2)r=0.$$
$`(33)`$
Since $`s^{}=0`$ at this point, from $`(21)^{}`$ we obtain

$$b_0s+b_1i+(b_2+\gamma )rb_0s^2(b_1+\beta +\lambda \epsilon _1)is(b_2\epsilon _2)sr=0.$$

and by using the relation $`s+i+r=1`$, we can write

$$b_1i+(b_2+\gamma )r+(b_0+\epsilon _1b_1\beta \lambda )is+(b_0+\epsilon _2b_2)sr=0.$$

Multiplying (3-3) by $`(s)`$ and adding it to the above expression, we get the following equality

$$b_1i+(b_2+\gamma )r+\lambda sr+(b_0+\epsilon _1+\alpha \lambda \beta _1)s=0.$$

But the left hand side is positive when $`R_01`$ and this is a contradiction. $`\mathrm{}`$

When $`R_0>1`$, the origin is a saddle point and it does not attract any point of $`D_1\{i=0\}`$. Thus the orbits with initial point in $`D_1\{i=0\}`$ must be attracted by some rest points in $`\underset{1}{\overset{}{D}}`$ by Theorem 2.1. These rest points belong to the line (3-3) and the conic section

$$r^{}=(\beta _2+\alpha )i(b_0+\epsilon _2+\gamma )r+(b_0+\epsilon _1b_1\beta )ir+(b_0+\epsilon _2b_2)r^2=0.$$

It follows that there are at most two rest points in $`\underset{1}{\overset{}{D}}`$. Notice that a nondegenerate rest point of the planar system is obtained by a transverse intersection of the line (3-3) and the above conic section. The following lemma has two immediate consequences which will be very helpful.

### Lemma 3.2.

The trace of the linearization of the system (3-1),(3-2) at a rest point in $`\underset{1}{\overset{}{D}}`$ is negative.

### Proof:

We compute the trace at a rest point.

$$\begin{array}{cc}\frac{i^{}}{i}\hfill & =(\lambda +\beta _1b_0\epsilon _1\alpha )+2(b_0+\epsilon _1\lambda b_1\beta )i+(b_0+\epsilon _2\lambda b_2)r,\hfill \\ \frac{r^{}}{r}\hfill & =(b_0+\epsilon _2+\gamma )+(b_0+\epsilon _1b_1\beta )i+2(b_0+\epsilon _2b_2)r.\hfill \end{array}$$

From $`i^{}=0`$ and $`r^{}=0`$ at a rest point, we can write

$$\frac{i^{}}{i}=(b_0+\epsilon _1\lambda b_1\beta )i\text{and}\frac{r^{}}{r}=(b_0+\epsilon _2b_2)r(\beta _2+\alpha )\frac{i}{r}.$$

Since $`s^{}=b_1i+(b_2+\gamma )r+(b_0+\epsilon _1\lambda b_1\beta )si+(b_0+\epsilon _2b_2)sr`$, it follows that $`(b_0+\epsilon _1\lambda b_1\beta )si+(b_0+\epsilon _1\lambda b_1\beta )sr<0`$ and hence $`{\displaystyle \frac{i^{}}{i}}+{\displaystyle \frac{r^{}}{r}}<0`$. $`\mathrm{}`$

### Corollary 3.3.

The system (3-1),(3-2) has no source point in $`\underset{1}{\overset{}{D}}`$

### Corollary 3.4.

Every nondegenerate rest point in $`\underset{1}{\overset{}{D}}`$ is hyperbolic.

### Theorem 3.5.

If $`R_0>1`$, then there exists a unique rest point $`(i^{},r^{})`$ in $`\underset{1}{\overset{}{D}}`$ which is hyperbolic and attracts $`D_1\{i=0\}`$.

### Proof.

When $`R_0>1`$, the origin is a saddle point with the unstable eigenvector

$$\left[\begin{array}{c}\lambda +\beta _1b_0\epsilon _1\alpha +b_0+\epsilon _2+\gamma \\ \beta _2+\alpha \end{array}\right].$$

Since $`R_0>1`$, we have $`\lambda +\beta _1b_0\epsilon _1\alpha >0`$ and hence this vector belongs to the first quadrant of the plane $`(i,r)`$. Since $`\underset{1}{\overset{}{D}}`$ is positively invariant it follows that a branch of the unstable manifold of the origin lies in $`\underset{1}{\overset{}{D}}`$. (See Figure 3.1.)

This helps us to find a piece-wise smooth Jordan curve $`C`$ on which our vector field $`X`$ is either tangent or inward. (See Figure 3.2.) The Poincaré index of such a Jordan curve is 1. (See , Lemma 5.1.)

We choose this Jordan curve so that it contains all rest points in $`\underset{1}{\overset{}{D}}`$. If there are two rest points in $`\underset{1}{\overset{}{D}}`$, they are obtain by a transverse intersection of the line (3-3) and the conic section $`r^{}=0`$ and hence both are nondegenerate. Thus their Poincaré index must be $`\pm 1`$ which contradicts $`I_X(C)=1`$. Therefore there is a unique rest point in $`\underset{1}{\overset{}{D}}`$ which attracts $`D_1\{i=0\}`$ by Theorem 2.1. It remains to prove that this rest point is hyperbolic. Suppose the contratry, then it must be nondegenerate by Corollary 3.4. Thus it is obtained by a tangent (non-transverse) intersection of the line (3-3) and the conic section $`r^{}=0`$. With a slight perturbation in $`\gamma `$, we will have either two transverse intersection in $`\underset{1}{\overset{}{D}}`$ or nothing. ($`\gamma `$ appears only in the coefficient of $`r`$ in $`r^{}=0`$ and does not appear in (3-3) and $`R_0`$). This is a contradiction with the uniqueness of the rest point in $`\underset{1}{\overset{}{D}}`$ proven above. $`\mathrm{}`$

Theorem 3.1. and Theorem 3.5. provide a complete global analysis of the planar system (3-1),(3-2) in $`\underset{1}{\overset{}{D}}`$. Since the dynamics of this system in $`\underset{1}{\overset{}{D}}`$ is equivalent to the dynamics of proportions system $`(21)^{}(23)^{}`$ in $`D`$, we have proved the following result which gives a complete global analysis of the proportions system in the feasibility region $`D`$.

### Theorem 3.6.

Consider the proportions system $`(21)^{}(23)^{}`$.

* If $`R_01`$, then the disease free equilibrium proportions $`(1,0,0)`$ is globally asymptotically stable in $`D`$.
* If $`R_0>1`$, then there is a unique rest point $`(s^{},i^{},r^{})`$ which is globally asymptotically stable in $`D\{i=0\}`$.

## 4 Analysis of the model equations

Consider the original model equation (2-1)-(2-3) and recall that the population equation is$`N^{}=b_0S+(b_1+\beta \epsilon _1)I+(b_2\epsilon _2)RdN.`$ Thus

$$\frac{N^{}}{N}=b_0s+(b_1+\beta \epsilon _1)i+(b_2\epsilon _2)rd.$$
$`(41)`$
If $`R_01`$, then $`(s,i,r)(1,0,0)`$ by Theorem 3.6., hence $`{\displaystyle \frac{N^{}}{N}}b_0d`$. Moreover if $`R_0>1`$ and $`I>0`$, then $`(s,i,r)(s^{},i^{},r^{})`$, i.e. the unique rest point in $`\underset{1}{\overset{}{D}}`$, and

$$\frac{N^{}}{N}b_0s^{}+(b_1+\beta \epsilon )i^{}+(b_2\epsilon _2)r^{}d.$$

We define the second threshold parameter which governs the total population as follows.

$$R_1=\{\begin{array}{cc}\frac{b_0}{d}\hfill & \text{if}R_01,\hfill \\ \frac{b_0s^{}+(b_1+\beta )i^{}+b_2r^{}}{d+\epsilon _1i^{}+\epsilon _2r^{}}\hfill & \text{if}R_0>1.\hfill \end{array}$$

Since $`I=0`$ is invariant with a linear equation, we may assume that $`I>0`$. Now from (2-2), we write

$$\frac{I^{}}{I}=(\beta _1d\epsilon _1\alpha )+\lambda s.$$
$`(42)`$
If $`R_01`$, then $`{\displaystyle \frac{I^{}}{I}}\beta _1d\epsilon _1\alpha +\lambda `$ and if $`R_0>1,{\displaystyle \frac{I^{}}{I}}\beta _1d\epsilon _1\alpha +\lambda s^{}`$. So we define the third threshold parameter which governs the total number of infective individuals.

$$R_2=\{\begin{array}{cc}\frac{\beta _1+\lambda }{d+\epsilon _1+\alpha }\hfill & \text{if}R_01,\hfill \\ \frac{\beta _1+\lambda s^{}}{d+\epsilon _1+\alpha }\hfill & \text{if}R_0>1.\hfill \end{array}$$

Notice that $`d`$ does not appear in the proportions system and hence $`(s^{},i^{},r^{})`$ is independent of $`d`$. The following results provide a rather complete global analysis of the model equations (2-1)-(2-3).

### Lemma 4.1.

If $`I(t)M`$ for every $`tt_0`$, then $`R(t)0`$ and if $`I(t)\mathrm{}`$, then $`R(t)\mathrm{}`$.

### Proof.

From (2-3), we have $`R^{}(t)=(\beta _2+\alpha )I(t)(d+\epsilon _2+\gamma )R(t)`$. If $`I(t)M`$ for $`tt_0`$, then $`R^{}(t)M(\beta _2+\alpha )(d+\epsilon _2+\gamma )R(t)`$ and by Granvell’s inequality , $`R(t)M(\beta _2+\alpha )e^{(d+\epsilon _2+\gamma )(tt_0)}`$ which follows that $`R(t)0`$. Now suppose that $`I(t)\mathrm{}`$. Then by (2-3), we have $`R^{}(t)+(d+\epsilon _2+\gamma )R(t)=(\beta _2+\alpha )I(t)`$ which implies that $`{\displaystyle \frac{d}{dt}}(R(t)e^{(d+\epsilon _2+\gamma )t})=(\beta _2+\alpha )I(t)e^{(d+\epsilon _2+\gamma )t}`$. For every $`M^+`$, there is a $`t_0`$ such that $`I(t)>M`$ for $`t>t_0`$. Thus

$$R(t)e^{(d+\epsilon _2+\gamma )t}R(t_0)e^{(d+\epsilon _2+\gamma )t_0}=_{t_0}^t(\beta _2+\alpha )I(t)e^{(d+\epsilon _2+\gamma )t}>M(\beta _2+\alpha )_{t_0}^te^{(d+\epsilon _2+\gamma )t}𝑑t.$$

$$R(t)>\left(R(t_0)\frac{M(\beta _2+\alpha )}{d+\epsilon _2+\gamma }\right)e^{(d+\epsilon _2+\gamma )(tt_0)}+\frac{M(\beta _2+\alpha )}{d+\epsilon _2+\gamma }.$$

Since $`e^{(d+\epsilon _2+\gamma )(tt_0)}`$ goes to zero as $`t\mathrm{}`$, $`R(t)>{\displaystyle \frac{M(\beta _2+\alpha )}{2(d+\epsilon _2+\gamma )}}`$ for large values of $`t`$. It means that $`R(t)\mathrm{}`$. $`\mathrm{}`$

### Theorem 4.2.

(i) If $`R_1>1`$, then $`N(t)\mathrm{}`$ and if $`R_1<1`$, then $`N(t)0`$.
(ii) If $`R_2>1`$, then $`(I(t),R(t))(\mathrm{},\mathrm{})`$ and if $`R_2<1`$, then $`(I(t),R(t))(0,0)`$.

### Proof:

First suppose that $`R_01`$ which implies that $`(s,i,r)(1,0,0)`$ by Theorem 3.6. Thus $`{\displaystyle \frac{N^{}}{N}}bd`$ by (4-1) and $`{\displaystyle \frac{I^{}}{I}}\beta _1d\epsilon _1\alpha +\lambda `$ by (4-2). If $`R_1<1`$, then $`bd<0`$, hence $`N(t)0`$. Similarly if $`R_1>1`$, then $`bd>0`$, hence $`N(t)0`$. Furthermore If $`R_2<1`$, then $`B_1d\epsilon _1\alpha +\lambda <0`$, hence $`I(t)0`$ and by the above lemma $`R(t)0`$. Similarly if $`R_2>1`$, then $`\beta _1d\epsilon _1\alpha +\lambda >0`$, hence $`I(t)\mathrm{}`$ and by the above lemma $`R(t)\mathrm{}`$.

Now suppose that $`R_0>1`$. We assumed that $`I(t)>0`$, hence $`i>0`$ and $`(s,i,r)(s^{},i^{},r^{})`$ by Theorem 3.6. Thus $`{\displaystyle \frac{N^{}}{N}}b_0s^{}+(b_1+\beta \epsilon _1)i^{}+(b_2\epsilon _2)r^{}d`$ which is positive if $`R_1>1`$, hence $`N(t)\mathrm{}`$ and negative if $`R_1<1`$, hence $`N(t)0`$. Moreover $`{\displaystyle \frac{I^{}}{I}}(\beta _1d\epsilon _1\alpha )+\lambda s^{}`$ which is positive if $`R_2>1`$, hence $`I(t)\mathrm{}`$ and then $`R(t)\mathrm{}`$ by the above lemma. Similarly $`\beta _1d\epsilon _1\alpha +\lambda s^{}<0`$ if $`R_2<1`$ and then $`(I(t),R(t))(0,0)`$. $`\mathrm{}`$

We summarize our results in the following table which is the same as Table 1. in .

| $`R_0`$ | $`R_1`$ | $`R_2`$ | $`N`$ | $`(s,i,r)`$ | $`(S,I,R)`$ |
| --- | --- | --- | --- | --- | --- |
| $`1`$ | $`<1`$ | $`<1^\alpha `$ | $`0`$ | $`(1,0,0)`$ | $`(0,0,0)`$ |
| $`>1`$ | $`<1`$ | $`<1^\alpha `$ | $`0`$ | $`(s^{},i^{},r^{})`$ | $`(0,0,0)`$ |
| $`1`$ | $`>1`$ | $`<1`$ | $`\mathrm{}`$ | $`(1,0,0)`$ | $`(\mathrm{},0,0)`$ |
| $`1`$ | $`>1`$ | $`>1`$ | $`\mathrm{}`$ | $`(1,0,0)`$ | $`(\mathrm{},\mathrm{},\mathrm{})`$ |
| $`>1`$ | $`>1`$ | $`>1^\alpha `$ | $`\mathrm{}`$ | $`(s^{},i^{},r)`$ | $`(\mathrm{},\mathrm{},\mathrm{})`$ |

<sup>α</sup> Given $`R_0`$ and $`R_1`$, this condition is automatically satisfied.

Conclusion Remarks. Here we will make some comparison between the result in  and the above results.

1. If we set $`b_0=b_1=b_2=b`$ and $`\beta _1=\beta _2=\beta =0`$ in the system (2-1)-(2-3), we obtain the system (2-1)-(2-3) in . Moreover by these assumptions, we get $`R_0={\displaystyle \frac{\lambda }{b+\epsilon _1+\alpha }}`$,

$$R_2=\{\begin{array}{cc}\frac{b}{d}\hfill & \text{if}R_01,\hfill \\ \frac{b}{d+\epsilon _1i^{}+\epsilon _2r^{}}\hfill & \text{if}R_0>1,\hfill \end{array}$$

$$R_3=\{\begin{array}{cc}\frac{\lambda }{d+\epsilon _1+\alpha }\hfill & \text{if}R_01,\hfill \\ \frac{\lambda s^{}}{d+\epsilon _1+\alpha }\hfill & \text{if}R_0>1.\hfill \end{array}$$

which are the same threshold parameters as in .

2. Comparing their threshold parameters with our ones, we see that the effect of $`b_0`$, $`b_1`$ and $`b_2`$ appears more clearly instead of $`b`$. For example when $`R_0>1`$, our $`R_1`$ is $`\frac{b_0s^{}+(b_1+\beta )i^{}+b_2r^{}}{d+\epsilon _1i^{}+\epsilon _2r^{}}`$, but they obtain $`R_1=\frac{b}{d+\epsilon _1i^{}+\epsilon _2r^{}}`$ in which the effects of $`b_0`$, $`b_1`$ and $`b_2`$ are hidden in $`b`$.

3. There are two vertical transmission parameter, $`\beta _1`$ and $`\beta _1`$, in our model. The effect of $`\beta _1`$ in $`R_0`$ and $`R_2`$ is crucial. Since $`\beta _1+\beta _2=\beta `$, we can decrease $`\beta _1`$ by increasing $`\beta _2`$. Therefore removing more infected newborns causes more safe situation.

### Acknowledgments.

The author would like to thank Institute for Studies in Theoretical Physics and Mathematics for supporting this research.


text:
# Hierarchical Clustering and Active Galaxies.

## 1 The Barcelona model

The model of Barcelona has been presented by E. Salvador-Solé in his oral contribution entitled “Modeling the density-morphology relation and the galaxy/AGN connection”. A detailed description is also given in Manrique et al., 2000, in preparation. Here, only a brief introduction will be made, to outline the major differences between the abovementioned model and previous works.

The big novelty of the model is the use of a modified Press - Schechter formalism of the hierarchical dark matter halo (DMH) formation. The internal structure of the dark matter halos is modeled and a distinction is made between dark matter halo merging and accretion. Next to the dark matter evolution the baryonic matter evolution is also followed, through the mechanisms of radiative cooling, star formation and re-heating, where all component (hot and cold gas, and stars) are taken into account.

The evolution of galaxies of different types is followed, as well as their localization within groups and/or clusters. The evolution of the central galaxy depends on the properties of the host halo and of the surrounding satellite galaxies that it captures. This captures depend on the orbits of the galaxies, in other words on the potential well and the initial orbital conditions. The capture of such satellite galaxies, which can be very numerous but whose sizes are small in comparison to the central galaxy, only produces minor effects. On the contrary, in the case of a capture of a galaxy with a mass comparable to the central galaxy, the latter’s disk (if it exists) can be destroyed, forming thus a new spheroid. This kind of capture is, therefore, crucial for the final configuration of the central galaxy and for the galactic halo gas, due to feedback mechanisms.

Generally speaking, a central galaxy is characterized by its total baryonic mass, the mass of its gaseous and stellar components and the respective metalicities, the star formation rate in the disk and the bulge, the surface density of the disk and the mass of the central black hole. The destiny of the black hole and the AGN activity are related to the history of the galaxy. This si what triggers the parallel study of “normal” galaxies and AGN.

## 2 Central Black Hole Evolution and Nourishing Mechanisms

After their formation black holes evolve nourished mainly by host galaxies interactions. Galactic bulges can collect material in three different ways: by cooling flows, by the direct infall of low momentum material from the galactic halo during a merger event between two galaxies of comparable sizes, and a mass transfer from the disk to the bulge through non-axisymmetric perturbations (e.g. spiral arms and bars).

### 2.1  Characteristic Timescales

A certain number of characteristic timescales are involved in the modeling of black holes and their connection to normal galaxies. Table 1 summarizes some of them, in increasing order. Not all of them appear in the present paper but all of them are used in the above described model. “AD” and “BLR” denote the accretion disk and broad line region, respectively.

### 2.2 Galaxy Mergers in the DMH Centers

When galaxies of comparable sizes merge, usually near the center of the dark matter halos, the disks that possibly exist are destroyed and a new spheroid is formed. Due to dynamical friction the two black holes (with masses $`M_{BH}^1`$ and $`M_{BH}^2`$) of the involved galaxies will soon find their way towards the central region and coalesce. A fraction $`ϵ`$ of the cool gas, $`M_{gaz}`$, of the spheroid will fall onto the galaxy center, nourishing the black hole, whose mass will now be:

$$M_{BH}=M_{BH}^1+M_{BH}^2+ϵ(f_{}M_{}+M_{gaz}),$$

where $`f_{}`$ denotes the fraction of the stelar mass accreted, $`M_{}`$, generally considered to be null.

The typical accretion and radiation time for a black hole is nowadays believed to be much shorter than the characteristic timescales of its host galaxy. However, it is this “short-term” evolution of the black hole that determines the light curve of an active galactic nucleus, and this why it should be modeled. In the case of a rapid growth of the bulge we suppose an exponential variation of the accretion rate, inspired by Dopita (1997), but asymmetrically bell-shaped:

$$\frac{dM_{BH}}{dt}\frac{\mathrm{\Delta }M_{BH}}{t_{acc}}\left[1exp\left(\frac{t}{t_{dyn}}\right)\right]exp\left(\frac{t}{t_{acc}}\right),$$

where $`\mathrm{\Delta }M_{BH}=ϵ(f_{}M_{}+M_{gaz})`$ and $`t_{acc}=3t_{dyn}`$, as explained in table 1.

The difference between this and the Dopita curves is the use of two different timescales: $`t_{dyn}`$ for the ascending part and $`t_{acc}`$ for the descending part. This accretion rate, in units of $`t_{dyn}`$, is presented schematically in figure 1. The same equation applies in the case of a black hole nourished through the process of cooling flows.

### 2.3 Slow Fueling

In the case of a slow but continual growth of the galactic bulge as spiral arms and bars transfer angular moment towards the outer regions and mass towards the center of the galaxy, the mass of the (new) black hole will be given by:

$$M_{BH}=M_{BH}^0+\stackrel{~}{ϵ}(f_{}M_{}^{DB}+M_{gaz}^{DB}),$$

where $`M_{}^{DB}`$ and $`M_{gaz}^{DB}`$ denote the stellar and gaseous mass, respectively, transfered from the disk to the bulge. In this case a constant accretion rate is adopted:

$$\frac{dM_{BH}}{dt}=\frac{\mathrm{\Delta }M_{BH}}{t_{acc}},$$

where $`\mathrm{\Delta }M_{BH}`$ is equal, this time, to $`\stackrel{~}{ϵ}(f_{}M_{}^{DB}+M_{gaz}^{DB})`$.

## 3 AGN Light Curves

According to the present model, the bolometric luminosity of an AGN of a given type is determined by its accretion rate, $`dM_{BH}/dt`$, the initial mass of the central black hole, $`M_{BH}`$, and the time elapsed since the beginning of the accretion process. The luminosity sustains the inflow rate and vice versa, in such a way as to respect the Eddington regime, since only an under-Eddington luminosity gives a stationary solution, as described in Manrique et al. For a given accretion rate the light emitted from a region of radius $`R_{acc}`$, and which in fact is the accretion disk, varies with time as:

$$L_{BH}(t)=ϵ_{Edd}L_{Edd}=\frac{L_{Edd}}{1+\frac{L_{Edd}R_{acc}}{GM_{BH}}}\times \frac{dM_{BH}}{dt},$$

where $`L_{Edd}`$ is the Eddington luminosity.

In the case of major mergers or cooling flows, this relation will give rise to a bell-shaped light curve with the same characteristic timescales as for the accretion rate, $`dM_{BH}/dt`$.

## 4 Conclusions and Perspectives

The model presented above will allow us to make predictions on several issues that relate quasars to normal galaxies and to the background cosmology. The problem of the existence of a BH within all galaxies and the time delay between the collapse of a dark matter halo and the formation of a black hole in its center will be studied first. The relations between the different types of AGN, their connection to normal galaxies, and the relations between the masses of the dark matter halos, the galactic disks and bulges, and the black holes can then be looked for. Issues like the role of AGN and obscured AGN in the reionization history of the Universe, the X-ray and UV backgrounds can be examined.

This model will give as an output a theoretical quasar luminosity function and its evolution with time. Its comparison with the observed luminosity function issued from recent or new large quasar samples (e.g. VIRMOS, 2dF, SDSS) will allow us to better adjust the values of our parameters and tests our assumptions.

## References