text:
# Human behavior as origin of traffic phases

\[

## Abstract

It is shown that the desire for smooth and comfortable driving is directly responsible for the occurrence of complex spatio-temporal structures (“synchronized traffic”) in highway traffic. This desire goes beyond the avoidance of accidents which so far has been the main focus of microscopic modeling and which is mainly responsible for the other two phases observed empirically, free flow and wide moving jams. These features have been incorporated into a microscopic model based on stochastic cellular automata and the results of computer simulations are compared with empirical data. The simple structure of the model allows for very fast implementations of realistic networks. The level of agreement with the empirical findings opens new perspectives for reliable traffic forecasts.

\]

The empirical observation of highway traffic has shown the existence of very complex spatio-temporal structures . Therefore it was questioned if it is possible to find a model which can reproduce the various observed phenomena and which is able to relate the origin of these phenomena to the human driving behavior. Compared to the analysis of a conventional liquid the difficulties in modeling traffic are already in describing the “microscopic” interactions between the vehicles. These do not follow the fundamental laws of nature but are rather complex decisions taken individually by each driver. Therefore any successful microscopic description has to extract the essentials of the human behavior in different traffic situations, which is clearly a difficult task. Most present modeling approaches concentrate on the fact that drivers want to avoid accidents. This has been implemented by adjusting the velocity according to differences in speed and/or the distances to the other cars. Traffic models based on this interaction only have been able to reproduce various observed phenomena but fail to give a complete description of the empirical results. Much progress has been made in recent years in understanding topological effects in highway traffic . It became clear that e.g. a change of the activity of an on-ramp may induce a phase transition on the highway. The observed phase transitions are between stable traffic states, i.e. the generated traffic states dissolve only if the external conditions change again. The robustness of the empirical observations give strong evidence that traffic states itself are a consequence of the human driving behavior rather than a response to different topological situations.

Below we will show that it is possible to overcome the problems in modeling traffic if one takes into account that people like to have a comfortable journey, i.e., they try to avoid strong accelerations and abrupt braking. This second essential demand of drivers has been incorporated in a recent traffic model, which introduces “brake-lights” for a timely adjustment of the velocity when approaching slow upstream traffic and “anticipation” by estimating the velocity of the leading vehicle. This approach goes far beyond the consideration of velocity differences since acceleration changes become visible and allows for an event driven anticipation of velocity reductions. These two features lead to a stabilization of the flow in dense traffic which is crucial to overcome the difficulties in describing the empirically observed phases and their transitions. Like any other theoretical approach traffic models have to be compared with empirical results. Here one has the difficulty that one cannot perform experiments but is restricted to pure observation of a given situation. Nevertheless much progress has been made in recent years by analyzing large data sets in different environments. So it is now widely believed that at least three traffic states exist, i.e., (i) free-flow (ii) synchronized traffic and (iii) wide moving jams . The characteristics of free-flow  and wide moving jams  are intuitively clear. For a long time it was believed that these states are the only stable traffic states. This commonly accepted picture was enhanced by establishing a second stable congested state, i.e., synchronized traffic. Synchronized traffic , which is typically observed at on- and off-ramps, is characterized by a large variance in flow and density measurements and a velocity which is significantly lower than in free-flow traffic. The origin of the notation “synchronized traffic” is the fact that the time-series of measurements on different lanes are highly correlated. But more important is the apparent absence of a functional flow-density form, i.e., the measurements of the flow, density and velocity of the traffic are distributed over a wide area . This observation has been confirmed quantitatively by means of vanishing cross-correlations between these two quantities. Synchronized traffic and wide moving jams differ also in their behavior at bottlenecks. If synchronized flow is generated at a bottleneck its downstream front is pinned at the bottleneck. In contrast the downstream front of wide moving jams propagates with constant velocity in upstream direction. Their velocity does not depend on the traffic state they cross and is even unchanged if they pass a bottleneck. The final characteristic feature are the phase transitions between the different states . In general they are of first order , e.g., the phase transition free-flow to synchronized traffic is characterized by a discontinuous change of the velocity.

Already the correct reproduction of synchronized traffic on a quantitative level is a difficult task, but the most puzzling point for any model is to reproduce stable synchronized traffic in coexistence with the other two traffic states. The discussed model is able to pass this most sensitive test.

The model  is based on the cellular automaton model of Nagel and Schreckenberg  and incorporates the desire for smooth and comfortable driving by anticipating the velocity of the leading vehicle and the introduction of ”brake lights”. For comparison with the empirical data we simulate a two-lane segment with an on-ramp (fig. 1). The additional input of cars triggers a dense traffic region behind the on-ramp which can be identified as synchronized traffic. Our simulations clarify how the synchronized state is related to the human factors in driving.

(i) Velocity anticipation: At the on-ramp the anticipation of the leaders velocity avoids abrupt braking of the traffic behind and therefore reduces the probability to form jams.

(ii) Retarded acceleration: Comfortable driving also implies that cars do not accelerate immediately in case of a larger gap ahead if they observe slow downstream traffic. On the one hand, this leads in some sense to a sub optimal gap usage, because the velocity is smaller than the headway allows. On the other hand, larger gaps in a dense region reduce the car-car interactions what may cut a chain of braking overreactions which is responsible for the formation of jams. These overreactions are a direct consequence of the delayed human behavior in adapting the velocity to the headway which can lead to an avalanche like amplification of the velocity fluctuations upstream and finally to the formation of jams.

(iii) Timely braking: Finally timely braking suppresses another mechanism of jam formation: When the velocity adjustment is only based on the distance to the next car ahead, jams often emerge in the layer between free-flow and synchronized traffic. In these models the jam formation arises from cars approaching a slow-moving cluster with high speed which leads to a compactified region. In contrast, our approach avoids this artificial mechanism to form a jam, the drivers adjust their speed to the vehicles ahead.

At this point we stress the fact that the observed phenomena are a consequence of the individual behavior of the drivers (See  for the calibration of the model). None of the microscopic model parameters has been changed throughout the simulations in order to optimize the accordance with the empirics. Already the single lane model on a periodic street without bottlenecks  shows the existence of synchronized traffic and wide moving jams. In this simulation the boundary conditions are used only to induce the transitions to these traffic states. By contrast, the excellent agreement with the empirics was obtained simply by applying the correct inflow at the upstream end and the on-ramp of the highway section. This side-steps another important question in traffic dynamics, i.e., the origin of phase transitions. Our simulation results support the view that the transitions are mostly induced by obstacles rather than by a spontaneous breakdown of the traffic stream.

The simulation protocol emulates a few hours of highway traffic including the realistic variations of the number of cars which are fed into the system. A busy on-ramp in combination with high flow generates in general synchronized flow on the highway segment. For the sake of simplicity we used only one type of cars in the simulations which leads to a smaller variance of the data points in the free flow regime compared with the empirical data.

The simulations show that we can recover the empirical results for the fundamental diagram quantitatively (see Fig.2). The agreement is not only for the average values but also concern the statistical properties of the results. This is mandatory for traffic forecasting, e.g., in order to calculate upper limits of individual travel times as well. But, as mentioned above, also the stability of the synchronized traffic state is described correctly. In order to verify this a jam was generated by an obstacle at the downstream end of the highway section. Figure 3 illustrates how the jam wave travels through the free flow region with constant velocity and also passes the section were the synchronized traffic is localized. This shows that we can superpose the different traffic states as observed empirically. The fact that the jam propagates with a constant velocity and passes the free-flow, on-ramp and synchronized regions without being disturbed is typical for wide moving jams.

Our simulation results have far reaching consequences for the theory of traffic flow as well as for practical applications. From a theoretical point of view, it is shown that the desire for smooth and comfortable driving is the origin of synchronized traffic and wide moving jams and is responsible for the stability of the different traffic phases. This stability allows for the application of phenomenological approaches. In particular the motion and formation of jam waves, which is most interesting for any traffic forecast, should be predictable within these approaches with high accuracy (see  for approaches of this kind).

From a practical point of view, our simulation results open the door the desire for smooth and comfortable driving allows for more realistic simulations and opens the door for a forecast of highway traffic which should outperform knowledge based approaches. Models of this type have already been used for fast microsimulations of large traffic networks. Here the simple structure of the model, i.e., its discreteness and the local rules, is of great importance.

Summarizing, we have shown that a rather simple cellular automaton model is able to reproduce the empirically observed phases of traffic flow even quantitatively. The model parameters can be related directly to the human behaviour, especially to the desire for smooth and comfortable driving. This desire turns out to be responsible for the occurance of complex spatio-temporal structures, like the so-called synchronized traffic phase.

Acknowledgments: We thank D. Helbing and B.S. Kerner for useful discussions. The authors are grateful to the Landesbetrieb Straßenbau NRW for the data support and to the Ministry of Economics and Midsize Businesses, Technology and Transport for the financial support. L. S. acknowledges support from the Deutsche Forschungsgemeinschaft under Grant No. SA864/1-2.


format: text

text:
# Statistical and Dynamical Study of Disease Propagation in a Small World Network

## 1 Introduction

To model disease propagation, it is necessary to define the corresponding social network connecting any two individuals in the world. The expected properties of such a network should be both the clustering (which excludes models of disorder like the random graphs ), and to allow a connection between any two individuals within a finite number of steps (which excludes the regular networks with only nearest neighbor connections). Indeed, for the latter feature Milgram showed in $`1967`$ that the average number of steps connecting any two individuals is six (called also six degrees of separation) . This behavior led recently Watts and Strogatz to propose the model of Small World Network ($`SWN`$) . They considered a low dimensional network with periodic boundary conditions for convenience (a ring for example) where they rewired some bonds with a probability $`\varphi `$ to a new site randomly chosen from the network. For small values of $`\varphi `$ this still corresponds to a regular network but with few long range connections called Short-Cuts ($`SC`$). A more recent work on the $`SWN`$ was proposed by Newman and Watts  where the number $`k`$ of Nearest Neighbors ($`NN`$) is conserved but instead of rewiring, they added an average density $`\varphi `$ of new bonds from each site $`i`$ to other nodes randomly chosen (except its nearest neighbors). A review of these models and their application to various fields and particularly epidemics can be found in the references . In these networks the percolation threshold was extensively investigated and its dependence to the $`NN`$ and $`SC`$ was found to satisfy the following equation

$$\varphi =\frac{(1p_c)^{k/2}}{p_c}$$
(1)
This threshold corresponds in epidemics to the smallest concentration of susceptibles leading to the outbreak . However, the statistical behavior of percolative $`SWN`$ networks is different from that of regular systems . In particular, at the percolation threshold, there is no diverging cluster for the $`SWN`$ because a $`SC`$ between the two ends of the system has a finite probability to occur for any non-vanishing number of these bonds. On the other hand, the characteristic length scale in such networks (which corresponds to the correlation length in regular lattices) behaves as $`\varphi ^{1/d}`$, $`d`$ being the euclidean dimension of the system . It is then obvious that this characteristic length does not diverge at the percolation threshold for such networks. Therefore, a further investigation of the cluster statistics around this new percolation threshold for such networks and the phase transition seems to be necessary.

Let us now consider the application of this model to epidemics which seems to be one of its main aims. From the large amount of works using $`SWN`$, there is no direct comparison with the existing data may be due to the complexity of the diseases features (incubation and latent periods, birth and death rates etc.). Furthermore, the parameters used (mainly $`\varphi `$) are very small and do not simulate the real connections between individuals. They use also commonly average values of the $`NN`$ and the $`SC`$ while these quantities strongly fluctuate in the real live (the number of contacts, friends, family members etc. varies from $`0`$ to few tens) which can influence sensitively the results on the density of susceptibles at the percolation threshold (epidemic outbreak). On the other hand, it is impossible in practice to measure the density of susceptibles systematically (it needs an extensive serological investigation in the epidemic phase). Generally, for large population samples epidemiologists measure the evolution with time of the number of cases for a given disease. It is then necessary to study the dynamical behavior of the propagation of the disease and relate it to the density of susceptibles. There are only few works which examined (only qualitatively) the dynamical behavior of the disease on social networks . The aim of these works was to show how the density of infected behaves in the endemic and epidemic phases.

In this article, we use a site percolation on a $`SWN`$ with parameters ($`k`$ and $`\varphi `$) representing a sample of school age children to study the effect of the fluctuations of $`NN`$ and $`SC`$ on the percolation threshold. Furthermore, in order to propose a formula for determining the serology of the sample from the rate of increase of the number of cases, we investigate also extensively the dynamical behavior of an infectious disease as well as its effect on the density of susceptibles below and above the percolation threshold. A new super-diffusion is found above the percolation threshold when the cluster is initially infected by one or a small number of infectious sites and its characteristic time dependence on the density of susceptibles is determined. We examined also the case where the infection probability of the $`NN`$ is different from that of $`SC`$, showing a double diffusion with two characteristic times. In the next section we describe the model and then present the results on the statistics of the clusters and the percolation threshold. The results on the dynamical behavior of the disease are presented in section 4.

## 2 Model description

We consider the one dimensional $`SWN`$ described by Newman and Watts  but $`\varphi `$ represents the total number of $`SC`$ generated for each site uniformly from all the other sites of the network. In the case where $`k`$ and $`\varphi `$ are not fixed they are generated randomly within a normal distribution centered at their average values with fluctuations $`\delta k`$ and $`\delta \varphi `$ respectively. The coordination number is the total number of bonds to a given site ($`z=k+\varphi `$). We study in this network a site percolation problem  by assuming each susceptible site $`j`$ (occupied) contracts the disease if it is connected with an ill site $`i`$ (occupied also). The occupied sites (susceptibles) are randomly generated with a concentration $`p`$ while the empty sites correspond to refracted individuals. For $`k`$ and $`\varphi `$ fixed, the percolation threshold $`p_c`$ is related to $`k`$ and $`\varphi `$ by Eq.(1) . This threshold corresponds to a transition from the endemic phase below $`p_c`$ to the epidemic one above this point . In $`SWN`$ networks, $`p_c`$ is the minimum concentration of occupied sites above which the average largest cluster size $`\xi `$ of the occupied sites becomes power-law increasing with the concentration ($`\xi `$ $`=(pp_c)^x`$ ), while it diverges in a regular network  (note here that the exponent $`x`$ is positive). By analogy with the regular lattices , we will check the universality of the exponent $`x`$.

We are interested in the application of such a model to a childhood disease like measles. In such diseases epidemiological investigations on school age children can be easily controlled and provide data with a minimum bias. We choose parameter values corresponding to such a disease by taking $`k=2`$ to be the average number of brothers, sisters and neighbors, while $`\varphi =30`$ represents the average number of children one can meet at the school. These parameters should correspond to a topology closer to that encountered in a real social network. Regarding the dynamical study, we assume the major contribution to the epidemics provided by the largest cluster. We restrict ourselves then to this cluster and start the infection with one or few infectious sites at time $`0`$. These sites will infect all the connected sites in the next step (after a time $`\mathrm{\Delta }t`$), which themselves infect their connected sites after $`2\mathrm{\Delta }t`$ and so on. We assume the latent and incubation periods smaller than $`\mathrm{\Delta }t`$ which is taken in the rest of this paper as a unit time. The number of infected sites in each step is averaged by varying the initial infectious site position through the whole cluster.

## 3 Percolation threshold and cluster distribution

In this section, we realize $`100`$ configurations of the network described in the previous section with a size fixed at $`100000`$ sites. We examine the effects of $`\varphi `$ and its fluctuations on the average cluster sizes, $`p_c`$ and $`x`$. Finally, we investigate the cluster size distribution around $`p_c`$ in order to determine the main contribution to the propagation of the disease.

In figure 1a we show the variation of the cluster size with the concentration of occupied sites for three different cases: $`\varphi =6,30`$ (fixed values) and for $`k`$ and $`\varphi `$ randomly generated with a normal distribution centered at $`2`$ and $`30`$ respectively, with a fluctuation of $`2`$ et $`15`$ respectively. We see clearly from this figure that in all cases the cluster sizes vary as a power law of ($`pp_c`$ ) above $`p_c`$. For fixed $`k`$ and $`\varphi `$ the value of $`p_c`$ is in a good agreement with the analytical predictions of Newman et al.  (Eq. 1). However, in the case of fluctuations of $`k`$ and $`\varphi `$, this threshold decreases sensitively (about $`50\%`$ in the case shown in Fig.1a) as the fluctuations increase. Therefore, the average values of $`k`$ and $`\varphi `$ are not sufficient to characterize an epidemic outbreak. The $`SC`$ fluctuations $`\delta \varphi `$ decrease $`p_c`$ as a power law with an exponent $`0.1`$ (as shown in figure 1b), indicating a sensitive participation of the larger values of $`\varphi `$ to built the largest cluster. Therefore, the percolation threshold behaves as

$$p_c\varphi ^1\delta \varphi ^{0.1}$$
(2)
From this behavior, we can estimate the percolation threshold in a real sample of school age children to be in the range $`2.3\%`$ to $`2.8\%`$.

Now let us restrict ourselves to the case of fixed $`k=2`$ and $`\varphi =6`$ in order to examine the statistical behavior of the clusters around $`p_c`$ (without loss of generality, these values are chosen only because $`p_c`$ is large enough to enable sufficient cluster statistics for such a sample size) . We found that the cluster size fluctuations are maximum at this threshold (see Fig.2a) implying a divergence of this quantity at $`p_c`$ which seems to be a signature of a phase transition. The cluster size distribution (see figure 2b) confirms this divergence since it decreases exponentially below $`p_c`$ while it is power-law decreasing at this threshold (this power law behavior is in agreement with the results of Castellano et al.  on other systems). Indeed, at $`p_c`$ this corresponds to a Lévy distribution  with an exponent of $`2.13`$ indicating the divergence of all its moments. We notice here that only the higher sizes (rare events) contribute to the outbreak at $`p_c`$ (as expected in such distributions). Above $`p_c`$ the small size clusters are absorbed by the largest one and we have again an exponentially decreasing distribution for small clusters while there is only one very large cluster (not shown in Fig.2b).

Since the cluster size does not diverge at $`p_c`$, it is obvious that $`x`$ is not universal (because it is not a critical exponent), but it is interesting to know how it depends on $`\varphi `$ in such lattices. In figure 2c, the exponent $`x`$ seems to vary only linearly for larger values of $`\varphi `$ but with a very small slope (about $`5.6`$ $`10^3`$). It is difficult to predict its behavior for very small values of $`\varphi `$ because in this case the network tends to a regular one and the cluster size becomes very large so that the sample sizes taken here do not allow us to measure this exponent accurately.

However, even if the parameters chosen in this model are close to those of a real social network, it seems impossible for epidemiologists to check these results. Indeed, as explained below, they cannot measure the density of susceptibles, except if they investigate systematically the serology of a sufficiently large sample of school age children (e.,e.g. for a city sample). Therefore, the behavior of $`p_c`$ should be checked for measurable quantities. In the case of disease propagation, the time dependence of the number of cases can be directly measured by epidemiological techniques. We will investigate this dynamical behavior in the following section.

## 4 Dynamical study of the propagation of a disease

In this section we restrict ourselves to the fixed values of $`k`$ and $`\varphi `$ ($`2`$ and $`30`$ respectively) to simulate a sample of school age children. From the results of Fig.2b, we assume that the main growing effect of the infection comes from the largest cluster and estimate the propagation time of the epidemics only from this cluster. We determine the evolution with time of the number of cases for both phases endemic ( $`p<p_c`$ ) and epidemic ( $`pp_c`$ ). As found in figure 2a the cluster size at $`p_c`$ strongly fluctuates and therefore, the time behavior of the number of cases also fluctuates. The variation of the number of cases with time is shown in figure 3a for three cases (just below $`p_c`$ , at $`p_c`$ and above $`p_c`$) with only one initial infectious site. In both cases, the number of cases increases up to a maximum and then decrease because the number of susceptibles decreases. In the endemic phase, the number of connections between occupied sites in the cluster is mostly $`1`$ and does not allow a significant increase of the number of cases (the behavior in this case is under estimated since all the clusters should contribute to this increase). For susceptible densities around $`p_c`$ this situation persists for a long time and the number of cases increases linearly with time showing a normal diffusion of the disease. In the epidemic phase the increase becomes exponential indicating a new kind of super-diffusion  of the disease, due mainly to the increasing number of connections in the cluster (as shown in figure 3b). This exponential growth is also observed for $`SIR`$ models  where the rate is proportional to the basic reproduction rate $`R_0`$ which correponds in our case to the average number of connections in the cluster. We have also performed a Monte-Carlo simulation to the measles propagation in a more realistic sample (births, deaths, incubation and latent periods etc.) where the average infections is $`2`$ for each infectious individual and found also an exponential growth of the infected cases . Therefore, this exponential growth does not seem to depend on the topology of the sample but the rate is sensitive to the geometry of the network. Note in the present work that in the case of more than one initial infectious site (see figure 3c) the exponential growth behavior does not change but the growing rate fluctuates due to the fluctuating number of connections. The average rate of the exponential growth $`\gamma `$ (corresponding to the characteristic time of the epidemics) is shown in figure 4 to increases as $`Ln(p)`$ above $`p_c`$ while the period of this epidemic behavior decreases. From this figure we can conclude that when the characteristic time decreases below $`5`$ (or $`\gamma `$ increases above $`0.2`$), the epidemic behavior takes place. This behavior seems to have a direct application in epidemiology since it provides a method for the estimation of the serological situation (density of susceptibles) from the characteristic time which is easily measurable. Therefore, this result stimulates a proposal for a serological examination for a given childhood disease in a sample of age school children, but during an epidemic period to compare a realistic behavior with that obtained in this paper.

Now let us examine the case of adding different infection probabilities to this system. We consider that a site $`i`$ infects another site $`j`$ with a probability $`p_n`$ if $`j`$ is a neighbor of $`i`$ and $`p_{sc}`$ if it is a short cut. The motivation of this investigation is that a susceptible child has a different probability to be infected by his brothers (or sisters) than by the other children meeting him at the school. We see clearly a double diffusion behavior in figure 5 (for $`p_n=0.1`$ and $`p_{sc}=0.9`$), where the number of infected starts growing exponentially up to the characteristic time ($`1/\gamma `$), then it increases as a power law up to a new characteristic time from which it grows again exponentially with the same rate. The slow diffusion is due to the small contact probability for the neighbors ($`p_n=0.1`$) and has been observed in other fields . This slow diffusion appears very short because the number of $`NN`$ is very small ($`k=2`$). It should be interesting to investigate this double diffusion for larger $`k`$ (which is the case of animal diseases).

## 5 Conclusion

We have investigated in this article the statistics of the cluster sizes in a one dimensional $`SWN`$ by taking into account the $`NN`$ and $`SC`$ fluctuations. We found that these fluctuations decrease $`p_c`$ as a power law with a small exponent leading to a new expression for the percolation threshold. We found also that cluster size fluctuations is the quantity governing the phase transition in such a network. On the other hand, in order to apply our results to the measured quantities in epidemiology, we have studied the dynamics of the disease propagation in such clusters. We found in epidemic phases a super-diffusive with an exponentially growing number of infected sites, while at $`p_c`$ this number increases as a power law. The growing characteristic time is larger than $`5`$ in the endemic situation and decreases linearly with $`Ln(p)`$ in the epidemic phase. This result provides a way to estimate the density of susceptibles in the epidemic phase. We propose then a serological investigation in epidemic situations to check this behavior. Finally, we examined the case where the infection probability is very small in the $`NN`$ compared to the $`SC`$. The dynamical behavior of infected cases shows a double diffusion with two characteristic times, and a power-law increase (deceleration) between them. We think that this effect is useful for samples with large $`NN`$ and shows a way to stop the propagation of the epidemic for other diseases.

ACKNOWLEDGEMENTS

One of the authors (NZ) would like to thank the Arab Fund for Economic and Social Development for the financial support and the staff of the I.U.S.T.I., Université de Provence - Marseille for hospitality during the progress of this work. We thank Professor M.Barthélémy for fruitful discussions and Professor A.M.Dykhne for drawing our attention to the double diffusion behavior.

Figure Captions

Figure 1 a)  Cluster size (number of sites in the cluster) versus concentration of the occupied sites for three cases: $`\varphi =6`$ (solid curve); $`\varphi =30`$ (dotted curve) and $`\varphi =30`$ with fluctuations $`\delta k=2`$, $`\delta \varphi =15`$ (dahs-dotted curve).

b)  $`p_c`$ versus $`\delta \varphi `$ (sites) for $`\varphi =30`$ sites. The solid line is a fit of the data.

Figure 2 a) Cluster size fluctuations (sites) versus the occupied sites concentration ($`\varphi =6`$ sites). The solid curve is a guide for the eyes.

b) Distribution of the cluster size (sites) for fixed ($`k=2`$, $`\varphi =6`$) in a semi-log plot
at $`p_c`$ (solid curve) and $`p=1\%`$ (dotted curve). The dotted line is a linear fit of the data below $`p_c`$. The insert is a log-log plot of the distribution at $`p_c`$ with a linear fit.

c) Variation of the exponent $`x`$ with the number of short cuts $`\varphi `$ (sites). The solid line is a linear fit of the data to $`5.6`$ $`10^3`$.

Figure 3 a) Number of cases versus time for three different cases; $`p=3.5\%`$ (solid curve), $`p=4.5\%`$ (dashed curve) and $`p=8\%`$ (dotted curve). Insert: log-log plot with a power law fit of $`p=4.5\%`$ and an exponential fit of $`p=8\%`$.

b) Distribution of the number of connections (acquaintances) in the largest cluster for $`p=3.5\%`$ (solid curve), $`p=5\%`$ (dashed curve) and $`p=10\%`$ (dotted curve).

c) The rate of the exponential growth (in arbitrary units) versus number of initial infectious sites. The horizontal line is the average rate.

Figure 4 The rate of the exponential growth versus $`p`$. The solid line is a fit of the curve linearly with $`Ln(p)`$.

Figure 5 The Number of infected cases versus time (in arbitrary units) for one initial infectious site and an infection probability one (solid curve), and the probabilities of infection: $`p_n=0.1`$ and $`p_{sc}=0.9`$ (dotted curve). The dashed curve is a power law fit of the second data in the region of the double diffusion.


text:
# The Onset of Mackey-Glass Leukemia at the Edge of Chaos

## 1 Introduction

Many years ago Mackey and Glass  proposed two different simple models to describe chronic and acute diseases, generated by the failure of physiological control systems, which may be identified by the altered periodicity of some observable. Recently Landa and Rosemblum  revisited one of the models proposed to describe the respiration control and introduced a modification in the original equation in order to take into account the activity of the brain on respiration. In this paper we revisit the model proposed to describe the fluctuations in the blood cell counts in normal individuals and the initial phase of chronic myelocytic leukemia (CML).

The human blood contains a remarkable variety of cells, each type precisely designed to its own vital function. All different types of blood cells are produced in the bone marrow and all of them are originated from earlier non-functional cells called stem cells. The stem cells replicate repeatedly and differentiate into various kind of secondary cells. The growth and development of cells in the bone marrow are carefully controlled to produce the correct number of each type of cell necessary to keep the body healthy. The blood-forming process (hematopoiesis) is subjected to active control by mechanisms acting to buffer the system against disturbances and to provide adaptive adjustment to non-stationary conditions . The hematopoiesis is regulated by feedback control systems and as other time-delayed processes  it exhibits oscillatory and cyclic behaviors.

Most of the models proposed to describe physiological control systems consider the differentiation steps of the secondary cells as different compartments described by differential equations. In general those models have a high-dimensional parameter space (\- ) which makes it quite difficult to compare theoretical results to the experimental ones. In contrast to those more realistic, but rather complicated models, Mackey and Glass  proposed a very simple model based on a single time-delayed differential equation, with parameters estimated from experimental data. The change of the time delay parameter leads to different kind of oscillations of the cell counts which may be associated to healthy and non-healthy behaviors when compared to experimental data obtained by Gatti et al . The onset of an abnormal dynamics, which may be associated to the onset of leukemia is then obtained in this model by the gradual tuning of a control (time-delay) parameter. Since its proposal this model has been mostly used as a paradigm for non-linear systems (\- ),and less attention was paid to its biological applications.

In this paper, we show that this simple model that generates very complex dynamical behaviors for cell counts, is more robust than thought to reproduce the onset of leukemia and it is also possible to use it to reproduce the two phases of leukemia and the effects of some treatments. The observed instabilities that occur at the onset of leukemia can be obtained for different values of the time delay parameter, corresponding to different possibilities of the system to achieve the edge of chaos. In other words, the instabilities resulting from complex interactions involved in the physiological control systems that may describe the oscillations observed in CML patients correspond to an state hovering between the incoherence of chaos and the spontaneous order found in periodic behaviors ().

In section 2 we present the biological motivation that inspired the model discussed here and other experimental data concerning the discussions we present along the paper. In section 3 we introduce the model. In section 4 we present the results obtained. We first obtain the bifurcation diagram in function of the time delay parameter for the Mackey-Glass model and show that the onset of leukemia is related to instabilities associated to transition regions between chaotic and periodic regimes, usually refereed in the literature as the edge of chaos. These transition regions occur for a set of different values of the time-delay parameter and not only for a unique value of this parameter. The bifurcation sequence observed in this diagram is quite uncommon and presents only some of the periods suggested by the authors . We also introduce a simple modification in the death rate parameter of this model in order to simulate the progressive increase of the minima counts in phase II of CML. The bifurcation diagram in function of the death rate parameter is also obtained and the effects of leukapheresis treatments is also discussed. Our conclusions are summarized in section 5.

## 2 The Biological Motivation

It has been observed spontaneous fluctuations with periodic behavior in the peripheral blood cell counts of healthy individuals, in chronic myelocytic leukemia (CML) patients  as well as in other pathological states . In healthy individuals constant or mild oscillations of the peripheral-blood leukocyte-counts occur with periods varying from 14 to 24 days , while in many patients with CML these oscillations have periods varying from 30 to 120 days (), with a cell accumulation process observed in the final stage (phase II) of the disease. Morley et al  suggested that in patients with cyclic leukocytosis and CML, the bone marrow is either unduly sensitive to feedback stimulus or being excited intermittently by an abnormal stimulus, generated outside the marrow.

From long term observations Gatti et al.  showed that the onset of leukemia (CML) is characterized by two distinct phases: phase I, a stable 15-month period, from time of diagnosis, during which leukocyte counts peaked every 72 days and then returned to the baselines level; and phase II, a 18-month period during which the minima levels of leukocyte counts became progressively higher indicating that those cells were accumulating in the peripheral blood. According to them, the increased marrow production in the disease does not represent a life-threatening situation, until cell accumulation begins. During phase I, the leukocyte counts return to baseline levels, indicating that the overproduction of cells in the peak of a given cycle could be cleared from the peripheral blood before the beginning of the next cycle. In contrast, in phase II, the leukocytes accumulate in the peripheral blood of the patient, suggesting primarily an increase in the cell production or a progressive loss of cell-clearing capacity with little or no changes in the production capacity .

Patients with CML were submitted to leukapheresis treatments in order to test out if this kind of treatment would be able to reverse the rising leukocyte levels providing a good therapy to maintain the patient healthy for years. The leukapheresis consists of 8 weeks of intense treatment to remove about $`3050\%`$ of the mature and premature leukocytes from marrow and peripheral blood. According to the results obtained , those treatments are able to retrieve the system from the behavior of phase II to the behavior observed in phase I, with low amplitude cycles and minima counts. However, after one year, subsequent cycles (of 72 days each) indicated re-accumulation of cells in the peripheral blood. In other words, these studies clearly demonstrated that leukapheresis is effective for the rapid reduction of elevated cell counts, but long-term treatment has not improved survival of the patients with CML, since this kind of treatment do not alter the mechanisms underlying the disease process.

Nowadays this kind of treatment is not used anymore, and have been replaced by other treatments like chemotherapy and radiotherapy that in some cases may lead to the patient cure. The mechanisms underlying the effectiveness (or not) of those modern treatments are not clear-they act by producing an intense immunosupression followed by an intense activation of the bone marrow in order to recover the system, involving different kinds of drugs- and it would be impossible to discuss its effects on this very simple model. However, it is possible to analyze and simulate the effects of leukapheresis in the MG model as will be shown in section 4.

## 3 The Model

The Mackey-Glass model  relates the number of precursor and mature cells resident on the bone marrow to the number of cells released into the blood stream. The time evolution of an homogeneous population of mature, circulating cells of density $`P`$ is then given by:

$$\frac{dP(t)}{dt}=\frac{\beta _0\theta ^nP_\tau }{\theta ^n+P_\tau ^n}\gamma P(t)$$
(1)
where the first term in the right hand side corresponds to the cellular production, which is described by a (Hill) function of $`P_\tau =P(t\tau )`$, where $`\tau `$ represents the delay between the production of the precursor cells and the release of the mature cells into the blood. The average (density independent) death rate of the population is given by $`\gamma `$. The constant values of $`\beta _0`$, $`\theta `$, $`n`$ and $`\gamma `$ are estimated from experimental data .

Like other theoretical models (-), the Mackey-Glass model associates the onset of cyclic leukocytosis and leukemia with instabilities in the dynamics of the model. They found that for $`\tau =6`$ days the solution has low amplitude oscillations with period of 20 days, which describes qualitatively the mild oscillations observed in the circulating levels of leukocytes in normal healthy adults . For $`\tau =20`$, the numerical solution reproduces the altered periodicity observed in phase I of the diagram of circulating leukocytes as a function of time of patients with diagnosed CML .

## 4  Results

By standard fourth-order Runge-Kutta integration techniques the numerical solution for the equation is obtained, where we adopt the same set of parameters of the original model . In order to characterize the dynamical instabilities associated to the onset of leukemia we have obtained the bifurcation diagram of the local maxima of the time evolution of $`P(t)`$ as a function of the time delay parameter $`\tau `$ as shown in figure 1. This diagram was built from the analysis of the cross section of the return map for different values of $`\tau `$. By increasing the time-delay parameter the oscillations become unstable and new cycles with different periods appear, leading to aperiodic and then to chaotic behavior. In contrast with the findings of Mackey-Glass, we found a quite uncommon bifurcation sequence: $`1`$, $`2`$, $`3`$, $`6`$, $`7`$, $`14`$, etc. with windows of period $`6`$ and $`9`$ in the midst of the aperiodic regime. We also note the appearance of a branch of sparse points just above $`\tau =17`$. This branch persists when we vary the initial conditions.

We have observed that the very same altered periodicity observed for $`\tau =20`$ (which was associated to the onset of the cyclic leukocytosis by Mackey and Glass ) could also be reproduced for other values of $`\tau `$, for instance, $`\tau =18.6`$ and $`21.2`$. The common feature among these values of the time delay parameter is that they correspond to the edge of the periodic windows which appear in the midst of aperiodic regime. This indicates the robustness of the model to describe the disturbances that would correspond to the onset of leukemia.

In this model, the regular behavior of the blood-forming process and the pathological states of disease correspond to different attractors of the dynamics. The standard methods of characterization of the attractors can be applied for small values of $`\tau `$ (periodic behavior), whose attractors have low embedding dimension, and for large values of this parameter (chaotic behavior), whose embedding dimension of the attractors can be approximated to infinity, for which simple asymptotic solutions exist . Nevertheless, for intermediate values of this parameter the attractors are embedded in a high-dimensional space and the characterization becomes quite difficult, specially in cases like this one, which involves time-delayed information. As far as we know up to now the available techniques do not allow a precise characterization of the dynamics for intermediate values of $`\tau `$. In figure 2 we show the power spectra for four values of $`\tau `$ corresponding to the three different dynamical behaviors described above. Our results indicate that the special values corresponding to the onset of leukemia would correspond to a quasi-periodic behavior: the power spectrum has numerous peaks and presents a slow decay to zero. This behavior should correspond to the zero value of the Lyapunov exponent. Nevertheless, because of the difficulties mentioned above, the calculations performed are still not conclusive.

As mentioned before, the phase II of the common pattern observed in CML patients may be associated to either a failure of the cell-clearing capacity of the system or a over-production of immature cells (the cells lose the ability to mature and differentiate) which leads to the accumulation of cells and to a progressive increase of the minima counts in the leukocyte population . In order to simulate the phase II, we introduce a simple modification to the death rate that would correspond to a failure in the cell-clearing capacity. Modifications in the cell production, in order to simulate the loss of the cell’s ability to mature and differentiate is not straightforward and will be discussed elsewhere. The different dynamical behaviors observed in the two distinct phases, will be related to changes in the death rate, since in phase II it will be no longer constant. The death rate $`\gamma `$ will be described during phase II by a step function that diminishes by $`10^4`$ every $`60`$ days. This simple modification simulates a system that slowly loses the ability to remove cells, and reproduces qualitatively the increase of the minima counts, observed by Gatti et al  in phase II of the disease, as shown in figure 3. We have also tried other possible functions, but this one was more appropriate to describe the accumulation of cells observed in phase II.

Figure 3 is obtained for the same set of parameters used in the original model. We have arbitrarily chosen to let the system evolve with $`\gamma `$ constant until $`t=4900`$ (which would correspond to the end of phase I) and after this point the death rate begins to decrease. The slope of the linear growth of the minima counts in our simulations has the same order of magnitude as the equivalent slope in the experimental data.

In figure 4 we present the bifurcation diagram as a function of $`\gamma `$, for $`\tau =20`$, in the range of interest. The diagram shows a bifurcation sequence which ends up in a limited region from $`0.08`$ to $`0.12`$ of aperiodic behavior, above which the system becomes periodic again. In the midst of the aperiodic region we also find windows of periodicity $`6`$.

We also simulate the effects of leukapheresis treatments during the phase II. We let the system evolve during some time ($`600`$ days) in phase II and then we simulate the first leukapheresis, in order to test if this cleaning procedure would revert the raising of the minima counts. The treatment is simulated maintaining the density of cells constant during 60 days (equivalent to the 8 weeks of real treatment) in a level that corresponds to $`60\%`$ of the value of the density in the beginning of the protocol . For $`\tau =20`$, after 60 days the system does not have memory of what happened before the treatment. After each treatment we let the system evolve under the same conditions, since the treatment only removes the excess of cells accumulated and do not interferes on the underlying mechanisms that could generate the disturbance. Just after the treatment the minima counts returned to the baseline levels, but after a while we have observed the same linear increase of the minimal counts of phase II, as in the experimental data. The numerical solutions obtained for two leukapheresis treatments with approximately one year interval ($`400`$ days) between each other is shown in figure 5. The results show a good qualitative agreement with the experimental results . Before the first treatment and between the two treatments $`\gamma `$ is practically constant (varies approximately from $`0.1`$ to $`0.099`$), maintaining itself in the aperiodic region of the $`\gamma `$-bifurcation diagram , with no change of attractor (see figure 4).

## 5 Conclusions

The Mackey-Glass model for the blood-forming process reproduces the main features observed in the peripheral blood-cell counts for normal individuals and pathologies like cyclic leukocytosis and the disorder that characterizes the accumulation of cells observed in a patient diagnosed with leukemia (CML) in phase I and II. We have shown that the onset of leukemia is related to special intermediate values of the time-delay parameter, which seem to correspond to quasi-periodic behaviors. From the bifurcation diagram as a function of the time-delay parameter $`\tau `$, we observed that the onset of leukemia is associated to the edge of periodic windows in the midst of aperiodic behavior. By varying the death rate parameter we could reproduce the phase II behavior observed in the experimental data, and with a simple cell-removing procedure we simulated the effects of leukapheresis. We have also obtained the bifurcation diagram as a function of the death rate $`\gamma `$.

## 6 Acknowledgements

We thank A.P. Serbeto, J. Stilck and A. Bernardes for the critical reading of this manuscript. This work was partially supported by the Brazilian Agencies Conselho Nacional para o Desenvolvimento Científico e Tecnológico - CNPq, CAPES and FINEP.


text:
# An all-cryogenic THz transmission spectrometer

## I Introduction

Many quantum systems of contemporary interest in condensed matter physics have energy levels in the meV range. These include, but are not limited to, two dimensional electron gas (2DEG) systems in high magnetic fields exhibiting the fractional and integer quantum hall effects, quantum nano-structures such as quantum dots and carbon nanotubes, and metallic single-electron transistors. In these systems, temperature dependent behavior in dc transport measurements is observed all the way down to tens of mK. The measurement of the far-infrared (THz) transmission through such systems at ultra-low temperatures could provide complementary information about the excitation spectrum of the system.footnotetext: a)Current address: Integrated Nanosystems Research Facility, University of California, Irvine 92697

The standard technique to measure THz transmission is with a Fourier Transform Infrared Spectrometer (FTIR) or a molecular gas laser which couples optically to a cold sample from room temperature. The disadvantage of this technique is that broadband, thermal blackbody radiation is also coupled to the sample; this is incompatible with ultra-low temperature (sub 100 mK) operation. Numerically, the (integrated) power density is 40 mW/cm2 at 300 K; typically this must be attenuated to of order pW in order to avoid sample heating in many ultra-low temperature systems; at the same time the radiation in the frequency band of interest must be detectable with reasonable signal to noise for the FTIR system to operate at all. In principle, this technique can be carried out. One must construct a set of filters with the desired pass functions from microwave to optical frequencies.

In this paper an all-cryogenic spectrometer is presented that bypasses the necessity to construct such filters. The spectrometer is based on magnetic field tunable cyclotron resonance (C.R.) emission from a high-mobility 2DEG in a GaAs/AlGaAs heterojunction, together with a broadband detector. Figure 1 illustrates the concept. Our demonstration experiment is carried out in a 4 K environment, but the technique is compatible with ultra-low temperature (dilution refrigerator) environments.

A transmission spectrometer using C.R. emission from bulk GaAs was proposed in 1980. Based on that proposal, a spectrometer was built by Knap in 1992 using C.R. emission from 2DEGs in GaAs, with a spectral resolution of 1.3 cm-1 using narrow and moderate band detectors. The spectrometer presented here is an improvement on that developed by Knap it two important ways. First, an ultra-broadband detector is used in place of a narrow band detector. This allows a much wider range of frequencies to be accessed in principle; in this particular case the range of 2-65 cm-1 are immediately accessible with a single sweep of the magnetic field. Second, we use a very low chopping frequency (1 Hz). Previous studies of C.R. emission typically used microsecond or millisecond pulses with low duty cycle for fear of overheating the sample and inadvertently broadening the emission linewidth. The experimental results presented here indicate that low chopping frequencies do not significantly degrade the emission linewidth. Low chopping frequencies will allow the use of state-of-the-art ultra-low-noise bolometers, which have exquisite sensitivity but time constants limited to typically tens of ms. In this paper, we also provide a quantitative estimate of the generated C.R. power and spectrometer noise performance, as well as quantitative estimates for the ultimate limits of this technique using current as well as future detector technology with the promise of single photon sensitivity in the THz band.

## II THz source

### A Principle of operation

In two-dimensional electron gas in the presence of a magnetic field perpendicular to the plane, the spectrum of states (in the absence of disorder and interactions) is given by a discrete set of Landau levels, with energies given by

$$E=\left(n+\frac{1}{2}\right)\mathrm{}\omega _c,$$
(1)
where

$$\mathrm{}\omega _c=eB/m^{}.$$
(2)
$`\omega _c`$ is referred to as the cyclotron frequency. At low temperatures, the lowest energy states are filled, and higher energy states are empty. If the electrons are heated to a higher temperature, then higher energy states are populated. The higher-energy electrons can decay to the lower energy levels via phonon emission or photon emission; in the latter case the photon frequency is given by $`hf=eB/m^{}`$. In our technique, we heat the electrons with a dc current; the emission frequency is then tunable by magnetic field. This cyclotron emission from 2DEGs has been studied by many groups. Here, we are mainly interested in characterizing the emission power and spectrum in order to use it as a THz source.

Our source consists of a 5x5 mm2 (cleaved) 2DEG formed in a GaAs/AlGaAs modulation-doped quantum well grown by molecular beam epitaxy. The mobility and density are $`1.2510^{11}cm^2`$ and $`600,000cm^2/Vs`$, respectively. Ohmic contacts of diffused Au/Ni/Ge are deposited along both edges of the sample to make low-resistance contacts. We have characterized the emission power using a broadband detector. To drive our source electrically, we apply a 0.5 Hz sinusoidal voltage to the terminals, typically with a 10 k$`\mathrm{\Omega }`$ resistor in series with the emitter. Since heating by the electrical current causes spontaneous emission, and since the power delivered to the emitter from the battery is proportional to $`V^2`$, the emission power is periodic with a frequency given by twice the electrical drive frequency. We thus detect the synchronous voltage on the detector at the second harmonic of the drive frequency, 1 Hz.

The detector is a commercially available broadband composite bolometer placed at one end of a 7 ” long 1/2 ” diameter evacuated (gold plated) light pipe. It consists Ge thermistor mounted on a thin-film (NiCr) coated sapphire absorber with dimensions 4x4 mm2. The emitter is placed at the other end of the light pipe, inside the bore of a superconducting magnet. The measured emission power at the detector for an electrical input power of 1 mW (electric field of roughly 20 V/cm) is shown in figure 2 as a function of the source magnetic field. The smooth variation of power as a function of magnetic field is due to the variation of the emitter resistance (and hence the power delivered from the battery) with the magnetic field. The periodic structures are interference fringes due to the Fabry-Perot etalon formed by the substrate. This effect was modeled and measured by Zinovev, and provides evidence that the source is indeed quasi-monochromatic. Three important figures of merit for our source are the optical beam pattern, emitted power, and linewidth. We discuss these each in turn below.

### B Optics

The emission from our sample is mostly into the substrate due to its high dielectric constant. We have mounted the sample with a thin layer of vacuum grease onto a highly-reflective (gold plated) mount. This reflects most of the light emitted into the substrate. We also tried evaporating a thin metallization layer on the back surface (Al), with no change in the output power. It is well-known that the emissivity from thin films is distributed very uniformly over the entire 2 $`\pi `$ steridian available solid angle. In our case, since the beam pattern is spread out even more after reflection and passing through the dielectric/air interface, the emission is essentially isotropically distributed as a function of solid angle. We use a light-pipe to guide the radiation to the detector; in principle a lens system can be used to focus it if that is desired. The radiation is circularly polarized for vertical emission, but multiple reflections off of the light pipe walls and the substrate/sample interface at varying angles served to randomize the polarization.

### C Power

The absolute power available is important, since this sets the sensitivity requirements for the detector in our arrangement. Absolute power calibrations at this frequency are typically no better than a factor of two, and our results for the measured signal are also reproducible to only a factor of two between cooldowns. This may be due to drift in the detector responsivity or to inherent changes in the sample between thermal cycles. We estimate the absolute power response of the detector using the method of dc substitution. In this method, we determine the detector responsivity (in V/W) to dc power, and assume the responsitivity is the same for THz power. This is generally the best agreed upon method for absolute power measurements at THz frequencies. This gives a responsivity of $`310^5V/W`$ for our detector; we find it also has a noise equivalent power (NEP) of $`0.5pW/\sqrt{Hz}`$ at a chopping frequency of 1 Hz. We estimate an optical efficiency from the source to the detector of order 10%; this is mainly due to the fact the we use a $`1/2\mathrm{"}`$ diameter light pipe, but the detector is only 4x4 mm2 in size. (No cone is used to concentrate the radiation onto the detector.) Thus, from the measured power we determine the absolute output power is roughly 100 pW for an electrical input power of 1 mW. In a dilution refrigerator environment, the emitter would have to be carefully heat sunk to a thermal stage that could handle the power load. This is possible without too much difficulty.

Our results for the emission power are consistent with those of Zinovev, who used microsecond electrical pulses and found roughly 50 pW of ac power out for roughly 1 mW in, but inconsistent with the results of Kawano, who find 50 pW out for 50 mW in. These variations in the output power within the literature may be related to the sample geometry or the optical coupling techniques used.

### D Linewidth

The spectral resolution of our spectrometer is limited by the linewidth of the source. In reality, the spectrum of states in a 2DEG in a magnetic field is broadened by disorder. We infer an upper limit on the linewidth of our source by measuring the transmission vs. magnetic field through a sample with a known resonant absorption behavior. This is discussed in further detail in section III. Based on these measurements, we infer an upper limit of 1.8 cm-1 (FWHM) on our source linewidth. These results are consistent with those of Knap and Komiyama, who studied C.R. emission from samples of similar mobility. The fact that we find such a narrow linewidth even when running the emitter quasi-continuous is important; this issue is discussed in more detail below.

If the current distribution and optical emissivity as a function of position on the surface of the sample were known, it should in principle be possible to model the interference fringes in figure 3 to determine the spectral linewidth. Zinovev modeled the interference fringes for the case of vertical emission only and found qualitative agreement with the measured fringe pattern. However, the emission angle is not purely vertical; in reality it is spread out over 2 $`\pi `$ steradians almost evenly (see above). Since we do not know the details of where in the sample the emission comes from, it is difficult to relate the contrast of the interference fringes to the linewidth of the radiation. Thus, while the presence of fringes indicates a component of the emission is monochromatic, it does not uniquely determine its linewidth, nor the amount of incoherent broadband radiation. The fringe pattern did not vary strongly with the mobility of the emitter sample ($`\mu =600,000cm^2/Vs`$ vs. $`2.710^6cm^2/Vs`$), suggesting that the “contrast” of the interference fringes is not a good measure of the emission linewidth.

### E Chopping technique

In contrast to most C.R. emission experiments, we run our sources quasi-continuous, and not pulsed, and still find reasonably narrow linewidth compared to the results found in the literature. This an important point if one wants to use this technique with ultimate state-of-the-art low-noise detectors. Fast detectors used for pulsed measurement typically have NEPs of $`10^{12}W/\sqrt{Hz}`$. With the output power of our source in the 10-100 pW range, this gives limited signal to noise. Recent progress on ultra-low noise detectors for measuring the cosmic microwave background have achieved NEPs of 2 10$`{}_{}{}^{18}W/\sqrt{Hz}`$ at an operating temperature of 100 mK. Although ultra-sensitive, they have time constants (typically 10 ms) that do not allow for pulsed operation. Thus, our quasi-continuous technique is important if those ultra-low noise detectors are to be used.

Initial work on fast low noise detectors based on lithographically fabricated hot-electron “microbolometers” began in the late 1980s. In a later variant of the idea, electrical NEPs of $`10^{17}W/\sqrt{Hz}`$ with a response time constant of 10 $`\mu `$s were measured. If these could be antenna coupled, they would make fast, sensitive detectors. This would obviate the need for the low chopping frequencies described in this paper. This general line of investigation is currently being pursued by several groups in order to demonstrate optical NEPs in that range or better with fast response time. Nonetheless, our results show that the development of fast, sensitive detectors is not a necessity for the technique described herein. Current state-of-the-are ultra-sensitive bolometer technology can be used without a degradation in the performance (i.e. linewidth) of the spectrometer due to the low chopping frequency.

## III THz Transmission Spectrometer

### A Principle of operation

In principle, our technique is simple and straightforward. The frequency of the source is tuned by its local magnetic field, and the light is guided through the sample and onto the detector. One can sweep the source frequency by sweeping its magnetic field; the power measured on the detector is a direct measure of the transmission through the sample at each frequency. An alternative mode of operation is to fix the source frequency (by fixing it’s local magnetic field) and sweep some sample parameter, to determine the transmission at a fixed frequency vs. the sample parameter. We discuss this in more detail in the next section.

### B Demonstration of operation

We have carried out a demonstration experiment to measure the transmission through another 2DEG at a fixed frequency as a function of magnetic field at the sample. For a 2DEG in a magnetic field, absorption of THz photons occurs only if photon energy matches the cyclotron energy, $`hf=eB/m^{}`$. Thus, as a function of frequency of the emitter (or magnetic field of the sample as realized here), the transmission is unity (neglecting the vacuum-dielectric mismatch) off resonance, and minimum on resonance. If the light is randomly polarized, the absorption coefficient depends on the sample mobility and density. The transmission minimum is limited theoretically to 50%, although weaker absorption is usually observed.

In our spectrometer, the light from the emitter passes through a $`1/2`$” diameter light pipe for about 3”, and is then passed through a 3 degree cone onto the transmission sample through a 1/8” hole, and then through another $`1/2`$” light pipe onto the detector. Instead of using two separate magnets, one to set the emitter frequency and one to set the sample magnetic field, we use a one-magnet setup indicated schematically in figure 3. By fixing the current through the superconducting magnet, we fix the magnetic field at the source and hence the emission frequency. (The emitter is located towards the bottom of the magnet at the 55 % field region.) We then monitor the detector voltage while sweeping the magnetic field at the sample; this is accomplished by physically moving the sample in the region where the field strength varies with position. The field profile is shown in figure 3; by moving the sample up and down two inches, we can vary the sample magnetic field from 25% to 75% of the value at the field center. With a 9 T superconducting magnet, we can achieve a maximum field of 5 T at the emitter, corresponding to a maximum emission frequency of 65 cm-1.

We plot in figure 4 the measured transmission coefficient as a function of sample magnetic field. The right axis is the measured power on the detector; the emitter frequency is fixed at 1.1 THz. We find a Lorentzian resonance with a FWHM linewidth of 0.26 Tesla (2.7 cm-1) for this sample ($`\mu =600,000cm^2/Vs`$). On another sample with somewhat higher mobility ($`2.710^6cm^2/Vs`$), we measure a linewidth of roughly 0.2 Tesla (1.8 cm-1). We conclude from this that the emitter linewidth is no broader than 1.8 cm-1; otherwise the sharp absorption feature in figure 4 would not be seen.

It is possible that there is a spectral component of the emitted power which is broadly distributed away from $`\omega _C`$. From the measurement described, we can set an upper limit on the amount of this component. Since the absorption for randomly polarized light is at most 50% and our absorption dip is about 10%, we conclude that at most 80% of the emitted power of the source is off resonance; and that at least 20% power is concentrated in the spectral region around $`\omega _c`$. However, this is only an upper limit on the out of band emission because we have not independently measured the transmission of our sample; if the true sample absorption is less than the theoretical limit of 50% on resonance, then the background component of the emission is less than 80%.

The linewidth of cyclotron resonance absorption in high mobility 2DEGs has been studied for many years by many groups. A very typical feature is that the linewidth varies periodically with filling factor $`\nu =nh/eB`$, and that the linewidth is an absolute maximum at filling factor 2. Our measurements presented here are near $`\nu =2`$, so that our measurement of the spectral resolution is probably an upper limit on the ultimate attainable spectral resolution using this technique.

In figure 6, we plot the measured linewidth vs. sample mobility at a filling factor of $`\nu =2`$, as well as several measured values from the literature. The temperature for the references varied between 1.3 K and 5 K. The linewidths we measure are comparable to those measured by other groups, suggesting that our emitter linewidth is narrow enough to measure roughly the correct value for the transmission linewidth.

We have also measured the transmission through the sample off of the resonant absorption as a function of source frequency. We find roughly unity transmission off resonance. We have also carried out this swept frequency transmission measurement though a blank sample (containing no 2DEG) of GaAs. In either case, we do not observe any Fabry-Perot fringes due to the finite (transmitting) sample thickness, even though the sample is not wedged. We conclude from this that the incoming radiation is equally distributed as a function of angle of incidence. This is an added advantage of this technique; the transmission samples do not need to be wedged.

Finally, even though the emitter must be placed in a magnetic field, many dilution refrigerators have compensating coils to minimize the magnetic field at the mixing chamber. By placing the transmission sample and detector near this compensated region, the technique described here can be applied to systems in low magnetic fields as well.

### C Ultimate sensitivity limits

Based on the noise of our detector, we can measure the power transmitted through the sample with a statistical uncertainty of roughly 0.5 pW in a one second integration time. In figure 4, the transmitted power is roughly 2 pW off resonance; we used an equivalent noise bandwidth of approximately 1 mHz for that measurement. The statistical error on each transmission measurement point is thus about 1 %. This noise performance is clearly marginal, and underscores the need for lower noise detector technology. If the spectrometer described in this paper were to be incorporated into a low temperature environment, low-noise state of the art bolometers with NEPs of as low as 2 10$`{}_{}{}^{18}W/\sqrt{Hz}`$ could be used. This would allow for a much better measurement of the transmission coefficient, with a statistical uncertainty of 10-5 %, given the power levels we use.

For a given statistical uncertainty in the measured transmission coefficient at each frequency, we can predict the statistical uncertainty on the measured linewidth and line position for a Lorentzian absorption profile. Based on numerical simulations of typical experimental parameters, we find a statistical uncertainty in the measured linewidth of roughly 25 Gauss (or 1 GHz in frequency) per 1% error in measured transmission coefficient, and roughly 2.5 Gauss (0.1 GHz) in the measured position per 1% error in measured transmission coefficient. With a detector with NEP of 2 10$`{}_{}{}^{18}W/\sqrt{Hz}`$, we could determine the linewidth and center frequency with a statistical uncertainty of order 1 kHz. We note that if this technique were used on samples with narrower absorption features than the emitter linewidth, it would not be possible to determine the shape of those absorption feature. However, it would be possible to determine the position (in frequency or magnetic field) of the absorption feature to a precision much better than the emitter linewidth.

Recent work on THz detectors has demonstrated single-photon sensitivity with unknown quantum-efficiency and proposals exist for single-photon sensitivity with excellent predicted quantum efficiency. These detectors must be operated in an ultra-low temperature environment which is very compatible with the spectrometer presented here. We now consider the ultimate limits of the spectrometer developed here if single photon THz detectors were to be used.

At low (dilution refrigerator) temperatures, the thermal background of THz photons is negligible. Therefore, the statistical uncertainty in the measured signal would be limited only by the “shot” noise of the incoming photons. Roughly $`10^9`$ photons/second are generated in a 1 pW beam. For classical statistics, the fluctuations are $`\sqrt{10^9}`$ photons/second, i.e. $`10^{17}`$ W; for quantum statistics the fluctuations can be lower or higher, corresponding to bunching or antibunching. Generally, the statistics of the radiation depends on the emission process. To our knowledge this has not been investigated (either theoretically or experimentally) in cyclotron resonance emission. These calculations suggest that perfect (noiseless) THz detectors would not significantly improve the noise performance of our spectrometer, as compared to state-of-the-art low noise detectors with NEPs of 2 10$`{}_{}{}^{18}W/\sqrt{Hz}`$, since the statistics of the signal dominate the noise performance. However, the ultimate limit of a few photons emitted per second with a single photon THz detector, where the photons interact with a quantum system of interest between the source and detector, may give experimentalists new tools to explore quantum information processing in condensed matter systems.

## IV Conclusions

Using a 2DEG as a cyclotron resonance source, we have developed an all-cryogenic THz transmission spectrometer for the spectral range of 2-65 cm-1, with a spectral resolution of at least 1.8 cm-1. Since a broadband detector is used, the entire frequency range can be accessed by the single sweep of the magnetic field at the emitter. Additionally, we have demonstrated that a low chopping frequency does not degrade the emission linewidth; this is important for future use with ultra-low noise detectors. The spectrometer was demonstrated by measuring the THz transmission through another 2DEG as a function of magnetic field of the 2DEG.

## V Acknowledgements

This work was supported by Sandia National Labs Grant No. DE-AC04-AL85000 and DOE Grant No. DE-FG03-99ER45766. One of the authors (P.J.B.) was supported in part by the Sherman Fairchild Foundation.


text:
# Varieties of Dirac equation and flavors of leptons and quarks

## I introduction

The standard model of the strong and electroweak interactions is used to describe the properties of the quarks and leptons. There are six known quarks  and six known leptons in nature. The leptons and quarks are spin $`1/2`$ fermions. They are grouped into three pairs or generations of quarks, and three generations of leptons. Quarks and leptons of same generation can form left-handed isospinors , and each member of a quark or lepton generation may be transformed into its partner via the charged-current weak interaction. Within the frame of the standard model, all these elementary fermion fields are assumed to satisfy the Dirac equation of the same form. The only difference between equations for different fermion fields are the fermions’ rest masses. Although the standard model has gained already great success, the above assumption is problematic. A quantum field is completely described by its equation and the quantization conditions. If all elementary fermion fields would satisfy the Dirac equation of the same form, then there is not a natural explanation for why a quark or lepton can be transformed only into another specific type of quark or lepton through the charged-current weak interaction, but not into other ones. In other words, there is no natural way to identify which two leptons or quarks should be grouped into a generation, because all leptons and quarks are assumed to satisfy equations of the same type. To be able to describe the behaviors of quarks and leptons under the charged-current weak interaction, we must introduce for each elementary fermion field a different equation. These equations must be grouped into pair in such a way that the vector-minus-axial-vector ($`VA`$) weak interaction current can be formed only from the fermion fields satisfying two equations of the same pair. All these equations are for spin $`1/2`$ massive fields, they must have similar properties as the original Dirac equation, and I call these equations as the varieties of the Dirac equation.

In this paper I will discuss in detail the varieties of the Dirac equation. Based on a careful analysis of the properties of the Dirac fields that satisfy these equations, I will try to answer to the questions: why there exist twelve elementary fermions , and why quarks and leptons appear in three generations.

## II Varieties of Dirac equation

A field of spin $`1/2`$ transforms under spatial rotation in the following way:

$$\psi (x)\psi ^{}(x^{})=e^{i\stackrel{}{s}\stackrel{}{\theta }}\psi (x),$$
(1)
where $`\stackrel{}{\theta }`$ is the rotation angle, and $`\stackrel{}{s}`$ the spin operator that satisfies the following relations:

$$[s_j,s_k]=iϵ_{jkl}s_l,s_j^2=\frac{1}{4},j,k,l=1,2,3.$$
(2)
For a four components spinor field, $`s_j`$ are $`4\times 4`$ matrices. There exist sixteen linear independent $`4\times 4`$ matrices, and they can be chosen as

$$I_4,\mathrm{\Sigma }_j^I,\mathrm{\Sigma }_j^{II},\mathrm{\Sigma }_j^I\mathrm{\Sigma }_k^{II},j,k=1,2,3,$$
(3)
where

$$I_4=\left(\begin{array}{cc}I_2& 0\\ 0& I_2\end{array}\right),\mathrm{\Sigma }_1^{II}=\left(\begin{array}{cc}0& I_2\\ I_2& 0\end{array}\right),\mathrm{\Sigma }_2^{II}=\left(\begin{array}{cc}0& iI_2\\ iI_2& 0\end{array}\right),\mathrm{\Sigma }_3^{II}=\left(\begin{array}{cc}I_2& 0\\ 0& I_2\end{array}\right),\mathrm{\Sigma }_j^I=\left(\begin{array}{cc}\sigma _j& 0\\ 0& \sigma _j\end{array}\right),j=1,2,3,$$
(4)
with

$$I_2=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\sigma _1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\sigma _2=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),\sigma _3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(5)
We have

$$[\mathrm{\Sigma }_j^I,\mathrm{\Sigma }_k^{II}]=0,[\mathrm{\Sigma }_j^I,\mathrm{\Sigma }_k^I]=2iϵ_{jkl}\mathrm{\Sigma }_l^I,[\mathrm{\Sigma }_j^{II},\mathrm{\Sigma }_k^{II}]=2iϵ_{jkl}\mathrm{\Sigma }_l^{II},j,k,l=1,2,3.$$
(6)
All $`4\times 4`$ matrices can be written as a linear combination of these sixteen matrices. It easy to verify that there exist only two sets of $`4\times 4`$ matrices which satisfy the relations (2), these are

$$s_j^I=\frac{1}{2}\mathrm{\Sigma }_j^I,s_j^{II}=\frac{1}{2}\mathrm{\Sigma }_j^{II},j=1,2,3.$$
(7)
The Dirac equation for fields of spin $`1/2`$ and rest mass $`m`$ is

$$\left[i\mathrm{}(\frac{}{t}+c\stackrel{}{\alpha })mc^2\beta \right]\psi (x)=0,$$
(8)
in which the matrices $`\stackrel{}{\alpha }`$ and $`\beta `$ satisfy the following conditions:

$$\begin{array}{ccc}\alpha _j^{}=\alpha _j,\hfill & \alpha _j^2=1,\hfill & [\alpha _j,\alpha _k]_+=0(jk)\hfill \\ \beta ^{}=\beta ,\hfill & \beta ^2=1,\hfill & [\alpha _j,\beta ]_+=0\hfill \end{array}\}j,k=1,2,3,$$
(9)
where the anticommutator is defined by

$$[A,B]_+AB+BA.$$
(10)
The matrices $`\stackrel{}{\alpha }`$ are closely related to the spin operator $`\stackrel{}{s}`$. The invariance of the Dirac equation under space rotation requests

$$[s_j,\alpha _k]=iϵ_{jkl}\alpha _l,[s_j,\beta ]=0j,k,l=1,2,3.$$
(11)
Conditions (9) and (11) are satisfied by matrices

$$\stackrel{}{\alpha }=\mathrm{\Sigma }_j^I\stackrel{}{\mathrm{\Sigma }}^{II},\beta =\mathrm{\Sigma }_k^I\text{and}\stackrel{}{\alpha }=\mathrm{\Sigma }_j^{II}\stackrel{}{\mathrm{\Sigma }}^I,\beta =\mathrm{\Sigma }_k^{II},kj,j,k=1,2,3.$$
(12)
Thus there exist totally twelve varieties of Dirac equation, which should be satisfied by twelve different Fermion fields.

It’s convenient to use the notation

$$\stackrel{}{\alpha }=2\lambda \stackrel{}{s}.$$
(13)
From relations (9) and (11) we get the following conditions for $`\lambda `$:

$$[\lambda ,\stackrel{}{s}]=0,[\lambda ,\beta ]_+=0,\lambda ^2=1.$$
(14)
According to Eqs. (12) we may have

$$\lambda =\mathrm{\Sigma }_j^{II},\beta =\mathrm{\Sigma }_k^{II}j,k=1,2,3,jk\text{for}\stackrel{}{s}=\stackrel{}{s}^I$$
(15)
and

$$\lambda =\mathrm{\Sigma }_j^I,\beta =\mathrm{\Sigma }_k^Ij,k=1,2,3,jk\text{for}\stackrel{}{s}=\stackrel{}{s}^{II}.$$
(16)
One may observe that $`\lambda =\mathrm{\Sigma }_1^{II},\stackrel{}{s}=\stackrel{}{s}^I,\beta =\mathrm{\Sigma }_3^{II}`$ just corresponds to the Dirac-Pauli representation. The Dirac equation (8) is invariant under Lorentz transformations. A Dirac field transforms under the Lorentz transformation in the following way:

$$\psi (x)\psi ^{}(x^{})=e^{\frac{1}{2}\stackrel{}{\alpha }\stackrel{}{\phi }}\psi (x^{})=e^{\lambda \stackrel{}{s}\stackrel{}{\phi }}\psi (x^{}).$$
(17)
By comparing the transformations (17) and (1), we find that the operator $`\lambda \stackrel{}{s}`$ plays a role in the Lorentz transformation similar to the role of $`\stackrel{}{s}`$ in the space rotation. Therefore we call $`\lambda `$ as the Lorentz multiplier. The Dirac equation is also invariant under space inversion, and we have under this transformation

$$\psi (x)\psi ^{}(x^{})=e^{i\varphi }\beta \psi (x).$$
(18)
Thus the matrix $`\beta `$ is nothing but the parity operator, and we find that the twelve varieties of Dirac equation differ among themselves by the spin operator, Lorentz multiplier and the parity operator.

## III flavors and generations of leptons and quarks

There are totally twelve varieties of Dirac equations. We can group these equations into six different pairs according to their Lorentz multipliers.

Let $`\psi _1(x)`$ and $`\psi _2(x)`$ be two different fermion fields with the same Lorentz multiplier that satisfy the following equations:

$$\left[i\mathrm{}(\frac{}{t}+c\mathrm{\Sigma }_1^{II}\stackrel{}{\mathrm{\Sigma }}^I)m_1c^2\mathrm{\Sigma }_2^{II}\right]\psi _1(x)=0.$$
(19)
and

$$\left[i\mathrm{}(\frac{}{t}+c\mathrm{\Sigma }_1^{II}\stackrel{}{\mathrm{\Sigma }}^I)m_2c^2\mathrm{\Sigma }_3^{II}\right]\psi _2(x)=0.$$
(20)
The 4-dimensional vector conservation current $`j(x)`$ of the field $`\psi _1(x)`$ is defined by

$$j_0(x)=\psi _1^{}(x)\psi _1(x),\stackrel{}{j}(x)=\psi _1^{}(x)\mathrm{\Sigma }_1^{II}\stackrel{}{\mathrm{\Sigma }}^I\psi _1(x).$$
(21)
Under space inversion $`\stackrel{}{j}(x)`$, the space components of the current $`j(x)`$, changes to $`\stackrel{}{j}(x^{})`$. We can also introduce the axial vector current $`j^A(x)`$ of the field $`\psi _1(x)`$ defined as

$$j_0^A(x)=\psi _1^{}(x)\mathrm{\Sigma }_1^{II}\psi _1(x),\stackrel{}{j}^A(x)=\psi _1^{}(x)\stackrel{}{\mathrm{\Sigma }}^I\psi _1(x).$$
(22)
The axial vector current $`j^A(x)`$ transforms as a 4-dimensional vector under continuous space-time transformations. But unlike the current $`j(x)`$, $`\stackrel{}{j}^A(x)`$ remains unchanged under the space inversion.

Because the fermion fields that satisfy Dirac equations with the same Lorentz multiplier transform in the same way under continuous space-time transformation, we may define two quasi-vector currents

$$j_0^q=\psi _1^{}(x)\psi _2(x),\stackrel{}{j}^q=\psi _1^{}(x)\mathrm{\Sigma }_1^{II}\stackrel{}{\mathrm{\Sigma }}^I\psi _2(x),$$
(23)
and

$$j_0^{qA}=\psi _1^{}(x)\mathrm{\Sigma }_1^{II}\psi _2(x),\stackrel{}{j}^{qA}=\psi _1^{}(x)\stackrel{}{\mathrm{\Sigma }}^I\psi _2(x),$$
(24)
It is easy to verify that the quasi-vector currents $`j^q(x)`$ and $`j^{qA}(x)`$ transform as 4-dimensional vectors under continuous space-time transformations. But under the space inversion, their behaviors are different from 4-dimensional vectors or axial vectors:

$`\stackrel{}{j}^q(x)\stackrel{}{j}^q(x^{})`$ $`=`$ $`e^{i(\varphi _2\varphi _1)}\psi _1^{}(x^{})\mathrm{\Sigma }_2^{II}\mathrm{\Sigma }_1^{II}\stackrel{}{\mathrm{\Sigma }^I}\mathrm{\Sigma }_3^{II}\psi _2(x^{})`$ (25)
$`=`$ $`ie^{i(\varphi _2\varphi _1)}\psi _1^{}(x^{})\mathrm{\Sigma }_3^{II}\stackrel{}{\mathrm{\Sigma }^I}\mathrm{\Sigma }_3^{II}\psi _2(x^{})`$ (26)
$`=`$ $`ie^{i(\varphi _2\varphi _1)}\stackrel{}{j}^{qA}(x)`$ (27)
$``$ $`\pm \stackrel{}{j}^q(x^{}).`$ (28)
The relation $`\mathrm{\Sigma }_2^{II}\mathrm{\Sigma }_3^{II}=i\mathrm{\Sigma }_1^{II}`$ was used in (25).

We may choose

$$e^{i(\varphi _2\varphi _1)}=i,$$
(29)
and we find that the current $`J_W(x)`$ defined by

$$J_W(x)=j^q(x)j^{qA}(x)$$
(30)
transforms just like the $`VA`$ current $`j(x)j^A(x)`$ under all space-time transformations. In the case of space inversion, we have

$`J_{W0}(x)=\psi _1^{}(x)(1\mathrm{\Sigma }_1^{II})\psi _2(x)`$ $``$ $`J_{W0}^{}(x^{})=\psi _1^{}(x^{})(1+\mathrm{\Sigma }_1^{II})\psi _2(x^{})`$ (31)
$`\stackrel{}{J}_W(x)=\psi _1^{}(x)(1\mathrm{\Sigma }_1^{II})\stackrel{}{\mathrm{\Sigma }}^I\psi _2(x^{})`$ $``$ $`\stackrel{}{J}_W^{}(x^{})=\psi _1^{}(x^{})(1+\mathrm{\Sigma }_1^{II})\stackrel{}{\mathrm{\Sigma }}^I\psi _2(x^{}).`$ (32)
In fact, $`J_W(x)`$ is just the $`VA`$ charged weak interaction currents. The possibility of forming the $`VA`$ charged weak interaction current from two different fermion fields depends on the explicit forms of the Dirac equations for these two fields. The $`VA`$ charged weak interaction currents can be formed only from Dirac fields with the same spin operator and the same Lorentz multiplier. By comparing this property of Dirac fields with the behavior of the leptons and the quarks under the charged-current weak interaction, we conclude that the leptons and the quarks of a same generation must satisfy Dirac equations with the same Lorentz multiplier and spin operator. We also conclude that different elementary fermion fields must satisfy different Dirac equations. Because there are twelve varieties of Dirac equations, there must be only twelve different favors of fermions.

We have shown so far that there exist twelve varieties of Dirac equations that define twelve different fermion fields. The matrices that define these twelve equations can transform among them under certain unitary transformations.

Let

$$U=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right),$$
(33)
then we have

$$\mathrm{\Sigma }_j^{II}=U^{}\mathrm{\Sigma }_j^IU,\mathrm{\Sigma }_j^I=U^{}\mathrm{\Sigma }_j^{II}U,j=1,2,3.$$
(34)
We have also

$$\frac{1}{\sqrt{2}}(1i\mathrm{\Sigma }_j^A)\mathrm{\Sigma }_l^A\frac{1}{\sqrt{2}}(1\pm i\mathrm{\Sigma }_j^A)=\pm ϵ_{jlk}\mathrm{\Sigma }_k^A,jl,j,l,k=1,2,3;A=I,II.$$
(35)
Thus we may transform any of these twelve Dirac equations into another by applying a unitary transformation. This unitary transformation can be either $`U`$, $`2^{1/2}(1\pm i\mathrm{\Sigma }_j^{I,II})`$, or a combination of them. Due to this symmetry in the forms of Dirac equations, one may choose any type of the Dirac equation to describe a certain lepton or quark field. But when more than one fermion fields are in consideration, one must respect the following rules: 1) fermion fields with different flavors satisfy Dirac equations of different types; 2) lepton or quark fields of the same generation have the same Lorentz multiplier and spin operator, but different parity operator.

For example, if we use the Dirac equation in the Dirac-Pauli representation for electron field $`\psi _e(x)`$

$$\left[i\mathrm{}(\frac{}{t}+c\mathrm{\Sigma }_1^{II}\stackrel{}{\mathrm{\Sigma }}^I)m_ec^2\mathrm{\Sigma }_3^{II}\right]\psi _e(x)=0,$$
(36)
then we must choose the following equation for the electronic neutrino field $`\psi _{\nu _e}`$:

$$\left[i\mathrm{}(\frac{}{t}+c\mathrm{\Sigma }_1^{II}\stackrel{}{\mathrm{\Sigma }}^I)m_{\nu _e}c^2\mathrm{\Sigma }_2^{II}\right]\psi _{\nu _e}(x)=0,$$
(37)
and perhaps the equation

$$\left[i\mathrm{}(\frac{}{t}+c\mathrm{\Sigma }_2^{II}\stackrel{}{\mathrm{\Sigma }}^I)m_\mu c^2\mathrm{\Sigma }_1^{II}\right]\psi _\mu (x)=0$$
(38)
for the muon field $`\psi _\mu (x)`$.

## IV conclusions

I have shown that there exist twelve independent Dirac equations for spin $`1/2`$ fermions. Each equation is defined by its spin operator, Lorentz multiplier and parity operator. There are two different spin operators, six different Lorentz multipliers and six different parity operators. These twelve varieties of Dirac equation can be classified into six groups according to the Lorentz multiplier. Each group contains two equations. These six groups of Dirac equations correspond to the three quark generations and the three lepton generations. According to the properties of the charged-current weak interaction, I conclude that different elementary fermion fields must satisfy different Dirac equations, and consequently there may not be more than twelve types of elementary fermions. Or in other words, there can not be elementary fermions other than the already known three generations of quarks and leptons. I have also shown that charged $`VA`$ currents can be formed only from fields with the same Lorentz multiplier. This property of Dirac fields implies that quarks and leptons may be transformed via the charged-current weak only into their partner within the same generation.


text:
# First Lasing of Volume FEL (VFEL) at Wavelength Range 𝜆∼4-6 mm

## 1 Introduction

The most essential feature of FEL and other types of generators is a feedback, which is formed by a system of mirrors, or distributed feedback based on diffraction in spatially periodic medium, when wave vectors of transmitted and reflected waves are colinear. The distinction of volume FEL (VFEL) is non-one-dimensional multi-wave volume distributed feedback (VDFB). VDFB performs two basic functions simultaneously:
1) it provides more effective interaction of an electromagnetic wave with an electron beam due to new dispersion law;
2) it forms volume distributed mirror retaining radiation in interaction region.

It is well known that each radiative system is defined by its eigenmodes and by the so-called dispersion equation, which in the case of small perturbations (linear regime) describes possible types of waves in system and relation between frequency and wave number of system eigenmodes. Thorough analysis of properties of FEL dispersion equation is done in . It is shown there that
1. dispersion equation for FEL in collective regime coincides with that for conventional travelling wave tube amplifier (TWTA) ;
2. FEL gain (increment of electron beam instability) is proportional to $`\rho ^{1/3}`$, where $`\rho `$ is the electron beam density.

For the first time the possibility of essential change of instability of an electron beam moving in a spatially periodic medium was indicated in . In this paper the dispersion equations are obtained and investigated for conditions of multi-wave diffraction. It is shown that there is a new law of electron beam instability in the points of diffraction roots degeneration. The amplification and generation gain of electromagnetic wave are sharply changed in these points. In  this statement is derived from the fact that the instability increment in the point with s-times degeneration is proportional to $`\rho ^{1/(3+s)}`$ (here $`\rho `$ is the electron beam density and it is supposed that the point in which s roots coincide is the point with (s-1)-times degeneration). This increment differs from the relevant increment for one-wave system which is proportional to $`\rho ^{1/3}`$. This result is also valid for electron beam which moves in vacuum close to the surface of spacially periodic medium  (or in a vacuum slot made inside a periodic medium). Explicit expressions for dependence of starting current $`j`$ on interaction length $`L`$ are obtained in degeneration points : $`j_{start}1/\{(kL)^3(k\chi _\tau L)^{2s}\}`$. The advantages of VFEL are exhibited in wide spectral range from microwaves to X-rays . In  the experimental simulation of electrodynamic processes in the volume diffraction grating was performed for a millimeter wavelength range. The possibility of obtaining of extremely high Q-factor for a system with two strongly coupled waves is experimentally confirmed.

In the present work the first lasing of VFEL in millimeter wavelength range is reported. The main parts of VFEL are two flat diffraction gratings with different spatial periods. The first diffraction grating is used for radiation generation due to Smith-Purcell radiating mechanism . The second grating provides distributed feedback using Bragg dynamical diffraction , the grooves of diffraction grating are oriented at the nonzero angle to the direction of electron beam velocity. It should be noted, that generation in nonrelativistic TWT devices in similar conditions is impossible at single harmonic, since the wavelength of emitted radiation considerably exceeds period of TWT spiral (or corrugation period). Lasing of considered VFEL type for the first time was described in , where the theoretical model of its operation was presented.

## 2 Experimental setup and results

Block-scheme of the experimental setup is illustrated in fig.1.

Volume resonator is the basic functional part of VFEL (see Fig.2). In considered setup it is formed by two diffraction gratings with different periods and two smooth side walls. Cross section of the resonator is rectangular and constant along all its length. But the distance between diffraction gratings can be varied in different experiments. Radiation is output from the resonator through end-wall.

The interaction of the first diffraction grating (exciting grating) with the electron beam arouses generation of Smith-Purcell radiation. The second (resonant) grating provides distributed feedback of generated radiation with electron beam. The distance between grooved surface of exciting grating and electron beam can be varied during experiment that allows to encrease efficiency of generation process. Rotation of resonant diffraction grating provides the possibility to change radiation frequency by varying the orientation of grating grooves with respect to electron beam velocity. Mechanical tuning of the radiation frequency is also provided by the adjustment of distance between gratings.

The sheet electron beam emitted by the thermal cathode (tungsten labilized by barium or calcium iridate) is formed in cylindrical Peirce gun. It is guided in the VFEL resonator by the magnetic field $`3KGs`$. Electrons are emitted in pulsed regime (unipolar pulse with sinusoidal shape and pulse duration $`10ms`$) in sequence of two or three voltage pulses (voltage amplitude can vary from 1 to 10 KeV).

Generated radiation is output through a radiotransparent window (plexiglas) to a radiation detector (thermistor detector M5-50 or power meter M3-22A) with a bandwidth of detected radiation 54-78 GHz.

Control and measurement system provides output of registed data to PC display. Also the following parameters are recorded on a hard disk drive:

* Cathode voltage;
* Total electron beam current;
* Current in a guiding magnetic system windings;
* Grating current;
* Power of microwave radiation.

The generation of microwave radiation was detected at electron energy $`5KeV`$ in millimeter wavelength range ($`\lambda 46`$ mm). Parameters of resonator are: the length of resonator is 100 mm;
the period of the exciting diffraction grating is 0.67 mm;
the period of the resonant diffraction grating is 3 mm.
The detected pulse power of generated radiation was about 3 - 4 W. Taking into account the fact, that only the part of the electron beam with cross size $`\delta \frac{\lambda u}{4\pi c}<0.1mm`$ interacts effectively with electromagnetic wave the efficiency of ”working” part of a beam can be evaluated as $`10\%`$. The oscillogram of VFEL lasing is shown in fig. 3.

The continuous generation of radiation was observed for electron beam currents higher than $`35mA`$ . To define the threshold conditions the beam current was changed by the varying of the cathode heating voltage. Threshold conditions of generation are illustrated in fig.4, where the dependence of grating current and emitted microwave radiation on cathode heating voltage is shown. The decreasing of heat voltage from 178 V to 158 V does not cause the decreasing of the grating current because of saturation mechanisn due to Child-Langmuir law. But at $`U_{heat}156V`$ the grating current decreases. It is clear from fig. 4 that reduction of the current leads to decreasing of radiation power. And the peaks on radiation power curve concentrates near electron beam energy $`5keV`$ when the grating current tends to the threshold value. The peaks on radiation curve appears at minimum value of grating current $`35mA`$. The lasing starts from voltage $`5kV`$. This value is defined by the diffraction gratings period. As a result at the electron beam energy $`5kV`$ the mode with the highest Q-factor is excited. For this reason when the grating current tends to the threshold value, the region of generation tends to a narrow band near to $`5kV`$. At higher current values the radiation appears in electron energy range 5 - 7.5 KeV that demonstrates the excitation of the other modes with smaller Q-factor, for which the working current exceeds the starting that.