text:
# 1 Introduction

## 1 Introduction

The requirements for reducing the practical execution time in GRACE have fostered our interest in parallelization of the multi-dimensional integration of GRACE. In addition, reducing the program size becomes more important to avoid the cash-miss, which increases the practical execution time, as the number of final particles becomes large. These motivations have lead to implement the parallelization of GRACE.

In GRACE, for the multi-dimensional integration, BASES/SPRING is used. BASES is a software package of Monte Carlo integration with an importance and stratified sampling method. For the parallelization of Monte Carlo integration, it is natural and efficient to distribute sampling points to processors. This parallelization is called the Data Parallel. When the integrand can be decomposed, the Function Parallel is applicable, that is, each part of the integrand is calculated in each different processor. Since 1992, we have investigated independently these two approaches to the parallel computation of the multi-dimensional integration in GRACE.

Based on above experiences, we have developed a new method, a hybrid use of Data Parallelism and Function Parallelism. In this method, we use MPI-1111The first version of MPI standard. It is standardized in May 1994. as the Message Passing Library which is standardized and is widely used for developing parallel code in both the distributed computing environment and MPP (Massively Parallel Processors) platform. A computing model adopted is SPMD (Single Program Multiple Data) computing model.

In this paper, in section 2 the details of the implementation of hybrid use of Data Parallelism and Function Parallelism is described. We also present the mechanism of the dynamic load balancing in section 3. The behavior of performance of the parallel computation is shown in section 4. Section 5 is devoted to a conclusion.

## 2 Hybrid Use of Data and Function Parallelism

In the hybrid method, we firstly distribute a group of hypercubes which are subspaces of the integral volume. In each hypercube, a definite number of sampling points are taken. The distribution of hypercubes corresponds to the distribution of sampling points. Secondly, together with distributing sampling points, we distribute the calculation of the scattering amplitudes, of which the integrand consists. The scattering amplitude is given as a sum of each matrix element corresponding to each Feynman diagram. Therefore the distribution of the calculation of the scattering amplitudes corresponds to distribution Feynman diagrams.

For the data-transfer among processors, we use a collective communication instead of the point-to-point communication. The collective communication is defined in a communicator which is one of the important concepts in MPI. It defines the communication space used for the communications among processors.

Fig. 1 shows the schematic view how the Data Parallel and the Function Parallel work together from the point of view of communicators, where hypercubes, as an example, are divided into two groups (DP1 and DP2) and then calculation of Feynman diagrams into three groups (FP1, FP2 and FP3). The communicator corresponding to the Data Parallelism is indicated as dpcom. On the other hand, the communicator corresponding to the Function Parallelism is indicated as fpcom. In fpcom, the global sum of the results of the calculation of the scattering amplitudes from each processor is carried out. Successively, in dpcom the global sum of several results needed for estimating the integral is carried out. In the figure, note that dpcom (hatched in the figure) is constructed by one of PE’s (Processing Elements) in each group, DP1 and DP2, to avoid making a global sum repeatedly.

Also note that the total number of parallel processors is given as the product of the degree of the Data Parallel and that of the Function Parallel. In this example, the number of parallel processors is $`2\times 3=6`$.

Summarizing, the parallel computation of the multi-dimensional integration is proceeded in the following way:

1. distribute hypercubes into parallel processors,
2. by using sampling points in hypercubes distributed, calculate a part of integrand (scattering amplitudes) in each parallel processor,
3. make a global sum of the scattering amplitudes for each sampling point in the communicator fpcom,
4. square the results obtained in (iii) and sum them up in the group, DP1 and DP2,
5. make a global sum in the communicator dpcom.

The above procedures are iterated until the multi-dimensional integration converges.

## 3 Dynamic Load Balancing

Once a physics process to be calculated is fixed, GRACE generates the Feynman diagrams according to defined physics model and defined order of the perturbation. Generated diagrams are numbered by GRACE for convenience. Assuming all these diagrams have the same numbers of vertices and internal lines for simplicity, the execution time needed to calculate each diagram is expected to be nearly identical. On this assumption, we distributed Feynman diagrams into parallel processors in the order conventionally numbered by GRACE.

However, indeed, the execution time required for calculating each Feynman diagram varies diagram by diagram because diagrams may have different kinds of couplings and internal lines. This small fluctuation causes the load imbalance among parallel processors and leads to the decrease of the performance.

We newly add the mechanism of the dynamic load balancing. The execution time of each diagram is automatically measured in the calculation of the scattering amplitudes. As the procedures 1. - 5. described in section 2 are iterated, diagrams are sorted in order of the height of the load. With these rearrangements of diagrams, the load on each parallel processor becomes well balanced. In Fig. 2, the way how the fluctuations of the load on each processor becomes small is shown when the degree of the Function Parallel is 8 and the number of Feynman diagrams is 232.

## 4 Performance Measurement

The performance of the parallel computation is often represented by the scalability which is the ratio of the resultant execution time to the execution time when the degree is 1.

Table 1 shows the scalability when we use up to 16 processors. As described in section 2, the number of parallel processors is given as the product of the degree of the Data Parallel and that of the Function Parallel though both can be set independently. The degree of the Function Parallel and the degree of the Data Parallel varies from 1 to 16 so as to keep the product of them equal to 16.

The physics process used in the measurement is :

* a physics process $`e^+e^{}b\overline{b}u\overline{d}\overline{\nu _\mu }\mu `$, and
* the total number of Feynman diagrams involved in this physics process is 232 at a tree level with the unitary gauge.

The performance measurement has been performed on AP3000 in Fujitsu Parallel Computing Research Center in Kawasaki, Japan. The AP3000 system consists of UltraSPARC -II 300MHz processors connected via AP-Net222AP-Net is a two-dimensional torus network. providing 200MB/s bandwidth per port. Of data-transfers among processors, the results gathered in fpcom are the biggest and the size of them is about 8.4MB per each transfer for the above physics process. In this measurement, the number of data-transfer is 8 times.

Table 1 clearly shows that the parallel computation gives an excellent effect to reduce the practical execution time of the multi-dimensional integration in GRACE.

## 5 Conclusion

We implemented a hybrid method of the Data Parallel and the Function Parallel in the multi-dimensional integration of GRACE. In the parallel computation of Monte Carlo integration, not only the sampling points but also the calculations of scattering amplitudes are distributed into parallel processors . We used MPI-1 as Message Passing Library. In the hybrid method, we can set the degree of the Data Parallel and that of the Function Parallel independently. The total number of parallel processors is given as the product of these two degrees. The computing model we used is SPMD model. To reduce the load imbalance among processors due to the small fluctuations of execution time for calculating each Feynman diagram, we implemented the mechanism of the dynamic load balancing for distributing Feynman diagrams to processors. With this mechanism, the effect to the reduction of the elapsed time in GRACE has been improved.

The reduction rate of the execution time has been measured on Fujitsu AP3000 system by using up to 16 processors for the physics process $`e^+e^{}b\overline{b}u\overline{d}\overline{\nu _\mu }\mu `$. In this study, we found when the program size is not large, the Data Parallel gives better performance than the Function Parallel. When, however, the program size is large, it is impossible to run the program on a single processor. In this case, we have to take the Function Parallel. Our study tells even the Function Parallel gives satisfactory results. Further, when the network speed is enough high, it is expected that the performance by the Function Parallel becomes high as the Data Parallel.

## Acknowledgments

We wish to thank the members of MINAMI-TATEYA collaboration for continuous discussions and many kinds of support. We are also grateful to express our sincere gratitude to Prof. Y.Shimizu for the valuable suggestions and continuous encouragements. Authors appreciate Prof. Y.Watase for the encouragements.


format: text

text:
# The Kerr–Schild ansatz for the Nariai spacetime and the generating conjecture

## The Kerr–Schild ansatz for the Nariai spacetime and the generating conjecture

MILAN ŠTEFANÍK, JAN HORSKÝ

Department of Theoretical Physics and Astrophysics, Faculty of Science,
Masaryk University, Kotlářska 2, 611 37 Brno, Czech Republic.
stefanik@physics.muni.cz, horsky@physics.muni.cz

PACS: 04.20.Jb, 04.40.Nr

Abstract

We will show that the Nariai metric, i.e. the static spherically symmetric vacuum spacetime with a cosmological constant, admits a conformally Kerr–Schild spacetime representation. We find the vacuum solutions of the Einstein–Maxwell equations for the Nariai metric using the Horský–Mitskievich generating conjecture.

### 1 Introduction

Some static spherically symmetric spacetimes can be written in the form

$$ds^2=e^{2\alpha }dt^2e^{2\beta }d\rho ^2\rho ^2\left(d\vartheta ^2+\mathrm{sin}^2\vartheta d\phi ^2\right).$$

Here $`(t,\rho ,\vartheta ,\phi )`$ is a coordinate system and $`\alpha ,\beta `$ are the functions of $`\rho `$ only. There exists the Nariai metric  which also describes the static and spherically symmetric solution of the Einstein equations with non zero cosmological constant $`\mathrm{\Lambda }`$. This spacetime has the form

$$ds^2=\left(1\mathrm{\Lambda }\rho ^2\right)dt^2\frac{1}{\left(1\mathrm{\Lambda }\rho ^2\right)}d\rho ^2\frac{1}{\mathrm{\Lambda }}\left(d\vartheta ^2+\mathrm{sin}^2\vartheta d\phi ^2\right).$$
(1)
In the second paragraph we will show that this static and spherically symmetric spacetime with a cosmological constant can also be transformed to the conformally Kerr–Schild form .

There are many techniques how to obtain the solutions of the Einstein equations and the Einstein–Maxwell equations of the general relativity, i.e. , , . For the Nariai metric (1) we use the Horský–Mitskievich generating conjecture originally formulated in . It is possible to find the exact solutions of the Einstein–Maxwell equations if a seed metric i.e. a spacetime without an electromagnetic field, has one Killing vector at least. This conjecture is a generalization of the fact that the timelike Killing covectors of the Schwarzschild and the Kerr solutions coincide (up to a constant factor) with the electromagnetic four-potentials of the Reissner–Nordström and the Kerr–Newman fields respectively . In the case of the Nariai metric (1) we will consider the Killing vectors $`_t`$ and $`_\phi `$.

### 2 The Kerr–Schild form of the Nariai metric

The Nariai metric (1) can be rewritten in the form

$$ds^2=\frac{1}{\mathrm{\Lambda }\rho ^2}\left[e^{2\alpha }dt^2e^{2\beta }d\rho ^2\rho ^2\left(d\vartheta ^2+\mathrm{sin}^2\vartheta d\phi ^2\right)\right],$$
(2)
where

$$e^{2\alpha }=\mathrm{\Lambda }\rho ^2\left(1\mathrm{\Lambda }\rho ^2\right),e^{2\beta }=\frac{\mathrm{\Lambda }\rho ^2}{\left(1\mathrm{\Lambda }\rho ^2\right)}.$$
(3)
Let us prove that (1) or (2) can be written in the conformally Kerr–Schild metric representation

$$ds^2=e^{2U}\left[ds_M^22H\left(d\tau \pm dr\right)^2\right],$$
(4)
where $`U=U(\tau ,r)`$, $`H=H(\tau ,r)`$ and $`ds_M^2`$ is the Minkowski metric (not necessarily in the Cartesian coordinates).

The calculations follow directly from the pattern presented in . Let us consider the transformation

$$dt=d\tau +f(\rho )d\rho ,\rho =\rho (r),$$
(5)
which transforms the spacetime (1) or (2) to the form (4). In the case of the Nariai metric the function $`f(\rho )`$ has the form

$$f(\rho )=\frac{1}{\mathrm{\Lambda }\rho ^2}\frac{\left(12\mathrm{\Lambda }\rho ^2\right)}{\left(1\mathrm{\Lambda }\rho ^2\right)},$$

and the term from the paper

$$ds^2=\frac{\rho ^2}{r^2}\left[ds_M^22H\left(d\tau \pm dr\right)^2\right],$$

takes the following form

$$ds^2=\frac{1}{\mathrm{\Lambda }r^2}\left[ds_M^22H\left(d\tau \pm dr\right)^2\right],$$
(6)
for the case of the Nariai metric. Here

$$2H=1\frac{r^2}{\rho ^2}e^{2\alpha },$$
(7)
and $`\alpha =\alpha (\rho (r))`$. The function $`\rho (r)`$ is found by integration of the relation

$$\frac{d\rho }{dr}=\frac{\rho ^2}{r^2}e^{\left(\alpha +\beta \right)},$$
(8)
where the functions $`\alpha `$ and $`\beta `$ are given by (3). Substituting for $`\alpha `$ and $`\beta `$ to (8) and after the integration we get

$$\rho =\frac{1}{\mathrm{\Lambda }r}.$$
(9)
Therefore the relationship (7) is equal to

$$2H=2\left(1\frac{\mathrm{\Lambda }r^2}{2}\right).$$
(10)
In such a way we have obtained the final form of (4)

$$ds^2=\frac{1}{\mathrm{\Lambda }r^2}\left[ds_M^22\left(1\frac{\mathrm{\Lambda }r^2}{2}\right)\left(d\tau dr\right)^2\right].$$
(11)
We have shown that the static and spherically symmetric spacetime with a cosmological constant (1) can also be transformed to the conformally Kerr–Schild form . We have completed, therefore, the assertion that it is possible to write every spherically symmetric solution of the Einstein equations in the conformally Kerr–Schild form (4).

### 3 The application of the Horský–Mitskievich conjecture to the Nariai metric

In 1989 the Horský–Mitskievich conjecture was formulated which states that  a test electromagnetic field having a potential proportional to the Killing vector of a seed vacuum gravitational field (up to a constant factor), is simultaneously an exact electromagnetic potential of the self consistent Einstein–Maxwell problem where the metric is stationary or static and it reduces to the mentioned seed metric when the parameter characterizing the charge of the electromagnetic field source vanishes.

This generating conjecture we will apply to the Nariai spacetime and to obtain a solution of the Einstein–Maxwell field equations.

The Nariai metric $`g_{mn}^{(\mathrm{seed})}`$ (1) has two Killing vectors

$$T^i=_t,\mathrm{\Phi }^i=_\phi .$$
(12)
The corresponding covector forms are

$`T_m`$ $`=`$ $`g_{mn}^{(\mathrm{seed})}T^n=\left(1\mathrm{\Lambda }\rho ^2\right),`$
$`\mathrm{\Phi }_p`$ $`=`$ $`g_{pq}^{(\mathrm{seed})}\mathrm{\Phi }^q={\displaystyle \frac{\mathrm{sin}^2\theta }{\mathrm{\Lambda }}},`$
and, consequently, the 4–potentials

$`A_m`$ $`=`$ $`\beta T_m=\beta \left(1\mathrm{\Lambda }\rho ^2\right),\text{for Killing vector}_t,`$ (13)
$`A_p`$ $`=`$ $`\beta \mathrm{\Phi }_p=\beta {\displaystyle \frac{\mathrm{sin}^2\theta }{\mathrm{\Lambda }}},\text{for Killing vector}_\phi .`$ (14)
Here $`\beta `$ is a parameter characterizing the charge of the electromagnetic field source.

The non-zero components of the electromagnetic field tensor are

$`F_{t\rho }`$ $`=`$ $`2\beta \mathrm{\Lambda }\rho ,\text{for}_t,`$
$`F_{\theta \phi }`$ $`=`$ $`{\displaystyle \frac{\beta }{\mathrm{\Lambda }}}\mathrm{sin}2\theta ,\text{for}_\phi .`$
Because the 4-potential (13) is a function of $`\rho `$-coordinate and the 4-potential (14) is a function of $`\theta `$-coordinate the charged metrics of the self consistent Einstein–Maxwell problem are supposed to have the form

$$ds_{new}^2=e^{2A(\rho )}dt^2e^{2B(\rho )}d\rho ^2e^{2C(\rho )}d\theta ^2e^{2D(\rho )}d\phi ^2,$$
(15)
for the 4-potential (13) and

$$ds_{new}^2=e^{2A(\theta )}dt^2e^{2B(\theta )}d\rho ^2e^{2C(\theta )}d\theta ^2e^{2D(\theta )}d\phi ^2,$$
(16)
for the 4-potential (14). From the Maxwell vacuum equations we obtain the following condition for the spacetime (15)

$$C^{}+D^{}A^{}B^{}+\frac{1}{\rho }=0,$$
(17)
and for the metric (16)

$$\dot{A}+\dot{B}\dot{C}\dot{D}+\frac{2\mathrm{cos}2\theta }{\mathrm{sin}2\theta }=0,$$
(18)
where a prime means the derivative with respect to $`\rho `$ and a dot means the derivative with respect to $`\theta `$.

Let us write the Einstein equations in the form

$$R_{ik}\frac{1}{2}Rg_{ik}\mathrm{\Lambda }g_{ik}=8\pi T_{ik}.$$

The Ricci tensor $`R_{ik}`$, the Ricci scalar $`R`$ and the energy–momentum tensor $`T_{ik}`$ are defined by . For the metric (15) these equations are

$`C^{\prime \prime }+C^2+D^{\prime \prime }+D^2B^{}\left(C^{}+D^{}\right)+D^{}C^{}+\mathrm{\Lambda }e^{2B}`$ $`=`$ $`4\left(\beta \mathrm{\Lambda }\rho \right)^2e^{2A},`$ (19)
$`C^{}A^{}+D^{}A^{}+D^{}C^{}+\mathrm{\Lambda }e^{2B}`$ $`=`$ $`4\left(\beta \mathrm{\Lambda }\rho \right)^2e^{2A},`$ (20)
$`A^{\prime \prime }+A^2+D^{\prime \prime }+D^2B^{}\left(A^{}+D^{}\right)+A^{}D^{}+\mathrm{\Lambda }e^{2B}`$ $`=`$ $`4\left(\beta \mathrm{\Lambda }\rho \right)^2e^{2A},`$ (21)
$`A^{\prime \prime }+A^2+C^{\prime \prime }+C^2B^{}\left(A^{}+C^{}\right)+A^{}C^{}+\mathrm{\Lambda }e^{2B}`$ $`=`$ $`4\left(\beta \mathrm{\Lambda }\rho \right)^2e^{2A}.`$ (22)
For the Einstein equations with respect to the spacetime (16) we obtain

$`\ddot{B}+\dot{B}^2+\ddot{D}+\dot{D}^2\dot{C}\left(\dot{B}+\dot{D}\right)+\dot{D}\dot{B}+\mathrm{\Lambda }e^{2C}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\beta }{\mathrm{\Lambda }}}\right)^2\left(1\mathrm{cos}4\theta \right)e^{2D},`$ (23)
$`\ddot{A}+\dot{A}^2+\ddot{D}+\dot{D}^2\dot{C}\left(\dot{A}+\dot{D}\right)+\dot{D}\dot{A}+\mathrm{\Lambda }e^{2C}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\beta }{\mathrm{\Lambda }}}\right)^2\left(1\mathrm{cos}4\theta \right)e^{2D},`$ (24)
$`\dot{B}\dot{A}+\dot{D}\dot{A}+\dot{D}\dot{B}+\mathrm{\Lambda }e^{2C}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\beta }{\mathrm{\Lambda }}}\right)^2\left(1\mathrm{cos}4\theta \right)e^{2D},`$ (25)
$`\ddot{B}+\dot{B}^2+\ddot{A}+\dot{A}^2\dot{C}\left(\dot{B}+\dot{A}\right)+\dot{A}\dot{B}+\mathrm{\Lambda }e^{2C}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\beta }{\mathrm{\Lambda }}}\right)^2\left(1\mathrm{cos}4\theta \right)e^{2D}.`$ (26)
The system of the equations (17) and (19 \- 22) should determine $`A\left(\rho \right),B\left(\rho \right),C\left(\rho \right)`$ and $`D\left(\rho \right)`$. We have thus to find the four functions from a system of the five equations. Fortunately, this system is not independent therefore we can easily obtain the solutions. The four unknown functions are

$`e^{2C\left(\rho \right)}`$ $`=`$ $`\left[c_2\left(\rho ^22c_1\right)\right]^2,`$ (27)
$`e^{2D\left(\rho \right)}`$ $`=`$ $`\left[c_2\left(\rho ^22c_1\right)\right]^2,`$ (28)
$`e^{2A\left(\rho \right)}`$ $`=`$ $`{\displaystyle \frac{\left(\beta \mathrm{\Lambda }\rho \right)^2\left(\rho ^624c_1^2\rho ^2+64c_1^3\right)}{\left(\rho ^22c_1\right)^2}}{\displaystyle \frac{c_3}{\left(\rho ^22c_1\right)^2}}{\displaystyle \frac{c_4\rho ^2\left(\rho ^46c_1\rho ^2+12c_1^2\right)}{\left(\rho ^22c_1\right)^2}},`$ (29)
$`\mathrm{\Lambda }e^{2B\left(\rho \right)}`$ $`=`$ $`\left[4\beta ^2\mathrm{\Lambda }^2\rho ^2e^{2A\left(\rho \right)}+2C^{}\left(\rho \right)A^{}\left(\rho \right)+C^{}\left(\rho \right)^2\right],`$ (30)
where $`c_1,c_2,c_3`$ and $`c_4`$ are integration constants.

For the system of the equations (18) and (23 \- 26) we get using the same way the solutions

$`e^{2A\left(\theta \right)}`$ $`=`$ $`\left[k_2\left(\mathrm{cos}2\theta +2k_1\right)\right]^2,`$ (31)
$`e^{2B\left(\theta \right)}`$ $`=`$ $`\left[k_2\left(\mathrm{cos}2\theta +2k_1\right)\right]^2,`$ (32)
$`e^{2D\left(\theta \right)}`$ $`=`$ $`k_4\mathrm{cos}^2\theta \left[1{\displaystyle \frac{1}{\left(\mathrm{cos}2\theta +2k_1\right)}}+{\displaystyle \frac{16k_1+8k_1^2+4k_1\mathrm{cos}^2\theta }{\left(\mathrm{cos}2\theta +2k_1\right)^2}}\right]`$ (33)
$`+k_3\mathrm{sin}^2\theta \left[1+{\displaystyle \frac{2\left(2k_1+1\right)}{\left(\mathrm{cos}2\theta +2k_1\right)}}+{\displaystyle \frac{2\mathrm{sin}^2\theta \left(1+2k_1\right)}{\left(\mathrm{cos}2\theta +2k_1\right)^2}}\right]`$
$`{\displaystyle \frac{2\beta ^2\mathrm{cos}^4\theta }{3\mathrm{\Lambda }^2}}\left[1+{\displaystyle \frac{2\left(2k_11\right)}{\left(\mathrm{cos}2\theta +2k_1\right)}}+{\displaystyle \frac{2\mathrm{cos}^4\theta }{\left(\mathrm{cos}2\theta +2k_1\right)^2}}\right],`$
$`\mathrm{\Lambda }e^{2C\left(\theta \right)}`$ $`=`$ $`{\displaystyle \frac{\beta ^2}{2\mathrm{\Lambda }^2}}\left(1\mathrm{cos}4\theta \right)e^{2D\left(\theta \right)}2\dot{A}\left(\theta \right)\dot{D}\left(\theta \right)\dot{A}\left(\theta \right)^2,`$ (34)
where $`k_1,k_2,k_3`$ and $`k_4`$ are integration constants. The invariants of the electromagnetic field of our solutions are

$`F_{ik}F^{ik}`$ $`=`$ $`8{\displaystyle \frac{\beta ^2\mathrm{\Lambda }^2\rho ^2}{e^{2A\left(\rho \right)+2B\left(\rho \right)}}}<0\text{for}_t,`$ (35)
$`F_{ik}F^{ik}`$ $`=`$ $`8{\displaystyle \frac{\beta ^2\mathrm{sin}^2\theta \mathrm{cos}^2\theta }{\mathrm{\Lambda }^2e^{2C\left(\theta \right)+2D\left(\theta \right)}}}>0\text{for}_\phi .`$ (36)
This means that the field is an electric type for the Killing vector $`_t`$ and a magnetic type for the Killing vector $`_\phi `$. Both solutions (2730) and (3134) admit three Killing vectors. In the first case they are $`_t,_\theta ,_\phi `$ and $`_t,_\rho ,_\phi `$ for the second one. This is differed from the seed spacetime (1) because we have assumed the higher symmetry of the charged solutions (15) and (16).

### 4 Conclusion

We have considered the static and spherically symmetric vacuum spacetime with the non zero cosmological constant which is called the Nariai metric (1). For this metric we have shown that this spacetime can be transformed to the conformally Kerr–Schild form (4). Therefore the statement that all static spherically symmetric spacetimes can be described by conformally Kerr–Schild metrics used in  really holds for all static spherically symmetric vacuum spacetimes with or without the cosmological constant.

For the Nariai metric (1) we have found the corresponding spacetime with an electromagnetic field generated by the Horský –Mitskievich conjecture. These solutions are really given explicitly by (2730) for the Killing vector $`T^i`$ and by (3134) for the Killing vector $`\mathrm{\Phi }^i`$. Each of these solutions already has three Killing vector fields, compared to two Killing vectors for the seed metric (1).

### Acknowledgements

We would like to thank Professor P. C. Aichelburg, Professor R. Beig and Professor H. Urbantke from the University in Vienna for their initial help. We are grateful to Professor N. V. Mitskievich for fruitful discussions.


text:
# 1 Introduction

## 1 Introduction

The immune system has become lately a classical example of emergence . It is an archetype for biological systems whose behavior cannot be directly deduced from their molecular and cellular properties .

As in many other biological systems, a large amount of microscopic data has been gathered. This data by itself is insufficient to explain the macroscopic behavior of the immune system.

Like any complex system, the immune one, requires classically a few complementary ways for describing it.

* The micro-biological level - molecules and cells..
* The generic properties that give the broad picture of what makes this system tick. For instance the ”clonal selection”: the random genetic changes of the antibodies and their selection by the affinity with appropriate antigens.
* The macroscopic mechanistic level. For example the aggregated action of groups of cell and cytokine: TH1 vs TH2, idiotypic vs anti idiotypic….
* The complex network of interactions. This represents the emergence of the immune response in terms of a highly connected (e.g neural) network. In particular, one uses neural networks to implement this view .
* The reaction-diffusion (partial) differential equations level that describes the IS as a nonlinear chemical system .

We are proposing here a methodology, which has the capacity to express much of the successes of the other methods while transcending their limitations and in particular the limitations which prevented the appropriate understanding of the central issue: the emergence of complex deterministic macroscopic functions out of simple stochastic microscopic interactions .

The methodology is based on recognizing the discreetness and spatial inhomogeneity of the biological units  and tailoring a specific model for each biological situation by representing the objects and interactions appropriate for each scale.

More precisely, we introduce microscopic simulation (MS)  methods, which on one hand represent directly the basic elements of the immune system but on the other hand allow the identification of the effective macroscopic objects relevant to the collective dynamics of the system.

In our MS we enact in the computer each reaction separately. More precisely we compute at each lattice point the probability of each reaction and then perform reactions randomly according to these probabilities.

We first show that very simple microscopic interactions may lead systematically (when studied with MS) to highly nontrivial effects which one would not have guessed using any of the classical methods.

We describe the resulting new phenomena and show how they are relevant for various biological applications (e.g. the generation of complex spatio-temporal localized self-organized adaptive structures).

In other words we show that systems that would be judged too simple to constitute an adaptive network and too linear to present non-homogeneous solutions in the differential equations context, do in fact emerge well defined/individualized macroscopic objects, which can move around, act in complex ways, and display adaptive behavior.

We present here our Microscopic Simulation and hybrid (discrete-continuum) models as a generic technology for representing in a direct intuitive way the microscopic knowledge on the immune system and for extracting non-trivial macroscopic information out of them (including the emergence of adaptive behavior and self-organization).

Our results are at variance with the predictions of the classical modeling methods , but are in encouraging agreement with many known IS features which were long waiting for an explanation.

In chap 2 we study the macroscopic effect of the proliferation of cells induced by a macroscopically homogeneous, but microscopically inhomogeneous signal distribution. We show that the proliferation rate predicted by the MS is many orders of magnitude larger than the one expected from ODE. In particular one obtains proliferation in conditions in which the naive continuum limit predicts decay.

In chap 3 we identify the source of the difference between MS and ODE and propose hybrid models (HM) that bridge the gap between the 2 limits.

In chap 4 we introduce additional reactions, which limit the growth of the proliferating entities. We show that in certain conditions the macroscopic inhomogeneity is even increased by the effect of competition.

Most of the effects analyzed in the four first chapters are due to the combination of auto-catalysis and of fluctuations.

In chap 5 We study the effects of discreetization in the presence of negative feedback loops.

We find that some effects predicted by the continuum simulations  (but not so frequent in nature) are not present when the system is analyzed using MS (e.g. the spurious oscillatory behavior of certain regulatory pairs).

In conclusion we are opening the way to a wide front re-evaluation of the dynamics that is supposed to emerge out of the known microscopic elementary interactions of the immune system. We provide the main effects to be re-scrutinized and the particular methods to do so. Paradoxically, the new methods that we propose are closer to the intuition of the practitioners in the immunological field. Indeed the objects which MS manipulates are the very same they encounter in their experimental work: cells, cytokines, epitopes (rather than abstract approximations such as differential equations, neural networks etc).

## 2 Dynamics of Discrete Proliferative Agents

Proliferating entities can be described as agents $`A`$ that duplicate whenever they meets a stimulating signal $`S`$. These agents have in general also a typical death rate $`d`$. Examples of such agents are : in immunology - immune cells reacting to an antigen  ; In ecology - animal proliferating whenever they find food…

We will deal with the dynamics of the stimulating signal in the next sections. In this section we will assume the stimulating agent $`S`$ is fixed in time.

This system can be described as a 2 reactions system (Figure 1). The first reaction is that whenever am $`A`$ agent meets a stimulating signal $`S`$ it multiplies with a probability $`d_1`$. We will note it as :

* $`A+SA+A+S`$. $`d_1`$

The second mechanism is the death of an agent with a probability $`d_2`$

* $`A\mathrm{\Phi }`$. $`d_2`$

In the homogeneous limit the growth rate of A is proportional to number of A and S pairs. Thus the number of new $`A`$ agents produced in a time interval $`\mathrm{\Delta }t`$ is:

$$\mathrm{\Delta }A=SAd_1\mathrm{\Delta }T$$
(1)
The number of A deaths is proportional to the total A population:

$$\mathrm{\Delta }A=d_2A\mathrm{\Delta }t$$
(2)
Thus the total change in A is:

$$\mathrm{\Delta }A=(Sd_1d_2)A\mathrm{\Delta }t$$
(3)
This equation becomes an ODE  when we take the limit $`\mathrm{\Delta }t0`$

$$\dot{A}=(Sd_1d_2)A$$
(4)
The solution of the ODE $`\dot{x}=\alpha x`$ is $`x=e^{\alpha t}`$ . In this case the solution is :

$$A=e^{(Sd_1d_2)t}$$
(5)
This solution has 3 regions of behavior (Figure 2):

a) If the average growth rate is higher than the death rate i.e. $`(Sd_1d_2)>0`$ then the $`A`$ concentration grows with an exponential growth rate of $`(Sd_1d_2)`$ This exponential growth will end when other mechanisms, discussed in the next section, become important.

b) If the average growth rate is lower than the decay rate i.e. $`(<S>d_1d_2)<0`$ the $`A`$ concentration decays to 0.

c) If the average growth rate is equal to the decay rate i.e. $`(<S>d_1d_2)=0`$ the $`A`$ concentration is constant in time.

This solution does not take into account the non-homogeneous distribution of $`S`$. In a real system the proliferation rate changes from one point to the other. Regions with a high $`S`$ concentration have a high proliferation rate, while regions with a low $`S`$ concentration have a low (or negative) proliferation rate.

In the ODE we assumed that the effective proliferation rate is equal to the average proliferation rate $`(<S>d_1d_2)`$. This assumption can be very mistaken in many cases. Indeed if the interaction between $`S`$ and $`A`$ are local then the local dynamics is described by the local population $`A_i`$ at each location i.

$$\dot{A_i}=(S_id_1d_2)A_i$$
(6)
The correct averaging of equation 6 is:

$$\dot{A}=d_1<S_iA_i>d_2A$$
(7)
Eq. 7 transforms into Eq. 4 if we can separate of $`<S_iA_i>`$ into $`<S_i><A_i>`$. Thus in order for Eq. 4 to hold S and A have to be independent variables . Eq. 4 is the very basis of any ODE approach. Yet it is obviously wrong in the case we present here since $`A_i`$ is high precisely in the sites where $`S_i`$ is high. Thus S and A are in fact strongly correlated.

Following this argument we expect to see differences between the microscopic representation of this system and its description by ODE.

Indeed if the signaling agent ($`S`$) is a random discrete agent its distribution is uniform over large scales but locally its density will contain some granularity. The difference between the large scale $`S`$ uniformity and its microscopic granularity is the origin of the macroscopic differences between the ODE and the MS results.

For example if $`(<S>d_1d_2)<0`$ (the average birth rate is lower than the average death rate) the ODE would yield a decaying solution but the MS yields an exponential growth (Figure 3). To understand this growth note that each point in space contains a different concentration of signaling agent ($`S`$). For most of the points $`i`$ one will have $`(S_id_1d_2)<0`$ and the local population will decay to zero. However for some points one may have $`(S_id_1d_2)>0`$ and the local population at these points will rise with an exponential rate. After a short time the contribution of the points with decaying population will be negligible and only the lattice points with an exponential growth contribute to the average population. In particular the population growth is dominated by the site with the highest $`S`$ number.

The previous description does not take into account the diffusion of the $`A`$’s. The $`A`$ diffusion will lead to the creation of large scale islands of high A population around each maxima of The $`S`$ density. These islands will keep growing until they fill all space. The diffusion of the $`A`$ does not change qualitatively the average population growth mentioned above for a very wide range of parameters . In order to avoid complex formulae we did not represent explicitly the diffusion term $`\mu \mathrm{\Delta }A`$. In all the simulations and analytic computations this term has been dully taken into account. We mention it in the sequel only when it has a major qualitative effect on the results.

## 3 How Well do Different Methods Deal with Discreetness

The MS helped us to uncover 2 basic features:

* The population increase in the MS is much higher than the one computed in the homogeneous case: The ODE is exponentially decaying while the MS is exponentially growing.
* Diffusion by itself is not enough to smooth $`A`$ concentration distribution.

The effect of the non-homogeneity can be observed in every experiment with a petri dish in which the proliferation of cells creates first local patterns of very high concentration even if the original distribution is homogeneous. The irrelevance of diffusion can be observed in many ecological systems, in which although the diffusion rate of each animal is much higher than its reproduction rate localization can still be observed.

### 3.1 Microscopic Simulation vs. PDE

One might think that using partial (rather than ordinary) differential equations will appropriately take into account the non-homogeneity in the a distribution. this is actually not the case: the fact that $`A`$ can diffuse will only increase the uniformity in the $`A`$ distribution, while the microscopic inhomogeneity of $`S`$ is lost upon expressing the $`S`$ distribution as a continuous function (recall that $`S`$ was spread with uniform probability across the space). Consequently there is no reason to expect that the partial differential equations acting on smooth distributions will do better than ODE. This is confirmed by rigorous analytical treatment.

In the next chapter we show that allowing $`S`$ to diffuse does not solve this problem either: the diffusion will only interchange one microscopically inhomogeneous configuration with another, equally inhomogeneous one. Their difference as well as their very inhomogeneity is lost once one expresses the distribution of discrete $`S`$ in terms of a continuous density function.

### 3.2 Discreetization Implies Fluctuations

The crucial role of microscopic discreetization in the emergence of the macroscopic inhomogeneity might be unfamiliar and even strange for the reader. After all, the microscopic corrections to the differential equations should morally lead to microscopic corrections in their solutions. That this is however not the case is known in the theory of electric conductivity since the late fifties. Indeed, by solving the Schroedinger equations for periodic potentials one obtains electronic wave functions that spread over the entire space (which implies electric conductivity) . However, at least in 2 dimensions, the slightest random uniform noise leads to spatially localized wave functions (insulator) . Similarly here the presence of microscopic discreetization of $`S`$ induces intrinsically microscopic spatial inhomogeneity in the $`S`$ distribution. The partial differential equation (PDE) for $`A`$ share the property that (at least in 2 dimensions) the slightest spatial inhomogeneity (of $`S`$) in the equations leads to macroscopic localization. We have checked that indeed the effect of $`S`$ discreetness can be fully taken into account by working with un-quantized (real valued) $`S`$’s with same distribution function as the quantized $`S`$’s (Figure 4). continuous

### 3.3 Hybrid Models

The previous section showed that the difference between the MS and ODE results are only due to the microscopic non homogeneity of the $`S`$ agents. We conceived a class of models containing elements from both ODE and MS. We call such models with the collective name of hybrid models (HM). The particular combination of ODE and MS techniques has to be tailored specifically for each biological situation.

In the present case we compute the local $`A`$ concentration changes in a continuous way, but use the $`S`$ distribution produced by the MS. However the spatial distribution of the $`A`$S is non-uniform, and in fact we don not assume that the $`A`$’s spatial distribution is .

In other words we replace the stochastic dynamics of the MS by one differential equation per lattice site, and we preserve thereby the inhomogeneity in the agents concentrations.

Other system require other HM that fit their specific aspects. For example a system containing regions with high interaction rates, and high concentrations, and regions with low interaction rates can be modeled using a HM that computes the interaction rate stochastically for low concentrations, and deterministically for high concentrations.

HM model has the advantage that when the number of A and S pairs per site is large, the local reaction rate can be computed directly rather than having to enact explicitly each and every $`A+S>A+A+S`$ reaction. Thus they have a much lower CPU cost than the MS.

intrinsically This very simple system shows that the results one obtains depends crucially on the appropriate use of models, and the recognition of the limits of applicability of the various approximations. It draws in rough lines the limits between the range of applicability of each type of model.

* ODE can be used when the diffusion rate is much higher than the maximum activation rate.
* MS are precise over a very wide range but have a very high cost in CPU time.
* 4D PDEs have all the severe problems related to the continuity assumption, and are not much cheaper than MS. However there is a small range of cases where they can provide analytical solutions .
* HM can reproduce most of the results produced in MS without the high CPU cost. They reproduce the effects of microscopic inhomogeneity.

The HM direction can be used in an adaptive way. In particular the collective objects emerging form the dynamics can be treated as the effective elementary agents for higher (coarser scale) level of analysis.

## 4 Mechanisms Limiting the Population Growth

The proliferation of the $`A`$ agents in the models of section 2 must be limited by some mechanism if we do not want them to take over the universe. We have shown that a high death rate is not enough in order to limit the population proliferation. The required mechanisms can be either an external regulator , or a limit on the resources available to the proliferating cells .

### 4.1 Local Competition

Many of the mechanisms limiting the A proliferation can be described as competition. Competition for a resource means that whenever 2 cells need the resource either to survive or to proliferate only one of them will get it and the other will die (or fail to proliferate). In ecology competition can be over food or water, while in the immune system competition can be over access to an antigen (Figure 5).

Such a system is sometimes called a Lotka Voltera system  and can be described using a local differential equation:

$$\dot{A}=d_1S_iA_id_2A_i\lambda A_i^2.$$
(8)
The system can exist in two regimes:

* If the proliferation rate is higher than the death rate, this system will proliferate until it reach the saturation level uniformly throughout the entire space. In such a case there is no difference between the MS and ODE results.
* If the proliferation rate is lower than the death rate than the ODE would eventually converge to a uniform vanishing concentration. By contrast the MS will reach a steady state of islands whose population density is at saturation level. These islands are embedded in a background empty of $`A`$. For high $`A`$ diffusion rate the MS yields non-zero concentration everywhere (Figure 6) in spite of the ODE prediction.

Paradoxically the conditions in which ODE and the MS show similar features are seldom realized in nature. They represent extreme cases in which a certain agent is capable to fill all the available space. In the animal world only a few species enjoy such a status (For example humans and their parasites). In immunology this regime corresponds to a cell type filling all the blood and immune system (systemic diseases). Except for these extreme cases the situation in which $`(<S>d_1d_2)<0`$ seems to be more appropriate for realistic systems.

### 4.2 Global Competition

The local competition mechanism presented in the previous section assumes that the limiting factor that creates the competition is homogeneous on the same scale as the $`A`$ agents concentration. (The scale of homogeneity is the spatial scale in which a cell type has a homogeneous concentration). Some reactions leading to competition have a larger homogeneity scale than the factor they are limiting. For example if cells proliferate and require for their proliferation a substance found in the blood. In other words, while the proliferation is a local mechanism, the factor leading to competition is a global one . In ecological systems the global character of the limiting factor may be related with the fact that the predators move much faster than their pray and that their population is proportional to the total prey population. Therefore their population constitutes an indirect interaction between distant prey locations. We will show that in such cases the effects of inhomogeneity are even more important than in the previous sections. As an aside let us note that this kind of reactions although formally local (A lion can eat a pray only if they meet) are very effective in inducing top down effects; that is the regulation of the local pray population according is enforced by their total population.

The overall homogeneity and uniformity assumption used in ODE has the benefit of merging into one formalism a large number of biologically different homogeneous systems. This advantage turns into a disadvantage when multiple spatial scales are involved in the system, as will be described in this section. MS on the other hand can describe in a natural way only local interactions. It is non-trivial to describe with MS factors acting on widely different scales. The solution in this case would be to create a new type of HM, which will allow in particular to express the interaction between agents and their averages over various ranges.

A simple example for such system is obtained by modifying the regulation mechanism of the system presented in the previous section (Figure 7). There the 2 mechanisms proliferation and competition were local interactions. We can now replace the local competition by a global competition. Even-though from the microscopic point of view the non-locality might look unnatural this is the only way that nonlinearities can arise at the macroscopic level. Indeed A term of the form $`<A>^2`$ can appear only as the average over terms of the type $`A_i<A>`$. Therefore the modification of Eq. 6 would be

$$\dot{A_i}=d_1S_iA_id_2A_i\lambda A_i<A_i>$$
(9)
That is we replaced the local competition in the system described in section 4.1 by a global competition term. Once again we find that passing from local to global reactions has crucial implications: By comparing the behavior of the system Eq. 9 with the system Eq. 8 one find that they are dramatically different in all the parameter regions. Even in the region in which the system Eq. 8 is described correctly by an ODE, the ODE are still incapable to describe properly the system Eq. 9. Indeed in the 9 system the entire space is empty of $`A`$s except for a single island centered around the point in which the $`S`$ concentration is maximal (Figure 8).

We can understand this system by simplifying it to :<sup>1</sup><sup>1</sup>1We rescaled the equation to $`d_1=1,\lambda =1`$ and added the death term $`d_2`$ into the $`S`$ density. This can be done with no loss of generality

$$\dot{A_i}=A_iS_iA_i<A_I>=A_i(S_i<A_i>)$$
(10)
One can see that in every point in which $`S_i`$ is lower than $`S_{max}`$ the $`A_i`$ density will vanish, while around the point $`i_{max}`$ in which $`S_i=S_{max}`$ a region of high $`A_i`$ density will appear. Suppose that one starts with all the $`A_i`$s very small; As long as the average concentration of the $`A`$’s ($`<A_i>`$) is lower than $`S_{max}`$, the local concentration $`A_i`$ at the point $`i_{max}`$ will rise (acording to Eq. 10). This in turn will raise the value of $`<A_i>`$. Eventually $`<A_i>`$ will reach the value of $`S_{max}`$. What will be the equilibrium value of the $`A_i`$s in every other point? In every point where $`S_i<S_{max}`$, one will have $`A_i(S_i<A_i>)=A_i(S_iS_{max})<0`$, which according to 10 means that the concentration at those points will decay to 0. It is clear therefore that the entire $`A`$ population responsible for $`<A_i>=S_{max}`$ is due to a very dense island situated around $`i_{max}`$. The equality $`<A_i>=S_{max}`$ is not affected by the diffusion, which only has the effect of spreading the distribution of the $`A_i`$’s in a limited region around $`i_{max}`$. One sees that the dynamics is dominated in this case by the most extreme stochastic local fluctuation in the $`S_i`$ distribution, and not by the global or average properties of the $`S_i`$s.

### 4.3 Dynamical $`S`$s

In the previous section we considered a static $`S_i`$ distribution. The intrinsic microscopic inhomogeneity of the $`S`$s produced large-scale inhomogeneity in the concentration of the $`A`$s. We will show that these effects survive the introduction of $`S`$ diffusion, as long as the maximal proliferation rate of the $`A`$’s is not lower than the diffusion rate of the $`S`$s. If the signaling agent diffusion rate of the $`S`$ is very high, the $`A`$s will “see” a blurred $`S`$ concentration and will react only to the average $`S`$ concentration, as in the ODE limit. If however the diffusion rate of the $`S`$ agents is not much higher than the $`A`$s proliferation rate then the $`A`$s do react to the local concentration of the $`S`$s. The $`A`$s produce peaks in their concentration at the current locations of the $`S`$ maxima. The locations of the $`S`$ maxima change since the $`S`$s diffuse. Consequently the regions of high $`A`$ concentration will follow them. This will lead to long term large (intermittent) fluctuations in the total $`A`$ concentration (Figure 9.).

In this system the difference between the ODE and MS is not due to any imposed inhomogeneity in the starting conditions. It is only due to the discreteness of the $`S`$s.

### 4.4 Emergence of Complexity

In the first chapters it has been shown that the discrete character of the elementary agents $`S`$ leads to unexpectedly rich and complex behavior of the $`A`$s population. The complex spatio-temporal emerging structure of the $`A`$ population cannot be predicted by an ODE approach. Moreover in this approach one cannot even ask the appropriate questions about the system’s dynamics. However a system which has only proliferation and death (linear) terms leads eventually to a divergent $`A`$ population. In the present chapter we introduced a few types of competition terms. Local competition limits the size and maximum concentration of each macroscopic $`A`$ island. However rather than interfering with the emergence of macroscopic features, the competition terms turned out to enhance the local and complex character of the $`A`$ islands. Global competition can endow these macroscopic islands with an emerging social (or sometime anti-social) behavior. The evolving macroscopic objects may now have apparent “goals”: seeking food, multiplying, dying or destroying the competing islands. The lifespan of these emerging islands is many order of magnitudes longer than the lifespan of their components (which of course remain non adaptive and totally mechanical during the entire history of the system).

## 5 Negative Feedback Loops

### 5.1 The Survival of the Weakest

The complex behavior observed in the preceding sections is the result of positive feedback loops (autocatalyis). One could think that the difference between the ODE and MS descriptions is specific to systems with autocatalytic components. We will show in this section that this is not the case. In fact we consider a simple negative feedback situation: The response of the immune system to a pathogen invasion. The high pathogen concentration leads to an increase in the immune cells concentration which in turn leads to a decrease in the pathogen concentration.

One can express the dynamics of this system using ODEs, in which $`x`$ is the concentration of the pathogen and $`y`$ is the concentration of the immune cells.

$`\dot{x}=axy`$ (11)
$`\dot{y}=bxycy=(bxc)y`$ (12)
Eq. 11 expresses the fact that upon meeting an immune cell the pathogen is destroyed. Eq. 12 expresses the fact that the immune cells proliferate when they meet a pathogen with a probability ($`b`$), and die with a constant rate ($`c`$). The solution of the ODE 11,12 consist in a fast rise of the immune cell population, and following it a decrease in the pathogen population, until the pathogen population is lower than $`\frac{c}{Oscillationsb}`$. After that point the immune cells population starts to decrease. The dynamics of the MS and the ODE are the same in this case, with a very notable difference in the final steady state of the system:

* Using ODE one gets 2 types of asymptotic steady states depending on the initial condition ($`x(0),y(0)`$). More precisely if $`bc(ln(bc)1)ay(0)+bx(0)cln(bx(0))<0`$ both the immune response and the pathogen decreases asymptotically to zero. Otherwise the immune response asymptotically decreases to zero, while the pathogen population approaches a non-zero asymptotic value. This can be interpreted as a constant latent infection.
* In the MS various spatial regions may converge either of the 2 steady states described above within a finite time. The regions containing a non zero concentration of pathogens (and no immune cells) constitute protected reservoirs of pathogens.

This can explain for example the fact that even people receiving potent anti HIV drugs for which the average concentration of HIV is undetectable for a long period reproduces the disease spontaneously when the drugs are stopped.

### 5.2 MS Eliminates Spurious Oscillations

It has been shown in the previous sections that ODE fail to represent emergent features of dynamical systems, but ODEs can also generate artifacts that do not appear in real systems. Such artifacts are exemplified by the following system consisting 3 cell types: An external antigen, idiotypic (Id) and suppressor cells (A-Id). The antigens stimulate the Idiotypic cells. The idiotypic cells in turn stimulate the anti idiotypic cells, which destroy idiotypic cells. We also introduce a saturation term that regulates the anti idiotypic cell population (Figure 12). The ODE describing this system is :

$`\dot{x}=axbxy`$ (13)
$`\dot{y}=cxydy`$ (14)
where x is the concentration of the Id cells and $`y`$ is the concentration of the a-Id cells. Such equations have an oscillatory solution around $`x=\frac{d}{c},y=\frac{a}{b}`$ with a period of $`\sqrt{da}`$ (Fgigure 13).

The same reactions simulated at the microscopic level do not produce oscillations . The MS produces regions containing a high concentration of Id cell with very low concentration of A-Id cells around them, or regions of A-Id cells that completely destroyed the Id cells around them dominate the concentration. The disappearance of the oscillations in the MS is due to the de-phasing of the different regions in space. Each regions variate independently so that the average yields a constant concentration. When one cluster has a maximum concentration of Id cells, other clusters can have a low concentration of Id cells. Therefore the oscillatory character of the ODE solution is related with the unrealistic assumption that the system is precisely homogeneous in space.

If the concentration of the Id and Anti Id are low, one can still have in the MS system quite significant fluctuations around the average due to the finite size of the system relative to the clusters size (Figure 14) (i.e the fact that the system contain at every moment a small number of clusters, so that the variation of every cluster influences the average behavior). This fluctuations involve the total disappearance of agents in certain regions of space. In such a situation it might take a macroscopic amount of time until this region is “colonized” again by agents immigrating from other clusters.

This simple system shows that the appearance of oscillations cannot occur spontaneously, unless there is an explicit mechanism phasing the whole system. This mechanism can be a very high diffusion and/or an external input forcing a unified timing over the whole system.Oscillations appear in many descriptions of systems using ODE . This simple system shows that these oscillations are many times only an artifact of the model continuity assumption, and not a realistic feature of the system.

## 6 Discussion

Dynamic Biological systems were traditionally described using ODE. ODE were invented and proved useful in simple physical systems that fullfiled certain restrictions:

* A very large density of interacting particles at every point in space.
* A high reaction rate.
* A diffusion rate higher than the local interaction rate.
* An immediate result of each interaction. i.e. no delay between the interaction and its result.

All of these assumptions are usually false in biological systems.

* The density of interacting particles is small. For example the number of interacting cells at every point is of order of tens.
* The reaction rate is low - The interaction between cells is a long process involving a long chain of molecular sub processes.
* The diffusion rate can in many cases be low. The motion of cells cannot be described as a diffusive process. The motion is best described as a slow directional motion directed by many chemical signals.
* Most biological processes require the production of new molecules or the division of cells. Each of these process can take hours - days. Thus there is a delay between an action and its result.

Moreover the dynamics of physical systems is dominated and in fact determined by conservation laws. In biology none of the classical quantities (mass,number,energy,momentum….) is conserved.

In this work we showed that even in simple cases ODEs fail to describe appropriately the dynamics. This is as well: Systems which were too simple to present complex emergent features in the framework of ODE turned out to present when studied appropriately a wide range of experimentally documented features . We used in addition to the straight-forward MS intermediate techniques (HM) that allow to take explicitly into account the intermediate relevant dynamical scales between the MS and ODEs.

### 6.1 Modeling Autocatalysis

One of the main sources of complexity is the combination of microscopic inhomogeneity and autocatalysis . In biological systems there is an inherent source of microscopic inhomogeneity- the discreetness of the entities composing biological systems. Moreover most biological systems involve proliferation which is the basic autocatalitic mechanism. In reality most situations are mixed containing both autocatalytic and regulatory mechanisms and can be modeled using custom made models describing globally the non catalytic parts and explicitly the microscopic details of the autocatlytic parts.

### 6.2 Inter-Scale Information Flow

Autocatalysis results in the flow of information from small scale to large-scale feature ; from the microscopic details of the system to its macroscopic behavior. The rapid evolution of small-scale features to a macroscopic scale affects the distribution of the macroscopic variables. In particular the macroscopic variables which are dominated by the sum of a small number of large values (instead of on the sum of a large number of small values) have a very specific (fractal-power law) space-time behavior. Such effects are enhanced when there is a direct effect of the macroscopic values on the microscopic mechanisms. Such a dialog between different scales is the obvious reality in biological systems .

Many biological systems evolve at time scales much longer than their fundamental microscopic time scales. In ODE the emergence of long time scales cannot occur naturally, while in MS they can emerge spontaneously. In MS, the agents self-organize  in collective structures with emergent adaptive properties and with dynamical space-time scales of their own. The effective dynamics of these emerging objects is much slower than the one of their components (and this fact constitutes their very definition ). Analyzing the system at the level of the emerging objects require automatically the transition to longer time scales and larger spatial scale. Each level corresponds to a new space-time scale in a hierarchy of complementary representations of the same complex macroscopic system.

We found in this paper that the core mechanism for the emergence of collective macroscopic objects is the interplay between autocatalysis and the discrete microscopic structure of the components of the system. This mechanism is totally missed by the mean field description, which predicts uniform extinction in conditions in which the real system emerges actually very lively collective adaptive objects.


text:
# 1 Hierarchical internal shock model

## 1 Hierarchical internal shock model

We propose that the internal shock scenario, which is the standard scenario proposed to explain the observed gamma–ray burst radiation (Rees & Mészáros 1994, Lazzati et al. 1999, Panaitescu et al. 1999), can work also for radio sources in general, and for blazar in particular. The central engine which gives origin to jets in radio sources may work intermittently, accelerating shells of plasma with slightly different mass, energy and velocity. Faster and later shells can then catch up slower earlier ones. In the resulting collisions a shock develops, converting some of the ordered bulk kinetic energy into magnetic field and random energy of the electrons which radiate.

The wind is treated as a sequence of $`N=t_w/t_v`$ shells, where $`t_w`$ is the duration time of the wind ejection from the central source and $`t_vt_w`$ is the average interval between consecutive ejections. Each shell is characterized by a mass $`M_j`$, a Lorenz factor $`\mathrm{\Gamma }_j`$ and an ejection time $`t_j`$. After setting the dynamics of the wind ejection, we calculate the radii where the shells collide, approximately given by:

$$R=\frac{2\alpha ^2}{\alpha ^21}\mathrm{\Gamma }_{\mathrm{min}}^2ct_v\alpha =\mathrm{\Gamma }_{\mathrm{max}}/\mathrm{\Gamma }_{\mathrm{min}}$$
(1)
where $`\mathrm{\Gamma }_{\mathrm{min}}\mathrm{\Gamma }_{\mathrm{max}}`$ is the range of bulk Lorentz factors of the shells. In the case of blazar jets $`\mathrm{\Gamma }10`$ and $`ct_v`$ is assumed to be of the order of the dimension of the central engine $`10^{14}`$ cm (for a BH with $`M10^9M_{}`$). Thus the collisions start at $`10^{16}10^{17}`$ cm, in agreement with the estimates of the source size given by the observed variability and by the requirement of negligible $`\gamma `$$`\gamma `$ collisions. The merged shells propagate on larger scales colliding again up to $`10^{20}`$ cm, where the differences among the shell velocities are completely smoothed out and the wind can be considered uniform.

For each collision we study the hydrodynamics, in order to determine the shock velocity, the compression ratio and the internal energy $`E_{sh}`$ of the shocked fluid. Assuming a given partition of $`E_{sh}`$ among protons, electrons ($`E_e=ϵ_eE_{sh}`$, with $`ϵ_e0.5`$) and magnetic field ($`E_B=ϵ_BE_{sh}`$, $`ϵ_B0.04`$), we calculate the relevant physical parameters in the shocked fluid and the emitted spectrum. The relativistic particles are assumed to have the same energy distribution (a broken power–law) throughout the entire emitting zone. This simplification does not allow to consider details of spectral changes on a timescale faster than the light crossing time of a single shell (few hours when $`R10^{17}`$ cm and a month on the parsec scale).

The relevant radiation processes are synchrotron, synchrotron self–Compton (SSC), and Compton scattering on the external radiation (EC), likely produced by the broad line region (BLR). The emission from each shocked region is assumed to last for the longest among the following timescales: the cooling time of the electrons, the time for the shock to cross the two colliding shells, and the time taken by photons to escape the source. We simulate the evolution of the total spectrum summing the locally produced spectra of those regions of the jet which are simultaneously active in the frame of the observer.

## 2 Results and Conclusions

In order to compare the results of the internal shock model with the observed spectrum and the temporal behavior of a powerful blazar such as 3C 279, we simulated a wind with an average kinetic luminosity $`L_w10^{48}`$ erg s-1. The injection lasts for 10 years with an average interval between two consecutive shells of 3 hours (30000 shells in total). The shell Lorentz factors are random between $`10<\mathrm{\Gamma }_j<30`$, the masses are also random with an average value of $`M_j4\times 10^3M_{}`$, and we require a constant luminosity for the wind. This implies that more energetic shells are followed by longer quiet times, during which the central engine “re-charges”. The BLR has a radius of $`R_{blr}=5\times 10^{17}`$cm and a luminosity of $`10^{45}`$ erg s-1.

Since the Lorentz factors of the colliding shells are only slightly different ($`\mathrm{\Gamma }_{\mathrm{max}}/\mathrm{\Gamma }_{\mathrm{min}}=3`$), the internal shock model is characterized by a low radiative efficiency: the fraction of kinetic energy converted into photons is less than 10%.

In Figure 1 we plot the average spectrum produced in different regions of the jet, in order to show the spectral evolution during the wind expansion. For the inner radii, $`RR_{blr}`$, the external Compton emission dominates: most of the radiation is emitted at $`10^{23}10^{24}`$Hz on time scales of the order of a few hours; the synchrotron spectrum peaks at $`10^{14}10^{15}`$ Hz, with a luminosity of $`10^{46}`$ erg s-1. At larger radii only the synchrotron and the SSC (first and second order) processes contribute. The shell–shell collisions become less efficient as the distance increases, and the synchrotron peak shifts to lower energy, reaching $``$10 GHz for collisions at $`10^{19}`$ cm. The time scale of the emission increases from few hours to few weeks, and the synchrotron dominates over the inverse Compton power at large distances.

The main observed properties explained by the internal shock model are:

1. Low efficiency — The radiative output of radio sources in general and blazars in particular must be a small fraction (say less than 10%) of the energy transported by the jet, since the extended radio structures require a power input much exceeding what is lost through radiation.
2. Minimum distance for dissipation — The bulk of the power we see must be produced at some distance ($`>10^{16}`$ cm) from the jet apex and from the accretion disk. Otherwise the high energy $`\gamma `$–rays are absorbed by a dense radiation field leading to electron–positron pairs and to a softer (especially X–ray) radiation, which is not observed.
3. Continuous energy deposition along the jet – Dissipation occurs all along the jet: the $`\gamma `$–ray flux is produced mainly in the inner part, and can vary rapidly, while the radio flux is produced by the parsec scale jet, with a variability time scale of the order of weeks–months.


text:
# Absence of Overscreened Kondo Effect in Ferromagnetic Host

## Abstract

We study the low temperature behavior of a boundary magnetic impurity $`S^{}=1/2`$ in an open ferromagnetic Takhatajian-Babujian spin-$`S`$ chain. For antiferromagnetic Kondo coupling, it is shown via Bethe ansatz solution that the impurity spin is always locked into the critical behavior of the bulk. At low temperature, a local composite of spin-($`S1/2`$) forms near the impurity site and its contribution to specific heat is of simple power law $`T^{\frac{1}{2}}`$. The absence of overscreened Kondo effect is due to the large correlation length of host spins which is divergent near the quantum critical point.

Multichannel Kondo problem was originally proposed to study the magnetic impurity behavior in real metals. At low temperature for antiferromagnetic Kondo coupling, the impurity spin $`S^{}`$ is completely compensated by the surrounding $`n`$-channel electrons when $`2S^{}=n`$; while it is undercompensated for $`2S^{}>n`$ or overcompensated for $`2S^{}<n`$. The later situation (so-called overscreened Kondo effect) bears no resemblance to a non-interacting gas of electrons and provides us a possible interpretation of non-Fermi liquid behavior observed in certain uraniun alloys as well as in certain tunneling problems in two-level systems. As this effect is immediately destroyed by weak anisotropy, it has been a challenge to deduce its signature relevant to other anomalous behavior in a variety of exotic metals, especially when the systems show reduced effective dimensionality, or undergo the proximity to quantum critical phase transitions. In this paper we will concentrate on multichannel Kondo problem in a typical quantum critical ferromagnet, with the number of channels $`n`$ being represented by the host spin, $`n=2S`$. Contrary to the problem in antiferromagnetic host chain, where overscreened Kondo effect appears when $`S>S^{}`$ , we will find that due to the divergent correlation length near the critical point, the overscreened Kondo effect does not exist in ferromagnetic host.

While the Kondo problems in paramagnetic or antiferromagnetic hosts have been intensively investigated over three decades, they attract little attentions in ferromagnetic host. The main reason is the difficulty in approaching this problem, i.e., usual perturbation techniques as well as bosonization fail in the critical region because of the absence of the “large energy scale”, i.e., the Fermi energy. Recent investigations have been shown that as the system closes to a quantum critical point, the impurity behavior depends strongly on the host properties and is non-universal. Another reason is due to the fact that the existence of Kondo effect in ferromagnetic host remains controversial. At a first glance, the traditional Kondo screening seems possible only at temperature higher than the Curie point where host magnetization vanishes. Otherwise, a local impurity could not change the total magnetization, opposite to what happened in dilute alloys. However, as far as the screening of a local moment is concerned, the issue should be how the impurity behavior changes as the temperature lowers down. The above picture does not necessary hold true even for the host goes away from the critical point where the long range order is destroyed by strong fluctuations. Yet there is no rigorous proof of the absence of overscreening in low temperature ferromagnet . In a earlier publication, Larkin and Mel’nikov first studied the Kondo effect in an almost ferromagnetic metal. With traditional perturbation theory they shown that the impurity susceptibility is almost Curie type with logarithmic corrections at intermediately low temperature and the usual Kondo screening is not effective. More recently, the authors studied a similar problem in 1D, i.e., a boundary Kondo impurity of spin $`S^{}`$ in $`S=1/2`$ ferromagnet. The model provides us a first exact soluble example for the Kondo problem near a quantum critical point. Via Bethe ansatz solution it has been shown that at sufficient low temperature, the impurity spin is completely screened by $`2S^{}`$ host spins. Thus a local composite forms near the impurity site and shows a simple power law dependence of specific heat on temperature. Therefore, these interest results strongly motivate us to know what will happen for a magnetic impurity in a multichannel ferromagnetic host at low temperature and if there exhibits the conventional overscreened Kondo effect.

In order to study the multichannel Kondo problem in ferromagnet via exact solution, we consider a magnetic impurity $`S^{}=1/2`$ coupled to an open spin-$`S`$ chain which is integrable by Bethe ansatz. A well known integrable extension of isotropic $`S=1/2`$ spin chain to arbitrary spin $`S`$ is

$`H_S=J_0{\displaystyle \underset{j=1}{\overset{N1}{}}}Q_{2S}(\stackrel{}{S}_j\stackrel{}{S}_{j+1}),`$ (1)
where $`Q_{2S}(x)`$ is a polynomial of degree $`2S`$ of $`SU(2)`$ invariant quantities $`x=\stackrel{}{S}_j\stackrel{}{S}_{j+1}`$. The general construction for $`Q_S(x)`$ is given in Ref., one recovers the standard Heisenberg model and Takhatajian-Babujian model for $`S=1/2`$ and $`S=1`$ respectively. Here we follow the definition of $`H_S`$ appeared in Ref. so that the single-magnon dispersion coincides with that of the standard Heisenberg ferromagnet. Thus for $`J_0<0`$, the model (1) is qualitatively equivalent to the standard spin-$`S`$ ferromagnetic Heisenberg chain. The boundary impurity spin is now coupled to the first host spin as

$`H_{imp}=J\stackrel{}{S}_1\stackrel{}{S^{}}.`$ (2)
The novelty of this construction is that the total Hamiltonian $`H_S+H_{imp}`$ is again integrable for any $`J_0`$ and $`J`$. The result is compared with the problem of a single impurity in periodic spin chains, where the integrability requires a fixed Kondo coupling $`J`$ and a fine turned impurity interaction instead of simple form (2). As long as spin dynamics is concerned, the integrable spin-$`S`$ chain represents the multichannel host, because Kondo screening of the impurity is equivalent to usual multichannel Kondo problem in antiferromagnetic case($`J_0>0`$): the impurity spin is completely screened for $`S^{}=S`$; partially screened for $`S^{}>S`$, with Schottky anomaly when an external magnetic field $`H`$ is applied; and overscreened for $`S^{}<S`$, giving rise to quantum critical behavior.

The model Hamiltonian is solved within the framework of Bethe ansatz approach to impurity problem. Eigenstates of the model are parameterized by a set of rapidities $`\{\lambda _\alpha \}_{\alpha =1}^M`$ ($`M`$ being the number of down spins), with total magnetization $`S_z=NS+1/2M`$. Each rapidity represents a magnon of energy $`J_0S/(\lambda _\alpha ^2+S^2)`$, they are solutions of Bethe ansatz equations

$`e_{imp}(\lambda _\alpha )e_{2S}^{2N}(\lambda _\alpha )={\displaystyle \underset{\pm }{}}{\displaystyle \underset{\beta \alpha }{\overset{M}{}}}e_2(\lambda _\alpha \pm \lambda _\beta )`$ (3)
where $`e_r(x)=\frac{x+ir/2}{xir/2}`$, $`e_{imp}(x)=e_{1+2c}(x)e_{12c}(x)`$ with $`c=\sqrt{(S+1/2)^2J_0/J}`$. In the following, we only consider the case $`J_0<0,J>0`$, so $`cS+1/2`$. In the thermodynamics limit the solutions are classified by $`n`$-strings or $`n`$-magnons, they are associated with functions $`\eta _n(\lambda )`$ ( with real rapidity $`\lambda `$ ), which in thermal equilibrium for finite temperature $`T`$ and field $`H`$ satisfy an infinite set of non-linearly coupled integral equations:

$`\mathrm{ln}\eta _n={\displaystyle \frac{\pi }{g}}{\displaystyle \frac{|J_0|}{T}}\delta _{n,2S}+𝐆\mathrm{ln}[(1+\eta _{n1})(1+\eta _{n+1})],`$ (4)
with the asymptotic condition $`lim_n\mathrm{}\frac{\mathrm{ln}\eta _n(\lambda )}{n}=H/T2x_0`$. Here, $`𝐆`$ is an integral operator with kernel $`g(\lambda )=1/2\mathrm{cosh}(\pi \lambda )`$. Notice that the free energy, expressed in terms of $`\eta _n(\lambda )`$, consists of three parts $`F=F_{bulk}+\frac{1}{N}(F_{imp}+F_{edge})`$, bulk part is always the dominate one, impurity and edge parts are both of order of $`1/N`$, edge part is only due to the open boundary condition and is irrelevant to spin dynamics. We will discuss bulk and impurity parts in the following.

Ferromagnetic host. The groundstate of ferromagnetic host has $`M=0`$ with finite magnetization $`S`$ per site. The pure states are those with spin aligned along an arbitrary axis. With finite $`T`$ and $`H`$ the elementary excitations are $`n`$-magnons, bound states carrying $`n`$ quanta of reduction of magnetization along the ordered axis. Each $`n`$-magnon has energy $`ϵ_n=2|J_0|a_{n,2S}(\lambda )`$ with crystal momentum $`K=4_{p=1}^{n,2S}\mathrm{arctan}\frac{n+2S+12p}{2|\lambda |}`$. In open boundary problem, $`K`$ is positive. The energy spectrum $`ϵ_n(K)`$ is a continuous function, increases monotonically as $`K`$ varies in $`[0,\pi min(n,2S)]`$. The long-wavelength limit, i.e., $`K0`$, $`ϵ_n(K)|J_0|SK^2/n`$, corresponds to $`|\lambda |\mathrm{}`$. The absence of gap for all $`n`$-magnons at boundaries of Brillouin zones $`\pi min(n,2S)`$ indicates the exact cancellation of Umklapp processes. Correspondingly in thermodynamics limit, the driving term in Eq.(4) dominates when $`T0`$, all $`\eta _n(\lambda )`$ functions are relevant. The situation is quite different to antiferromagnetic host ($`J_0>0`$), where only $`\eta _{2S}(\lambda )`$ is relevant as $`T0`$. So the solutions of Eq.(4) are more complicated for $`J_0<0`$ . Here two limiting cases are important: (a) Weak coupling limit: $`T\mathrm{}`$. The driving term disappears, there is only one parameter in equations, $`x_0=H/T`$, the solutions are

$`\eta _n={\displaystyle \frac{\mathrm{sinh}^2[(n+1)x_0]}{\mathrm{sinh}^2x_0}}1`$ (5)
for all $`\lambda `$ and $`n`$, this is the free spin limit. (b) Strong coupling limit: $`T0`$. In this case the driving term diverges, $`\mathrm{ln}(1+\eta _n)\mathrm{ln}\eta _n`$, we obtain

$`\eta _n(\lambda )=\mathrm{exp}[(nH+\pi a_{n,2S}(\lambda )|J_0|)/T].`$ (6)
At zero temperature, $`\eta _n=\mathrm{}`$ for all $`n`$, i.e., no magnons, the ferromagnetic groundstate with $`S_z/N=S`$ and $`E/N=H`$ is reproduced. At low but finite temperature, $`\eta _n(\lambda )`$ as function of $`\lambda `$ and $`n`$ show crossovers between the two asymptotic solutions (5) and (6). In principle, they can be solved numerically for fixed $`T`$ and $`H`$. To study the critical behavior of $`S=1/2`$ ferromagnetic Heisenberg chain, Schlottmann proposed an analytic method to solve these equations, based on an elegant but simple correlation-length approximation, and the result coincides with the numerical one very well. His idea is now developed for arbitrary $`S`$ as follows.

Let us first notice that for sufficient large $`|\lambda |`$, the driving term becomes negligible and one reproduces the weak coupling solution Eq.(5) for all $`n`$, up to $`𝒪(e^{\pi |\lambda |}|J_0|/T)`$. Similarly, for sufficient large string index $`n`$, $`\eta _n(\lambda )`$ is again close to Eq.(5), because the effect of the driving term is also negligible. On the other hand, for small $`|\lambda |`$ and $`n`$, the driving term dominates, $`\eta _n(\lambda )`$ approach the strong coupling solution (6). For intermediate $`\lambda `$ and $`n`$ we have a crossover region. By equating two solutions (5),(6) for small $`\lambda `$ and $`H`$, we define a crossover scale $`n_c(T)`$ as $`2S|J_0|n_c=T\mathrm{ln}[n_c(n_c+1)]/2`$, or solving it iteratively for $`T<<|J_0|`$,

$`n_c(T){\displaystyle \frac{2S|J_0|}{T}}[{\displaystyle \frac{1}{\mathrm{ln}(\frac{2S|J_0|}{T})}}+{\displaystyle \frac{\mathrm{ln}\mathrm{ln}(\frac{2S|J_0|}{T})}{\mathrm{ln}^2(\frac{2S|J_0|}{T})}}+\mathrm{}].`$ (7)
The meaning of $`n_c(T)`$ is clear: for $`n>n_c`$, $`\eta _n`$ is close to weak coupling solution while for $`n<n_c`$ it is close to strong coupling one. Therefore $`n_c(T)`$ serves as the correlation length of host, because it is the average number of the correlated spins in thermal equilibrium. Similarly, for each $`n`$-magnon, we define a crossover scale for rapidity,

$`\lambda _n^c(T)\sqrt{{\displaystyle \frac{n}{\mathrm{ln}(n+1)}}{\displaystyle \frac{2S|J_0|}{T}}}.`$ (8)
Obviously $`\lambda _n^c(T)`$ is the correlation length in momentum space: for $`|\lambda |>\lambda _n^c(T)`$ (long wave-length or low energy), $`\eta _n(\lambda )`$ is close to weak coupling solution, while for $`|\lambda |<\lambda _n^c(T)`$ (and $`n<n_c`$) it is close to strong coupling one. Thus we adopt the following strategy to solve $`\eta _n(\lambda )`$, or more conveniently, to calculate $`\zeta _n(\lambda )=\mathrm{ln}[1+\eta _n(\lambda )][\pi a_{n,2S}(\lambda )|J_0|+nH]/T`$: (i) As a zero order approximation, we assume $`\eta _n(\lambda )=\mathrm{}`$ for $`n<n_c`$ and $`|\lambda |<\lambda _n^c`$ , and elsewhere is given by Eq.(5). (ii) Substituting this approximation into the right hand side of Eq.(4), we get the first order approximation for $`\zeta _n(\lambda )`$. The result is

$`\zeta _n(\lambda ){\displaystyle \underset{m=1}{\overset{n_c}{}}}\{\mathrm{ln}[1+{\displaystyle \frac{1}{m(m+2)}}]{\displaystyle \frac{2}{3}}x_0^2\}B_{mn}(\lambda _m^c\lambda )`$ (9)
$`+2n\mathrm{ln}{\displaystyle \frac{\mathrm{sinh}(1+n_c)x_0}{\mathrm{sinh}n_cx_0}},`$ (10)
with $`B_{mn}(\lambda _m^c\lambda )=(_{\mathrm{}}^{\lambda _m^c}+_{\lambda _m^c}^{\mathrm{}})A_{mn}(\lambda \lambda ^{})d\lambda ^{}`$. The leading contribution to bulk free energy $`F_{bulk}=T𝑑\lambda g(\lambda )\zeta _{2S}(\lambda )`$ for small $`T`$ and $`H/T`$ is estimated as

$`F_{bulk}=1.1({\displaystyle \frac{2S}{|J_0|}})^{1/2}T^{3/2}`$ (11)
$`0.42|J_0|{\displaystyle \frac{(2SH)^2}{T^2}}[{\displaystyle \frac{1}{\mathrm{ln}(\frac{2S|J_0|}{T})}}+{\displaystyle \frac{\mathrm{ln}\mathrm{ln}(\frac{2S|J_0|}{T})}{\mathrm{ln}^2(\frac{2S|J_0|}{T})}}+\mathrm{}]`$ (12)
The zero field dependence of the free energy, as given above, is due to the $`n`$-magnon contributions from those within the cut-off, $`n_c`$; while the contributions beyond the cut-off is of higher order, $`T^2`$. So this part eventually approach the exact one as $`n_c\mathrm{}`$ or $`T0`$. The finite-field contribution of the free energy, which comes from all kinds of $`n`$magnons, is a result of competition of both parts, and it is approximately proportional to the correlation length $`n_c`$, giving rise logarithmic corrections to the bulk susceptibility. Notice that by rescaling $`2SJ_0J_0`$, the bulk free energy is proportional to that of the standard spin$`1/2`$ ferromagnetic Heisenberg chain, so one obtains $`C_{bulk}=2SC_{1/2}T^{1/2}`$ and $`\chi _{bulk}=2S\chi _{1/2}T^2\mathrm{ln}^1T`$.

Boundary impurity. When the impurity is coupled to the ferromagnetic host via antiferromagnetic $`J>0`$, Eq.(3) has a pure imaginary mode $`\lambda =i(c1/2)`$. This mode contributes a negative energy $`\epsilon _{imp}=|J_0|S/[S^2(c1/2)^2]<0`$ to the system, indicating the formation of a bound state, i.e., a local composite made of the impurity and its neighboring spin S. One notices that $`\lambda =i(c1/2)`$ is the only possible boundary mode for $`S^{}=1/2`$ and no boundary string is allowed because $`\lambda =i(c+1/2)`$ is forbidden as shown in Eq.(3).

Therefore, by taking into account the pure imaginary mode, the BAE for the bulk rapidities read

$`e_{2S}^{2N}(\lambda _\alpha ){\displaystyle \frac{e_{2c3}(\lambda _\alpha )}{e_{2c1}(\lambda _\alpha )}}={\displaystyle \underset{\pm }{}}{\displaystyle \underset{\beta \alpha }{\overset{(M1)}{}}}e_1(\lambda _\alpha \pm \lambda _\beta ).`$ (13)
When $`2c`$ is an integer, the impurity free energy is the difference of contribution from two effective ”ghost” spins $`S_+=c1/2`$ and $`S_{}=c3/2`$

$`F_{imp}={\displaystyle \frac{T}{2}}{\displaystyle 𝑑\lambda g(\lambda )[\zeta _{2c1}(\lambda )\zeta _{2c3}(\lambda )]}.`$ (14)
One finds that the leading contributions from the impurity is the same as that of $`\frac{1}{2}`$spin(rescaling of $`J_0`$ is used), but is always negative, i.e., $`C_{imp}=C_{1/2}T^{1/2}`$, $`\chi _{imp}=\chi _{1/2}(T^2\mathrm{ln}T)^1`$, indicating the freezing of some degrees of freedom of the bulk. Notice that $`c`$-dependence arises only in the next orders. When $`2c`$ is not an integer, the situation is somewhat more complicated. In this case, there are additional contributions coming from the non-zero residue $`\{2c\}=2c[2c]`$ ($`[2c]`$ is the integer part of $`2c`$), which renormalize the effective Kondo temperature $`T_k1/\mathrm{cos}(\pi \{2c\}/2)`$ as well as the physical quantities. Explicit calculation shows that they only change the subleading behavior.

The effect of the local composite can be analyzed in the same way. This local composite is made of the impurity spin and the neighboring host spin. In contrast to the single channel problem where the local composite is always a spin singlet, the present one is not a spin singlet, but with a residual spin $`S1/2`$. Its feature is encoded in BAE as follows

$`e_S^{2(N1)}(\lambda _\alpha )=\mathrm{exp}\{i\varphi (\lambda )\}{\displaystyle \underset{\pm }{}}{\displaystyle \underset{\beta \alpha }{\overset{(M1)}{}}}e_1(\lambda _\alpha \pm \lambda _\beta ).`$ (15)
where $`\varphi (\lambda )=i\mathrm{ln}[e_{2S}^2(\lambda )\frac{e_{2c3}(\lambda )}{e_{2c1}(\lambda )}]`$ represents the phase shift of a bulk spin wave scattering off the local composite. In the limit $`J\mathrm{}`$, one finds a momentum dependent phase shift $`\varphi (\lambda )0`$ indicating no complete compensation of the impurity spin in multichannel ferromagnetic host. (except for $`S=1/2`$ where $`\varphi (\lambda )=0`$). Similarly, For $`2c=`$integer, we estimate the local composite contribution to the free energy

$`F_{loc}={\displaystyle \frac{T}{2}}{\displaystyle 𝑑\lambda g(\lambda )[2\zeta _{2S}(\lambda )+\zeta _{2c3}(\lambda )\zeta _{2c1}(\lambda )]}.`$ (16)
One immediately obtains $`C_{loc}=(S1/2)C_{1/2}T^{1/2}`$, $`\chi _{loc}=(S1/2)\chi _{1/2}(T^2\mathrm{ln}T)^1`$. So the leading order behavior of the local composite is similar to that of ($`S1/2`$)-spin, indicating that the impurity is neither completely screened as in one channel problem nor overscreened as in antiferromagnetic host, but is always locked into the critical behavior of the bulk.

Discussions and conclusions. In this paper, we studied a spin-1/2 boundary impurity coupled to the ferromagnetic spin-$`S`$ chain. The results were obtained with the assumption that the crossover from strong coupling to weak one is abrupt. A smooth crossover region should not modify the low-$`T`$ and small $`H/T`$ dependences of the free energies obtained. Moreover, if the correlation lengthens $`n_c`$ and $`\lambda _n^c`$ are scaled, the low-$`T`$ dependence of our results remain unchanged, only the amplitudes are rescaled. The leading zero-field dependence of the bulk free energy is a rigorous result, even the amplitude agrees well with those obtained numerically for $`S=1/2`$ , or analytically for $`S>1/2`$ via spin-wave theory. However, since all $`n`$strings contribute, the situation is more complicated for the susceptibility, logarithmic corrections in susceptibility arise even without Kondo impurity, indicating the existence of a marginal variable which does not exist in classical spin wave solution.

The impurity effects are controlled by the bulk properties and the impurity coupling. In the case of antiferromagnetic host, see Ref., a boundary impurity will shows the usual Kondo effect when the impurity spin $`S^{}`$ is less than or equal to the host spin $`S`$, i.e., $`S^{}>S`$, otherwise, it shows a typical overscreened Kondo effect, i.e., $`S>S^{}`$. Depending on the Kondo coupling $`J`$, or $`c`$, there exhibits some novel quantum critical phenomenon in both cases. In the case of ferromagnetic host, the single channel problem ($`S=1/2,S^{}1/2`$) shows a novel Kondo diamagnetic effect, the impurity spin is completely screened by surrounding $`2S^{}`$ host spins for antiferromagnetic $`J>0`$. Other magnons see a local composite formed around the impurity site and interact with it weakly.

The situation for multichannel ferromagnetic is dramatically changed. Due to the presence of antifrromagnetic coupled impurity, there is always a pure imaginary mode in addition to the bulk $`n`$-magnons. The magnetization is suppressed due to the local composite which is a spin $`S1/2`$ object. But the impurity spin is not overscreened, in contrast to that in antiferromagnetic host. Rather, it is locked into the critical behavior of the bulk, with impurity specific heat showing a simple power law $`C_{imp}T^{\frac{1}{2}}`$. The local composite now effectively coupled ferromagnetically to the other host spins via quantum fluctuation, showing the same power law specific heat $`C_{loc}T^{\frac{1}{2}}`$. The absence of multichannel Kondo effect is due to the large correlation length $`n_c`$ which is divergent when $`T0`$. Though these results are limited to 1D , similar phenomena may exist in higher dimensions. In multichannel ferromagnet, the antiferromagnetic Kondo coupling always leads to a local composite of a smaller spin, and leads to ferrimagnetic state when there is a finite density distribution of the impurities.

The authors are grateful to Profs. N. Anderi, I.N. Karnaukhov and Yu Lu for helpful discussions. We acknowledge the financial support of NSF of China.


text:
# Microcanonical Lattice Gas Model for Nuclear Disassembly

## I INTRODUCTION

The Lattice Gas Model (LGM) is frequently used to compute observables in heavy ion collisions. The applications are numerous; see  for references. The calculations are normally done for a fixed temperature so that we will also call the the usual model canonical lattice gas model CLGM. The temperature is fixed so there will be fluctuations in energy. It is often argued that it might be more appropriate to keep energy fixed. In the following we set up a scheme for doing calculations with fixed energy. We will call this MLGM.

## II Calculational Procedure

We assume the reader is familiar with the usual LGM model for nuclear disassembly.

In CLGM, $`N`$ neutrons and $`Z`$ protons are put in $`N_s`$ lattice sites using a Metropolis algorithm. Because of bonds between nearest neighbours ($`ϵ_{np}=5.33`$MeV, $`ϵ_{nn}=ϵ_{pp}=0`$, ) and coulomb interaction between protons, there is a potential energy which we denote by $`E_{pot}`$. In Metropolis method, switch is attempted between occupied sites and unoccupied sites and also between occupied neutrons and protons. If, in the switch, the energy goes down the move is accepted. If the energy goes up, the move is accepted but only with a probability $`\mathrm{exp}(\mathrm{\Delta }E/T)`$. After many such switches an event is chosen. Once an event is picked, momenta are assigned from Monte-Carlo sampling of a Maxwell-Boltzmann distribution at temperature $`T`$. We denote the kinetic energy by $`E_{kin}`$. The total energy of the system is $`E_{tot}=E_{pot}+E_{kin}`$ which will fluctuate from event to event.

The change to MLGM can be made in the following way. We start from a given configuration, hence a given $`E_{pot}`$. The total energy which will be kept fixed is $`E_{tot}`$, so the kinetic energy is $`E_{kin}=E_{tot}E_{pot}`$. The available phase space of $`A=N+Z`$ nucleons having this kinetic energy is known analytically:

$`{\displaystyle \delta (E_{kin}\underset{1}{\overset{A}{}}p_i^2/2m)\mathrm{\Pi }d^3p_i}={\displaystyle \frac{2(\sqrt{\pi })^{3A}}{\mathrm{\Gamma }(3A/2)}}(2mE_{kin})^{3A/21/2}`$ (1)
We will call the value of the integral $`\mathrm{\Omega }_1(E_{kin})`$. We now attempt a switch in the configuration space. The potential energy will change to $`E_{pot}^{}`$. To conserve energy the kinetic energy of the system should be fixed at $`E_{tot}E_{pot}^{}=E_{kin}^{}`$. If $`\mathrm{\Omega }_1(E_{kin}^{})/\mathrm{\Omega }_1(E_{kin})>1`$, the move is accepted. If it is less, the probablity is given by the ratio. Since all configurations have identical weights, this satisfies the principle of detailed balance. After many such switches an event is accepted. We now have to assign momenta to the nucleons so that the total kinetic energy of the $`A`$ nucleons add up to correct kinetic energy which we denote by $`\stackrel{~}{E}_{kin}`$. The correct way to do this such that the sole criterion is phase space is the following. Choose a sphere of radius $`P`$. Do a Monte-Carlo sampling on $`A`$ particles for uniform distribution in this sphere. This means fixing $`p,\theta _p,\varphi _p`$ for each particle from $`p=P(x_1)^{1/3}`$, cos$`\theta _p=12x_2`$ and $`\varphi _p=2\pi x_3`$ where $`x_1,x_2`$ and $`x_3`$ are random numbers. Finally normalise $`P`$ so that the total energy equals $`\stackrel{~}{E}_{kin}`$. We are now ready to calculate all relevant quantities including cluster decomposition.

In microcanonical simulation observables can be calculated without having to invoke a temperature. But it is useful to extract a temperature in the model. The “temperature” will vary from event to event. The event temperature is taken from $`\frac{3}{2}(A1)T\frac{3}{2}AT=E_{kin}`$. The ensemble average gives the average temperature for the given microcanonical total energy. This is obviously attractive from an experimental point of view but we can also justify it from more basic principles. Let us write

$`\mathrm{\Omega }_T(E_{tot})={\displaystyle \mathrm{\Omega }_1(E_{kin})\mathrm{\Omega }_2(E_{tot}E_{kin})}`$ (2)
Although we are talking of one system only which has both kinetic and potential energy, formally, the right hand side is the same as two systems in “thermal” contact whose total energy is fixed but each one’s individual energy can vary. This is very standard statistical mechanics . For large systems, the sum is dominated by the maximum in the product $`\mathrm{\Omega }_1(E_{kin})\mathrm{\Omega }_2(E_{tot}E_{kin})`$ which is obtained when $`\frac{ln\mathrm{\Omega }_1(E_{kin})}{E_{kin}}=\frac{ln\mathrm{\Omega }_2(E_{pot})}{E_{pot}}`$ which then defines the inverse of the average temperature. This of course leads again to the same identification as above and is consistent with eq.(2.1). The sharpness of the maximum in the sum depends on the size of the system. For a small system we will then expect the fluctuation in temperature $`<T^2><T>^2`$ to be larger for the same value of $`<T>`$.

Finally although the starting point of Metropolis algorithm can be fairly arbitrary, it is helpful computationally to start from close to equilibrium. The method of generatiing the starting point followed the scheme given in  which was modified to take into account that there are two kinds of bonds and the coulomb force between protons.

Below we consider a few applications. Many more can be made.

## III The Caloric Curve

As the first application, we show in Fig.1, the caloric curve, $`E^{}/A`$ against temperature $`T`$, for 84Kr (an intermediate mass disintegrating system) and for 197Au (a heavy system). For the microcanonical model, $`<T>`$ is used for $`T`$ and for the canonical one, $`<E^{}/A>`$ is used for $`E^{}/A`$. Even though the fluctuation $`<T^2><T>^2`$ is more than four times in the case of Kr (6.59 MeV2 at $`<T>=5.0`$ MeV as compared to 1.1 MeV2 for Au at the same average temperature) there is almost no difference in the $`E^{}/A`$ vs. $`T`$ curves. We should point out, all calculations use two kinds of bonds and include the coulomb interaction. For 84Kr we use a cubic lattice of size $`6^3`$; for 197Au, we used $`8^3`$.

## IV IMF emission probability

One might expect that while average quantities will be nearly the same in both the models, there would be larger differences in fluctuations of observables. With this in view, we have investigated Intermediate Mass Fragment ($`IMF,Z`$ between 3 and 20) emissions. We show in Fig. 2, the plots of $`\sigma ^2<n_{IMF}^2><n_{IMF}>^2`$ (upper panel) and $`\sigma ^2/<n_{IMF}>`$ (lower panel) at different temperatures. Not a great deal of difference is found between CLGM and MLGM calculations for intermediate or heavy masses. Since the fluctuation in temperature for the intermediate case is much larger (more than a factor of four) this calls for an explanation. We think this is the reason. Referring to eq. (2.2) we see that when the temperaure is low, $`E_{pot}`$ is high. That means in the particular configuration the number of attractive bonds is less. However, the probability of an attractive bond being able to bind two nucleons is much higher since the temperature is low (the energy of relative motion $`p_r^2/2\mu `$ has a lesser chance of exceeding the bond energy $`ϵ_{np}`$, see ). Similar arguments hold when the temperature is high. This means $`E_{pot}`$ is low, so that there are many more attractive bonds. However, higher temperature will be able to break these bonds more easily. The two effects seem to cancel each other quite efficiently.

In Fig. 3 we compare $`P(n_{IMF})`$, the probability of emitting $`n`$ IMF’s, in the two models for the case of 197Au for four different temperatures. Again, the results are quite close.

We verified that in Fig.3, for temperatures 4.5 MeV and 5 MeV both the microcanonical and canonical probabilities are fitted quite well by a binomial distribution: $`P(n)=\frac{m!}{n!(mn)!}p^n(1p)^{mn}`$ where $`m`$ and $`p`$ are obtained from $`<n>=mp`$ and $`\sigma ^2/<n>=1p`$. This is a topic that has been discussed thoroughly in recent times . This parametrisation does not work for temperatures 3.5 MeV and 4 MeV where $`\frac{\sigma ^2}{<n>}>1`$. Here our calculated points are fitted quite well by a negative binomial distribution: $`P(n)=\frac{\mathrm{\Gamma }(N+n)}{\mathrm{\Gamma }(N)n!}p^n(1p)^N`$. The increase of $`\sigma ^2`$ at 4 MeV is not unnatural in LGM. This is happening because at this temperature the system is crossing the co-existence curve  with the accompaniment of a maximum in the fluctuation.

## V Summary

To summarise, microcanonical calculations in LGM are no harder to do than canonical calculations. The lowest mass number that we used was 84. Down to this size at least there are no serious departures from canonical results.

## VI Acknowledgments

This work was supported in part by the Natural Sciences and Engineering Council of Canada and by le Fonds pour la Formation de chercheurs et l’Aide à la Recherche du Québec.