format: text

text:
# Application of Coarse Integration to Bacterial Chemotaxis

## I Introduction

A recurring bottleneck in computational modeling of physical processes is the existence of multiple space and/or time scales. Important examples include patterns in fluids, defect dynamics in solids, and molecular dynamics of macromolecules. For example, in the latter case, the evolution of the physical configuration from initial to final state typically occurs on the time scale of milliseconds, while the interactions between constituent components must be resolved on the picosecond time scale. Furthermore, it is increasingly important in simulating biological systems to have a hybrid computational framework for interacting stochastic and deterministic processes, which occur on different time scales. Multiscale computational methods combine processing on fine space/time scales according to the governing microscopic description with macroscopic changes on a coarse grid (see for example, ref:Shenoy ; ref:Garcia ; ref:Kevghkrt ).

Here we present a computational scheme for “coarse projective” time integration of the macroscropic dynamics of stochastic processes that involve multiple, widely-separated time scales. To illustrate the method we apply it to a microscopic model that generates the biased random walk describing the macroscopic motion of a bacterium such as E. coli in an attractant gradient. Details of this Monte Carlo model of bacterial chemotaxis that integrates signal transduction with the motor response are given in Section III. The focus of our analysis in this paper is the computational scheme, and as such, our results address its plausibility rather than the phenomenology of chemotaxis.

In Section II we discuss previous work on coarse time integration of microscopic dynamics and the new features of the method developed here, and we motivate its application to chemotaxis. The Monte Carlo description is presented in Section III. We outline the coarse integration scheme in Section IV, and present numerical results and analysis of errors in Section V. We conclude with future extensions and applications of this work.

## II Previous work

The traditional approach to studying long-term dynamics of multiscale processes involves (a) the derivation of a “coarse-grained” set of evolution equations, followed by, (b) the analytical and/or computational study of these reduced equations using established, continuum numerical analysis tools. Recently, a so-called “equation-free” approach to the study of the coarse-grained behavior of such problems has been proposed which circumvents the first step (ref:TQK ; ref:GKT ; ref:LKGK ). This computational approach is based on the “coarse,” or macroscopic time-stepper, a map from the coarse variables at time $`t=0`$ to those at time $`t=T`$, where $`T`$ is typically much larger than characteristic microscopic time scales in the system. This map is not obtained directly from the macroscopic evolution equations, which we may not know, but rather through short time evolution intervals of appropriately initialized microscopic simulations. The initial macroscopic variables are lifted to microscopic variables to initialize the microscopic simulation. At the completion of a burst of microscopic simulation, the microscopic variables are restricted back to macroscopic variables, providing an approximation to the macroscopic time step. This provides a chord of the macroscopic solution, which is an approximation to the time derivative of the macroscopic solution. This value can then be used in any conventional continuum numerical method for the macroscopic equations. This approach has been applied in several microscopic contexts, and the results appear to be promising ref:MMK ; ref:MMPK ; ref:SGK ; ref:HK .

In this work we apply the coarse timestepper in a projective integration study of a spatially-distributed kinetic Monte Carlo simulation of a biased random walk. When, as in this case, the microscopic equations are stochastic, an effective algorithm must reduce the variance inherent in individual realizations of the stochastic process to a level that can be tolerated by the continuum numerical algorithm applied to the coarse system. This can be done by lifting to multiple copies, or, as is done here, by estimating the derivative from a least-squares fit to a large number of microscopic time steps. The projective integration method applied here uses a derivative estimate for the projective step, so that the stochastic noise is amplified by the inverse of the effective step length used in the derivative estimate. Hence, variance reduction is very important. Augmenting the number of copies of the simulation (which, for our noninteracting particle model corresponds also to a simulation with a larger number of cells) is the most direct approach to variance reduction; other variance reduction schemes are discussed in ref:Oettinger1 ; ref:Oettinger2 .

In the context of bacterial chemotaxis, the computational scheme developed here can be viewed as a direct method to compute the macroscopic evolution of the cell density in space and time without actually deriving these equations. An alternate approach begins with the transport equation for the velocity jump process, in which discontinuous changes in the direction (or speed) of an individual cell are generated by a Poisson process. It can be shown rigorously that this reduces to the chemotaxis equation under suitable scaling of space and time ref:Othmer:DLT:2002

$$\frac{\mu _0}{t}=(D\mu _0\mu _0\chi (S)S),$$
(1)
In the above, $`\mu _0(x,t)`$ represents the density of particles at spatial position $`x`$ at time $`t`$, $`D`$ is the diffusion constant, $`S`$ is the concentration of the chemotactic attractant/repellent, and $`\chi (S)`$ is the chemotactic sensitivity. However, there is as yet no analytical procedure available for determining the diffusion constant and chemotactic sensitivity when dependence on internal state variables determining the cell’s response to the external signal is explicitly included the transport equation. We show that a macroscopic time-stepper can be used here in the time evolution of the spatial density, even though the macroscopic evolution equations are not known.

## III Bacterial chemotaxis

### III.1 Signal transduction

The ability to sense and respond to environmental cues is necessary for the survival of most organisms. Escherichia Coli (E. coli) is a common and well-known single cell organism, with roughly $`4000`$ genes, whose chemotaxis network has emerged as a prototype for understanding signal transduction networks in general ref:BactRev . Its genome is known, the crystal structures of many proteins have been obtained, and a large number of mutant strains exist, allowing detailed behavioral studies.

For each cell, chemotactic behavior begins when attractant or repellent molecules bind to membrane receptors, triggering a cascade of chemical reactions inside the cell that culminates in the production of the phosphorylated form of a response regulator protein (CheY-P) which controls the direction of rotation of the flagellar motor. The series of reactions that converts the extracellular signal (attractrant/repellent) into cellular response is referred to as the signal transduction pathway. Flagella possess an inherent chirality, such that counter-clockwise (CCW) rotation results in bundling-up of the six to eight flagella per cell, allowing them to act as a single propeller and leading to smooth swimming motion of the bacterium. When the flagella rotate clockwise (CW), the bundle flies apart and the bacterium tumbles. As conditions become increasingly favorable due to increase in chemoattractant concentration, a bacterium extends its run-length; otherwise, it tumbles, and the subsequent direction of motion is randomly chosen, allowing a more favorable direction to be discovered. The resulting motion is a biased random walk toward favorable conditions and away from less favorable ones.

The sequence of biochemical reactions that take place inside a bacterium, starting with the binding of an attractant or repellent molecule to receptors on the cell surface and leading to the change in concentration of the response regulator species, CheY-P, has been extensively studied. Both deterministic models of these biochemical pathways, based on the law of mass action and Michaelis-Menten kinetics ref:Barkai\_Leibler ; ref:Spiro\_Othmer , as well as fully stochastic models exist ref:Bray\_98 . Measurement of concentration change is achieved through a temporal rather than a spatial comparison: fast sampling of the present external concentration is compared with memory of that concentration some time ago. Order-of-magnitude analyses for why measurement of concentration changes as spatial gradients across the cell length is not physically feasible have been given ref:Berg\_Purcell .

Memory in the signal transduction network is achieved through the existence of fast and slow reaction time scales. The fast reactions are receptor-ligand binding and phosphorylation kinetics; the slow reactions are methylation and demethylation. Ligand binding reduces the autophosphorylation rate of the corresponding membrane-bound receptor, in turn decreasing the rate of transfer of phosphoryl groups to CheY, and resulting in a (fast) drop in \[CheY-P\]. Addition of methyl groups to ligand-bound receptors restores the autophosphorylation rate, resulting in (slow) increase in \[CheY-P\]. The rate of demethylation becomes significant once a high methylation level is achieved. Hence, \[CheY-P\] reflects the balance between the fraction of ligand-bound receptors and methylated receptors. Perfect adaptation refers to the return of \[CheY-P\] to the same steady state level, regardless of the constant concentration level of the external stimulus. This value falls within the fixed operational range of the motor response to CheY-P.

A minimal model representing the signal transduction process includes fast excitation and slow adaptation to an external stimulus and is given in ref:Othmer\_Schaap :

$`{\displaystyle \frac{du_1}{dt}}`$ $`=`$ $`{\displaystyle \frac{\left(f(S)u_2\right)u_1}{\tau _e}},`$ (2)
$`{\displaystyle \frac{du_2}{dt}}`$ $`=`$ $`{\displaystyle \frac{f(S)u_2}{\tau _a}},`$ (3)
where $`\tau _e`$ and $`\tau _a`$ are the excitation and adaptation times, respectively, with $`\tau _e\tau _a`$. We identify $`u_1`$ with the deviation of \[CheY-P\] from its steady state value, and $`u_2`$ as the number of methylated receptors per unit volume. $`f(S)=f(S(t))`$ is a function of the external stimulus; for example, it is the number of bound receptors per unit volume in the presence of an external signal concentration, $`S`$:

$$f=N_TS/(K_L+S).$$
(4)
$`N_T`$ is the total number of receptors per unit volume, and $`K_L`$ is the dissociation constant for the ligand-binding process.

The formal solution to Eqs. (2)–(3) can be easily obtained for each cell as

$`u_2(t)`$ $`=`$ $`{\displaystyle \frac{e^{t/\tau _a}}{\tau _a}}{\displaystyle ^t}f\left(S(t^{})\right)e^{t^{}/\tau _a}𝑑t^{},`$ (5)
$`u_1(t)`$ $`=`$ $`{\displaystyle \frac{e^{t/\tau _e}}{\tau _e}}\left({\displaystyle ^t}f\left(S(t^{})\right)e^{t^{}/\tau _e}𝑑t^{}{\displaystyle \frac{1}{\tau _a}}{\displaystyle ^t}e^{t^{}\left(1/\tau _a1/\tau _e\right)}{\displaystyle ^t^{}}f\left(S(t^{\prime \prime })\right)e^{t^{\prime \prime }/\tau _a}𝑑t^{\prime \prime }𝑑t^{}\right).`$ (6)
However, because $`S(t)=S(x(t),t)`$, where $`x(t)`$ is a random variable that represents the cell’s position at time $`t`$ in a given external concentration field, the integration takes place along the cell trajectory, which is a biased random walk: $`x(t^{})`$ for $`t^{}<t`$ biases the probability of $`x(t)`$. Hence, $`S(t)`$ is a stochastic variable, and the integration must be carried out computationally.

### III.2 Motor and cell response

We adopt a stochastic approach to modeling the response of the flagellar motor to the regulator species, CheY-P. The motor is assumed to be in one of two states, CW or CCW, corresponding to clockwise or counterclockwise rotation ref:Block\_83 ; ref:Macnab . The rates of transition between these states are functions of \[CheY-P\]:

$$CCW_{i}^{}{}_{\stackrel{}{k_{}}}{}^{\stackrel{k_+}{}}CW_ii=1,\mathrm{},N_f,$$
(7)
where $`N_f`$ is the number of flagella per cell. For any concentration of CheY-P, the motor has a nonzero probability of being in either state. The probability of waiting at least a time $`T`$ for a transition from CCW to CW, $`\mathrm{\Phi }_+`$, or a transition from CW to CCW, $`\mathrm{\Phi }_{}`$, is given by:

$$\mathrm{\Phi }_\pm =e^{_0^Tk_\pm (t)𝑑t}.$$
(8)
If the motor is probed at a time $`T\mathrm{\Delta }t`$, such that $`k_\pm \mathrm{\Delta }t1`$, then the probability that a change in direction will occur in $`(T,T+\mathrm{\Delta }t)`$ can be approximated as

$$1\mathrm{\Phi }_\pm k_\pm (0)\mathrm{\Delta }t.$$
(9)
Systematic studies have been undertaken that show simulation results are qualitatively independent of the choice of $`\mathrm{\Delta }t`$ once the above restrictions are met.

In this work, we use equilibrium transition rates, $`k_\pm `$, consistent with recent experimental measurements ref:Cluzel  of the single flagellum clockwise bias, $`b_{CW}`$, and the reversal frequency, $`w`$, as functions of \[CheY-P\]. The rates are related to these measured quantities according to ref:Scharf\_Turner :

$`b_{CW}`$ $`=`$ $`{\displaystyle \frac{k_+}{k_++k_{}}},`$ (10)
$`w`$ $`=`$ $`{\displaystyle \frac{2k_+k_{}}{k_++k_{}}},`$ (11)
where the reversal frequency is the geometric mean of the transition rates. We take $`b_{CW}`$ to be described by a Hill function, $`b_{CW}=y_{p}^{}{}_{}{}^{h}/(K_{M}^{}{}_{}{}^{h}+y_{p}^{}{}_{}{}^{h})`$, where $`y_p=[\mathrm{CheY}\mathrm{P}]`$, with experimentally measured values for the Hill coefficient and dissociation constant, equal to $`h=10.3`$ and $`K_M=3.1\mu \mathrm{M}`$. In these experiments, it was noted that the data were consistent with $`w=(2\mu \mathrm{M}s^1)b_{CW}/y_p`$. Based on these results, we use

$`k_+`$ $`=`$ $`{\displaystyle \frac{hy_{p}^{}{}_{}{}^{h1}}{K_{M}^{}{}_{}{}^{h}+y_{p}^{}{}_{}{}^{h}}},`$ (12)
$`k_{}`$ $`=`$ $`{\displaystyle \frac{1}{y_p}}{\displaystyle \frac{hK_{M}^{}{}_{}{}^{h}}{K_{M}^{}{}_{}{}^{h}+y_{p}^{}{}_{}{}^{h}}}.`$ (13)
In a constant chemoattractant field, ref:Ishihara\_83 ; ref:Spiro\_Othmer  show that the “Voting Hypothesis” is successful in producing the correct running bias of the cell, $`B_{CCW}`$, from the individual flagellar bias, $`b_{CCW}`$:

$$B_{CCW}=\underset{j=\nu }{\overset{N_f}{}}\left(\stackrel{N_f}{j}\right)b_{CCW}^j(1b_{CCW})^{N_fj},$$
(14)
where $`\nu `$ is the minimum number of flagella required to be in the CCW state for the cell to run. Although cooperativity among the flagella is not yet fully understood, we similarly adopt the “majority rules” algorithm to determine whether each cell will run or tumble in the presence of a chemoattractant gradient:

$$\nu =\underset{i=1}{\overset{N_f}{}}s^{(i)}.$$
(15)
where $`s^{(i)}=0,1`$ is the state of the $`i^{th}`$ flagellum, corresponding to CW and CCW directions of rotation, respectively. Based on the experimentally measured gain of a single flagellum ref:Cluzel , a steady state value of $`\overline{y}_p=2.95\mu \mathrm{M}`$ yields $`b_{CCW}0.64`$ to within experimental error. The resulting running bias for the cell, according to Eq. 14, is $`B_{CCW}0.86`$ for $`N_f=6`$ and $`\nu =3`$, in agreement with experimental results. Hence, if three or more flagella are determined to be in the CCW state, the cell in our simulation runs; otherwise, it tumbles.

### III.3 Monte Carlo scheme

For each cell, the Monte Carlo scheme can be outlined as follows. At $`t_n`$:

* Each flagellum has a state $`s_n`$, where $`s_n=0`$ if CW, and $`s_n=1`$ if CCW. To determine $`s_{n+1}`$, for each flagellum, we draw a uniformly distributed random number, $`\zeta [0,1]`$.
* We determine an approximation to the probability for changing direction of rotation, given by Eq. 9. If the values of $`k_\pm `$ are such that this approximation is not valid, $`\mathrm{\Delta }t`$ is reduced until it is. This probability is a function of \[CheY-P\], which depends on the trajectory that the cell has taken along the external chemoattractant gradient.

+ If the flagellum is in the CW state ($`s_n=0`$) and $`\zeta <k_{}\mathrm{\Delta }t`$, it continues in the CW state ($`s_{n+1}=0`$); else, it switches to the CCW state ($`s_{n+1}=1`$).
    + If the flagellum is in the CCW state ($`s_n=1`$) and $`\zeta <k_+\mathrm{\Delta }t`$, it continues in the CCW state ($`s_{n+1}=1`$); else, it switches to the CW state ($`s_{n+1}=0`$).
* If three or more flagella are now in the CCW state, $`\nu _{n+1}3`$ in Eq. 15, and the cell runs; else, it tumbles. If the cell is determined to run and $`\nu _n<3`$, the direction of motion is chosen to be left or right with equal probability. Otherwise, it continues to run in the same direction.
* The position of the $`j^{th}`$ cell, $`x_n(j)`$, is accordingly updated to $`x_{n+1}(j)`$ (using the accepted time step $`\mathrm{\Delta }t`$), with the cell speed, $`v_{cell}`$, assumed to be constant and independent of position.
* The signal transduction variables $`(u_1,u_2)`$ are integrated in time, using the forward Euler scheme. Their time history is a function of each cell’s trajectory.

### III.4 Combining signal transduction and motor response

Although we use a toy excitation-adaptation model to describe signal transduction, it is coupled to an experimentally realistic description of the flagellar switch. To do so requires introducing (i) a shift in the steady state value of $`u_1`$, and (ii) an amplification or gain factor, $`g_0`$, for the signal transduction network, required so that the output of this network spans the dynamic range of the motor:

$$y_p=\overline{y}_pg_0u_1$$
(16)
The issue of gain in chemotaxis continues to receive much attention ref:amplification . Given the level of simplicity of our signal transduction model, we do not imply this treatment of network gain to be physically realistic. An obvious drawback is that, like other parameters in the toy model (excitation and adaptation times), it is static, whereas in the actual biological system, these parameters depend on the external signal.

In Table 1, we report the numerical values for the MC model parameters that operationally interface the toy signal transduction model with the realistic model of the flagellar switch. Hence, although the choice of toy model parameters was “physically motivated”, we do not attempt to make direct correspondence with experimental values. We caution the reader against making quantitative connection between our numerical results at these values and experimental results.

## IV Coarse integration

For the present problem of simulating the chemotactic response of a population of independent bacteria, the microscopic phase space is a directproduct of the phase space for each cell, consisting of its

* signal transduction variables: $`\stackrel{}{u}=(u_1,u_2)`$,
* flagellar state: $`\stackrel{}{s}`$, where $`s^{(i)}=0,1`$ for $`i=1\mathrm{}N_f`$,
* position: $`x`$,
* direction of motion, $`d=R,L,T`$, corresponding to running right, running left and tumbling, respectively.

Rather than evolving all microscopic degrees of freedom, a coarse integration scheme identifies suitable reduced variable(s) to be integrated. A clearly relevant reduced variable in population studies is the spatial density distribution of cell positions, $`\mu _0(x,t)`$, obtained from the set of individual cell positions, $`x_j`$, which are a subset of the full microscopic phase space. However, it is possible that the unknown equations of the coarse description use additional variables, for example, the densities of the right moving, left moving, and tumbling cells, $`\rho _i`$ where $`i=R,L,`$ or $`T`$. In this case, $`\mu _0(x,t)=_i\rho _i(x,t)`$. Indeed, an optimal reduced representation balances accuracy and efficiency of modeling the physical process. Here, we propose to retain only the spatial density distribution $`\mu _0(x,t)`$ as the relevant macroscopic variable in coarse modeling of chemotaxis since, as we will demonstrate numerically, the density of cells in each state rapidly approaches a functional of the total density.

A systematic approach to testing the adequacy of a particular coarse description is discussed in ref:SGK . It is based on locating the same fixed point at different levels of coarse description, and examining the slow eigenvalues and corresponding eigenvectors of the linearization of the coarse timestepper at these fixed points (see also the discussion in ref:MMK ).

The separation of characteristic time scales describing the underlying macroscopic and microscopic dynamics – for example, the (long) time scale on which the spatial density distribution of a population of bacteria moving in a chemoattractant profile changes versus the (short) mean runtime of a single cell – allows taking time steps in evolving the macroscopic state that are “long” in comparison with the “fast” microscopic time scales in the problem, resulting in a computationally - efficient time evolution of the macroscopic state. We demonstrate this assumption of the separation of time scales explicitly in Section V. Figure 1 shows a cartoon sketch of the coarse integration (CI) procedure, where the solid trajectories denote the restriction of the full dynamics onto a suitable low-dimensional subspace (see Section IV.1 for how a low dimensional representation of the macroscopic spatial distribution, $`\mu _0(x,t)`$, is constructed). At each CI step, the coarse-integrated solution is lifted from the lower dimensional space into the full microscopic state space and evolved according to the Monte Carlo scheme, allowing the error incurred during the coarse time step to relax to the slow manifold parameterized by the cell density. Because the internal state $`\{\stackrel{}{u},\stackrel{}{s}\}`$ of cells, and their directions of motion, $`d`$, are ignored in obtaining $`\mu _0(x,t)`$, they must be suitably reinitialized to construct the initial condition for the next MC step (see Section IV.2).

### IV.1 SVD-based restriction scheme

Below, we outline our coarse integration scheme which operates on data from singular value decomposition of the distribution of cell positions from Monte Carlo time evolution results. In the following, the solution at a “reporting step” refers to the sorted cell positions at specified time intervals, $`T_{step}`$. Note that this time interval should be distinguished from a single Monte Carlo iteration, of which it contains a large number.

1. The Monte Carlo description is simulated for $`k`$ reporting steps, corresponding to a total time interval equal to $`T_{settle}`$. Step by step results are accumulated for a further $`m`$ reporting steps, corresponding to a total time interval equal to $`T_{fit}`$.
2. At each reporting step, we sort the cell positions, $`x_i(j)`$, to obtain $`X_i(j)`$, $`i=k+1,\mathrm{},k+m`$, and $`j=1,\mathrm{},N_{\mathrm{cells}}`$. The matrix $`A_{\mathrm{}}`$ is constructed (columnwise): $`A_{\mathrm{}}=\{\stackrel{}{X}_{k+1},\stackrel{}{X}_{k+2},\mathrm{},\stackrel{}{X}_{k+m}\}`$. Figures 2 and 3 show the chemoattractant profile, and representative time sequences of sorted cell positions and corresponding histograms of the spatial density distributions, respectively. However, we find it more convenient to work with the cumulative probability distribution function, defined as

$$P(X,t)=\frac{1}{L}_{x_{min}}^X\mu _0(x^{},t)𝑑x^{},$$
(17)
which is related to the sorted cell positions at that time according to

$$j/N_{cells}=P(X_i(j),t_i),$$
(18)
where $`Lx_{max}x_{min}`$ is the spatial domain. We will use approximations to $`P(X,t)`$, or equivalently the sorted cell positions, $`\stackrel{}{X}_i`$, as the macroscopic variable ref:Gear\_distribution .
3. A low dimensional representation of the sorted cell positions is constructed in terms of orthonormal numerical basis functions, $`\{\stackrel{}{u}^{(r)}\}`$,

$$\alpha ^{(r)}(t_i)=\stackrel{}{u}^{(r)}\stackrel{}{X_i}=\underset{j=1}{\overset{N_{cells}}{}}u^{(r)}(j)X_i(j),$$
(19)
where the basis functions are obtained from singular value decomposition of the MC data as described below.

We distinguish between two approaches to obtaining the numerical basis set.

* Local basis functions: The singular value decomposition of $`A_{\mathrm{}}`$, constructed in Step 2, is computed:

    $$A_{\mathrm{}}=U_{\mathrm{}}W_{\mathrm{}}V_{\mathrm{}}^{}$$
    (20)
    The columns of $`U_{\mathrm{}}`$, given by $`\{\stackrel{}{u_{\mathrm{}}}^{(r)}\}`$, $`r=1,\mathrm{},m`$, are the numerical basis functions. We assume that they remain valid basis functions during the projected step. $`W_{\mathrm{}}`$ is a diagonal matrix, where $`\{w_{\mathrm{}}^{}{}_{1}{}^{},w_{\mathrm{}}^{}{}_{2}{}^{},\mathrm{},w_{\mathrm{}}^{}{}_{m}{}^{}\}`$ are the singular values.
    * Global basis functions: Using a single MC evolution from initial conditions to steady state, $`T_f`$, the matrix $`A_g`$ is constructed from the full data according to Step 2: $`A_g=\{\stackrel{}{X}_1,\stackrel{}{X}_2,\mathrm{},\stackrel{}{X}_M\}`$, where $`M=T_f/T_{step}`$. Singular value decomposition of $`A_g`$ gives the the global basis functions, $`\{\stackrel{}{u_g}^{(r)}\}`$. (In practice, initial data from $`t=0`$ to $`t=T_{settle}`$ are not included.)
4. We perform linear least squares extrapolation of $`\{\alpha ^{(r)}(t_i)\}`$ using $`i=k+1,\mathrm{},k+m+1`$ to obtain $`\{\alpha ^{(r)}(t_{k+m+1+p})\}`$, corresponding to a projected time equal to $`T_{proj}`$. The projected solution is given by:

$$\stackrel{}{Y}_{k+1+m+p}=\underset{r}{}\alpha ^{(r)}(t_{k+m+1+p})\stackrel{}{u}^{(r)}.$$
(21)
5. In practice, we use only the numerical basis functions associated with the dominant singular values, $`r=1,\mathrm{},r_{max}`$, where the truncated number of basis functions, $`r_{max}`$, depends on the spectrum of singular values.

The use of empirical orthogonal eigenfunctions, also referred to as the Karhunen-Loève expansion or the Proper Orthogonal Decomposition, obtained through singular value decomposition of experimental or simulation data, for model reduction in systems with spatiotemporal dynamics originates with Lorenz in the context of weather prediction ref:Lorenz  and has found widespread use in the dynamical systems context since the mid-80’s (see the monograph ref:Holmes  and references therein).

The SVD process is computationally expensive, so if the computational cost of the microscopic cell simulations were small, using global basis functions generated in a preliminary run would be less costly than using local basis functions. However, in a simulation of a large number of cells, the total microscopic simulation cost for the cell population may be much larger that for the SVD process. For example, in the numerical experiments reported here the CPU time for one simulation step of an individual cell was 2.07 microseconds, while the time for one SVD step was 22.4 milliseconds. However, in a simulation of $`10^4`$ cells using local basis functions, the microscopic integration routine was called $`4\times 10^9`$ times versus 16 calls on the SVD routine so that 0.0041% of the time was spent on SVD.

### IV.2 Reinitialization of internal variables

The question of consistent schemes for combining microscopic and macroscropic descriptions of a physical process is central to multiscale modeling. Here, to alternate each CI step with MC, in principle allowing relaxation of the numerical error in the projective step, we must define an appropriate reinitialization procedure: We know the position of each cell, $`x`$, after the CI step, but we discard all information about its internal state, $`\{\stackrel{}{u},\stackrel{}{s}\}`$, given by the values of the signal transduction variables and the flagellar states, and its direction of motion, $`d_j`$.

Our reinitialization protocol is empirically motivated: the signal transduction variables are set equal to their local equilibrium values, $`\stackrel{}{u}_j=\{0,f\left(S(x_j)\right)\}`$, where $`j`$ refers to the $`j^{th}`$ cell. This allows the subsequent response of each cell to be within the most sensitive range of the motor gain. For simplicity, all (six) flagella are restarted in the CW state, $`\stackrel{}{s}_j=\{0,0,\mathrm{},0\}`$. In Figure 4, we show relaxation of the ratios of the numbers of right, left and tumbling cells at a set location along the chemoattractant profile. These ratios at steady state for a flat chemoattractant profile are

$$\frac{\rho _R+\rho _L}{\mu _0}=B,$$
(22)
where the right and left moving ratios are equal. For the present choice of numerical parameters, $`B0.86`$ and $`\rho _R=\rho _L0.43`$. Other flagellar reinitialization schemes lead to similarly rapid relaxation rates.

## V Numerical results

### V.1 Low-dimensional representation

In Figures 5 and 6, we show the first four dominant global basis functions and singular values obtained from MC evolution of $`N_{cells}=10^4`$ from initial conditions to steady state. Figure 7 shows the mean coefficients over $`N_{\mathrm{MC}}=10`$ MC realizations

$$\overline{\alpha }^{(r)}(t)=\frac{1}{N_{\mathrm{MC}}}\underset{k=1}{\overset{N_{\mathrm{MC}}}{}}\alpha _k^{(r)}(t),$$
(23)
and the root-mean-square of the distribution of coefficients, given by

$$\sigma _\alpha ^{(r)}(t)=\left\{\frac{1}{N_{\mathrm{MC}}1}\underset{k=1}{\overset{N_{\mathrm{MC}}}{}}\left[\alpha _k^{(r)}(t)\overline{\alpha }^{(r)}(t)\right]^2\right\}^{1/2}.$$
(24)
Hence, the error bars in this figure denote the expected statistical variation in the values of these coefficients as a function of time.

In Figure 8, we show the low-dimensional representation of the coarse-integrated solution for $`T_{settle}=T_{fit}=5\times 10^2`$ and $`T_{proj}=5\times 10^3`$. For each CI step, after evolving the Monte Carlo for $`T_{settle}`$, we use a linear fit to $`\{\alpha ^{(1)}(t),\alpha ^{(2)}(t),\alpha ^{(3)}(t),\alpha ^{(4)}(t)\}`$ during the interval $`T_{fit}`$ to project the solution forward in time by $`T_{proj}`$. The points correspond to averages over $`N_{\mathrm{MC}}=5`$ realizations, and the error bars give the rms of the distribution of these coefficients. For reference, we have included the average MC coefficients in this figure. These results indicate that for the higher order coefficients, $`\{\alpha ^{(3)}(t),\alpha ^{(4)}(t)\}`$, whose coarse dynamics are described by a shorter characteristic time-scale, a higher order time integration scheme would be more effective than the explicit Euler method used here ref:RM . Consequently, although for intermediate times the difference between the CI and MC results the straightforward CI scheme has captured the macroscopic dynamics of the solution.

These results also illustrate the implicit assumption of separation of time scales: The coefficients $`\left\{\alpha ^{(r)}(t)\right\}`$, governing the macroscopic behavior of the system, vary on a time scale of $`𝒪(10^4)`$ units, while the longest microscopic time scale (adaptation time of the signal transduction model) is $`𝒪(10^2)`$.

### V.2 Analysis of Errors

For comparison of solutions at different values of the CI parameters, $`(T_{settle},T_{fit},T_{proj})`$, with the MC, we construct the following measure of relative error:

$$\epsilon (t)=\left\{\frac{_{r=1}^4\left[\alpha ^{(r)}(t)\overline{\alpha }^{(r)}(t)\right]^2}{_{r=1}^4\left[\overline{\alpha }^{(r)}(t)\right]^2}\right\}^{1/2}.$$
(25)
The coefficients $`\{\alpha ^{(1)}(t),\alpha ^{(2)}(t),\alpha ^{(3)}(t),\alpha ^{(4)}(t)\}`$, used in computing the relative error of the CI solution are obtained as inner products of the solution with the global basis set, regardless (a) of whether global or local basis sets were used in the restricting/lifting step, and (b) the dimensionality of the restricted space.

In Figure 9, we first plot the average relative error for $`N_{\mathrm{MC}}=10`$ MC runs, given by:

$$\overline{\epsilon }(t)=\frac{1}{N_{\mathrm{MC}}}\underset{k=1}{\overset{N_{\mathrm{MC}}}{}}\epsilon _k(t).$$
(26)
The error bars denote the error on the mean, given by $`\sigma _{\overline{\epsilon }}=\sigma _\epsilon /\sqrt{N_{\mathrm{MC}}}`$, where $`\sigma _\epsilon `$ is the root-mean-square of the distribution of the error

$$\sigma _\epsilon (t)=\left\{\frac{1}{N_{\mathrm{MC}}1}\underset{k=1}{\overset{N_{\mathrm{MC}}}{}}\left[\epsilon _k(t)\overline{\epsilon }(t)\right]^2\right\}^{1/2}.$$
(27)
We similarly compute the relative error for CI solutions, obtained using either (i) global or (ii) local SVD basis functions in each step. Figure 10 shows the mean errors for solutions with CI parameters given by $`(T_{settle},T_{fit},T_{proj})=10^3\times (\frac{1}{2},\frac{1}{2},5),10^3\times (1,1,10)`$ and $`10^3\times (2,2,20)`$. Note that for these parameter values, the efficiency of the CI scheme remains the same. In this figure, solid lines show results obtained using global SVD basis functions, $`\{\stackrel{}{u}_g^{(1)},\stackrel{}{u}_g^{(2)},\stackrel{}{u}_g^{(3)},\stackrel{}{u}_g^{(4)}\}`$, in the CI steps. Dashed lines show these results obtained using local SVD basis functions. In the latter case, since singular value decomposition is applied to a short segment of the MC, we found the higher order local basis functions to be relatively “noisy”; hence, $`\{\stackrel{}{u}_l^{(1)},\stackrel{}{u}_l^{(2)}\}`$ were used in each CI step. Note that the rms of the distribution of errors obtained for CI solutions using local basis function is larger; in this case, the basis functions themselves are also subject to statistical variations. We find that the relative errors using global SVD basis functions are generally smaller than those corresponding to CI solutions obtained using local SVD basis functions. Using either global or local basis functions in the CI step, the smallest relative error is achieved with the first two sets CI of parameters, corresponding to projected time intervals short enough to resolve the macroscopic dynamics of the higher order coefficients of the low-dimensional representation of the solution.

As a benchmark, in Figure 11 we show the error obtained for $`(T_{settle},T_{fit},T_{proj})=10^3\times (1,1,0)`$, using global SVD functions $`\{\stackrel{}{u}_g^{(r)}\}`$, with $`r=1,\mathrm{},4`$ (dotted line) and $`r=1,\mathrm{},8`$ (solid line). These results show the error incurred in restricting the spatial distribution of cell positions to the low dimensional representation given by $`\{\alpha ^{(r)}\}`$, without the contribution from the projective step. First, we note that the maximum error obtained from the restriction to $`\{\alpha _g^{(1)},\mathrm{},\alpha _g^{(4)}\}`$ and $`\{\alpha _g^{(1)},\mathrm{},\alpha _g^{(8)}\}`$ are comparable, indicating that using a higher dimensional restriction is not significantly advantageous here. Secondly, this error is comparable to the total error in Figure 10(b), which includes projective integration. Given that the total error increases with a longer projective time step, $`T_{proj}`$, as shown in Figure 10(c), these results point out that the “optimal” $`T_{proj}`$, for which the errors due to the projective time step and restriction are separately comparable, is achieved in Figure 10(b).

Finally, in Figure 12, we compare CI errors using two different reinitialization schemes and find them to be equivalent.

## VI Concluding remarks

We have demonstrated, through a coarse projective integration algorithm, that short bursts of appropriately initialized microscopic simulations can be used to simulate the macroscopic evolution of the cell density distribution in bacterial chemotaxis. The outer level of our computational structure was a traditional, continuum Euler integration scheme; the time-derivatives of the cell density field it employed, however, were estimated from short time evolution intervals of the Monte-Carlo description of this process, and not evaluated from a known, macroscopic expression. This approach leads to a general framework for the computer-assisted analysis of systems whose dynamics are given at a microscopic/stochastic “fine” level, but for which we require averaged, macroscopic information.

It is interesting to discuss the benefits and shortcomings of such a procedure. If accurate closed chemotactic equations exist, one should use them instead of the two-tier modeling we propose here. In addition to the insight gained from exact or approximate analytical solutions, computer-aided time-evolution or bifurcation analysis using explicit equations is more efficient than kinetic Monte Carlo simulations. However, if such model equations are not available, our hybrid computational approach can become more economical than long-time direct Monte Carlo simulation. Furthermore, when one is interested in qualitative transitions or bifurcations of the macroscopic behavior, marginally stable or unstable stationary states may be difficult to identify through direct microscopic simulations, while coarse timestepping holds promise when combined with traditional bifurcation algorithms ref:MMK ; ref:SGK ; ref:RM .

It appears, therefore, that a modeler would ultimately gain in obtaining quantitative computational answers efficiently, but perhaps lose in the fundamental understanding of a physical process that macroscopic model equations would offer. Therefore a promising research direction is to use such algorithms to test the validity of explicit closures that assume slaving of certain higher order moments of the evolving distribution to lower order ones. In the chemotaxis context, when signal transduction and corresponding motor response of the cell to an external signal are taken into account as we do here, macroscopic model equations derived systematically from the microscopic description do not exist. The assumption implicit our the coarse integration algorithm was that the macroscopic description closes at the level of the spatial density distribution. In one spatial direction, for example, it would be interesting to learn when this assumption breaks down, requiring separate evolution of right-going, left-going and tumbling cell density distributions.

More generally, we believe that the approach illustrated here for coarse time integration of a Monte Carlo description of bacterial chemotaxis, is broadly applicable with different types of microscopic simulators describing, for example, Brownian, Lattice-Boltzmann, or molecular dynamics, and leading to emergent macroscopic dynamic behavior. Such a hybrid scheme allows efficient simulation of the macroscopic behavior, and may provide insight into macroscopic model equations.

## ACKNOWLEDGEMENTS

The work reported here was supported by the Princeton Council on Science and Technology Fellowship (SS), NIH Grant #GM29123 and NSF Grant #DMS 0074043 (HGO), and AFSOR (Dynamics and Control) and NSF-ITR grants (CWG and IGK).


text:
# 1 Introduction

## 1 Introduction

The properties of strange mesons within nuclear matter has been subject to numerous investigations. Especially nuclear collisions at energies near or below the production threshold of nucleon-nucleon collisions should be very sensitive to the kaon properties in matter . Early theoretical approaches based on effective chiral Lagrangians  predicted an attractive scalar potential which together with an isovector potential leads to a strong attractive K- and a moderately repulsive K+ potential. These potentials depend only weakly on the momentum. Indeed kaonic atoms require a strong attractive K- potential and also the large K- rates observed in heavy-ion collisions, carried out by the FOPI  and KaoS  collaborations, seemed to support these predictions. However, more sophisticated theoretical investigations  demonstrated a strong momentum dependence of the potentials which even became repulsive at large density and momentum for both K+ and K- mesons. These potentials do not comply with the early analyses of the measurements of K- production. However one has to keep in mind that the elementary cross sections for K- production used are not very well known from experiment.

The knowledge of the production mechanism of kaons is a necessary condition to probe the theoretical predictions for the potentials. In this respect it is useful to study also nucleon-nucleus collisions as a further source of information.

In nucleus-nucleus collision antikaons have been observed at ion bombarding energy of 1.5 GeV per nucleon . This energy is far below the threshold energy of 2.5 GeV for nucleon-nucleon collisions. This nucleon-nucleon production channel is unimportant in heavy-ion collisions due to the smallness of the kaon pair production cross section  even if Fermi motion would help to overcome the threshold. The antikaons can however be produced by multistep-scattering processes. In the sequence of the collisions combined with Fermi motion strangeness transfer via $`\pi `$Y $``$ NK- reactions can take place with their comparably large cross sections.

Only a few data for antikaon production in nucleon-nucleus collision near threshold are available. We compare our calculations with the KaoS data . A comparison to additional existing antikaon data from FHS  below and KEK-PS  above threshold seems not to be reasonable, because of their very strong kinematical constraint. In nucleon-nucleus collisions the above mentioned $`\pi `$Y channel is a rather improbable three step process, because the single incoming proton alone has to produce both reaction partners before. Also the second chance collision $`\pi `$N $``$ NK+K- channel has the small cross section of the kaon pair production. Therefore, in a second chance collision the antikaon could mainly be produced via the NY $``$ NNK- channel which is open for an incident proton energy larger than 1.73 GeV. Additional Fermi motion may lead to a further reduction of the threshold.

Therefore we draw our attention to the nucleon-hyperon channel NY $``$ NNK- cross section. This cross section was already calculated early in the one-pion exchange model  and later in one-boson exchange approximation  and was found to be unimportant for heavy-ion collisions. Here we reevaluate the cross section within a different approach avoiding the uncertainties arising from the badly known formfactors when applying an effective perturbation theory. Our results roughly agree with the cross sections obtained in ref. . In the relevant kinetic energy region the cross section reaches nearly one mb and is about 50 times larger than the related cross sections of the reactions NN$``$ NYK+. As the ratio of hyperons as well as of kaons to the participating nucleons is about 10-4 this would lead to a K- to K+ ratio of about 10-2 which is nearly the magnitude of the measured ratio at 2.5 GeV beam energy. Thus, this channel should compete with the NN and $`\pi `$Y production channels.

In addition we mention that the in-medium cross sections may considerably differ from their vacuum values. This was pointed out e.g. in ref.  where a considerable enhancement of the pion-hyperon channels have been predicted. These results are based on coupled channel calculations and are connected with a shift of the masses of the $`\mathrm{\Lambda }(1405)`$ and $`\mathrm{\Sigma }(1385)`$ resonances in nuclear matter. This effect will not be considered here.

## 2 Elementary cross sections

To estimate the NY $``$ NNK- cross sections we start with the Feynman diagrams shown in Fig. 1 which are similar for NY $``$ NNK- and NN $``$ NYK+ processes. The K± mesons are generated by the subprocess where the intermediate meson (here a pion) interacts with the second baryon. This subprocess cannot be calculated for certainty because many resonances contribute with unknown coupling constants. The results of such one-boson exchange reactions can be found in ref. . Here we will approach the study of these cross sections differently by making the assumption that the most important meson exchange is that of a pion. The $`\pi `$B $``$ K± cross sections are known experimentally from which we can extract the square of the transition matrix elements $`T_{\pi \mathrm{B}}`$ illustrated by the hatched areas in Fig. 1. Then, these values are used to calculate the cross sections in accordance with the diagram in Fig. 1. We calculate both the kaon and the antikaon production in order to check the method since the cross sections for NN $``$ NYK+  has been calculated and adjusted to the partially known pp cross sections .

Thus we consider the process N1+B2 $``$ N3+B4 \+ K($`\overline{\mathrm{K}}`$), where the symbol N denotes a nucleon with isospin $`I_1=I_3=1/2`$ and one of the symbols B2 or B4 stands for a hyperon. The pion-nucleon coupling in the left hand vertex in Fig. 1 is described by the lagrangian

$$_{\pi NN}=g\overline{\psi }\gamma ^5\stackrel{}{\tau }\psi \stackrel{}{\pi }$$
(1)
with $`\psi `$ denoting the nucleon field, $`\stackrel{}{\pi }`$ the pion field, $`\stackrel{}{\tau }`$ the isospin Pauli matrix, $`\gamma ^5`$ is a Dirac matrix, and g = 13.6 fixes the pion-nucleon coupling constant. The spin and isospin averaged cross section for the process at the center-of-mass energy $`\sqrt{s}`$ reads

$`\sigma _\mathrm{K}`$ $`=`$ $`{\displaystyle \frac{g^2}{2\lambda ^{1/2}(m_1^2,m_2^2,s)}}{\displaystyle \frac{1}{(2\pi )^5}}{\displaystyle \frac{1}{8(2I_2+1)}}`$
$`\times `$ $`{\displaystyle \underset{spin,isospin}{}}|<I_3|\stackrel{}{\tau }|I_1>|^2{\displaystyle d^4p_\pi \frac{d^3𝐩_4}{2p_4^0}\frac{d^3𝐩_K}{2p_K^0}(\overline{u}_3\gamma ^5u_1)^2}`$
$`\times `$ $`\delta ((p_1p_\pi )^2m_3^2)\delta ^4(p_\pi +p_2p_3p_K){\displaystyle \frac{1}{p_\pi ^2m_\pi ^2}}T_{\pi B_2:B_4K}^2,`$
where the indices of the momenta $`p`$ and masses $`m`$ are those used in Fig. 1, $`u`$ denotes the nucleon spinor, and $`\lambda (a,b,c)=(sab)^2(2ab)^2`$ is the triangle function, and $`I_2`$ stands for the isospin of particle B2. The integral over the outgoing momentum $`p_3`$ of particle N3 has been substituted by the pion momentum $`p_\pi =p_1p_3`$. The symbol $`T_{\pi B_2:B_4K}`$ represents the encircled part in Fig. 1 which determines the kaon (antikaon) production in pion-baryon collisions the cross section of which is given by

$`\sigma _{\pi B_2}`$ $`=`$ $`{\displaystyle \frac{1}{2\lambda ^{1/2}(m_\pi ^2,m_2^2,s_{\pi B_2})}}{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle \frac{1}{6(2I_2+1)}}`$ (3)
$`\times `$ $`{\displaystyle \underset{spin,isospin}{}}{\displaystyle \frac{d^3𝐩_4}{2p_4^0}\frac{d^3𝐩_K}{2p_K^0}\delta ^4(p_\pi +p_2p_3p_K)T_{\pi B_2:B_4K}^2}.`$
This cross section depends on the square of the center-of-mass energy $`s_{\pi B_2}=(p_\pi +p_2)^2`$ with the on-shell condition $`p_\pi ^0=\sqrt{m_\pi ^2+𝐩_\pi ^2}`$. Notice that in Eq.(2) the pion momentum $`p_\pi `$ is off-shell.

Inserting Eq.(3) into Eq.(2), summing over the spin quantum numbers of the nucleons and integrating out the time component $`p_\pi ^0`$ of the pion momentum we obtain the cross section

$`\sigma _\mathrm{K}`$ $`=`$ $`{\displaystyle \frac{g^2}{\lambda ^{1/2}(m_1^2,m_2^2,s)}}{\displaystyle \frac{1}{(2\pi )^3}}{\displaystyle \frac{d^3𝐩_\pi }{\sqrt{m_3^2+𝐩_\pi ^2}}\frac{p_1p_3m_1m_3}{(p_\pi ^2m_\pi ^2)^2}}`$ (4)
$`\times `$ $`f^4(p_\pi )\lambda ^{1/2}(m_\pi ^2,m_2^2,s_{\pi B_2}){\displaystyle \underset{I_\pi }{}}\sigma _{\pi B_2}(s_{\pi B_2}),`$
where the formfactor

$$f(p_\pi )=\frac{\mathrm{\Lambda }^2m_\pi ^2}{\mathrm{\Lambda }^2p_\pi ^2}$$
(5)
with $`\mathrm{\Lambda }`$ = 1.6 GeV  has been introduced.

The cross sections $`\sigma _{\pi B_2}`$ needed in Eq.(4) for the reactions $`\pi `$Y $``$ NK- can be derived from the measured inverse reactions K-p $``$ $`\mathrm{\Sigma }^+\pi ^{}`$, $`\mathrm{\Lambda }\pi ^0`$, $`\mathrm{\Sigma }^{}\pi ^+`$ and K-n $``$ $`\mathrm{\Lambda }\pi ^{}`$, $`\mathrm{\Sigma }^{}\pi ^0`$ which are given in ref. . The K0N cross sections are derived by isospin reflection.

The results for the antikaon production in N$`\mathrm{\Sigma }`$ and N$`\mathrm{\Lambda }`$ collisions are represented by the open symbols in Fig. 2. The lines through the symbols show fits with the standard parametrization $`\sigma (ss_0)^a/s^b`$. The K+ cross section is compared to the parametrization of ref.  which gives in the low energy region nearly the same cross section but underestimates the cross section above 100 MeV. The parametrization in ref.  is based on a model which uses diagrams which have the same structure as those of Fig. 1. In that investigation an exchange of $`\pi `$, $`\eta `$ or $`\rho `$ mesons is included, and it is shown that the pion gives the main contribution, a fact that supports our assumption of the dominance of the pion exchange.

In Fig. 3 we compare the curves in Fig. 2 to the results of previous calculations . There are only small deviations from the calculations in the one-boson exchange approximation  which also shows the dominance of the N$`\mathrm{\Sigma }`$ channel. The result  has a weaker energy dependence but provides comparable cross section in the relevant energy region of about 100 MeV above threshold. In this investigation it was also shown that this channel could contribute about 10% to the K- production in heavy-ion collisions. In our following calculations we find that even half of the antikaons can stem from the NY channel in proton-nucleus collisions.

## 3 Comparison to data

It is our aim to study the role of the NY $``$ NNK- reaction in pA collisions. Usually these channels are not included in standard analyses of those reactions which consider only the elementary BB and $`\pi `$B collisions. To study this question we additionally incorporate the NY channels into a transport model calculation which is based on the Boltzmann-Ühling-Uhlenbeck equation . Furthermore, the production rate of the K- mesons also depends sensitively on the attractive K- potential. Therefore, the inclusion of different production channels will effect predictions on the size of this potential when derived from comparison with data.

In Fig. 4 we present the differential cross sections at a laboratory angle of 40 for collisions of protons with 12C and 197Au at 2.5 GeV beam energy. The dotted lines are calculated without using potentials for the kaons and antikaons. These calculations underestimate clearly the measured data . The attractive antikaon potential,

$`V_{\overline{K}}=0.08GeV{\displaystyle \frac{n}{n_0}}`$ (6)
in addition with the NY channels leads to an increase of the cross section as shown by the full lines. Such a potential improves the agreement with the data. The dashed curves are results where the NY$``$K- channels has been excluded. Disregarding these channels the cross section diminishes by about 50% in p+Au collisions. This shows the importance of the NY$``$K- channels when one intends to determine the K- potential. For the light C target the influence is much smaller as the hyperons have a smaller chance to collide with further nucleons before leaving the reaction zone.

Finally, we compare in Fig. 5 our calculations for K+ and K- production with data obtained by the KaoS collaboration  for proton-nucleus collisions at bombarding energies of 2.5 GeV and 3.5 GeV on C and Au targets and K meson emission angles of $`40^{}`$ and $`56^{}`$. The kinetic beam energy $`T_{\mathrm{kin}}=2.5`$ GeV was close to the production threshold in nucleon-nucleon collisions. In the calculations we have used the parametrizations of ref.  and obtained kaon cross sections which are slightly smaller than the data for the gold target for both angles.

The kaon and antikaon production yields are related in the threshold region since nearly all of the antikaons are created in collisions of pions and nucleons with hyperons, the number of which equals to the kaon number because of strangeness conservation. This fact also holds if a chemical equilibrium between K- and Y is reached via the K-N $``$ $`\pi `$Y reaction . Therefore it is interesting to compare the ratio of the angle integrated cross sections for K- to those of K+. In Tab. 1 we display the calculated values and the experimental results obtained within the measured phase space by the KaoS collaboration . Our calculations overestimate this ratio for 2.5 GeV but underestimate it for the higher energy of 3.5 GeV.

Tab. 1 Comparison of calculated and experimentally obtained antikaon to kaon cross-section ratios for proton collisions on carbon $`R(C)`$ and gold $`R(Au)`$. Experimental values with errors of about 20% are taken from .
$`T_{kin}`$(GeV) $`R_{exp}(C)`$ $`R_{calc}(C)`$ $`R_{exp}(Au)`$ $`R_{calc}(Au)`$    2.5 0.0085 0.013 0.0074 0.010  3.5 0.028 0.024 0.027 0.024

## 4 Conclusions

We have calculated the cross sections of antikaon-production in near-threshold proton collisions on carbon and gold targets. Comparison of our calculations with both kaon and antikaon data from the KaoS collaboration  were made at beam energies of 2.5 and 3.5 GeV for laboratory angles of $`40^{}`$ and $`56^{}`$. Including the NY $``$ NNK- channels nearly doubles the antikaon yield for collisions of protons on heavy targets like Au. No significant influence was found for the light carbon target. Calculated ratios of K- to K+ cross sections came also reasonably close to the data. We conclude that the NY $``$ NNK- channels have to be included in realistic calculations, especially for heavy targets, in order to study the properties of strange mesons in nuclear matter at normal nuclear density.

### Acknowledgments

Valuable discussions with Werner Scheinast and Burkhard Kämpfer are acknowledged. We thank Christoph Hartnack for informing us on similar results in his study. This work was supported in part by the BMBF grant 06DR121.


text:
# Classifications of the Host Galaxies of Supernovae, Set II

## 1 Introduction

Archaeologists derive much of their knowledge of ancient civilizations from digging in ancient cemeteries. By the same token astronomers can learn a great deal about the evolutionary history of galaxies from the frequency, and distribution of, supernova (SN) outbursts which signal the occurrence of stellar deaths.

Until quite recently the discovery of a supernova was a rather rare and haphazard event. However, the advent of modern systematic and automated searches now enables us to study the occurrence of supernovae (SNe) in a much more systematic fashion. The largest such automated search program is currently the Lick Observatory Supernova Search (LOSS) with the 0.75-m Katzman Automatic Imaging Telescope (KAIT). LOSS, which started in 1997 and recently expanded to the Lick Observatory and Tenagra Observatory Supernova Searches (LOTOSS; Schwartz et al. 2000), has been described by Treffers et al. (1997), Li et al. (2000), Filippenko et al. (2001), and Filippenko (2003). In order to derive the full benefit from LOSS, one needs to know both the spectroscopic types of all of these automatically discovered SNe and the morphological classifications of the SN host galaxies. By combining the new morphological classifications of 231 host galaxies with those of 177 hosts in a previous paper, we now have available morphological types for the host galaxies for a total of 408 SNe that were discovered (or independently rediscovered) during the course of the LOSS/LOTOSS surveys. This large and homogeneous data set will be discussed below.

## 2 New Morphological Classifications

The present investigation represents a continuation of the work of van den Bergh, Li, & Filippenko (2002, hereafter Paper I), in which the Hubble types of 177 SNe discovered during the course of LOSS were classified on the DDO System (van den Bergh 1960a,b,c). In the DDO system galaxies are assigned a Hubble type \[E0…E7, S0, Sa, Sb, Sc, Ir\], a luminosity class \[I…V\], and to a form family \[S, S(B) or SB\]. An “n” denotes smooth-armed spirals, a “\*” denotes patchy-armed spirals, a “t” denotes tidally distorted objects, and a “:” is used for uncertain values.

New classifications of the host galaxies of SNe that were discovered (or rediscovered) during the course of LOSS/LOTOSS are listed in Table I. For brevity, in the table N = NGC, I = IC, U = UGC, and M = MCG. In addition to the KAIT CCD images themselves, we also had available either the red images of the Palomar Sky Survey (POSS), or the SERC-J survey in blue-green light. For southern galaxies red IIIaF + RG610 images were also available from the UK Schmidt Equatorial Red Survey or from the UK Schmidt Second Epoch red survey. For additional information on the observational database the reader is referred to Paper I. For some galaxies only images in red light were available. This posed no problem for elliptical galaxies, but required some “mental extrapolation” for the spirals that had no $`B`$ images available. All morphological classifications were made by one of the authors (SvdB). Because of the limited dynamic range of the images it was often not possible to distinguish with confidence between E and S0 galaxies. Although the majority of galaxies with redshift $`z0.0`$ fit comfortably within the Hubble classification system, there are some that do not. No attempt was made to shoe-horn peculiar objects into the Hubble scheme.

## 3 Classifications of Galaxy Morphology and Supernova Type

### 3.1 Luminosity Classifications

Table 1 contains classifications on the DDO system for 44 galaxies that had also been classified previously by van den Bergh (1960c). One of these objects (NGC 4653) was apparently misidentified in that paper and has been excluded from the comparisons given below. For the remaining 24 galaxies that had luminosity classifications the mean difference (in the sense old minus new) was $`0.02\pm 0.09`$ luminosity classes. The root-mean-square dispersion in the difference between the old and the new classifications was 0.43 luminosity classes, so that the intrinsic dispersion of a single luminosity classification on the DDO system is $``$0.3 luminosity class. This close agreement is particularly pleasing because our comparison includes objects with uncertain luminosity classifications (that had been marked “:”).

### 3.2 Classifications of Hubble Types

Of the host galaxies of SNe that are listed in Table 1, 36 were also assigned Hubble types by van den Bergh (1960c). Thirty galaxies (83%) were independently assigned the same Hubble type in 1960 and in 2003, 4 (11%) were classified half a Hubble class later in the 2003 investigation, and 2 (6%) were assigned half a Hubble class earlier in the 2003 classification. For no object did the 1960 and the 2003 classifications differ by 1.0 or more Hubble classes. It is remarkable (and satisfying!) to see that morphological classifications made independently (but by the same person) with a time separation of over four decades exhibit no measurable systematic differences.

### 3.3 Spectral Classifications

The classifications of SN spectra that are given in Table 1 were drawn from the IAU Circulars. Supernovae of Type Ia were divided into “normal” and “peculiar” categories on the basis of careful inspection of the spectroscopic information in the IAU Circulars: objects that showed the strong Si II $`\lambda `$5970 feature or Ti II absorption lines near 4200 Å (which are evidence for a subluminous SN 1991bg-like event; Filippenko 1997), or weak Si II $`\lambda `$6350 absorption or strong Fe III absorption (which indicates a possibly overluminous, SN 1991T-like event; Filippenko 1997) were classified as “peculiar” SNe Ia.

## 4 Discussion

### 4.1 Possible Massive Supernovae in E and S0 Galaxies

It is generally believed (e.g., van den Bergh & Tammann 1991) that only SNe of Type Ia occur among the old stellar populations that inhabit elliptical galaxies. The greatly enlarged modern database of accurate DDO morphological classifications and uniform SN classifications provided in the present paper, and in Paper I, allows us to test this hypothesis.

Table 2 lists the five SNe, which are not of Type Ia, that have been discovered in host galaxies that appear to be of types E or S0. It should be emphasized that late-type galaxies of unusually high surface brightness may be misclassified as being of early type. A good example of this effect is provided by the high surface brightness face-on late-type spiral NGC 3928 = Markarian 190 (van den Bergh 1980) which was misclassified as an E0 in the Second Reference Catalogue (de Vaucouleurs, de Vaucouleurs, & Corwin 1976). It would be of interest to obtain large-scale images of the galaxies that are listed in Table 2 to see if any of these galaxies were misclassified, or show evidence for a recent merger with an object of later type that might have contained massive stars. Alternatively, a few early-type galaxies with late-type companions might have picked up some gas during tidal encounters. In this connection it is noted that Sandage & Bedke (1994) refer to very subtle dust patches in NGC 2768. Massive stars and SN progenitors might, of course, also have formed recently in the gas that was so acquired. Finally, some of the SNe II that appear to be located in E or S0 galaxies might actually be situated in dim late-type companions to these objects.

On the basis of the data listed in Table 2 it appears that not more than $``$3% of all SNe with massive progenitors occur in E and S0 galaxies. However, note that some of our E/S0 classifications might be erroneous because the SN host galaxy was a late-type object with an unusually high surface brightness. Such a high surface brightness might, for example, result from an intense burst of star formation triggered by a rather recent tidal encounter. As listed in the notes of Table 2, out of the 5 galaxies mentioned, two (UGC 2836 and NGC 3720) are possible Sa galaxies, and two (IC 2461 and IC 3203) are viewed edge-on and have somewhat uncertain classifications (note, for example, NED111NED is the NASA/IPAC Extragalactic Database, http://nedwww.ipac.caltech.edu. lists the two galaxies as Sb and Sbc, respectively). Even NGC 2768 might not be a convincing early-type galaxy, being classified as E3/Sa; moreover, it is a LINER and may host a nuclear starburst. It is thus inconclusive whether any SNe with massive progenitors have been discovered in normal early-type (E/S0) galaxies from this study.

### 4.2 Relative Supernova Numbers in Galaxies of Differing Types

The number of SNe of various types in Table 1 is listed in Table 3. In Table 4 these new data are combined with those that we had previously published in Paper I. Inspection of these data reveals the following trends:

1. A Kolmogorov-Smirnov test shows that the distribution of 117 SNe Ia and 119 SNe II (excluding SNe IIn) over host galaxy type is different at 99.7% confidence. The observed difference is in the expected sense, with the SNe II (which are believed to have massive core-collapse progenitors) being most common in late-type galaxies of Hubble types Sb–Sc, whereas SNe Ia are observed to occur in hosts of all Hubble types. A small number of possible cases in which SNe II appear to have occurred in environments that seem to be overwhelmingly composed of old stars is given in Table 2. Note, however, that such objects account for only $``$3% of our total sample of core-collapse SNe (i.e., SNe II, Ib, Ib/c, Ic).
2. A comparison between the distributions of 119 SNe II and of 15.5 SNe IIn shows no hint of any difference in the distribution of host types for these two kinds of SNe, but this result suffers from small-number statistics (an object classified as S0/Sa was counted as 0.5 S0 and 0.5 Sa). In this connection it is noted that Hamuy et al. (2003) have recently speculated that some SNe IIn might be SNe Ia that are interacting with circumstellar material. If our conclusion is correct, most SNe IIn are probably not SNe Ia masquerading as SNe II. Only if SNe IIn come from the subset of SNe Ia having the most massive possible progenitors could the two observed distributions be consistent.
3. A comparison between the distributions of normal and peculiar SNe Ia hints at the possibility that SNe Ia-pec might occur preferentially in early-type hosts. A Kolmogorov-Smirnov test shows that there is only a 7% probability that they are drawn from the same parent population of host galaxies. Since some normal SNe Ia and SNe Ia-pec occurred in S0 galaxies, which lie outside the linear E-Sa-Sab-Sb-Sbc-Sc-Ir scheme, the sample can be enlarged by dividing galaxies into (1) elliptical and S0 galaxies, and (2) spiral and irregular galaxies. Such a division has the advantage that it allows us to include in the statistics spirals classified as “S” but that could not be assigned to types “Sa,” “Sb,” or “Sc.” Table 5 shows a compilation of the frequency distribution of normal SNe Ia and SNe Ia-pec in E + S0 galaxies and in S + Ir galaxies. For the data in this table, $`\chi ^2=4.7`$. A $`\chi ^2`$ test with one degree of freedom therefore yields only a 3% chance that SNe Ia-pec have the same probability distribution over Hubble types as do normal SNe Ia.
4. There is not even a hint of a difference in the distribution of host-galaxy Hubble types between 45.5 SNe Ibc (defined here to include SNe Ib, Ib/c, and Ic) and 119 SNe II. Moreover, although the currently available data sample is rather small, there is no evidence that SNe Ib, SNe Ib/c, and SNe Ic have different distributions over host-galaxy types.
5. A Kolmogorov-Smirnov test of the data for 117 SNe Ia and 46.5 SNe Ibc shows that there is only a 0.1% probability that the host galaxies of these two types of objects are drawn from the same distribution of Hubble types.

These results are consistent with the widely held suspicion that SNe Ib, SNe Ib/c, SNe Ic, and SNe II all have massive progenitors, whereas those of SNe Ia do not.

### 4.3 Peculiar Type Ia Supernovae

The data in Table 1 show that the majority of the 16 SNe Ia with spectra resembling that of the peculiar, subluminous SN 1991bg occurred in early-type galaxies, whereas all three of those with spectra like the peculiar, possibly overluminous SN 1991T erupted in spirals of types Sb or Sc. A Kolmogorov-Smirnov test shows that there is only a 4% probability that objects like SN 1991bg and those resembling SN 1991T are drawn from the same distribution of host galaxy classification types. Our results suggest that luminous SNe Ia resembling SN 1991T probably had younger (or more massive) progenitors than do those of the SN 1991bg variety. The present data confirm and strengthen a similar conclusion by Howell (2001).

### 4.4 Supernovae of Types II-L and II-P

Ever since Minkowski (1941), it has been customary to classify SNe primarily according to their spectra. However, it is also possible to classify them from their light-curve morphology. This enables one to assign SNe II to “L” (linear) and “P” (plateau) subtypes (e.g., Doggett & Branch 1985). Schlegel (1996) and Filippenko (1997) have suggested that SNe II-P and SNe II-L can be spectroscopically distinguished on the basis of the presence or absence of H$`\alpha `$ absorption. However, this speculation still needs to be confirmed with additional data.

A recent listing of SNe II that have been assigned to the linear and plateau subtypes is given in the Asiago Supernova Catalog (http://merlino.pd.astro.it/$``$supern/). For the majority of these objects, high quality and homogeneous morphological classifications of the host galaxies, based on inspection of plates obtained with large reflectors, are given by Sandage & Tammann (1981). These data show a very similar distribution of SNe II-L and of SNe II-P over Hubble types. It is found that 17 out of 25 (68%) of all SNe II-P occur in galaxies of type Sc. This does not differ significantly from the distribution of SNe II-L, for which 8 out of 11 (73%) are located in spirals of type Sc. Luminosity classes are available for only 34 of the host galaxies of SNe assigned to types II-P and II-L. For this small sample there is a hint that SNe II-L might occur in more luminous (metal-rich?) galaxies than do SNe II-P. However, a Kolmogorov-Smirnov test of the data shows that this possible difference does not have a respectable level of statistical significance. It is therefore concluded that the currently available database does not provide evidence for a significant difference between the host populations of SNe II-P and SNe II-L.

## 5 Summary

Morphological types on the DDO system are now available for a homogeneous sample of the host galaxies of 408 SNe that were discovered (or rediscovered independently) during the course of the LOSS/LOTOSS surveys. Based on these data, we come to the following conclusions.

1. The present (2003) classifications of Hubble types and luminosity classes by van den Bergh are, to within statistical errors, found to be on exactly the same system as his DDO classifications that were made in 1959 and 1960.
2. Five SNe of Types Ib/c, II, and IIn (which are thought to have massive progenitors) are found to have occurred in host galaxies nominally classified as types E or S0. However, in each case the galaxy classification is uncertain, or newly inspected images show evidence suggesting a later classification. Thus, we find no clear evidence of SNe with massive progenitors occurring in E or S0 galaxies. It would clearly be important to obtain new, large-scale images of the five early-type galaxies to improve the classifications and to search for possible evidence of young population components.
3. The distribution of SNe Ia and SNe Ia-pec over host-galaxy Hubble types is found to differ from that for SNe II + SNe IIn at 99.7% confidence.
4. No statistically significant differences are found between the distributions over Hubble type of the host galaxies of SNe II and SNe IIn. Hence, it appears unlikely that most SNe IIn are actually SNe Ia embedded in circumstellar material.
5. Compared with normal SNe Ia, those objects classified as SNe Ia-pec are found to be more common in early-type galaxies than they are in spirals of later type. A Kolmogorov-Smirnov test shows only a 1.6% probability that the 117 SNe Ia and the 29.5 SNe Ia-pec in Table 4 are drawn from the same parent population of Hubble types. This suggests that the progenitors of normal SNe Ia may differ in age or metallicity from the ancestors of SNe Ia-pec.
6. Among SNe Ia-pec, those that resemble SN 1991T appear to be younger (or more massive) than those like SN 1991bg — i.e., all three SN 1991T-like objects occurred in spirals of types Sb to Sc, whereas 14 out of 16 objects resembling SN 1991bg were hosted by galaxies of type E, S0, or Sa. A Kolmogorov-Smirnov test shows only a 1.4% probability that the three SN 1991T-like and the sixteen SN 1991bg-like objects are drawn from the same population of host galaxy Hubble types.
7. There is only a 0.1% probability that the SNe Ia and SNe Ibc (defined here to include SNe Ib, Ib/c, and Ic) in the present sample are drawn from the same parent population. SNe Ia are found to occur in all kinds of galaxies, whereas the (presumably more massive) progenitors of SNe Ibc mostly occur in spirals of types Sb–Sc.
8. No significant differences are found between the distributions over Hubble type of the small numbers of SNe Ib, SNe Ib/c, and SNe Ic contained in the present sample.
9. SNe II-P and SNe II-L appear to occur in galaxies having similar stellar populations.

All of the data discussed above refer to the relative rates of SNe of different types. The implications of the present morphological classifications for absolute frequency determinations of SNe of different types will be the subject of a future investigation.

We are indebted the the referee (Harold G. Corwin, Jr.) for a number of helpful suggestions and for providing useful information on various galaxies. We are also grateful to Allan Sandage for information on the classification of NGC 3719 and NGC 3720, to David Branch for consultations on SN spectra, to Ryan Chornock for helpful discussions, and to Avishay Gal-Yam for valuable information that affected our discussion of the core-collapse SNe that may have occurred in early-type host galaxies. The work of A.V.F.’s group at U. C. Berkeley is supported by National Science Foundation grant AST-0307894, as well as by the Sylvia and Jim Katzman Foundation. KAIT was made possible by generous donations from Sun Microsystems, Inc., the Hewlett-Packard Company, AutoScope Corporation, Lick Observatory, the National Science Foundation, the University of California, and the Katzman Foundation.


text:
# Computation of the spectra of the quasars.

## 1 Introduction

The credibility of the standard explanations of the spectra of the quasars appears low: The Lyman lines of the atomic hydrogen appear with a lot of different redshifts generally considered as produced by a Doppler effect or by an expansion of the Universe. Both hypothesis have flaws because they require clouds of hot atomic hydrogen whose speed or stability cannot be explained using regular physics.

Studying the propagation of light in low pressure gases, one must take into account the “Coherent Raman Effect on time-Incoherent Light”(CREIL) which transfers energy by frequency shifts from the hot modes of light to the thermal modes without any blur of the images and of the spectra; the relative frequency shifts $`\mathrm{\Delta }\nu /\nu `$ due to CREIL are nearly constant.

In a previous paper , we explained that the spectra can be obtained assuming that a nearly homogeneous cloud of Lyman pumped, atomic hydrogen is perturbed by a variable magnetic field: where the field is very low, there is nearly no redshift, the absorption (or emission) lines are written visibly into the spectrum; else, the redshift simultaneous with absorption (or emission) blurs the lines which become invisible. The spectrum results not from a modulation of the absorption, but from a modulation of the redshifting power of the gas. We showed the existence of a non-linearity able to increase the contrast of the lines.

However our previous explanation requires an important variation of the magnetic field for each line, thus, the existence of a large number of magnetised satellites. Worse, the spectra show a pseudo-stochastic repartition of the redshifts of the lines, the difference of two redshifts being the product of a constant $`z_f=0.062`$ by an integer . It seems very difficult to explain such a spectroscopic regularity of the “intrinsic redshift”by the presence of objects on the line of sight; the non-linearity we had introduced provides a very simple solution.

Section 2 reminds the key spectroscopic property which provides the intrinsic redshifts, the “Coherent Raman Effect on time-Incoherent Light”(CREIL) which is an avatar of the well known “Impulsive Stimulated Raman Scattering”(ISRS).

Section 3 applies the CREIL to a halo of atomic hydrogen; to get the required resonances, H must have a non-zero orbital quantum number, so that the necessary Lyman pumping introduces a non-linearity.

## 2 The “Coherent Raman Effect on time-Incoherent Light”as a limit case of the “Impulsive Stimulated Raman Scattering”.

The “Impulsive Stimulated Raman Scattering”(ISRS) was discovered in 1968  and is commonly used to study the matter . ISRS is a light-matter interaction which transfers energy from a pulsed laser beam to a colder pulsed laser beam, by frequency shifts of the beams, the temperature of the beams being deduced from Planck’s law. ISRS is space-coherent, that is the wave surfaces are left unchanged. ISRS is a parametric effect which uses matter as a catalyst, leaving it unchanged; it may be considered as a combination of two nearly simultaneous Raman scatterings giving virtual, opposite transitions. The relaxation times of the matter used as a catalyst must be longer that the length of the light pulses; these relaxation times are the period of (at least) a Raman-allowed transition, and, in a gas, the collisional time.

Replacing the ultrashort, powerful pulses of lasers by the pulses which make the ordinary incoherent light is a large change of the power and of the length of the pulses .

The change of power has a qualitative effect: In ISRS, the Raman scattered amplitude which, by interference with the exciting amplitude, produces the frequency shifts, is an increasing function, nearly quadratic, of the exciting amplitude. At low power, all scatterings are linear: If the coherence is broken by frequent collisions (ordinary Raman), the scattered intensity is proportional to the exciting intensity; else the coherently scattered amplitude is proportional to the exciting amplitude (Refraction, that is coherent Rayleigh scattering, Bragg scattering, coherent Raman scattering and photon echoes in microwaves…). This behaviour which makes the frequency shifts independant on the intensities justifies a different name: “Coherent Raman Effect on time-Incoherent Light”(CREIL). ISRS and CREIL may be computed from the tensors of polarisability of the molecules which do not depend much on the exciting frequencies, so that, in a first approximation, the relative frequency shift $`\mathrm{\Delta }\nu /\nu `$ does not depend on the exciting frequency $`\nu `$.

As the length of the pulses, of the order of 5 nanoseconds, is larger than using lasers, the collisional time must be longer than 5 ns, so that the pressure must be lower than some Pascals. The Raman resonance period must be longer than 5 ns too, so that, to be active in CREIL, the molecules must have Raman resonances in the MHz range. Both conditions decrease the intensity of the frequency shift : therefore its observation requires a very long path which is not available in an usual laboratory.

In astrophysics, the thermal radiation, at least the 2.7K radiation, provides low temperature exciting beams whose blueshift corresponds to an amplification or a heating; the required long path is evidently available; it remains finding molecules having convenient Raman resonances. As hydrogen is the most common gas, two molecules may be considered:

\- H$`{}_{}{}^{+}{}_{2}{}^{}`$ has Raman resonances close to 30 MHz, but it has probably not been observed although it should be produced by an ionisation of H2 by the ultraviolet radiation. The reason is very simple: as H$`{}_{}{}^{+}{}_{2}{}^{}`$ is destroyed by the collisions with H2 , its life-time is long only where the pressure is low enough to allow CREIL. A simultaneous absorption and CREIL makes absorption lines as wide as the redshift; therefore the almost forbidden lines of H$`{}_{}{}^{+}{}_{2}{}^{}`$ are so wide and weak that they cannot be seen.

\- atomic hydrogen may have convenient transitions if it has been pumped by a Lyman transition so that hyperfine structures producing two photons resonances in the MHz range appear.

## 3 Propagation of light in atomic hydrogen.

Consider the propagation of light having a continuous spectrum (constant intensity $`I`$, in particular in the Lyman region), in an homogeneous atmosphere of low pressure atomic hydrogen. In the fundamental state (principal quantum number n=1), the distance between the hyperfine levels (1420 MHz) is too large. In the other states, hyperfine transitions have convenient frequencies for the Raman allowed selection rule $`\mathrm{\Delta }F=1`$, for instance : 178 MHz in $`2s_{1/2}`$, 59 MHz in $`2p_{1/2}`$ and 24 MHz in $`2p_{3/2}`$.

Set $`\mathrm{\Delta }L`$ the length of path for which the redshift is equal to the linewidth $`\delta \nu `$ of the Lyman $`\alpha `$ line, and assume that the atoms which are active in CREIL are mostly pumped by the Lyman $`\alpha `$ transition.

Set $`\mathrm{\Delta }\nu `$ the redshift along the path $`\mathrm{\Delta }L`$, which would result from a complete Lyman $`\alpha `$ absorption of the intensity $`I`$ , and suppose that, in a first approximation, the whole redshift results from the Lyman $`\alpha `$ absorption.

\- case a: If $`\mathrm{\Delta }\nu `$ is larger than the Lyman $`\alpha `$ linewidth $`\delta \nu `$ , that is if $`I`$ is large enough, $`I`$ is not fully absorbed, only the constant intensity $`\mathrm{\Delta }I`$ which produces the redshift $`\delta \nu `$ is subtracted from $`I`$, that is from the spectrum while, by the redshift, the Lyman line crosses it. Thus, the contrast of lines which have been written into the spectrum is increased.

\- case b: If, on the contrary, $`\mathrm{\Delta }\nu `$ is lower than $`\delta \nu `$, the first approximation fails, a part of the redshift must result from other Lyman absorptions or other active atoms. Assuming a low redshifting power for these effects, a long path $`\mathrm{\Delta }L`$ is necessary to get the redshift $`\delta \nu `$, so that the absorption of all lines is strong.

If the intensity $`I`$ is constant and high, except for a single absorption line, the redshift and the absorption are constant (case a), except at a coincidence of the line with a Lyman line; at this coincidence, the redshifting power decreases (strongly if case b is reached), so that the absorption of all lines of the gas is increased; similarly, a written emission line increases the redshifting power, so that the decrease of absorption appears as an emission; the coincidence by redshift of a line already written in the spectrum with a Lyman line writes the whole spectral pattern of the gas into the spectrum.

## 4 The “forest”is made of “trees”.

Suppose that a single Lyman pattern is written in the spectrum. The coincidence of the written, redshifted Lyman $`\beta `$ (resp. Lyman $`\gamma `$) line with the Lyman $`\alpha `$ line of the gas writes the Lyman pattern into the gas. Both patterns differ by the shift of frequencies $`\nu _{(\beta \mathrm{resp}.\gamma )}\nu _\alpha `$ of the $`\alpha `$ and $`\beta `$ (resp. $`\gamma `$) lines. As in the standard computations the lines are considered as Lyman $`\alpha `$, the frequency shift is relative to the Lyman $`\alpha `$ frequency:

$$z_{(\beta \mathrm{resp}.\gamma ),\alpha }=\frac{\nu _{(\beta \mathrm{resp}.\gamma )}\nu _\alpha }{\nu _\alpha }\frac{11/(3^2\mathrm{resp}.4^2)(11/2^2)}{11/2^2}$$
(1)
$$z_{(\beta ,\alpha )}5/270.185230.0617;z_{(\gamma ,\alpha )}=1/4=0.25=40.0625.$$
(2)
Similarly $`z_{(\gamma ,\beta )}7/1080.065.`$ The redshifts appear, with a good approximation as the products of $`z_f=0.062`$ by an integer $`q`$.

The intensities of the Lyman lines are decreasing functions of the final principal quantum number $`n`$, so that the inscription of a pattern is better for $`q=3`$ than for $`q=4`$ and a fortiori for $`q=1`$.

Iterating, the coincidences of the shifted lines frequencies with the Lyman $`\beta `$ or $`\alpha `$ frequencies build a tree, final values of $`q`$ being sums of the basic values 4, 3 and 1. Each step being characterised by the value of q, a generation of successive lines is characterised by successive values of $`q:q_1,q_2\mathrm{}`$ As the final redshift is $`q_Fz_f=(q_1+q_2+\mathrm{})z_f`$, the addition $`q_F=q_1+q_2+\mathrm{}`$ is both a symbolic representation of the successive elementary processes, and the result of these processes.

The name “tree”is not very good because branches of the tree may be sticked by coincidences of frequencies. A remarkable coincidence happens for $`q=10`$, this number being obtained by the effective coincidences deduced from:

$$10=3+3+4=3+4+3=4+3+3=3+3+3+1=\mathrm{}$$
(3)
$`q=10`$ is so remarkable that $`z=0.62`$ may seem experimentally a value of $`z`$ more fundamental than $`z_f`$. In these computations, the levels for a value of the principal quantum number $`n`$ larger than 4 are not taken into account, assuming that the corresponding transitions are too weak.

Is the forest made of a single or several trees?

## 5 Conclusion

The present computation is a quantitative explanation of the frequencies observed in the spectra of the quasars. It requires only elementary physics, no unknown matter. The computation should be improved by an evaluation of the intensities.

The CREIL appears as a key in the study of the quasars. It should probably be helpful for other studies: for instance, some astrophysicists think that the dark matter needed to get the gravitational stability could be simply molecular hydrogen; if this is true, it exists some H$`{}_{}{}^{+}{}_{2}{}^{}`$, and a contribution of the CREIL not only to the “intrinsic”redshifts, but to the “cosmological”redshift.