text:
# Production, Supply, and Traffic Systems: A Unified Description

## 1 Modeling production systems as supply networks

Our production model assumes $`u`$ production units or suppliers $`j`$ which deliver $`d_{ij}`$ products of kind $`i\{1,\mathrm{},p\}`$ per production cycle to other suppliers and consume $`c_{kj}`$ goods of kind $`k`$ per production cycle. The coefficients $`c_{kj}`$ and $`d_{ij}`$ are determined by the respective production process, and the number of production cycles per unit time (e.g. per day) is given by the production speed $`Q_j(t)`$. That is, supplier $`j`$ requires an average time interval of $`1/Q_j(t)`$ to produce or deliver $`d_{ij}`$ units of good $`i`$. The temporal change in the number $`N_i(t)`$ of goods of kind $`i`$ available in the system is given by the difference between the inflow

$$Q_i^{\mathrm{in}}(t)=\underset{j=0}{\overset{u}{}}d_{ij}Q_j(t)$$
(1)
and the outflow

$$Q_i^{\mathrm{out}}(t)=\underset{j=1}{\overset{u+1}{}}c_{ij}Q_j(t).$$
(2)
In other words, it is determined by the overall production rates $`d_{ij}Q_j(t)`$ of all suppliers $`j`$ minus their overall consumption rates $`c_{ij}Q_j(t)`$:

$$\frac{dN_i}{dt}=Q_i^{\mathrm{in}}(t)Q_i^{\mathrm{out}}(t)=\underset{j=1}{\overset{u}{}}(d_{ij}c_{ij})Q_j(t)Y_i(t).$$
(3)
Herein, the quantity

$$Y_i(t)=\underset{\mathrm{consumption}\mathrm{and}\mathrm{losses}}{\underset{}{c_{i,u+1}Q_{u+1}(t)}}\underset{\mathrm{inflow}\mathrm{of}\mathrm{resources}}{\underset{}{d_{i0}Q_0(t)}}$$
(4)
comprises the consumption rate of goods $`i`$, losses, and waste (the “export” of material), minus the inflows into the considered system (the “imports”). In the following, we will assume that the quantities are measured in a way that $`0c_{ij},d_{ij}1`$ (for $`1ip`$, $`1ju`$) and the “normalization conditions”

$$d_{i0}=1\underset{j=1}{\overset{u}{}}d_{ij}0,c_{i,u+1}=1\underset{j=1}{\overset{u}{}}c_{ij}0$$
(5)
are fulfilled. Equations (3) can then be interpreted as conservation equations for the flows of goods.

### 1.1 Adaptation of production speeds

The production speeds $`Q_j(t)`$ may be changed in time in order to adapt to varying consumption rates $`Y_i(t)`$. We will assume that it takes a typical time period $`T`$ to adapt the actual production speeds $`Q_j(t)`$ to the desired production ones, $`W_j`$:

$$\frac{dQ_j}{dt}=\frac{1}{T}\left[W_j(\{N_i(t)\},\{dN_i/dt\},\{Q_k(t)\})Q_j(t)\right].$$
(6)
Here, the curly brackets indicate that the so-called management or control function $`W_j(\mathrm{})`$ may depend on all inventories $`N_i(t)`$ with $`i\{1,\mathrm{},p\}`$, their derivatives $`dN_i/dt`$, and/or all production speeds $`Q_k(t)`$ with $`k\{1,\mathrm{},u\}`$. The resulting dynamics of the supply network can be investigated by means of a linear stability analysis, which has been carried out for linear supply chains in Ref.  and for supply networks with $`d_{ij}=\delta _{ij}`$ in Ref. . Some of the main results will be discussed in the following.

### 1.2 Modelling one-dimensional supply chains

For simplicity, let us investigate a model of one-dimensional supply chains, here, which corresponds to $`d_{ij}=\delta _{ij}`$ and $`c_{ij}=\delta _{i+1,j}`$, where $`\delta _{ij}=1`$ for $`i=j`$ and $`\delta _{ij}=0`$ otherwise. This implies $`Q_i^{\mathrm{in}}(t)=Q_i(t)`$ and $`Q_i^{\mathrm{out}}(t)=Q_{i+1}(t)=Q_{i+1}^{\mathrm{in}}(t)`$. The assumed model consists of a series of $`u`$ suppliers $`i`$, which receive products of kind $`i1`$ from the next “upstream” supplier $`i1`$ and generate products of kind $`i`$ for the next “downstream” supplier $`i+1`$ at a rate $`Q_i(t)`$ . The final products are delivered at the rate $`Q_u(t)`$ and removed from the system with the consumption rate $`Y_u(t)=Q_{u+1}(t)`$. The consumption rate is typically subject to perturbations, which may cause a “bullwhip effect” , i.e. growing variations in the stock levels and deliveries of upstream suppliers. This is due to delays in the adaptation of their delivery rates.

The inventory of goods of kind $`i`$ changes in time $`t`$ according to

$$\frac{dN_i}{dt}=Q_i^{\mathrm{in}}(t)Q_i^{\mathrm{out}}(t)=Q_i(t)Q_{i+1}(t),$$
(7)
while the temporal change of the delivery rate will be assumed as

$$\frac{dQ_i}{dt}=\frac{1}{T}\left[W_i(N_i(t),dN_i/dt,Q_i(t),Q_{i+1}(t))Q_i(t)\right].$$
(8)
The management strategy

$$\frac{dQ_i}{dt}=\frac{1}{T}\left\{\frac{N_i^0N_i(t)}{\tau }\beta \frac{dN_i}{dt}+ϵ[Q_i^0Q_i(t)]\right\}$$
(9)
appears to be appropriate to keep the inventories $`N_i(t)`$ stationary, to maintain a certain optimal inventory $`N_i^0`$ (in order to cope with stochastic variations due to machine breakdowns etc.), and to operate with the equilibrium production rates $`Q_i^0=Y_u^0`$, where $`Y_u^0`$ denotes the average consumption rate. The supplier-independent parameter values $`\tau `$, $`\beta `$, and $`ϵ`$ can be justified with suitable scaling arguments .

### 1.3 Dynamic solution and resonance effects

In the vicinity of the stationary state characterized by $`N_i(t)=N_i^0`$ and $`Q_i(t)=Q_i^0`$, it is possible to calculate the dynamic solution of the one-dimensional supply chain model . For this, let $`n_i(t)=N_i(t)N_i^0`$ be the deviation of the inventory from the stationary one, and $`q_i(t)=Q_i(t)Q_i^0`$ the deviation of the delivery rate. The linearized model equations read

$$\frac{dn_i}{dt}=q_i(t)q_{i+1}(t)$$
(10)
and

$$\frac{dq_i}{dt}=\frac{1}{T}\left(\frac{n_i}{\tau }+\beta \frac{dn_i}{dt}+ϵq_i\right).$$
(11)
Deriving Eq. (11) with respect to $`t`$ and inserting Eq. (10) results in the following set of second-order differential equations:

$$\frac{d^2q_i}{dt^2}+\underset{=2\gamma }{\underset{}{\frac{(\beta +ϵ)}{T}}}\frac{dq_i}{dt}+\underset{=\omega ^2}{\underset{}{\frac{1}{T\tau }}}q_i(t)=\underset{=f_i(t)}{\underset{}{\frac{1}{T}\left[\frac{q_{i+1}(t)}{\tau }+\beta \frac{dq_{i+1}}{dt}\right]}}.$$
(12)
This corresponds to the differential equation for the damped harmonic oscillator with damping constant $`\gamma `$, eigenfrequency $`\omega `$, and driving term $`f_i(t)`$. The two eigenvalues of this system of equations are

$$\lambda _{1,2}=\gamma \pm \sqrt{\gamma ^2\omega ^2}=\frac{1}{2T}\left[(\beta +ϵ)\sqrt{(\beta +ϵ)^24T/\tau }\right],$$
(13)
i.e. for $`(\beta +ϵ)>0`$ their real parts are always negative, corresponding to a stable behavior in time. Nevertheless, we will identify a convective instability below, i.e. the oscillation amplitude can grow from one supplier to the next one upstream.

The set of equations (12) can be solved successively, starting with $`i=u`$ and progressing to lower values of $`i`$. For example, assuming periodic oscillations of the form $`f_u(t)=f_u^0\mathrm{cos}(\alpha t)`$, after a transient time much longer than $`1/\gamma `$ we find

$$q_u(t)=f_u^0F\mathrm{cos}(\alpha t+\phi )$$
(14)
with

$$\mathrm{tan}\phi =\frac{2\gamma \alpha }{\alpha ^2\omega ^2}=\frac{\alpha (\beta +ϵ)}{\alpha ^2T1/\tau }$$
(15)
and

$$F=\frac{1}{\sqrt{(\alpha ^2\omega ^2)^2+4\gamma ^2\alpha ^2}}=\frac{T}{\sqrt{[\alpha ^2T1/\tau ]^2+\alpha ^2(\beta +ϵ)^2}},$$
(16)
where the dependence on the eigenfrequency $`\omega `$ is important to understand the occuring resonance effect. Equations (12) and (14) imply

$$f_{u1}(t)=\frac{1}{T}\left[\frac{q_u(t)}{\tau }+\beta \frac{dq_u}{dt}\right]=f_{u1}^0\mathrm{cos}(\alpha t+\phi +\delta )$$
(17)
with

$$\mathrm{tan}\delta =\alpha \beta \tau \text{and}f_{u1}^0=f_u^0\frac{F}{T}\sqrt{(1/\tau )^2+(\alpha \beta )^2}.$$
(18)
### 1.4 “Bull-whip effect” and stop-and-go traffic

The oscillation amplitude increases from one supplier to the next upstream one, if

$$\frac{f_{u1}^0}{f_u^0}=\left\{1+\frac{\alpha ^2[ϵ(ϵ+2\beta )2T/\tau ]+\alpha ^4T^2}{(1/\tau )^2+(\alpha \beta )^2}\right\}^{1/2}>1.$$
(19)
One can see that this resonance effect can occur for $`0<\alpha ^2<2/(T\tau )ϵ(ϵ+2\beta )/T^2`$. Therefore, variations in the consumption rate are magnified under the instability condition

$$T>ϵ\tau \left(\beta +ϵ/2\right).$$
(20)
Supply chains show this “bullwhip effect” (which corresponds to the phenomenon of convective, i.e. upstream moving instability), if the adaptation time $`T`$ is too large, if there is no adaptation to the equilibrium production speed $`Q_i^0`$, corresponding to $`ϵ=0`$, or if the management reacts too strong to deviations of the actual stock level $`N_i`$ from the desired one $`N_i^0`$, corresponding to a small value of $`\tau `$. The latter is very surprising, as it implies that the strategy

$$\frac{dQ_i}{dt}=\frac{1}{T\tau }[N_i^0N_i(t)],$$
(21)
which tries to maintain a constant work in progress $`N_i(t)=N_i^0`$, would ultimately lead to an undesireable bullwhip effect. In contrast, the management strategy

$$\frac{dQ_i}{dt}=\frac{1}{T}\left\{\beta \frac{dN_i}{dt}+ϵ[Q_i^0Q_i(t)]\right\}$$
(22)
would avoid this problem, but it would not maintain a constant work in progress. The control strategy (9) with a sufficiently large value of $`\tau `$ would fulfill both requirements.

The “bullwhip effect” has, for example, been reported for beer distribution , but similar dynamical effects are also known for other distribution or transportation chains with significant adaptation times $`T`$. It has, for example, some analogy to stop-and-go traffic , where delayed adaptation also leads to an unstable behavior. In the case $`\beta =0`$ and $`ϵ=1`$, the stability condition (20) agrees exactly with the one of the optimal velocity model , which is a particular microscopic traffic model. This car-following model assumes an acceleration equation of the form

$$\frac{dv_i(t)}{dt}=\frac{V_{\mathrm{opt}}(d_i(t))v_i(t)}{T}$$
(23)
and the complementary equation

$$\frac{dd_i(t)}{dt}=[v_i(t)v_{i+1}(t)].$$
(24)
In contrast to the above supply chain model, however, the index $`i`$ represents single vehicles, $`v_i(t)`$ is their actual velocity of motion, $`V_{\mathrm{opt}}`$ the so-called optimal (safe) velocity, which depends on the distance $`d_i(t)`$ to the next vehicle ahead. $`T`$ denotes again an adaptation time. Comparing this equation with Eq. (9), the velocities $`v_i`$ would correspond to the delivery rates $`Q_i`$, the optimal velocity $`V_{\mathrm{opt}}`$ to the desired delivery rate $`W_i(N_i)=Q_i^0+(N_i^0N_i)/\tau `$, and the inverse vehicle distance $`1/d_i`$ would approximately correspond to the stock level $`N_i`$ (apart from a proportionality factor). This shows that the analogy between supply chain and traffic models concerns only their mathematical structure, but not their interpretation, although both relate to transport processes. Nevertheless, this mathematical relationship can give us hints, how methods, which have been successfully applied to the investigation of traffic models before, can be generalized for the study of supply networks. Compared to traffic dynamics, supply networks and production systems have some interesting new features: Instead of a continuous space, we have discrete production units $`j`$, and the management strategy (6) is generally different from the velocity adaptation (23) in traffic. With suitable strategies, in particular with large values of $`\tau `$ and $`\beta `$, the oscillations can be mitigated or even suppressed. Moreover, production systems are frequently supply networks with complex topologies rather than one-dimensional supply chains, i.e. they have additional features compared to (more or less) one-dimensional freeway traffic. They are more comparable to street networks of cities .

### 1.5 Calculation of the cycle times

Apart from the productivity or throughput $`Q_i`$ of a production unit, production managers are highly interested in the cycle time $`T_i`$, i.e. the time interval between the beginning of the generation of a product and its completion. Let us assume that the queue length $`L_i(t)`$ of products waiting to be processed by production unit $`i`$ is given by the inventory $`N_{i1}(t)`$ of product $`i1`$. The change of the queue length $`L_i(t)`$ in time is then determined by the difference between the arrival rate $`Q_i^{\mathrm{arr}}(t)`$, which corresponds to the inflow $`Q_{i1}^{\mathrm{in}}=Q_{i1}(t)`$ from production unit $`i1`$, and the departure rate $`Q_i^{\mathrm{dep}}(t)`$, which corresponds to the outflow $`Q_{i1}^{\mathrm{out}}(t)=Q_i(t)`$ to production unit $`i`$:

$$\frac{dL_i}{dt}=\frac{dN_{i1}}{dt}=Q_{i1}(t)Q_i(t)=Q_i^{\mathrm{arr}}(t)Q_i^{\mathrm{dep}}(t).$$
(25)
On the other hand, the waiting products move forward $`Q_i^{\mathrm{dep}}(t)=Q_i(t)`$ steps per unit time, as $`Q_i(t)`$ is the processing rate (production speed). For this reason, the overall time $`T_i(t)`$ until having been processed is given by the implicit equation

$$N_{i1}(t)=L_i(t)=\underset{t}{\overset{t+T_i(t)}{}}𝑑t^{}Q_i^{\mathrm{dep}}(t^{})=\underset{\mathrm{}}{\overset{t+T_i(t)}{}}𝑑t^{}Q_i^{\mathrm{dep}}(t^{})\underset{\mathrm{}}{\overset{t}{}}𝑑t^{}Q_i^{\mathrm{dep}}(t^{}),$$
(26)
if the queue of length $`L_i(t)`$ was joined at time $`t`$. Accordingly, high inventories imply long cycle times, which favours just-in-time production with small stock levels. From Eqs. (25) and (26), one can finally derive a delay-differential equation for the waiting time under varying production conditions :

$$\frac{dT_i}{dt}=\frac{Q_i^{\mathrm{arr}}(t)}{Q_i^{\mathrm{dep}}(t+T_i(t))}1=\frac{Q_{i1}(t)}{Q_i(t+T_i(t))}1.$$
(27)
This equation can be solved numerically as a function of the production rates $`Q_i(t)`$, since the production initially starts with a cycle time of $`T_i(0)=1/Q_i(0)`$, corresponding to the average processing time $`1/Q_i(0)`$ when the production unit $`i`$ is started. In this way, it is possible to determine the process cycle times $`T_i`$ and the overall production time as the sum of the cycle times of all single production steps, taking into account the respective time delays $`T_i`$: When a specific product enters the queue before production unit $`i`$ at time $`t=t_i`$, it enters the queue before production unit $`i+1`$ at time $`t_{i+1}=t_i+T_i(t_i)`$. The time of delivery to the customer is given by $`t_{u+1}=t_u+T_u(t_u)`$.

### 1.6 Modeling of discrete units

The above “fluid-dynamic” model equations can not only be used to represent approximate mean values of large numbers of products. They can also be transfered to the treatment of discrete units (such as single units of a product), if their dynamics is sufficiently deterministic. For example, one could represent the time interval $`\mathrm{\Delta }T`$, during which a discrete unit occupies a certain production unit $`i`$ by a step function. However, as step functions are not differentiable everywhere, we will replace them by smooth functions.

One possible specification uses a Fourier approximation of the step functions. To lowest order, one may take

$$Q_i(t)=A\underset{k}{}\{1\mathrm{cos}[\pi (tt_i^k)/\mathrm{\Delta }T]\}$$
(28)
where $`t_i^k`$ denotes the starting time of occupation by object $`k`$. The cosinus function is set to zero for $`t<t_i^k`$ and $`t>t_i^k+\mathrm{\Delta }T`$. The prefactor $`A`$ is determined in a way that satisfies the normalization condition,

$$\underset{t_i^k}{\overset{t_i^k+\mathrm{\Delta }T}{}}𝑑t^{}Q_i(t^{})=1,$$
(29)
which implies $`A=1/\mathrm{\Delta }T`$. Under this condition, the integral

$$N_i^{\mathrm{tot}}(t)=_0^t𝑑t^{}Q_i(t^{})$$
(30)
counts the number of units that have passed the cross section or production unit under consideration.

Another possible specification assumes

$$Q_i(t)=B\underset{k}{}(tt_i^k)^2(tt_i^k\mathrm{\Delta }T)^2,$$
(31)
where we set the fourth order polynomial to zero for $`t<t_i^k`$ and $`t>t_i^k+\mathrm{\Delta }T`$. The normalization condition (29) implies $`B=30/(\mathrm{\Delta }T)^5`$.

Together with (28) or (31), the relationship (30) can also be applied to situations where several units occupy a production unit at the same time. As before, the management function can be chosen as a function of the inventories $`N_i(t)`$, but in many cases, it would be reasonable to replace a reaction to temporal changes $`dN_i/dt`$ in the inventories by exponentially smoothed values.

### 1.7 Re-entrant production

Semiconductor production  relies on some extremely expensive (lithographic) production units, which are therefore used for many similar subsequent production steps $`i`$, namely the production of the various layers of a chip. This is the reason for re-entrant production  (see Fig. 1), posing particular control problems not only for the overall arrival rate $`Q^{\mathrm{arr}}(t)`$ in the re-entrant production area, but also concerning the fractions $`p_i(t)`$ of products in the different production stages $`i`$ fed into it. If $`N_i(t)`$ denotes the stock level of chips after the $`i`$th entry (production step), $`Q_i^{\mathrm{in}}(t)`$ the respective inflow, and $`Q_i^{\mathrm{out}}(t)`$ the outflow, the related balance equation is again

$$\frac{dN_i}{dt}=Q_i^{\mathrm{in}}(t)Q_i^{\mathrm{out}}(t).$$
(32)
Assuming also a buffer before the first and a buffer after the last ($`u`$th) production step, this equation applies to $`i\{0,1,\mathrm{},u\}`$, where $`Q_u^{\mathrm{out}}`$ corresponds to the removal rate of products that have completed the re-entrant (lithographic) production steps. The outflow $`Q_{i1}^{\mathrm{out}}(t)`$ determines the arrival rate $`Q_i^{\mathrm{arr}}(t)`$ for the $`i`$th re-entrant production step, while the departure rate $`Q_i^{\mathrm{dep}}(t)`$ determines the inflow $`Q_i^{\mathrm{in}}(t)`$ into the subsequent buffer. The overall stock level $`N(t)`$ of chips at various production stages in the re-entrant production area changes in time according to

$$\frac{dN}{dt}=\underset{i=1}{\overset{u}{}}[Q_{i1}^{\mathrm{out}}(t)Q_i^{\mathrm{in}}(t)]=\underset{i=1}{\overset{u}{}}[Q_i^{\mathrm{arr}}(t)Q_i^{\mathrm{dep}}(t)]=Q^{\mathrm{arr}}(t)Q^{\mathrm{dep}}(t).$$
(33)
The overall departure rate $`Q^{\mathrm{dep}}(t)`$ from the re-entrant production area depends on the number $`N(t)`$ of products processed at the same time. Moreover, if $`p_i(t)`$ denotes the fraction of chips that have entered the re-entrant area for the $`i`$th production step at time $`t`$ and $`T^0(t)`$ is the related cycle time of these chips in the re-entrant area, their departure rate is given by

$$Q_i^{\mathrm{dep}}(t)=p_i(tT^0(t))Q^{\mathrm{dep}}(N(t))=Q_i^{\mathrm{in}}(t).$$
(34)
As before, the delay-differential equation for the temporal change of the cycle time is

$$\frac{dT^0}{dt}=\frac{Q^{\mathrm{arr}}(t)}{Q^{\mathrm{dep}}(t+T^0(t))}1.$$
(35)
The overall arrival rate may be adapted according to

$$\frac{dQ^{\mathrm{arr}}}{dt}=\frac{1}{T}[W(N(t),dN/dt,Q^{\mathrm{arr}}(t),Q^{\mathrm{dep}}(t),N_u(t),dN_u/dt,Q_u^{\mathrm{out}}(t))Q^{\mathrm{arr}}(t)],$$
(36)
where the management function $`W(\mathrm{})`$ depends not only on the variables $`N(t)`$ and $`Q^{\mathrm{dep}}(t)`$ characterizing the re-entrant production area, but also on the (desired) removal rate $`Q_u^{\mathrm{out}}(t)`$ and the stock level $`N_u(t)`$ of the final buffer. Finally, the specific arrival rates of chips for the $`(i+1)`$st re-entrant production step are

$$Q_{i+1}^{\mathrm{arr}}(t)=p_{i+1}(t)Q^{\mathrm{arr}}(t)=Q_i^{\mathrm{out}}(t),$$
(37)
where their relative percentages

$$p_{i+1}(t)=\frac{W_i(N_i(t),dN_i/dt,Q_i^{\mathrm{in}}(t),Q_i^{\mathrm{out}}(t),Q_u^{\mathrm{out}}(t))}{_jW_j(N_j(t),dN_j/dt,Q_j^{\mathrm{in}}(t),Q_j^{\mathrm{out}}(t),Q_u^{\mathrm{out}}(t))}$$
(38)
are controlled by another management function $`W_i(\mathrm{})`$, which depends on the stock levels $`N_i(t)`$ of the respectively preceding buffers and on the removal rate $`Q_u^{\mathrm{out}}(t)`$. The management functions $`W_i(\mathrm{})`$ allow one to specify different priorities such as push or pull strategies. Specifications of the management function, which can avoid bullwhip (resonance) effects and long cycle times will be presented in a forth-coming paper. Generalizations to the simultaneous production of various products are possible, requiring the consideration of an individual number of re-entrant steps with potentially different production speeds and cycle times (“overtaking”) in the re-entrant area. The problem is similar to the treatment of heterogeneous multi-lane traffic . In the following, we will present a simpler traffic model for uniform vehicles, which shows how to transfer the above approach from production processes to transportation processes in street networks.

## 2 Queueing model of vehicle traffic in freeway networks

In the following, we will propose a traffic model , which was inspired by the above model of supply networks. Here, the elements to be served are the vehicles. The formula for the queue length determines the density in a road section, and the cycle time corresponds to the travel time.

When we now specify road traffic as a queueing system, we will take into account essential traffic characteristics such as the flow-density relation or the properties of extended congestion patterns at bottlenecks. In fact, traffic congestion is usually triggered by spatial inhomogeneities of the road network , and queueing effects are normally not observed along sections of low capacity, but upstream of the beginning of a bottleneck. Therefore, we will subdivide roads into sections $`i`$ of homogeneous capacity and length $`l_i^{\mathrm{max}}`$, which start at place $`x_i`$ and end with some kind of inhomogeneity (i.e. an increase or decrease of capacity) at place $`x_{i+1}=x_i+l_i^{\mathrm{max}}`$. In other words, the end of a road section $`i`$ is, for example, determined by the location of an on- or off-ramp, a change in the number of lanes, or the beginning or end of a gradient.

The model can be derived from the “fluid-dynamic” continuity equation

$$\frac{\rho (x,t)}{t}+\frac{Q_i(x,t)}{x}=\text{source terms}$$
(39)
describing the conservation of the number of vehicles. Here, $`\rho (x,t)`$ denotes the vehicle density per lane at place $`x`$ and time $`t`$, and $`Q_i(x,t)`$ the traffic flow per lane. The source terms originate from ramp flows $`Q_i^{\mathrm{ramp}}(t)`$, which enter the road at place $`x_{i+1}`$. Let us define the arrival rate at the upstream end of road section $`i`$ (the “inflow”) by $`Q_i^{\mathrm{arr}}(t)=I_iQ_i(x_i+dx,t)`$, where $`dx`$ is a differential space interval and $`I_i`$ the number of lanes of road section $`i`$. Analogously, the departure rate from the downstream end of this section is defined by $`Q_i^{\mathrm{dep}}(t)=I_iQ_i(x_{i+1}dx,t)`$. (Note that this “outflow” from section $`i`$ is to be distinguished from the outflow $`Q_{\mathrm{out}}`$ per lane from congested traffic). The conservation of the number of vehicles implies that the departure rate plus the ramp flow determine the arrival rate in the next downstream section $`i+1`$:

$$Q_{i+1}^{\mathrm{arr}}(t)=Q_i^{\mathrm{dep}}(t)+Q_i^{\mathrm{ramp}}(t).$$
(40)
In order to guarantee non-negative flows, we will demand for the ramp flows that the consistency condition $`Q_i^{\mathrm{dep}}(t)Q_i^{\mathrm{ramp}}(t)Q_{i+1}^{\mathrm{arr}}(t)`$ is always met.

Integrating the continuity equation over $`x`$ with $`x_i<x<x_{i+1}`$ provides a conservation equation for the number $`N_i(t)=_{x_i}^{x_{i+1}}𝑑xI_i\rho (x,t)`$ of vehicles in road section $`i`$. It changes according to

$$\frac{dN_i(t)}{dt}=Q_i^{\mathrm{arr}}(t)Q_i^{\mathrm{dep}}(t)=Q_{i1}^{\mathrm{dep}}(t)+Q_{i1}^{\mathrm{ramp}}(t)Q_i^{\mathrm{dep}}(t).$$
(41)
In the terminology of queueing theory, equation (41) reflects the change in the number $`N_i(t)`$ of vehicles waiting to be served by the downstream end of the section, and the number $`I_i`$ of lanes corresponds to the number of “channels” serving in parallel. Moreover, $`N_i(t)`$ corresponds to the “queue length”. However, this does not necessarily mean that we have congested traffic as, in the terminology of traffic theory, “queue length” refers to something else, namely the spatial extension $`l_i(t)>0`$ of a congested road section. In the following, we will try to express the traffic dynamics and the travel times only through the flows at the cross sections $`x_i`$, taking into account the features of traffic flow in a simplified way.

For free flow, i.e. below some critical vehicle density $`\rho _{\mathrm{cr}}`$ per lane, the relation between the traffic flow $`Q_i`$ per lane and the density $`\rho `$ per lane can be approximated by an increasing linear relationship, while above it, a falling linear relationship is consistent with congested flow-density data (in particular, if the average time gap $`T`$ is treated as a time-dependent, fluctuating variable) . This implies $`Q_i(x,t)Q_i(\rho (x,t))`$ with

$$Q_i(\rho )=\{\begin{array}{cc}Q_i^{\mathrm{free}}(\rho )=\rho V_i^0\hfill & \text{if }\rho <\rho _{\mathrm{cr}}\hfill \\ Q_i^{\mathrm{cong}}(\rho )=(1\rho /\rho _{\mathrm{jam}})/T\hfill & \text{otherwise.}\hfill \end{array}$$
(42)
Here, $`V_i^0`$ denotes the average free velocity, $`T`$ the average time gap, and $`\rho _{\mathrm{jam}}`$ the density per lane inside of traffic jams. Moreover, we define the free and congested densities by

$$\rho _i^{\mathrm{free}}(Q_i)=Q_i/V_i^0\text{and}\rho _i^{\mathrm{cong}}(Q_i)=(1TQ_i)\rho _{\mathrm{jam}}.$$
(43)
The quantity $`Q_{\mathrm{out}}=(1\rho _{\mathrm{cr}}/\rho _{\mathrm{jam}})/T`$ corresponds to the outflow per lane from congested traffic . Depending on the parameter specification, the model describes a continuous flow-density relation (for $`\rho _{\mathrm{cr}}V_i^0=Q_{\mathrm{out}}`$) or a capacity drop at the critical density $`\rho _{\mathrm{cr}}`$ and high-flow states immediately before (if $`\rho _{\mathrm{cr}}V_i^0>Q_{\mathrm{out}}`$).

According to shock wave theory , density variations at place $`x`$ propagate with velocity $`C(t)=[Q_i(x+dx,t)Q_i(xdx,t)]/[\rho (x+dx,t)\rho (xdx,t)]`$. Accordingly, the propagation velocity is $`C=V_i^0`$ in free traffic, and $`C=c=1/(T\rho _{\mathrm{max}})`$ in congested traffic. Therefore, it takes the time period $`T_i^{\mathrm{free}}=l_i^{\mathrm{max}}/V_i^0`$ for a perturbation to travel through free traffic, while it takes the time period $`T_i^{\mathrm{cong}}=l_i^{\mathrm{max}}/c`$, when the entire road section $`i`$ is congested.

Now, remember that congestion in section $`i`$ starts to form upstream of a bottleneck, i.e. at place $`x_{i+1}`$. Let $`l_i(t)`$ denote the length of the congested area and $`x(t)=x_{i+1}l_i(t)=x_i+l_i^{\mathrm{max}}l_i(t)`$ the location of its upstream front. Then, we have free traffic between $`x_i`$ and $`x_i+l_i^{\mathrm{max}}l_i(t)`$, i.e. $`Q_i(xdx,t)=Q_i^{\mathrm{arr}}(t(xx_i)/V_i^0)`$ (considering $`dx0`$), and congested traffic downstream of $`x(t)`$, i.e. $`Q_i(x+dx,t)=Q_i^{\mathrm{dep}}(t(x_{i+1}x)/c)`$. With $`dx/dt=dl_i/dt=C(t)`$ and Eq. (43) we find

$$\frac{dl_i}{dt}=\frac{Q_i^{\mathrm{dep}}(tl_i(t)/c)/I_iQ_i^{\mathrm{arr}}(t[l_i^{\mathrm{max}}l_i(t)]/V_i^0)/I_i}{\rho _i^{\mathrm{cong}}(Q_i^{\mathrm{dep}}(tl_i(t)/c)/I_i)\rho _i^{\mathrm{free}}(Q_i^{\mathrm{arr}}(t[l_i^{\mathrm{max}}l_i(t)]/V_i^0)/I_i)}.$$
(44)
The capacity of a congested road section $`i`$ is approximated as the outflow $`Q_{\mathrm{out}}=(1\rho _{\mathrm{cr}}/\rho _{\mathrm{jam}})/T`$ from congested traffic per lane times the number $`I_i`$ of lanes, minus the maximum bottleneck strength at the end of this section. This may be given by an on-ramp flow $`Q_i^{\mathrm{ramp}}(t)>0`$ or analogously by $`(I_iI_{i+1})Q_{\mathrm{out}}`$ in case of a reduction $`I_{i+1}I_i<0`$ in the number of lanes, or in general by some time-dependent value $`\mathrm{\Delta }Q_i(t)`$ in case of another bottleneck such as a gradient:

$$Q_i^{\mathrm{cap}}(t)=I_iQ_{\mathrm{out}}\mathrm{max}[Q_i^{\mathrm{ramp}}(t),(I_iI_{i+1})Q_{\mathrm{out}},\mathrm{\Delta }Q_i(t),0].$$
(45)
Analogously, the maximum capacity $`Q_i^{\mathrm{max}}(t)`$ of the road section $`i`$ under free flow conditions is given by the maximum flow $`I_i\rho _{\mathrm{cr}}V_i^0`$ minus the reduction by bottleneck effects:

$$Q_i^{\mathrm{max}}(t)=I_i\rho _{\mathrm{cr}}V_i^0\mathrm{max}[Q_i^{\mathrm{ramp}}(t),(I_iI_{i+1})\rho _{\mathrm{cr}}V_i^0,\mathrm{\Delta }Q_i(t),0].$$
(46)
Moreover, we have to specify the departure rate $`Q_i^{\mathrm{dep}}(t)`$ as a function of the respective traffic situation. Focussing on the cross section at location $`x_{i+1}`$ and considering the directions of information flow (i.e. the propagation direction of density variations), we can distinguish three different cases:

* If we have free traffic in the upstream section $`i`$ and free or partially congested traffic in the downstream section $`i+1`$, density variations propagate downstream and the departure rate $`Q_i^{\mathrm{dep}}(t)`$ at time $`t`$ is given as the arrival rate $`Q_i^{\mathrm{arr}}(tT_i^{\mathrm{free}})=Q_{i1}^{\mathrm{dep}}(tT_i^{\mathrm{free}})+Q_{i1}^{\mathrm{ramp}}(tT_i^{\mathrm{free}})`$, since the vehicles entering section $`i`$ at time $`tT_i^{\mathrm{free}}`$ leave the section after an average travel time $`T_i`$ of $`T_i^{\mathrm{free}}`$.
* In the case of partially or completely congested traffic upstream and free or partially congested traffic downstream, the departure rate $`Q_i^{\mathrm{dep}}(t)`$ is given by the capacity $`Q_i^{\mathrm{cap}}(t)`$ of the congested road section $`i`$.
* In the case of congested traffic on the entire downstream road section $`i+1`$, the departure rate $`Q_i^{\mathrm{dep}}(t)`$ is given by the departure rate $`Q_{i+1}^{\mathrm{dep}}(tT_{i+1}^{\mathrm{cong}})`$ from the downstream section at time $`tT_{i+1}^{\mathrm{cong}}`$ minus the ramp flow $`Q_i^{\mathrm{ramp}}(t)`$ entering at location $`x_{i+1}`$.

Summarizing this, we have

$$Q_i^{\mathrm{dep}}(t)=\{\begin{array}{cc}Q_i^{\mathrm{arr}}(tT_i^{\mathrm{free}})\hfill & \text{if }l_{i+1}(t)<l_{i+1}^{\mathrm{max}}\text{ and }l_i(t)=0\hfill \\ Q_i^{\mathrm{cap}}(t)\hfill & \text{if }l_{i+1}(t)<l_{i+1}^{\mathrm{max}}\text{ and }l_i(t)>0\hfill \\ Q_{i+1}^{\mathrm{dep}}(tT_{i+1}^{\mathrm{cong}})Q_i^{\mathrm{ramp}}(t)\hfill & \text{if }l_{i+1}(t)=l_{i+1}^{\mathrm{max}}.\hfill \end{array}$$
(47)
A numerical solution of the above defined section-based queueing-theoretical traffic model is carried out as follows: First, calculate the new arrival and departure rates by means of Eqs. (40) and (47), taking into account the boundary conditions for the flows at the open ends of the road network. Second, determine the queue lengths $`l_i(t)`$ in all road sections $`i`$: (i) If traffic in the road section flows freely ($`l_i(t)=0`$) and the maximum capacity $`Q_i^{\mathrm{max}}(t)`$ is not reached, i.e. $`Q_i^{\mathrm{arr}}(tT_i^{\mathrm{free}})<Q_i^{\mathrm{max}}(t)`$, we have $`dl_i(t)/dt=0`$ and the traffic flow in the road section remains free. (ii) If the road section is completely congested ($`l_i(t)=l_i^{\mathrm{max}}`$) and the arrival rate $`Q_i^{\mathrm{arr}}(t)`$ is not below the departure rate at time $`tT_i^{\mathrm{cong}}`$, i.e. $`Q_i^{\mathrm{dep}}(tT_i^{\mathrm{cong}})Q_i^{\mathrm{arr}}(t)`$, the road section $`i`$ stays fully congested and $`dl_i/dt=0`$. (iii) In other cases, we have partially congested traffic in road section $`i`$ and the length $`l_i(t)`$ of the congested area changes according to Eq. (44). Next, one continues with the first step for the new time $`t+dt`$, and so on. It is obvious, that this numerical solution is significantly more simple and robust than the numerical solution of the Lighthill-Whitham model, as shock waves (i.e. the interfaces between free and congested traffic) are treated analytically and the propagation velocities of perturbations within the free and congested regions are constant.

The travel time $`T_i(t)`$ of a vehicle that enters road section $`i`$ at time $`t`$ can be calculated analogously to Eq. (27) , i.e. via the delay-differential equation

$$\frac{dT_i(t)}{dt}=\frac{Q_i^{\mathrm{arr}}(t)}{Q_i^{\mathrm{dep}}(t+T_i(t))}1=\frac{Q_{i1}^{\mathrm{dep}}(t)+Q_{i1}^{\mathrm{ramp}}(t)}{Q_i^{\mathrm{dep}}(t+T_i(t))}1.$$
(48)
According to this, the travel time $`T_i(t)`$ increases with time, when the arrival rate $`Q_i^{\mathrm{arr}}`$ at the time $`t`$ of entry exceeds the departure rate $`Q_i^{\mathrm{dep}}`$ at the leaving time $`t+T_i(t)`$, while it decreases when it is lower. It is remarkable that this formula does not explicitly depend on the velocities on the road section, but only on the arrival and departure rates. The calculation of the travel time based on the velocity is considerably more complicated: Let $`v(t)`$ be the velocity and $`x(t)=x_i+_{t_0}^t𝑑t^{}v(t^{})`$ the location of a vehicle at time $`t`$, when it enters section $`i`$ at time $`t_0`$. Its travel time $`T_i(t_0)`$ on section $`i`$ is given by the implicit equation $`x(t+T_i(t_0))=x_{i+1}`$, which says that this vehicle reaches place $`x_{i+1}=x_i+l_i^{\mathrm{max}}`$ at time $`t_0+T_i(t_0)`$. The vehicle speed $`v(t)`$ is also difficult to determine, as it depends on the (free or congested) traffic state at its respective location $`x(t)`$: It is $`v(t)=V_i^0`$ in free flow, i.e. for $`V_i^0(tt_0)<l_i^{\mathrm{max}}l_i(t)`$. In congested flow, i.e. for $`t_0+[l_i^{\mathrm{max}}l_i(t)]/V_i^0tt_0+T_i(t_0)`$, it is determined via $`v(t)=Q_i(x(t),t)/\rho (x(t),t)=[1/\rho (x(t),t)1/\rho _{\mathrm{max}}]/T`$ with $`\rho (x(t),t)=\rho (x_{i+1},t[x_{i+1}x(t)]/c)=\{1TQ_i^{\mathrm{dep}}(t[x_{i+1}x(t)]/c)/I_i\}\rho _{\mathrm{jam}}`$.

## 3 Summary and conclusions

In this contribution, traffic and production systems have been treated in a uniform way as dynamic queueing networks, since it does not matter whether one treats the transport of goods from one production unit to the next one or of vehicles from one cross section of the road network to the next. Consequently, the stability conditions for supply chains looked similar to the ones of some specific traffic models. Moreover, we have presented a delay-differential equation for the determination of cycle (production) times, which can be also applied to calculate the travel times of vehicles. Interestingly enough, this formula did not require to calculate the vehicle speed on freeway sections, but only the in- and outflows at certain cross sections of the street network, namely, where the road capacity changed.

The derived instability conditions allow one to choose appropriate management strategies which can avoid the well-known bullwhip effect. This effect describes an amplification of variations in the delivery rate and inventory from one supplier to the next one upstream. Such a convective instability can occur despite of a stable behavior in time because of the possibility of resonance effects. Apart from this, we have sketched the treatment of re-entrant production processes and of discrete units.

The advantage of “fluid-dynamic” models of traffic, supply, and production networks is their great numerical efficiency and their consideration of non-linear interaction effects. Therefore, in contrast to most classical queueing theoretical approaches and to event-driven (Monte Carlo) simulations, they are suitable for on-line control and the treatment of variations in the consumption rate or in the production program.

Present research focusses on the effect of the topology of supply networks on their dynamics . The dynamics of production processes can even be chaotic . Other studies concentrate on the subject of optimal control (including chaos control) , which is particularly challenging for re-entrant production . Finally, the presented traffic model is now being applied to the simulation of city networks with adaptive traffic light control and to dynamic assignment problems.

### Acknowledgements

The author would like to thank Stefan Lämmer, Martin Treiber, and Thomas Seidel for their valuable comments. This work was partially supported by the German Research Foundation (DFG), grant no. He 2789/5-1.


format: text

text:
# Dark Energy, Gravitation and Electromagnetism

## 1 Introduction

Some years ago the author had put forward a model in which dark energy in the form of a background Zero Point Field would lead to an accelerating and expanding universe . Subsequently observations of distant Type 1A Supernovae by Perlmutter, Schmidt and others confirmed that the universe was indeed accelerating and expanding . However final confirmation for the dark energy contribution to this observed acceleration and expansion was awaited. This confirming evidence has now been obtained through the observations of the Wilkinson Microwave Anisotropy Probe and the Sloan Digital Sky Survey . In the light of the above observational evidence, we will now try to reexamine and characterize gravitation and electromagnetism.

## 2 Gravitation

Our starting point is the 1967 analysis of Sakharov . He starts with the space time and curvature action integral of General Relativity

$$S(R)=\frac{1}{16\pi G}(dx)\sqrt{g}R$$
(1)
He then observes that the action (1) implies a metrical elasticity of space, that is generalized forces opposing the curving of space. He then goes on to the suggestion of Zeldovich , “… the gravitational interactions could lead to a “small” dusturbance of this equilibrium and thus to a finite value of Einstein’s cosmological constant, in agreement with the recent interpretation of the astrophysical data.”
Expanding the Lagrangian density function in (1) in a power series in $`R`$, the scalar curvature, we get

$$L(R)=L(0)+Ak𝑑kR+B\frac{dk}{k}R^2+\mathrm{}$$
(2)
where $`A`$ and $`B`$ are of the order of $`1`$.
Sakharov then identifies the universal constant of gravitation from the second term in (2) as

$$G=\frac{1}{16\pi Ak𝑑k},A1$$
(3)
There is a divergence in the denominator of (3). To circumvent this Sakharov postulates a minimum length cut off at the Planck scale, that is a cut off at the Planck mass, $`10^5gm`$.
If we implement this cut off in (3) we get

$$Gm_P^21$$
(4)
The importance of the above considerations is that they show that, “…the magnitude of the gravitational interaction is determined by the masses and equations of motion of free particles, and also, probably, by the momentum cut off.” More explicitly the above analysis demotes gravitation to the analogue of elasticity in Chemical Physics, that is a mere statistical measure of residual energies (Cf.ref.). (We will further justify this in this sequel). On the other hand it must also be borne in mind that the Planck mass particle is not an elementary particle. It would at the scale of Planck dimensions decay in $`10^{42}`$seconds. Moreover Rosen has shown that it would constitute a mini universe in itself .
Nevertheless it is possible to extend the above arguments and get results in the observable realm of elementary particles by considering the background Zero Point Field (ZPF) or Dark Energy alluded to earlier (Cf.ref..
We note that as is well known, such a background ZPF can explain the Quantum Mechanical spin half as also the anomalous $`g=2`$ factor for an otherwise purely classical electron . The key point here is (Cf.ref.) that the classical angular momentum $`\stackrel{}{r}\times m\stackrel{}{v}`$ does not satisfy the Quantum Mechanical commutation rule for the angular momentum $`\stackrel{}{J}`$. However when we introduce the background Zero Point Field, the momentum now becomes

$$\stackrel{}{J}=\stackrel{}{r}\times m=\stackrel{}{v}+(e/2c)\stackrel{}{r}\times (\stackrel{}{B}\times \stackrel{}{r})+(e/c)\stackrel{}{r}\times \stackrel{}{A}^0,$$
(5)
where $`\stackrel{}{A}^0`$ is the vector potential associated with the ZPF and $`\stackrel{}{B}`$ is an external magnetic field introduced merely for convenience, and which can be made vanishingly small.
It can be shown that $`\stackrel{}{J}`$ in (5) satisfies the Quantum Mechanical commutation relation for $`\stackrel{}{J}\times \stackrel{}{J}`$. At the same time we can deduce from (5)

$$J_z=\frac{1}{2}\mathrm{}\omega _0/|\omega _0|$$
(6)
Relation (6) gives the correct Quantum Mechanical results referred to above.
From (5) we can also deduce that

$$l=r^2^{\frac{1}{2}}=\left(\frac{\mathrm{}}{mc}\right)$$
(7)
(7) shows that the mean dimension of the region in which the fluctuation contributes is of the order of the Compton wavelength of the electron.
As a further confirmation, we can similarly deduce that there is a mean time interval of the order of the electrons’ Compton time. For this we note that the energy of the ZPF is given by

$$\frac{1}{8\pi }E^2+B^2=\frac{1}{8\pi }\underset{\lambda =1}{\overset{2}{}}d^2kh^2(\omega _{\mathrm{}})$$

$$=_{k=0}^{\mathrm{}}𝑑kk^2h^2(\omega _{\mathrm{}})=_{\omega =0}^{\mathrm{}}𝑑\omega _{c^3}^{\omega ^2}h^2(\omega ),$$
(8)
where the spectral energy density is given by

$$\rho (\omega )=\frac{\mathrm{}}{2\pi ^2}\frac{\omega ^3}{c^3}$$
(9)
If now we denote the impulsive change of momentum in the interval $`\tau `$ due to the buffetting of the particle by the ZPF by

$$\mathrm{\Delta }^2^{\frac{1}{2}}$$

then we can easily deduce that by (8) and (9),

$$\mathrm{\Delta }^2=\left(r\mathrm{\Gamma }\pi ^4c^4\tau /5\omega ^2\right)\rho ^2(\omega ,T)$$
(10)
where

$$\rho (\omega ,T)=\left(\frac{\omega ^2}{2\pi ^2c^3}\right)mv^2$$
(11)
In (11) $`v`$ is the root mean square velocity of the rapidly vibrating particle, which we take to be $`c`$, the velocity of light. Then we can deduce from (10), on using the fact that the magnitude of the impulsive change of the momentum would be $`mc`$, that

$$\tau \frac{\mathrm{}}{mc^2}$$
(12)
(12) shows that $`\tau `$ of the order of the Compton time.
If we use the fact that there is the Compton time uncertainty $`\tau `$ in (3) and replace

$$G=G\left(t\right)$$

by

$$G\left(t\tau \right)$$

and further use the fact to be elaborated that

$$\dot{G}=G/t$$
(13)
We can easily see that the factor $`A`$ acquires an extra term $`\frac{At}{\tau }`$. This immediately leads to the length cut off being

$$l=l_P\sqrt{\frac{t}{\tau }}$$
(14)
(14) shows that over and above the cut off in (3) and the Planck mass in (4), there is due to the Compton scale average a further contribution at the cut off length $`l`$ which is the Compton wavelength of an elementary particle, the corresponding mass being the mass of an elementary particle.
(13) is well known . Its significance in the above context will be seen below. If we now use the extra contribution to (4) with (14) we deduce the well known otherwise, empirical electromagnetism - gravitation coupling strengths ratio

$$\frac{Gm^2}{e^2}10^{40}$$
(15)
In other words the extra contribution to the action integrals in (1) or (2) from the Compton scale ZPF effects leads to a characterization of a gravitational constant in terms of the electric charge. Looking at it another way, this points to a unified description of gravitation and electromagnetism, because (15) is no longer ad hoc or empirical, but rather is a consequence of the theory.
Let us now look at how relation (13) can be arrived at from the point of view of fluctuations in the number of particles $`N10^{80}`$ in the large scale universe. Due to this fluctuation which is of the order $`\sqrt{N}`$, the excess electromagnetic energy is given by the expression

$$\frac{e^2\sqrt{N}}{R}$$

If we now use the well known relation (Cf.ref.)

$$R=GM/c^2=\frac{GmN}{c^2}$$

in the above, the excess energy comes out as

$$\text{Excess}=\frac{e^2c^2}{Gm\sqrt{N}}$$
(16)
$`M`$ is the mass of the universe which equals $`Nm,m`$ being the mass of a typical elementary particle. (16) is the contribution due to the fluctuation in the number of particles. The question is can this energy due to the particle number fluctuation be identified with the ZPF energy within the Compton scale? Infact this latter energy is given by (Cf.ref.), independently,

$$\frac{1}{8\pi }E^2+B^2=_\omega \frac{\mathrm{}\omega ^3}{c^3}𝑑\omega \frac{\mathrm{}\omega ^4}{c^3}$$
(17)
(Cf.ref. for details). The ZPF energy (of oscillators) over an extension $`l`$ is given by the expression

$$\frac{\mathrm{}c}{l^4}$$

This is seen to be the same as (17) because,

$$\omega \frac{1}{\tau }\text{and}c\tau l$$

So the energy in a volume is $`l^3`$ is from the above

$$\text{Excess}=\frac{\mathrm{}c}{l}$$
(18)
which can be seen to be of the same order as (16) or its predecessor.
In other words the fluctuating energy of the ZPF at the Compton scale gives the entire energy of the elementary particle, on the one hand. On the other hand this equality of (16) and (18) shows that

$$G=\frac{lc^2}{m\sqrt{N}}=\frac{l^2c\tau }{mt}$$
(19)
where $`t`$ is the age of the universe. (19) immediately leads to (13).

## 3 Discussion

1. We have seen that the ZPF can be associated with the fluctuation in the number of particles in the considerations leading to (16). Another way of looking at Sakharov’s formulation in (2) and the following considerations is that if we write the Lagrangian in (2) for the $`\sqrt{N}`$ particles, then $`A`$ gets multiplied by this factor and instead of (4), we will get, as before, (15). Indeed the fact that the two approaches lead to identical conclusions is due to the fact that $`t=\sqrt{N}\tau `$, a relation which infact can be deduced from the theory (Cf.ref.), and the relation which is equivalent to the well known Eddington formula $`r=\sqrt{N}l`$.
2. It is interesting that the ZPF effects as discussed in equations (5), 7), (8), (12) and (17) bear a strong resemblance to similar zitterbewegung effects at the same scale in the Quantum Mechanical treatment of an electron, and can be identified with the latter. This is similar to the case of the ZPF at the scale of the atomic radius giving rise to the Lamb Shift.
3. The considerations leading to (7) or (12) via (8) and (10), as also the equation (17) point to the fact that it is possible to consider the mass of an elementary particle being generated by the ZPF at the Compton scale. This has been discussed in detail in (Cf.ref.).
4. The significance of the ZPF considerations vis-a-vis the curvature of space is the followintg: It is well known that the fluctuation in the scalar curvature in a region of extent $`\lambda `$ is given by

$$\mathrm{\Delta }R\frac{l_P}{\lambda ^3}$$
(20)
where $`l_P`$ is the Planck length $`10^{33}cms`$ (Cf.ref.). Usually this quantity is miniscule. When $`\lambda =l`$ the Compton wavelength of for example is electron, then it can be seen from (20) that $`\mathrm{\Delta }R`$ is not only not small but infact is of the order $`1`$. Thus we have a huge fluctuation at the Compton scale.
5. It may be mentioned that the fluctuation $`\sqrt{N}`$ referred to above, has been shown to be related via the order parameter to a Landau-Ginsburg type Phase transition (from the Quantum vacuum) .
6. In the context of the minimum Compton scale it would be interesting to consider a modified Coulomb potential proportional to $`\frac{1}{xl}`$, where $`x`$ is the radial distance from the origin. In this case we consider as in the usual theory, the Schrodinger equation which now becomes

$$u^{\prime \prime }+\left[k^2U(x)\right]u=0,$$
(21)
where $`U(r)=\frac{\beta }{xl}`$
For $`x<<l,U`$ in (21) has the expansion

$$\frac{\beta }{l}\left(1+\frac{x}{l}\right)=const\frac{\beta x}{l^2}$$

This potential agrees with the confining quark potential .
Furthermore the wave function solution of (21) is now expressible in terms of the well known Airy function, and infact in this region can be approximated by

$$u\mathrm{exp}\left(\gamma x^{3/2}\right),$$

$$\gamma \text{a \hspace{0.17em} constant}$$

This shows that the wave function and therefore the particle is confined to a region roughly of the order of the Compton wavelength $`l`$ itself. We thus get a rationale for quark confinement. For the region $`x`$ much greater than $`l`$, the potential in (21) behaves like the normal Coulumb potential.
7. Following the above remarks, it is interesting to note the following. Very much in this spirit, it was argued that  a neutrino has a small electrical charge

$$e_\nu 10^6e,$$

$`e`$ is the electron charge.
This work infact predicted that the neutrino would have a mass, a hundred millionth that of the electron, as was subsequently confirmed by the superkamiokande experiment (Cf. also. ). If we now consider a tight neutrino anti-neutrino bound state, as was also argued for the pion, we get

$$e_\nu ^2/r^2=m_\nu c^2/r$$

which leads to the correct quark Compton wavelength and mass which $`10^4`$ times the electron mass.
8. In the context of the minimum Compton length considerations, it was argued that this would lead to a non-integrable derivative

$$\frac{a^\mu }{x^\nu }\frac{a^\mu }{x^\nu }\mathrm{\Gamma }_{\lambda \nu }^\mu a^\lambda $$

$$=\frac{a^\mu }{x^\nu }\mathrm{\Gamma }_{\lambda \nu }^\mu \delta _\rho ^\lambda a^\rho $$
(22)
If further we have

$$a^\mu =^\mu Q$$

from here we can show that the net result is that there is a covariant derivative given by

$$_\nu ^\mu _\nu ^\mu \mathrm{\Gamma }_{\lambda \nu }^\mu ^\lambda $$

whence,

$$^\mu ^\mu \mathrm{\Gamma }_\nu ^{\mu \nu }$$
(23)
The equation (23) shows the emergence of the electromagnetic potential, as has been discussed in detail (Cf.ref.).
9. A problem which has vexed Physicists is the matter anti-matter asymmetry. It may be pointed out that in the light of the above considerations, it is possible to give an argument from probability theory itself. Let us suppose that the probability for the creation of particles and anti-particles is the same, as for example the probability for giving birth to a boy or a girl. However in this latter case, the probability that any family with several children has exactly half the number of boys is infact less than half. A similar argument would apply to matter and anti-matter.


text:
# Solving the Cooling Flow Problem of Galaxy Clusters by Dark Matter Neutralino Annihilation

(received Jan. 03, ’04 and accepted Mar. 25, ’04 to Physical Review Letters)

## Abstract

Recent X-ray observations revealed that strong cooling flow of intracluster gas is not present in galaxy clusters, even though predicted theoretically if there is no additional heating source. I show that relativistic particles produced by dark matter neutralino annihilation in cluster cores provide a sufficient heating source to suppress the cooling flow, under reasonable astrophysical circumstances including adiabatic growth of central density profile, with appropriate particle physics parameters for dark matter neutralinos. In contrast to other astrophysical heat sources such as AGNs, this process is a steady and stable feedback over cosmological time scales after turned on.

Diffuse thermal X-ray emission by bremsstrahlung of intracluster gas at a temperature of $``$ 10 keV has been observed from galaxy clusters for many tens of years. More than half of clusters are called cooling flow (CF) clusters, since cooling time of central cores is less than the Hubble time, and theorists predicted the existence of strong CF in such systems with a rate of $`100M_{}\mathrm{yr}^1`$. However, recent X-ray observations failed to reveal evidence of CFs, requiring that, somewhat ironically, there must be some heat source to suppress CFs in “cooling flow clusters”cf-obs . Required amount of heating is $`10^{45}`$ erg/s over a time scale of the cluster age ($`10^{10}`$ yr).

Thermal conduction is probably playing a role, especially by preventing the gas from becoming thermally unstable, but a fine tuning is necessary and thermal conduction alone does not successfully explain all clustersconduction-problem . Another heating source popularly discussed is active galactic nuclei (AGNs)AGN , but efficiency must be very high ($``$ 10% of the black hole rest mass energy)Fabian-L . Generally AGNs are intermittent activity, and accretion rate is likely determined by dynamics of small region around supermassive black holes (SMBHs). Hence it might be somewhat surprising if all clusters are kept stable over $``$ 100 kpc scale by feedback of central AGNs.

Clusters are gravitationally dominated by the cold dark matter, for which the leading candidate is the lightest supersymmetric (SUSY) particles, plausibly the neutralino $`\chi `$. The neutralino mass is limited in the range 30 GeV $`m_\chi `$ 10 TeV, and thermally averaged annihilation cross section is related to the relic density as $`\sigma \upsilon 3\times 10^{27}/(\mathrm{\Omega }_\chi h^2)\mathrm{cm}^3\mathrm{s}^1`$ (see, e.g., DM-SUSY-rev  for a review). Detectability of annihilation products from high density regions such as the Galactic center (GC) has been widely discussed (DM-SUSY-rev ; Bergstrom  and references therein). Here I consider a possibility that annihilation products may contribute to heating of intracluster gas. Exotic particle dark matter interacting with baryons has been proposed to solve the CF problemQin , but our scenario is based on theoretically better-motivated neutralino dark matter. Neutralino annihilation in galaxy clusters has been considered by ref.Cola , to explain diffuse radio halos observed in some clusters. Correlation of radio halos with merging clustersBuote , however, indicates that such halos are formed by cosmic-ray electrons produced by merger shocks. The change of central density profile by SMBH was not taken into account in ref.Cola , without which the annihilation luminosity is too small to solve the CF problem.

Since we are interested in relatively central region of a cluster, I use the following form of the dark matter density profile at $`rr_0`$: $`\rho =\rho _0(r/r_0)^\gamma `$ with $`r_0=`$ 0.5 Mpc and $`\rho _0=10^{25}\mathrm{g}\mathrm{cm}^3`$, for a typical rich cluster of $`M_{15}M_{\mathrm{cl}}/(10^{15}M_{})1`$. In latest numerical simulations of cosmological structure formation, the density profile around the center generally becomes cusps with $`\gamma `$ 1–1.5density . Here I use $`\gamma =1`$, since $`\gamma =1.5`$ is not supported by a recent X-ray observationLewis . It is easy to see that neutralino annihilation has only negligible effects on cluster energetics, when simply this density profile is applied. The contribution from the cusp to annihilation rate ($`\rho ^2`$) is convergent with $`r0`$ and unimportant when $`\gamma <1.5`$. Therefore it is a good approximation to estimate the total annihilation luminosity, $`L_{\chi \chi }`$, using a mean cluster density, $`\rho _0`$. It becomes $`L_{\chi \chi }2\sigma \upsilon \rho _0M_{\mathrm{cl}}c^2/m_\chi 2\times 10^{40}\sigma \upsilon _{26}m_2^1M_{15}\mathrm{erg}/\mathrm{s}`$, where $`m_2=m_\chi /(100\mathrm{G}\mathrm{e}\mathrm{V})`$ and $`\sigma \upsilon _{26}=\sigma \upsilon /(10^{26}\mathrm{cm}^3\mathrm{s}^1)`$. This is about five orders of magnitude lower than that required to suppress the CF.

However, the situation drastically changes if there is a density “spike” associated with a central SMBH. Adiabatic growth of a SMBH at the center of a preexisting halo produces a spike in the density profile within $`rr_s`$ from the center, which is even steeper than original cuspsGS99 . The power-law index of spike density profile, $`\gamma _s`$, is related to that of the original cusp as: $`\gamma _s=(92\gamma )/(4\gamma )=`$ 2.33–2.4 for $`\gamma =`$ 1–1.5. The spike radius $`r_s`$ is a radius within which the enclosed mass of the halo is the same as the central SMBH. Since $`\gamma _s>1.5`$, the volume integration of the annihilation rate is divergent with $`r0`$, and hence a huge enhancement of $`L_{\chi \chi }`$ is possible from the very central part of the halo.

If this enhancement is happening in the center of our Galaxy, much stronger flux of annihilation gamma-rays and cosmic-rays are expected than previously thought, and already existing experimental/observational limits exclude a considerable SUSY parameter spaceGS99 . However, several authors have argued that such a process is rather unlikely to occur in the GCUllio-spike ; Merritt . The initial SMBH mass before adiabatic growth should be much smaller than the final one, and it must be placed to the dynamical center of the halo, but the dynamical friction time for that may be larger than the hubble time. The mass density of the GC is dominated by baryons rather than by dark matter, and violent processes such as star formation and supernova explosions might disrupt cold orbits of dark matter particles required for the spike formation. Mergers between halos containing SMBHs may lead to a flatter density profile than original cusps.

However, the center of galaxy clusters, especially those having cooling cores within the Hubble time, appear to be the best site for the spike growth to happen. All CF clusters seem dynamically well evolved systems and a single giant cD galaxy is placed on the gravitational center where X-ray surface brightness peaksJones-Forman . Recent Chandra observations (e.g., Lewis ) have shown that the dark matter dominates the mass density with a profile consistent with $`\gamma =1`$ down to $``$ kpc, and probably further down to the SMBH scale. If strong CF with the theoretically predicted accretion rate ($`10^2M_{}\mathrm{yr}^1`$) occurred at some epoch in early cluster evolution and this accretion is used to feed the SMBH of a cD galaxy, the SMBH mass could grow to $`M_{,10}M_{}/(10^{10}\mathrm{M}_{})1`$ within $`10^8`$yr depending on the efficiency of mass accretion on to the SMBH. This is much shorter than typical cluster age, while the orbital period at the spike radius $`r_s=1.5M_{,10}^{1/2}`$ kpc is $`5.7\times 10^7M_{,10}^{1/4}`$ yr, which is shorter than the SMBH growth time scale and hence it satisfies a requirement for the adiabatic growth. The final SMBH mass is much larger than initial mass that is probably similar to typical SMBH mass in normal galaxies ($`10^{69}M_{}`$), and hence another condition of adiabatic growth mentioned by Ullio-spike  is satisfied.

Thus, it seems reasonable to suppose that the adiabatic spike growth occurred in the past, within $``$ kpc of the cluster center where even the latest X-ray satellites cannot resolve the density profile. Dominant energy production by annihilation occurs in the central core, where the core density $`\rho _c`$ is limited by annihilation itself over typical cluster age, $`t_{\mathrm{cl}}10^{10}t_{10}`$ yr, as: $`\rho _c\sigma \upsilon t_{\mathrm{cl}}/m_\chi 1`$. Equating this core density and the spike density profile, $`\rho =\rho _0(r/r_s)^{\gamma _s}(r_s/r_0)^\gamma `$, I find the core radius $`r_c=0.17M_{,10}^{2/7}m_2^{3/7}\sigma \upsilon _{26}^{3/7}t_{10}^{3/7}`$ pc. (It is much larger than the Schwarzschild radius of the SMBH, $`0.95M_{,10}`$ mpc.) The annihilation luminosity within this core radius is given by:

$`L_{\chi \chi }`$ $`=`$ $`2m_\chi c^2\sigma \upsilon \left({\displaystyle \frac{\rho _c}{m_\chi }}\right)^2\left({\displaystyle \frac{4\pi }{3}}r_c^3\right)`$
$`=`$ $`1.9\times 10^{44}M_{,10}^{6/7}m_2^{2/7}\sigma \upsilon _{26}^{2/7}t_{10}^{5/7}\mathrm{erg}\mathrm{s}^1.`$
This is a rate within $`r_c`$, and adding integration at $`r>r_c`$ increases the rate by a factor of 2.8. The above equation suggests a modest time evolution of $`L_{\chi \chi }t^{5/7}`$, and time average over the cluster age is increased by a factor of 7/2. Then I finally obtain $`L_{\chi \chi }10^{45}`$ erg/s, which is very close to the number required to heat the cooling cluster cores. In contrast to other astrophysical heating sources, this process is stable with only mild time evolution once it is turned on, unless the spike is disrupted by violent events such as major mergers. When such events destroy the spike, the cooling core may also be destroyed, changing a CF cluster into the other category of non-CF ones.

It must be examined whether this energy production is efficiently converted to the heat of intracluster gas. The annihilation products are eventually converted to stable particles. I used the DarkSUSY packageDS  to calculate the amount and energy spectrum of these annihilation yields, and find that, with only weak dependence on SUSY parameters, about 1/4, 1/6 and 1/15 of the total annihilation energy goes to continuum gamma-rays, $`e^\pm `$’s, and $`p\overline{p}`$’s, respectively. The other energy goes to neutrinos, which are not useful for heating. The spectral energy distribution of particles per logarithmic interval, $`ϵ^2dN/dϵ`$, peaks at about 0.05, 0.05, and 0.1 times $`m_\chi c^2`$, in the same order. Therefore most of the annihilation energy will be carried away by particles of about 5 GeV, for example, for $`m_\chi `$ 100 GeV.

First I consider the fate of electrons/positrons. For simplicity, I assume that a fraction $`f_\pm 1`$ of the total annihilation energy is given to $`e^\pm `$ and all $`e^\pm `$’s have the same energy of $`ϵ_0ϵ_\pm /(1\mathrm{GeV})1`$. They are produced at very dense environment, and their relativistic motion would result in much higher pressure than the environment. Therefore $`e^\pm `$’s would expand until their pressure becomes comparable with the intracluster pressure. Buoyancy may result in intermittent formation of bubbles of relativistic particles. The $`e^\pm `$ density in the bubbles can be written by external pressure, as: $`n_\pm P_{\mathrm{ext}}/ϵ_\pm 6.3\times 10^7P_9ϵ_0^1\mathrm{cm}^3`$, where $`P_9P_{\mathrm{ext}}/(10^9\mathrm{erg}\mathrm{cm}^3)1`$. This value and gas density $`n_1n/(0.1\mathrm{cm}^3)1`$ are taken from observed values within $`r`$ 1–10 kpcLewis .

Heating of intracluster gas by cosmic ray electrons produced by AGNs has been discussed in literatureLea-Holman ; Scott , and similar treatments can be applied to estimate the energy loss time scale of $`e^\pm `$’s, $`\tau _\pm `$. The upper limit on $`\tau _\pm `$ is given by the ordinary Coulomb collisions with background gas: $`\tau _{\pm ,\mathrm{cc}}=5.1\times 10^8n_1^1ϵ_0`$ yr. This should be compared with radiative energy loss time scale by inverse-Compton scattering (ICS) of the cosmic microwave background (CMB) photons: $`\tau _{\pm ,\mathrm{ic}}=1.2\times 10^9ϵ_0^1`$ yr, and synchrotron radiation: $`\tau _{\pm ,\mathrm{sync}}=\tau _{\pm ,\mathrm{ic}}(B/B_{\mathrm{CMB}})^2`$, where $`B_{\mathrm{CMB}}=3.30\mu \mathrm{G}`$. I also found that stellar radiation energy density in cD galaxies is comparable to that of CMB, estimating it by $`U_{\mathrm{st}}L_{\mathrm{cD}}/(4\pi r_{\mathrm{cD}}^2c)`$, where $`L_{\mathrm{cD}}`$ and $`r_{\mathrm{cD}}`$ are typical observed stellar luminosity and spatial size, respectivelyLewis . Therefore ICS of stellar photons should also have a comparable effect. Bremsstrahlung loss time scale is $`\tau _{\pm ,\mathrm{br}}5.7\times 10^8n_1^1`$ yr and annihilation time scale of positrons is $`1.0\times 10^{10}n_1^1ϵ_0`$ yr. Comparing these time scales, it can be seen that a considerable part of $`e^\pm `$ energy can be converted into thermal energy.

Furthermore, even more efficient energy loss is possible when collective effects of plasma, such as the relativistic two-stream instability, are importantScott . Injection of huge amount of relativistic $`e^\pm `$’s within the small core radius $`r_c`$ would lead to strong wind of these particles, making their momentum distribution strongly anisotropic, which is necessary for the collective effect. Using formulae given in these references, I found $`\tau _{\pm ,\mathrm{tsi}}=3.1\times 10^3(n_\pm /n)^2(n/\mathrm{cm}^3)^{0.5}ϵ_0^2\mathrm{s}=8.3\times 10^6P_9^2n_1^{1.5}ϵ_0^4`$ yr. Dependence on $`ϵ_\pm `$ is large, and almost all energy could be quickly converted to thermal energy when $`ϵ_\pm `$ 1 GeV.

The annihilation energy given to $`p\overline{p}`$’s would also be partially converted into thermal energy by at least the rate of Coulomb collisions that is similar to $`\tau _{\pm ,\mathrm{cc}}`$. Some fraction of energy may be lost by inelastic interaction with ambient thermal protons, with a time scale of $`\tau _{\mathrm{pp}}(n\sigma _{\mathrm{pp}}c)^1=3.3\times 10^8n_1^1`$ yr, and secondary $`e^\pm `$’s produced by this interaction would be again used for heating. It is likely that primary annihilation gamma-rays do not contribute much to heat the cooling gas, because the optical depth to Compton scattering is $`10^31`$. However, if there are dense clouds having large optical depth around the annihilation core, the Compton heating may also have considerable effect, as considered for heating by AGNs (e.g., Ciotti ).

What is the preferred value of the neutralino mass in this context? The peak of $`ϵ^2dN/dϵ`$ should be less than $``$ a few GeV, otherwise the energy loss time scale to heating ($`\tau _{\pm ,\mathrm{cc}}`$ or $`\tau _{\pm ,\mathrm{tsi}}`$) becomes longer than that for radiative loss, leading to inefficient heating. Combined with the relation $`L_{\chi \chi }m_\chi ^{2/7}`$, the neutralino mass should be $``$ 100 GeV for the proposed process to efficiently work.

Now I discuss the observability of any signature of the neutralino annihilation. Typically about $``$30 continuum gamma-rays are produced at $`ϵ_\gamma >`$ 100 MeV per annihilationDS , and expected gamma-ray flux becomes $`F_\gamma (>100\mathrm{M}\mathrm{e}\mathrm{V})7\times 10^8L_{45}m_2^1d_2^2\mathrm{cm}^2\mathrm{s}^1`$, where $`L_{45}L_{\chi \chi }/(10^{45}\mathrm{erg}\mathrm{s}^1)`$ and $`d_2=d/(100\mathrm{M}\mathrm{p}\mathrm{c})`$ is a typical distance to nearby rich clusters. This is just almost the same as the EGRET sensitivity limit, and hence the prediction is marginally consistent with no reported gamma-ray detection from nearby galaxy clustersReimer . In fact, there are many positional coincidences between known galaxy clusters and unidentified EGRET sources, and the detection from clusters has not yet been claimed because of low statistical significance. The next generation gamma-ray satellite, GLAST will very likely detect the continuum gamma-ray flux as steady and point sources at cluster centers.

Even if the continuum gamma-rays may be discriminated from other astrophysical sources by spectrum, variability, and/or extension, a conclusive evidence would come from line gamma-rays. Line flux (in $`\mathrm{cm}^2\mathrm{s}^1`$) is expected to be smaller than the continuum flux by a factor of about $`15\times 10^3`$, where the former factor (30/2) comes from the ratio of photon number produced per annihilation and the latter comes from the branching ratio into $`\chi \chi \gamma \gamma `$ or $`Z\gamma `$ modesBergstrom . Following the line sensitivity estimate given in this reference, I found that the line flux expected at $``$ 100 GeV may produce several photons from a cluster center for five-year operation of GLAST, compared with negligibly small background rate of $`10^3`$ events within angular and energy resolution. Co-added analysis of many CF clusters would even increase the sensitivity. Future air Cerenkov telescopes may have even better sensitivity for the line flux, but the threshold energy must be lower than $``$ 100 GeV since we expect $`m_\chi `$ 100 GeV.

Since the flux enhancement by the density spike is so drastic, the best target to search neutralino annihilation may be nearby CF clusters, rather than the GC or nearby galaxies. On the other hand, $`L_{\chi \chi }`$ without the density spike ($`10^{40}`$ erg/s) is below the sensitivity limit even for the next generation gamma-ray telescopes.

Equally intriguing is the detectability of synchrotron radiation in the radio bands, by $`e^\pm `$ pairs produced by annihilation. It is known that central cD galaxies in CF clusters have a higher probability of becoming radio sources than those in non-CF clusters, and there seems a correlation between radio power and CF strength radio , suggesting that the annihilation $`e^\pm `$ may be responsible for radio emission. Typical radio luminosity from CF cluster cores is $`10^{4142}`$ erg/sradio , which is about $`10^3`$ times lower than the cluster X-ray luminosity or $`L_{\chi \chi }`$. This factor can be explained by a few effects as follows. The magnetic field strength is poorly known, which may be weakened by relativistic bubble formation. Hard X-ray tails observed in some cluster spectra are popularly interpreted as ICS of CMB photons by cosmic-ray electrons, and their luminosity is a few hundreds times larger than diffuse radio halo luminosity, indicating that synchrotron is inefficient process compared with ICSICS-vs-sync . The energy loss by heating can be much more efficient than radiative loss by the collective effect, and then the total radio flux will be further reduced. Characteristic synchrotron frequency is $`\nu _{\mathrm{sync}}0.011ϵ_0^2(B/\mu \mathrm{G})`$ GHz. Since the injection of relativistic $`e^\pm `$ occurs into relatively narrow energy range around $`0.05m_\chi c^2`$ compared with broad power-law spectrum of cosmic-ray electrons, most of radio emission may occur at $``$ GHz, which is out of typical observing frequencies.

Here I give a size estimate of the radio emitting region. It is expected that the size of $`e^\pm `$ bubble, $`r_b`$, is determined by the pressure balance and the lifetime of these particles, as: $`r_b=(3\tau _\pm \dot{N}_\pm /4\pi n_\pm )^{1/3}`$, where $`\dot{N}_\pm =f_\pm L_{\chi \chi }/ϵ_\pm `$ is the $`e^\pm `$ production rate. Using $`\tau _{\pm ,\mathrm{tsi}}`$ for the lifetime of $`e^\pm `$’s, I found $`r_b=13n_1^{1/2}ϵ_0^{4/3}P_9^1f_\pm ^{1/3}L_{45}^{1/3}`$ kpc. Typical extension of the radio emission associated with cD galaxies is $``$ 5–10 kpcradio . Large scale radio halos ($``$ Mpc) are found only in non-CF clusters like Coma, but CF clusters often have less extended “minihalos” ($``$ 100 kpc) in which a strong radio cD galaxy is centered radio-minihalo . It may be speculated that the $`e^\pm `$ bubble corresponds to the central strong radio source, while radio minihalos are made by leaking $`e^\pm `$ from the bubble, perhaps by diffusion. The morphology of such radio galaxies and minihalos in CF clusters is poorly collimated and more spherical, compared with clear bipolar jet-like structures generally found in radio galaxiesradio . Such a trend may be difficult to explain if CF is suppressed by AGN jets, while it is naturally understood if the steady, isotropic energy production by neutralino annihilation is responsible. If this interpretation is true, it means that we have already observed, though indirectly, the dark matter for many tens of years!

Annihilation $`e^\pm `$ pairs would also produce a similar luminosity to synchrotron radiation in the X-ray and MeV bands by ICS of CMB and stellar photons, respectively, but they are difficult to detect because of the strong thermal X-ray emission and difficulty of MeV observations.

The simple estimate of $`L_{\chi \chi }`$ without the density spike, $`10^{40}`$ erg/s, is comparable to the power of diffuse radio halos in merging clusters and hence neutralinos may be the origin, as pointed out by ref.Cola . A necessary condition for this is that the synchrotron must be the dominant energy loss process. However, as argued above (and refs. ICS-vs-sync ), observed hard X-ray tails in cluster spectra indicate that synchrotron is rather inefficient process, and diffuse synchrotron halo by annihilation $`e^\pm `$’s is likely still under the current detection limit.

It is expected that stars in cD galaxies are also affected by the adiabatic growth, and it might be interesting to seek for any signature in the central luminosity density profile of cD galaxies by high resolution optical observations. In fact, adiabatic growth by SMBHs is one of the proposed explanations of cusps seen in central surface brightness profile of elliptical galaxiesvandenMarel . However, there is no reason to believe that stellar density profile should be the same as that of dark matter, and large elliptical galaxies generally have flat stellar density coresvandenMarel . A stellar density spike should be weak and may not be detectable, if it is formed from such a flat density core.

###### Acknowledgements.

This work has partially been supported by a Grant-in-Aid for the COE program for physics at Kyoto University.


text:
# Some aspects of radiative corrections and non-decoupling effects of heavy Higgs bosons in two Higgs Doublet Model

## 1 Introduction

Although the two-Higgs-doublet extension of the standard model (THDM) was invented about 30 years ago , it still belongs among viable candidates for a theory beyond the electroweak standard model (SM). Despite its simplicity it is quite popular, namely because of its capability to include various aspects of “new physics” like for example the additional sources of CP violation (see e.g. , ). Moreover, its two Higgs doublet structure mimics many features of the Higgs sector of perhaps the most popular SM extension, the minimal supersymmetric standard model (MSSM). On the other hand, since the Higgs sector of THDM is less constrained, it can lead to various effects which are not present in MSSM, in particular to the non-decoupling behaviour of the heavy Higgs boson contributions in the electroweak scattering amplitudes.

As in the MSSM, the presence of the additional doublet leads to five physical Higgs states: 2 CP even Higgs scalars $`h^0`$ and $`H^0`$, a CP odd pseudoscalar $`A^0`$ and a charged pair of $`H^\pm `$. The lightest scalar $`h^0`$ is quite similar to the SM Higgs boson $`\eta `$ i.e. the mass of $`h^0`$ should be close to the weak scale. On the other hand the typical mass scale of the other Higgses ($`M_H`$) is not so constrained in general, the unitarity bounds  permit $`M_H`$ around one TeV (if there is no new physics in the game at this scale). Therefore a natural question arises as to whether these additional Higgs bosons tend to decouple from the weak-scale amplitudes. As we shall see in the next section, the answer is ’not in general’.

## 2 Non-decoupling of heavy Higgs bosons in THDM

The reason why the heavy Higgs bosons need not decouple from the weak-scale physics in the THDM, but they do so within MSSM  is roughly the following. Since the Higgs self-couplings are driven by SUSY, the only way to make the four additional Higgs bosons ($`H^0`$, $`A^0`$ and $`H^\pm `$) sufficiently heavy in MSSM is to adjust the $`SU(2)_LU(1)_Y`$ singlet mass parameters in the Higgs potential; in such case these masses have to decouple in accord with the famous Appelquist-Carazzone theorem . In the THDM case one can do the job also by a convenient choice of the Higgs couplings $`\lambda _i`$ and the SSB parameter $`\mathrm{tan}\beta `$, keeping at the same time the singlet mass parameters small. Notice that even the violation of the simple unitarity bounds could be fully compatible with the requirement of perturbativity of the Higgs sector ($`\lambda _i<1`$) provided one can choose a sufficiently large value of $`\mathrm{tan}\beta `$. As an illustration, consider the following tree-level THDM Higgs mass relations:

$`m_{h^0}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1\kappa )M^2+\left[B_2\mathrm{sin}^2\beta A_1\mathrm{cos}^2\beta +{\displaystyle \frac{1}{4}}C(1+\mathrm{cos}2\alpha \mathrm{cos}2\beta )\right]{\displaystyle \frac{v^2}{\mathrm{cos}2\alpha }}`$
$`m_{H^0}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1+\kappa )M^2+\left[A_2\mathrm{cos}^2\beta B_1\mathrm{sin}^2\beta {\displaystyle \frac{1}{4}}C(1\mathrm{cos}2\alpha \mathrm{cos}2\beta )\right]{\displaystyle \frac{v^2}{\mathrm{cos}2\alpha }}`$
$`m_{A^0}^2`$ $`=`$ $`M^2{\displaystyle \frac{1}{2}}\left[2\lambda _5^R+\lambda _6^R\mathrm{cot}\beta +\lambda _7^R\mathrm{tan}\beta \right]v^2`$ (1)
$`m_{H^\pm }^2`$ $`=`$ $`M^2{\displaystyle \frac{1}{2}}\left[\lambda _4+\lambda _5^R+\lambda _6^R\mathrm{cot}\beta +\lambda _7^R\mathrm{tan}\beta \right]v^2`$
where (using the superscript $`R`$ to denote the real part of a quantity)

$`M^2{\displaystyle \frac{m_{12}^{2}{}_{}{}^{R}}{\mathrm{sin}\beta \mathrm{cos}\beta }}`$ $`\kappa {\displaystyle \frac{\mathrm{cos}2\beta }{\mathrm{cos}2\alpha }}`$
$`A_1\lambda _1\mathrm{sin}^2\alpha \lambda _7^R\mathrm{tan}\beta \mathrm{cos}^2\alpha `$ $`B_1\lambda _2\mathrm{sin}^2\alpha \lambda _6^R\mathrm{cot}\beta \mathrm{cos}^2\alpha `$
$`A_2\lambda _1\mathrm{cos}^2\alpha \lambda _7^R\mathrm{tan}\beta \mathrm{sin}^2\alpha `$ $`B_2\lambda _2\mathrm{cos}^2\alpha \lambda _6^R\mathrm{cot}\beta \mathrm{sin}^2\alpha `$
$`C\lambda _7^R\mathrm{tan}\beta \lambda _6^R\mathrm{cot}\beta `$ $`D\lambda _7^R\mathrm{tan}\beta +\lambda _6^R\mathrm{cot}\beta `$
The model and notation are those used in . Notice that in the case $`\lambda _6=\lambda _7=0`$ one recovers the relations obtained in . Moreover, using

$$\mathrm{cos}^2(\alpha \beta )=\frac{m_{h^0}^2m_L^2}{m_{H^0}^2m_{h^0}^2}\mathrm{and}m_L^2\lambda _1\mathrm{cos}^4\beta +\lambda _2\mathrm{sin}^4\beta +\frac{1}{2}\left(\lambda +2D\mathrm{sin}^2\beta \right)v^2$$
(2)
it is easy to se that if the weak-scale contributions (square-brackets in (2)) are small compared to $`M^2`$, the requirement of having $`h^0`$ light and the others much heavier forces the heavy multiplet to be almost degenerate with masses proportional to $`M`$, which is the signature of the so-called decoupling regime . Therefore, it is the distortion of the heavy Higgs spectrum which matters concerning the possible nondecoupling effects of the additional Higgses in THDM.

In this work I would like to demonstrate these issues at the particular case of the amplitude of the proces $`e^+e^{}W^+W^{}`$ at one-loop level in THDM in comparison with the well-known one-loop SM result . Note that there are already earlier papers on this topic in the literature , but these usually make use of some specific approximations (in particular, the equivalence theorem for longitudinal vector bosons ) which we would like to avoid.

## 3 The process $`e^+e^{}W^+W^{}`$

For the considered process, the central quantity of our interest is the deviation of the differential cross-section, calculated within THDM, from its SM value; this is defined by

$$\delta \mathrm{d}\sigma ^{THDM}/\mathrm{d}\sigma ^{SM}1$$
(3)
Expanding the THDM amplitude around the SM value and keeping just the leading terms, one gets

$$\delta 2\mathrm{R}\mathrm{e}\frac{\mathrm{\Delta }[\mathrm{\Delta }\mathrm{\Gamma }^{TGV}]_{1loop}}{_{tree}^{SM}}$$
(4)
Here $`\mathrm{\Delta }_{1loop}`$ stands for the difference of the THDM and SM 1-loop amplitudes, which descends primarily (the leading term) from the differences of the triple gauge vertex corrections $`\mathrm{\Delta }\mathrm{\Gamma }^{TGV}`$:

$$\mathrm{\Delta }[\mathrm{\Delta }\mathrm{\Gamma }^{TGV}]_{1loop}=\left[\text{}\right]_{THDM}\left[\text{}\right]_{SM}$$

Since most of the technical aspects of the calculation are covered in , let us emphasize only several salient points. i) We have chosen to work in the on-shell renormalization scheme. There are two main reasons for that: the overall number of diagrams to be calculated is reduced with respect to other schemes and the mass-parameters we are playing with are the true physical masses. The only disadvantage is the need of treating carefully the finite parts of the counterterms which must be computed by means of Ward identities. On the other hand, the cancellation of UV-divergences provides a non-trivial consistency check. ii) There is also a simple consistency check for the finite parts of $`\mathrm{\Delta }\mathrm{\Gamma }^{TGV}`$: they should tend to vanish in the decoupling regime, i.e. in the case where the masses of heavy Higgs bosons are large and almost degenerate.

## 4 Summary of results and conclusion

Due to the large number of diagrams contributing to $`\mathrm{\Delta }\mathrm{\Gamma }^{TGV}`$ it is hard to get an analytic expression even for the leading terms in $`\mathrm{\Delta }_{1loop}`$. The numerical analysis shows that the formfactors $`\mathrm{\Delta }\mathrm{\Pi }^{\gamma WW}`$ and $`\mathrm{\Delta }\mathrm{\Pi }^{ZWW}`$ defined in  behave in accordance with the consistency conditions mentioned above. For example, let us look at $`|\mathrm{\Delta }\mathrm{\Pi }_1^{\gamma WW}|`$ as a function of the mass of the $`A^0`$, (fig.1): since the other Higgs masses are kept close to the weak scale, the heavy Higgs spectrum distortion grows with $`m_{A^0}`$ and the non-decoupling effect in the formfactor as well. Concerning $`\delta `$, one naturally expects a similar behaviour because it is linear in the formfactors (at the leading order, see ). Let us take the particular case: $`e_L^+e_R^{}W_L^+W_L^{}`$ (in this setup the leading term turns out to be $`\mathrm{cos}\theta ^{}`$-independent which allows us to draw simpler pictures). As can be seen in fig.2, for large distortions of the heavy Higgs spectrum one can get an effect of several percent. At least in principle, the nondecoupling effects of relatively heavy additional Higgs degrees of freedom can be used in an indirect exploration of the EW Higgs sector at future colliders.

Acknowledgements: The work was partially supported by ”Centre for Particle Physics”, project No. LN00A006 of the Czech Ministry of Education. I would like to thank Prof. Jiří Hořejší for useful discussions. I am grateful to S.I.S.S.A. for the financial support and the organizers for the possibility to take part at this nice event.


text:
# Exploration of Scale-Free Networks

## 1 Introduction

In recent years networks have become one of the most promising frameworks to describe systems as diverse as the Internet and the WWW, email and social communities, distribution systems, food-webs, protein interaction, genetic and metabolic networks AB02 . The collected data have allowed the discovery of many important properties: in particular two of them have become prominent, namely the small-world WS98  and scale-free features BA99 . Small-world implies that the average distance between nodes of the network increases at most logarithmically with the number of nodes, and formalizes the concept of ”six degrees of separation” typical in social contexts. Scale-free refers to the lack of an intrinsic scale in some of the properties of the network. In particular, the quantitiy that has been most thoroughly studied is the degree (or connectivity) distribution: the degree $`k`$ of a node is the number of other nodes it has links to (here we do not distinguish between directed and undirected links), and the degree distribution $`P(k)`$ is simply the histogram of the number of nodes with a given degree $`k`$. Scale-free networks exhibit a power-law behavior of the distribution $`P(k)k^\gamma `$, with $`\gamma `$ values often between $`2`$ and $`3`$ AB02 . The small-world and scale-free properties turned out being quite ubiquitous and some general, qualitative, models to explain their appearence have been put forward. At the same time various versions of these models have been also proposed in order to capture also the detailed values of some quantities, such as the exponent $`\gamma `$. Yet, as this new field is slowly coming of age, and as a consequence it is also becomeing more quantitative, an analysis of the data, and of their reliability, is due. The main problem that should be addressed is whether the data we are using have been skewed somehow by the detection method. In the lack of such analysis on real data and methods, we propose to work on syntetic models and data to explore their robustness in some simple test case. In the next section we address the tree-like exploration technique and discuss how it can bias the measurements, and in the third section we show that a random graph can be distorted by the exploration so to look like a SF one: in this case the exponent $`\gamma `$ is completely spurious.

## 2 Tree-like Exploration of Scale-Free Networks

Scale-free (SF) networks can be explored in many different ways. One of the most popular methods, that has been extensively used for example for the Internet, is a sort of tree-like exploration implemented by the recursive use of the traceroute command. In short, traceroute finds a path (usually a short one, but not necessarily the shortest) from the node where the command is executed to another given node. By repeating the procedure asking traceroute to find paths to all other possible nodes (addressed by their IP number), one ends up with a representation of the Internet that shows just a small amount of loops. This is due to the fact that traceroute mostly uses the same paths: if a node $`D`$ can be reached from $`A`$ through both $`B`$ and $`C`$, traceroute most of the times detects only one of them. Actually, chances are that traceroute can find more than a single path if traffic over an already discovered one is so high that it becomes more convenient to switch to a different path. Data collected with this technique have shown that degrees in the Internet are distributed according to a power-law with exponent $`\gamma 2.2\pm 0.1`$ MCP . In order to analyze the effects of a tree-building exploration algorithm on SF networks, we have synthesized our own networks according to two different models: the Barabasi-Albert (BA) model BA99 , and the hidden variable model CCDM02 .

The BA model describes the growth of a network as new nodes are added at a constant rate, and they connect to older nodes in the network according to the preferential attachment rule. Preferential attachment means that an old node has a probability proportional to its degree of aquiring a connection from a new one. It is useful to recall a simple derivation of the degree distribution starting from these two simple rules, growth and preferential attachment. The rate of change of the degree $`k_i`$ of node $`i`$ is

$$\frac{dk_i(t)}{dt}=m\frac{k_i(t)}{2mt}$$
(1)
where $`m`$ is the number of connections that a new node establishes with older ones, and the denominator in the right hand side of Eq.1 represents the sum over all the degrees of the network. Eq.1 has the simple solution $`k_i(t)=m(t/\tau _i)^{1/2}`$, where $`\tau _i`$ is the time at which node $`i`$ entered the network. Since the relation between $`k_i`$ and $`\tau _i`$ is monotonous we can classify nodes according either to their degree $`k`$ or to their age $`\tau `$. As a consequence we can apply the usual formula to transform probability distributions: $`P(k)dk=\rho (\tau )d\tau `$. Since nodes enter the network at a constant rate, we have $`\rho (\tau )=const`$ and therefore $`P(k)k^3`$. As mentionent above for the Internet, the exponent $`\gamma `$ is in general not equal to the BA prediction $`\gamma =3`$; yet it has been shown that, as long the attachment rate in (1) is asymptotically linear in $`k`$, the distribution $`P(k)k^\gamma `$ with $`2<\gamma <\mathrm{}`$ depending on the pre-asymptotic behavior KRL . Another important feature of the BA model, that we are going to exploit for our analytical approach, is the lack of correlations between the degree of a node and the degrees of its neighbors. This is best represented through the average neighbor degree $`k_{nn}(k)`$, that is the average degree of the neighbors of a node of degree $`k`$: this quantitiy is essentially constant for the BA model.

The first step in our analysis is to build a BA network, whose degree distribution is shown in Fig.1. Then, starting from a node (we choose a highly connected node) we begin our exploration procedure: 1) each edge connecting that node to its neighbors is followed with probability $`p`$; edges that are lost at this stage are lost forever 2) from each of the reached nodes, repeat step 1) until no new nodes are reachable. In this procedure, edges to nodes that have already been reached are not followed. The result of this algorithm is a network that has fewer nodes than the original one, and, on the average, fewer links per node, and that is, topologically, a tree. The intuitive result would be that every node sees just a fraction $`p`$ of its edges, so that all degrees should be reduced of a factor $`p`$, without consequences on the power-law behavior of $`P(k)`$. Actually the effects of the probability $`p`$ are much more dramatic. As it can be seen from our simulations (Fig.1), the measured exponent $`\gamma _m`$ actually changes. For a network grown with $`m=1`$ and explored with $`p=0.5`$ the measured exponent is close to $`\gamma _m=2.5`$. We can therefore wonder whether this is a crossover effect, and the correct exponent is recovered for very large networks, or whether this change of exponent is real. A simple analytical argument in favour of this second interpretation can be formulated using the lack of correlations in the BA model. Indeed, since there is no correlation between the degree of a site and the degrees of its neighbors, due to (1) there is no correlation between the age of a node and the age of its neighbors. This allows us to look at exploration during the growth of the network. In particular we can say that, in a growing network formalism, any time a new node is added to the network, we label it as reachable if it connects to at least a reachable site through a followed connection (with probability $`p`$). We assume that the first site is reachable. Then, the density of reachable nodes at time $`t`$ is given by

$$\frac{dN(t)}{dt}=1\left(1p_0^t\frac{dN(t^{})}{dt^{}}q(t^{})𝑑t^{}\right)^m$$
(2)
where $`q(t^{})`$ is the probability to choose a node introduced in the network between $`t^{}`$ and $`t^{}+dt^{}`$: the preferential attachment rule translates to $`q(t^{})=1/[2(tt^{})^{1/2}]`$ (this trick is similar to assigning to each node a hidden variable corresponding to the time $`t^{}`$ at which it entered the network, with a connection probability that depends on the hidden variables of both the new and old nodes; for more details see below BPS03 ; CBS03 ). Since $`N(t)`$ can grow at most linearly, we make the assumption that $`dN(t)/dtt^\alpha `$ with $`\alpha `$ expected to be negative. After some algebra, and keeping only the leading terms, we find $`\alpha =(mp1)/2`$: as long as $`mp<1`$ the density of reachable nodes decreases in time. Then, the measured degree distribution can be again obtained from the relation $`P_m(k)dk=\rho (\tau )d\tau `$, with $`\rho (\tau )\tau ^\alpha `$, from which we obtain $`P_m(k)k^{\gamma _m}`$, with $`\gamma _m=2+mp`$ For $`m=1`$, $`p=0.5`$ we have $`\gamma _m=2.5`$, in agreement with simulations. We expect therefore that, as long as $`mp<1`$, the measured exponent could be different from the real one.

To check whether the distortion of the exponent is a feature only of BA networks, we have also studied networks generated according to the hidden variables model. Hidden variable networks are characterized by a quantity $`x`$ (the ”fitness”) assigned to every node and taken from some probability distribution $`p(x)`$; every pair of nodes $`i`$ and $`j`$ is connected then with a probability $`q(x_i,x_j)`$. As a consequence the average degree of a node of fitness $`x`$ is

$$k(x)=N_Uq(x,x^{})p(x^{})𝑑x^{}$$
(3)
where $`N`$ is the number of nodes in the network and $`U`$ is the support of $`p(x)`$. Eq.3 gives a relation between $`k`$ and $`x`$ that can be used in $`P(k)dk=p(x)dx`$ to obtain $`P(k)`$. Suitable choices of $`p(x)`$ and of $`q(x,x^{})`$ give SF networks. Here we use the same examples provided in CCDM02 , namely Zipf and exponential fitness distributions. Zipf distributed fitnesses are inspired by the idea that many quantities such as personal wealth, company size, city population and others are power-law distributed Zipf . In this case a connection probability $`q(x,x^{})xx^{}`$ ensures that the resulting network is SF. Using $`p(x)x^3`$ we obtain $`P(k)k^3`$ and in Fig.2 we show that also this network’s exponent has changed following the exploration, with a value $`\gamma _m2.7`$.

Next we look at SF networks obtaned using an exponential fitness distribution $`p(x)=exp(x)`$ and a connection probability $`q(x,x^{})=\theta (x+x^{}x_c)`$, that is, a link between two nodes of fitness $`x`$ and $`x^{}`$ is present only if the sum $`x+x^{}>x_c`$. Eq.3 yelds $`k(x)=Ne^{x_c}exp(x)`$ and as a consequence $`P(k)k^2`$ (see Fig.3CCDM02 . The tree-like exploration of this network shows that, also in this case, the measured exponent can change, $`\gamma _m1.5`$. We do not have at the moment an analytical derivation of the measured exponents. Indeed, degree-degree correlations between nearest neighbors and the lack of an explicit time evolution hinder the formulation of some equations similar to (2).

Interestingly, in all cases we have analysed, the measured exponent $`\gamma _m<\gamma `$, an indication that the exploration process penalizes nodes with small degree with respect to nodes with large degree. This is reasonable, since a node with few connections has fewer paths reaching it (and some bottlenecks, since all these paths have ultimately to flow through its few connections) than a high degree node. The final result is therefore that high degree nodes are fairly well represented in the final distribution, whereas the number of nodes with few connections is underestimated. This intuitive picture rationalizes our numerical and analytical finding that the measured exponent is smaller than the real one.

## 3 Edge-Picking Exploration of Erdös-Rényi and Complete Graphs

The last example is also a good example of how other link detection techniques can change the apparent topological properties of networks in an even more dramatic way, in particular if the probability to detect a link depends on some intrinsic properties of the nodes it connects. Indeed, we can interpret the generation of SF networks from hidden variables as the result of the exploration of complete or Erdös-Rényi ER  (ER) networks. Starting from a complete or ER graph, we assign to each node $`i`$ a variable $`x_i`$ taken from a probability distribution $`p(x)=exp(x)`$. Then, we prune the graph by discarding (that is, they are not detected) all those edges that join nodes $`i`$ and $`j`$ such that $`x_i+x_j<x_c`$, with $`x_c`$ a properly chosen threshold. This is tantamount to say that, during the exploration of the graph, the edge between nodes $`i`$ and $`j`$ is detected with probability $`q(x_i,x_j)=\theta (x_i+x_jx_c)`$, which brings us back to the formalism used above, that shows that the resulting network is scale-free with degree distribution $`P(k)k^2`$ (data shown in Fig.3, black symbols, refer to the procedure over a complete graph). Fig.4 (main panel) shows the results for a starting graph that is an Erdös-Rényi network of $`12800`$ vertices with a connection probability $`p=0.025`$, corresponding to a graph with average degree $`<k>=320`$; the threshold $`x_c=13`$. The resulting degrees are power-law distributed with exponent $`\gamma 2`$ (black thick line). The peak for large values is what is left of the Poisson distribution of the underlying ER network: whenever a node has a variable $`x>x_c`$, all of its connections are detected, and its degree is not distorted. Since in SF networks a special role is played by hubs, that is, nodes with a very large degree, we checked whether the hubs of the explored network fall into the Poisson cutoff or in the power-law distributed part. The result, shown in the inset of Fig.4 clearly shows that the distribution of the maximum degree nicely obeys a Frechet distribution $`P_{max}(k)=(\gamma 1)k^\gamma exp(k^{\gamma +1})`$ MAA02 , with $`\gamma =2.1(1)`$. Such a Frechet distribution is indeed the expected distribution of the maxima of set of variables taken from a power-law distribution $`k^\gamma `$. This implies that, apart from a $`1\%`$ of the networks (we show the distribution of $`k_{max}`$ over $`10000`$ networks), the maximum degree is almost always drawn from the power-law part of the degree distribution, an indication that typical networks can be considered genuinely scale-free.

## 4 Conclusions

The increase in the amount of real network data is prompting the community to study networks in more detail and to elaborate models able to predict qualitatively, but also quantitatively, the measured properties. Yet, before taking these data by face value, a thorough investigation of the measurement techniques is necessary to ascertain if and what kind of data distortion they could introduce. This is customary in physics, where systematic errors have always to be taken into account and possibly to be corrected, and the same kind of attention should be paid also to data from different disciplines. In this work we have shown, through simple examples, that tree-like explorations, inspired by the traceroute command, can indeed skew the data so that the measured exponent of the degree distribution of a scale-free network can change with respect to the real one. In the simplest case (BA networks), simulations and simple analytical arguments agree with each other. We have also shown that a recently proposed model of SF networks based on hidden variables can be interpreted as an exploration technique leading to the appearence of power-law degree distribution where the underlying network is topologically much simpler.

This work has been supported by the FET Open Project IST-2001-33555 COSIN, and by the OFES-Bern (CH).