text:
# Optimization of Green-Times at an Isolated Urban Crossroads

## I Introduction

Over the past decade, the vehicular traffic problems has been intensively investigated within the context of statistical physics ( for a review see Ref.  ). Modeled as a system of interacting  particles driven far from equilibrium, vehicular traffic presents the possibility to study various aspects of truly non-equilibrium systems which are of current interest in statistical physics . The majority of these studies have been allocated to the highway traffic . The other area in vehicular traffic is urban traffic which have also been studied by statistical physicists . The simulation of traffic flow in a large-sized city is a formidable task and many degrees of freedom have to be involved (see e.g. ). In practice, streets of a city form a network of junctions that are linked together. Each crossroads receives demands (vehicles attempting to pass the cross ) and at each crossroads there exists a traffic light which, with some certain programming, controls the transportation. The first model introduced by statistical physicists for the description of the city network, known as the BML model in the literature, uses a deterministic cellular automata framework  and predicts a sort of phase transition from free-flow to a jammed state. In the model, each site of a square lattice represents the crossing of a single-lane east-west street and a single-lane north-south street (no turning of vehicles are allowed). The state of east-bound vehicles are updated synchronously at every odd discrete time steps whereas those of the north-bound vehicles are updated in parallel at every even time step following a rule which is simple extension of totally asymmetric simple exclusion process (TASEP) .

The BML model has been generalized to take into account several realistic features of traffic in cities such as asymmetric distribution of cars , faulty traffic lights , independent turning of the vehicles  and green-wave synchronization . At first sight, the BML model and the above-mentioned extensions seem unrealistic because the vehicles hop from one crossing to the next. In a more realistic generalization, each connecting bond of neighbouring junctions were replaced by a decorated bond  . Later, a more serious model of city network was introduced by Chowdhury and Schadschneider (CS model)  which developed a more detailed fine-grained  description of city traffic. This model combines the BML mode together with the Nagel-Schreckenberg model of highway traffic . In the CS model, the signals are synchronized in such a way that all the signals remain green for the east-bound vehicles for a time interval $`T`$ and then, simultaneously turn red for the east-bound vehicles. To the best of our knowledge, The CS model and its generalization  is the most realistic model introduced by the statistical physicists for city network but nevertheless, despite its nice formulation, there are still a lot of simplifications which prevent it from being an effective  and applicable model to a city network.

In spite of developing models for describing the city network, no detailed description of an isolated crossroads has yet been explored. We strongly believe that in order to have a better insight to the problem of city network, one must have a clear picture at single crossroads which would definitely be of great importance for an optimal programming of the entire network. The empirical mechanism by which the traffic lights are controlled is generally divided into two distinct methods : fixed time and real-time. In the fixed-time method, a fixed value of time is allocated to the traffic light as well as its sub-phase times. In the real-time method, which is becoming increasingly popular in great cities, the ensemble of crossroads are intelligently controlled by a central controller. The control mechanism is usually based on the concept of producing a kind of green waves  between the crossroads. These waves interact with each other, and if the passing-demands are of high values, the green waves may have destructive effects on each other, and hence, the concept of green wave may fail to be the final solution for optimizing the overall flow. Knowing the local optimum behaviour of single intersections could give us an appropriate criterion to control the set of crossroads in an adaptive manner. Additionally, not all the intersections of a city are highly affected by the neighbouring intersections. To a good approximation, marginal intersections could be regarded as isolated crossroads and therefore single-crossroads optimization strategies need to be investigated. In this paper we aim to analyze a single crossroads in detail in order to find a better insight to the problem of optimizing the total flow in cities.

The organization of this paper is as follows: In section two, we introduce the model, state our strategy for optimizing the traffic flow at the one-way to one-way intersection and finally obtain the optimum green-times of the corresponding phases. In section three, we extent our model to allow for a three-phase one-way to two-way  crossroads. Section four is devoted to some empirical data on two of the Tehran crossroads which are under the control of an intelligent traffic controller system. We compare the averaged green-time proposed by system to the optimum green-time of our fixed time theory. Our theory input-parameters are obtained via the empirical data. Finally we give our concluding remarks in section five.

## II  Formulation of the Model

### A One-Way to One-Way Crossroads

Let us consider a single crossroads which is the result of the intersection of two perpendicular streets. In their simplest structure, these streets can each direct a one-way traffic flow. With no loss of generality, we take them as one-way South to North (S-N) and West to East (W-E) streets. Cars arrive at the south and the west entrances of the crossroads. In our model, we assume that the arrival rates of the cars, i.e., the number of cars reaching the crossroads per second, are constant in time. Although everyday driving experiences in cities indicates that these rates have inevitable fluctuations in the course of time, yet in definite time intervals, the assumption of the constant arrival rates could be justified at least on an average level. As will be seen in what follows, this assumption leads to great simplifications. We take the arrival rates to be $`\alpha _1`$ (for S-N cars ) and $`\alpha _2`$ (for the W-E cars) respectively. Also we denote the passing-rate of cars ( number of cars passing the crossroads in the unit of time during the green-phase ) by $`\beta _1`$ and $`\beta _2`$. The period of the traffic lights is taken to be a definite value $`T`$ which is assumed to remain constant. The starting time of each cycle of the traffic light is the moment at which the light turns green for the S-N street. The S-N light remains green for $`T_1`$ seconds. At $`T_1`$ the traffic lights turn red for the S-N street and simultaneously changes to green for the W-E street. This is the beginning of the second phase which continues from $`T_1`$ to $`T`$ (end of the cycle). During Phase I ( $`0tT_1`$ ), the S-N cars can pass the crossroads northwards and W-E cars are stopped for the red light. Over Phase II ( $`T_1tT`$ ), the S-N cars must wait behind the red light whilst the W-E cars are eastwards passing the crossroads.

Now the basic question is ” how should traffic engineers adjust the value of $`T_1`$ in order to optimize the traffic flow through the intersection?”. With the assumption that the number of passengers in each car takes an equal average value for each direction, the optimization task is realized by minimizing the total waiting time of the cars per cycle of the traffic light. For this purpose, we introduce two quantities $`N_1`$ and $`N_2`$ which represent the number of cars stopping (queues) at the red lights in the red phases of the S-N and W-E streets respectively. Clearly $`N_1`$ and $`N_2`$ are functions of time and, in general, are divided into different lanes on the streets. The dynamics of $`N_1`$ and $`N_2`$ are read from the following equations in which $`n`$ denotes the cycle number of the traffic light:

$$N_1(nT+T_1)=[N_1(nT)+(\alpha _1\beta _1)T_1]\theta $$
(1)
$$N_2(nT+T_1)=N_2(nT)+\alpha _2T_1$$
(2)
$$N_1((n+1)T)=N_1(nT+T_1)+\alpha _1(TT_1)$$
(3)
$$N_2((n+1)T)=[N_2(nT+T_1)+(\alpha _2\beta _2)(TT_1)]\theta $$
(4)
The $`\theta `$ symbols ensure the positiveness of the quantities in the brackets, i.e., the value of $`\theta `$ is one if the quantity in the bracket is positive and zero elsewhere. This limitation is dictated to us since, by definition, the quantities $`N_1`$ and $`N_2`$ can only take positive values ( they denote queues’ lengths). Let us investigate the different situations in more details. For instance in the equation $`(1)`$, the case $`\theta =1`$ which corresponds to $`N_1(nT+T_1)>0`$ describes the situation that after the S-N lights goes red (in the $`n`$-th cycle), the whole queue of vehicles has not pass the cross and only a part of the queue has managed to pass during the green-phase. The other case, i.e., $`\theta =0`$ which corresponds to $`N_1(nT+T_1)=0`$ indicates that the $`n`$-th queue waiting in S-N direction has completely passed the cross during the time interval $`nTtnT+T_1`$. The same arguments apply to eq. (4).

We now define the total waiting time (TWT) of the vehicles per cycle of the traffic light. It is the total time wasted by the vehicles during their stop in the red phases. The TWT is the sum of the sub-waiting times of each direction. Denoting the TWT and the sub waiting times by $`T^w`$, $`T^{(w,1)}`$ and $`T^{(w,2)}`$ respectively, the following equations could be written for the $`n`$-th cycle.

$$T_{nn+1}^{(w,1)}=N_1(nT+T_1)(TT_1)+\frac{1}{2}\alpha _1(TT_1)^2$$
(5)
$$T_{nn+1}^{(w,2)}=N_2(nT)T_1+\frac{1}{2}\alpha _2T_1^2$$
(6)
Let us consider eq.(5). The waiting time of the $`n+1`$-th S-N queue is divided into two part. The first part is related to the initial length of the queue: $`N_1(nT+T_1)`$. If this initial part has a non-zero length, then the time wasted by the initial vehicles is simply their number times the total period of the red phase. This leads to the first term of eq.(5). The second part is related to the contribution given by the new oncoming vehicles arriving at the S-N direction of the crossroads during the red period which lasts for $`TT_1`$ seconds. Since we have assumed the vehicles arrive at a constant rate $`\alpha _1`$, the time wasted by the vehicles arriving in the infinitesimal interval $`[t,t+dt]`$ of the red interval is simply their number ($`\alpha _1dt`$) times the remaining time to the green signal ( $`TT_1t`$). Therefore the total contribution is given by integrating over the red period.

$$_0^{TT_1}\alpha _1𝑑t(TT_1t)=\frac{1}{2}\alpha _1(TT_1)^2$$

Which is the second term on the right hand side of eq.(5). Similar arguments lead the eq.(6).

As clearly can be seen, the analytical expression of the TWT strongly depends on the positiveness of the queues’ lengths just after the lights go red. In what follows, we show that different traffic status can be identified according to the behaviour of the quantities $`N_1(nT+T_1)`$ and $`N_2(nT)`$. Let us look at the first cycle, i.e., $`n=0`$. If $`N_2(0)=0`$ (a complete passing of the previous W-E queue) then it is easily seen that in order to have a complete passing of the next W-E queue, one should have the condition $`\alpha _2T\beta _2(TT_1)0`$. In this case, one has $`N_2(T)=0`$. It can be easily verified that that provided that the above stability condition holds, we have no W-E queues after the W-E light goes red for the general $`n`$-th cycle. This characterizes a light traffic state in which the whole queue can pass the crossroads during one green time. Similar arguments for the S-N direction shows that provided the stability condition $`\alpha _1T\beta _1T_10`$ holds, we have a stable, light traffic state in the S-N direction and that $`N_1(nT+T_1)=0`$ for general $`n`$. The case which traffic condition is light in both directions ( State I ) is characterized by the following stability conditions :

$$\alpha _2T\beta _2(TT_1)0;\alpha _1T\beta _1T_10$$
(7)
Which result in the following relations:

$$N_1(nT+T_1)=N_2(nT)=0,n=1,2,\mathrm{}$$
(8)
In this state, The TWT is independent of the cycle number $`n`$. Putting the above conditions into eqs.(5,6) and

minimizing the TWT with respect to $`T_1`$ leads to the following equation:

$$T_1=\frac{\alpha _1}{\alpha _1+\alpha _2}T$$
(9)
Inserting the above answer in the stability conditions (7) yields to the following constraints among the rates.

$$\beta _2\alpha _1+\alpha _2,\beta _1\alpha _1+\alpha _2$$
(10)
We now investigate a totally different situation, i.e., a crowded crossroads in both directions. Let us again look at the first cycle. Supposing the first cycle is characterized by the conditions $`N_2(0),N_1(T_1)>0`$. one could easily verify that provided the following relations hold :

$$\alpha _1>\beta _1,\alpha _2T>\beta _2(TT_1)$$
(11)
We have a stable condition in the next cycles.

$$N_1(nT+T_1),N_2(nT)>0$$
(12)
In sharp contrast to the State I, in this state which is characterized by the equation (11) and referred to as the state II, the values of $`N_1`$ and $`N_2`$ are functions of the cycle number. This is easily seen by the following relations:

$$N_1(nT)=N_1(0)+n(\alpha _1T\beta _1T_1)$$
(13)
$$N_2(nT)=N_2(0)+n[\alpha _2T\beta _2(TT_1)]$$
(14)
$$N_1(nT+T_1)=N_1(0)+n(\alpha _1T\beta _1T_1)+(\alpha _1\beta _1)T_1$$
(15)
$$N_2(nT+T_1)=N_2(0)+n[\alpha _2T\beta _2(TT_1)]+\alpha _2T_1$$
(16)
The above relations show that queues’ lengths grow linearly with time and a complete passing of a queue in one cycle is not possible. Drivers should wait more than one cycle in order to pass the crossroads. It can be shown that in the large cycle-number limit, the sub-waiting times are as follows:

$$T_{nn+1}^{(w,1)}n(TT_1)(\alpha _1T\beta _1T_1)$$
(17)
and

$$T_{nn+1}^{(w,2)}nT_1(\alpha _2T\beta _2(TT_1))$$
(18)
We now minimize the TWT with respect to $`T_1`$ which leads to the following equation:

$$T_1=\frac{\beta _1+\beta _2+\alpha _1\alpha _2}{2(\beta _1+\beta _2)}T$$
(19)
the consistency of this solution with stability conditions (11) yields the following constraints:

$$\alpha _1>\beta _1,\alpha _1\beta _2+2\alpha _2\beta _1+\alpha _2\beta _2\beta _1\beta _2+\beta _2^2$$
(20)
Also positiveness of $`T_1`$ itself imposes the extra restriction $`\alpha _2<\beta _1+\beta _2+\alpha _1`$.

Next we consider the situation (state III) where in the first cycle of the traffic light, one has $`N_1(T_1)=0`$ but $`N_2(0)>0`$. This corresponds to the situation in which the S-N street has a light traffic flow while the W-E street has a heavy one. The conditions for a stable pattern are:

$$\alpha _1T\beta _1T_10,\alpha _2T\beta _2(TT_1)0$$
(21)
The above stability conditions ensures the following relations:

$$N_1(nT+T_1)=0,N_2(nT)>0$$
(22)
In the large cycle-number limit, minimizing of the TWT leads to the value $`T_1=\frac{\beta _2\alpha _2}{2\beta _2}T`$. It could be easily checked that the above solution is inconsistent with the stability conditions, and hence, is not acceptable as an optimum signalization of traffic light. In the region determined by the stability conditions (21), the TWT is positive definite and therefore its minimum coincides with the upper limit of the inequalities (21). Therefore one finds:

$$T_1=max(\frac{\alpha _1}{\beta _1},\frac{\beta _2\alpha _2}{\beta _2})T$$
(23)
A similar argument applies to the final state ( state IV ) which is characterized by $`N_1(nT+T_1)>0`$ and $`N_2(nT)=0`$ ( it is sufficient to interchange the indices one and two ).

## III One-Way to Two-Way Crossroads

At this stage, we consider another frequent type of a crossroads. Here we let the vehicles move in both S-N as well as N-S directions but still the vehicles in the W-E street are restricted to move eastwards. This situation describes a one-way to two-way urban intersection. Consequently each cycle consists of three phases. In the first phase which lasts for $`0tT_1`$, the traffic light is green for the S-N cars and red for the other two directions. During the second phase which starts at $`T_1`$ and finishes at $`T_2`$, the traffic light is green for the N-S cars and red for the other two directions. In the final phase, which lasts for $`T_2tT`$, the traffic light remains green for the W-E cars and red for the N-S as well as S-N directions. The entrance rates are taken to be $`\alpha _1,\alpha _2`$ and $`\alpha _3`$ and we denote the passing rate by $`\beta _1,\beta _2`$ and $`\beta _3`$ for each direction respectively. The starting time of the cycles is chosen to be the moment at which the traffic light turns green for the S-N direction. Similar equations for the queues’ lengths $`N_1,N_2`$ and $`N_3`$ could be written down and in principle one can evaluate the TWT in terms of these quantities. The exact form of the TWT strongly depends on the positiveness of the queues’ lengths just after the traffic lights goes red for the respective direction. In the case under consideration, eight different possibilities are identified due to the traffic conditions (two possibilities for each movement direction). Now we write the explicit expression for the TWT during the $`nn+1`$ cycle.

$$T_{nn+1}^{(w,1)}=N_1(nT+T_1)(T_2T_1)+\frac{1}{2}\alpha _1(TT_1)^2+N_1(nT+T_2)(TT_2)$$
(24)
$$T_{nn+1}^{(w,2)}=N_2(nT)T_1+\frac{1}{2}\alpha _2T_1^2+N_2(nT+T_2)(TT_2)+\frac{1}{2}\alpha _2(TT_2)^2+\alpha _2T_1(TT_2)$$
(25)
$$T_{nn+1}^{(w,3)}=N_3(nT)T_1+\frac{1}{2}\alpha _3T_2^2+N_3(nT+T_1)(T_2T_1)$$
(26)
Let us only discuss the most probable one which corresponds to light traffic conditions in all directions:

$$N_1(nT+T_1)=N_2(nT+T_2)=N_3(nT)=0$$
(27)
It could be easily verified that the triple stability condition for the validity of the above assumptions are as follows:

$$\alpha _1T\beta _1T_10,\alpha _2T\beta _2(T_2T_1)0\text{and}\alpha _3T\beta _3(TT_2)$$
(28)
Inserting eq.(27) into eqs.(24-26) and minimizing the TWT with respect to $`T_1`$ and $`T_2`$ leads to the following fixation of $`T_1,T_2.`$

$$T_1=T\frac{\alpha _1(\alpha _2+\alpha _3)\alpha _2\alpha _3}{\alpha _1(\alpha _2+\alpha _3)+\alpha _2\alpha _3}$$
(29)
$$T_2=T\frac{2\alpha _1\alpha _2}{\alpha _1(\alpha _2+\alpha _3)+\alpha _2\alpha _3}$$
(30)
One directly observes that in the symmetric case of equal arrival rates ($`\alpha _1=\alpha _2=\alpha _3`$), $`T_1`$ and $`T_2`$ take the expected values $`\frac{T}{3}`$ and $`\frac{2T}{3}`$ respectively. The other point which must be mentioned is that in this traffic state, $`T_1`$ and $`T_2`$ do not depend on the passing rates, and are solely determined by the arrival rates.
Another extreme is the situation when all of the three directions are carrying a heavy traffic flow. It could be anticipated that the stability conditions for the positiveness of $`N_1(nT+T_1),N_2(nT+T_2)`$ and $`N_3(nT)`$ are $`\alpha _1>\beta _1`$, $`\alpha _2T_2\beta _2(T_2T_1)>0`$ and $`\alpha _3T\beta _3(TT_2)>0`$. In this case, all the three sub-waiting times depend on the cycle number $`n`$ and it can be shown that in the large $`n`$ limit, the minimization of the TWT give rises to the following values for $`T_1`$ and $`T_2`$.

$$T_2=\frac{T}{2}\frac{2\beta _1\beta _2+\beta _2\alpha _1+\beta _1\alpha _2\alpha _3(\beta _2+\beta _1)+\beta _3(\beta _1+\beta _2)}{\beta _1\beta _2+\beta _2\beta _3+\beta _3\beta _1}$$
(31)
$$T_1=\frac{T}{2}\frac{\beta _1\beta _2+\beta _3\beta _1+\beta _2\alpha _1+\beta _3\alpha _1\alpha _2\beta _3\beta _2\alpha _3}{\beta _1\beta _2+\beta _2\beta _3+\beta _3\beta _1}$$
(32)
Here in contrast to (29) and (30) , the passing rates appear in the expressions of $`T_1`$ and $`T_2`$ and it can be seen that in the fully symmetric condition, one again obtains that $`T_1`$ and $`T_2`$ are one-third and two-thirds of $`T`$ respectively.

## IV  Empirical Data

For comparison of our model to the empirical data, a time-series analysis on two of Tehran intersections were carried out. The central part of Tehran is under control of the SCATS (Sydney Co-ordinated Adaptive Traffic System)  that is an intelligent traffic controller. The strategy followed by most of urban traffic control systems is based on establishing green-waves along the major streets of cities. One popular strategy consists of dividing the city intersections into different sets. Each set has a leading mother crossroads (the prime crossroads in the set) and a lot of offspring crossroads. Each set is linked to the neibouring ones. The signalization of the crossroads network is determined by the implementation of green wave between adjacent sets. The allocated green times at a crossroads is proportional to the number of vehicles (traffic volume) approaching to the crossroads. Since in general, there is a natural fluctuation in the traffic volume, the amount of green times and hence the complete cycle time of traffic lights are variables in an intelligent control method.

We considered two different intersections. The first one which is located in Tehran downtown connects Valiasr St. to Takhtejamshid St. Both of these are one-way and major streets. The other crossroads connects Abbasabad St.(one-way) to Mahnaz St.(one-way) . In contrast to the previous case, here the first street is a major while the second street is a minor one. The data set is provided by magnetic counting loops which are installed just before the pedestrian-lines of each crossroads. The data were collected on second of July 2000. The data consist of the sub-phase green times and the numbers of vehicles passed during the green times for each cycle of the traffic lights (traffic volumes). In our major-to-minor crossroads and for each cycle of the traffic light, we evaluated the ratio of the major street green-time to the total cycle time. We call this quantity the time-cycle-ratio . Similarly can can consider the number of passed vehicles, i.e., traffic volume in a cycle and introduce the volume-cycle-ratio  which for each cycle is obtained by dividing the number of passed vehicles during a sub-phase green time to the whole number of passed vehicles during a complete cycle. Figures 1-3 belong to the major-to-minor crossroads. As seen from the graphs, the time-cycle-ratio allocated to the major street strongly fluctuates due to demand fluctuations received by the crossroads. In order to have a rough estimation of our model parameters, we considered a two-hour time interval between 12:30-14:30 during which we have the least fluctuations in the traffic volume. It is empirically observed that in this two-hour period, the traffic state is light and the queues are cleared during one cycle. Therefore, the optimal green times should be evaluated from eq.(9). According to this equation, we only need to know the the ratio of $`\frac{\alpha _1}{\alpha _2}`$.

We should mention that due to the position of the magnetic counting loops, they are unable to measure the upstream fluxes which are directly related to the parameters $`\alpha _1`$ and $`\alpha _2`$. For a better estimation of the in-flow parameters, one should install another set of magnetic loops a few meters upwards the pedestrian-lines. We approximated the ratio of $`\frac{\alpha _1}{\alpha _2}`$ by the ratio of the traffic volume which has passed through the major street to the traffic volume of the minor one during the two-hour period. The two-hour traffic volumes are simply obtained by adding the cycle volumes together. This yields the value $`\frac{\alpha _1}{\alpha _2}=0.35`$. Putting this value into equation (9) yields $`\frac{T_1}{T}=0.74`$. On the other hand, the averaged value of the empirical time-cycle-ratio of the major street over the two-hour period leads to the result $`\frac{T_1}{T}=0.64`$ which differ by ten percent from the value predicted by the theory. Figures 4-6 belong to our major-to-major crossroads. Here we focused on the interval 13:30 - 15:30. Analogous to the major-to-minor, the traffic state is light. The empirical averaged value of $`\frac{T_1}{T}`$ is 0.52 while the same amount evaluated from eq. (9) is 0.56 ($`T_1`$ refers to the green time of Takht. Street). Here we observe that the difference between the fixed time method (theory) and real-time method (intelligent control) is less than one obtained in the major-to-minor crossroads. As depicted from the diagrams, in the major-minor crossroads, we observe more fluctuations in the time-cycle-ratio in comparison with the major-major crossroads. These fluctuations are enhanced in the volume-cycle-ratio. The least fluctuation belongs to the time-cycle-ratio of the major-major crossroads. Comparing the averaged values of these time-cycle-ratio values over certain intervals leads to a better understanding of the traffic state.

## V Summary and Conclusion

In conclusion, we have developed a prescription for the traffic-light programming at a single urban crossroads. The method is based on minimizing the total waiting-time of cars stopping in the red phases of the traffic light. In our model the total period of the cycle is assumed to be a fixed value. We have also assumed that vehicles arrive at the crossroads with constant time rates. This is equivalent to a constant time-headway between cars. In reality we have a fluctuating time-headway due to the natural fluctuation in the traffic volume. As a first stage of an analytical treatment, we have taken the arrival rates to be constant in time. However, for a more realistic description, one should remove this restriction and assume that the time headway satisfies a random distribution function. Work along this assumption is in progress. The other point concerns the passing-rates. One should note that throughout the paper, the passing-rate of cars from the crossroads are taken to be constants. This is valid only if the green-phase time is not so long such that the time-headways between the cars exceed certain values. The value of $`T`$ should be so tuned that during the green phases, time-headway is less than a certain value. In fact, in the model, the values of the passing rate refer to the maximum capacity of cross ( maximum number of cars passing the cross in the unit of time ), which is plausible if the $`T`$ is appropriately adjusted with the congestion of the crossroads. The empirical value of the passing rates are determined by the crossroads characteristics such as road conditions, number of lanes, speed limits etc. Our model is more appropriate for rush ours during which the time-headways are minimum and the crossroads are operating with their maximum capacities. In the following table, we have summarized our theoretical optimized green times including stability conditions in four different traffic states of the one-way to one-way crossroads.

Table

| Traffic state | Optimized green time | Stability conditions |  |
| --- | --- | --- | --- |
| State I (S-N light, W-E light) | $`T_1=\frac{\alpha _1}{\alpha _1+\alpha _2}T`$ | $`\beta _1,\beta _2\alpha _1+\alpha _2,`$ |  |
| State II (S-N heavy, W-E heavy) | $`T_1=\frac{\beta _1+\beta _2+\alpha _1\alpha _2}{2(\beta _1+\beta _2)}T`$ | $`\alpha _1\beta _1,\alpha _1\beta _2+2\alpha _2\beta _1+\alpha _2\beta _2\beta _1\beta _2+\beta _2^2`$ |  |
| State III (S-N light, W-E heavy) | $`T_1=max(\frac{\alpha _1}{\beta _1},\frac{\beta _2\alpha _2}{\beta _2})T`$ | $`\alpha _1\beta _1,\beta _2\alpha _2`$ |  |
| State IV (S-N heavy, W-E light) | $`T_1=max(\frac{\alpha _2}{\beta _2},\frac{\beta _1\alpha _1}{\beta _1})T`$ | $`\alpha _2\beta _2,\beta _1\alpha _1`$ |  |

Optimizing the traffic at each crossroads is the stating point of the more comprehensive problem of city network. Nevertheless, our model best suits those marginal intersections of cities where the effect of the other crossroads is suppressed.

## VI  Acknowledgments

M.E.F is grateful to Tehran Traffic Control Company  and in particular to F. Zolfaghari for the data support. We highly appreciate V. Karimipour for fruitful comments and critically reading the manuscript and would like to express our gratitude to R. Sorfleet, R. Gerami, M. Arab salmani, M. Salmani and M. Khoshechin for useful helps.


format: text

text:
# On the Masses of the Leptons, Bosons and Quarks

## 1 Introduction

From a postulated model system for the elementary particles , we have found a simple mass equation that relates the observed rest mass of the massive leptons to the mass of the W boson. The relationship in logarithmic form is:

$$\mathrm{ln}m=\frac{1}{2\pi \alpha }\left(\frac{1}{2}\right)^\omega +\mathrm{ln}\left(\frac{m_W}{\pi }\right)$$
(1)
where $`m`$ is the lepton mass in MeV, $`\alpha `$ is the fine structure constant (1/137), $`\omega `$ is the winding number of the particle gauge field and equals 1, 2, or 3 for the electron, muon, and tau respectively; $`m_W`$ = 80,430 MeV (the observed mass of the W boson). A plot of this relationship is shown in Fig. 1 where the points are the experimental rest masses for the three generations of massive leptons and the line is given by Eqn. (1.1). While the fit shown in Fig. 1 looks good, one wonders if it is really significant or perhaps just coincidence. Having only three points on the line is worrisome. We address this worry here by showing how Eqn. (1.1) can be extended to include the massive gauge bosons and the quarks. The resulting equation, when interpreted physically in model terms, leads to surprising conclusions concerning the origin of particle mass, and why the fermion masses (excepting the top quark) are so small compared to the boson masses.

Our model system for the elementary particles  will be summarized briefly in order to define the quantum numbers and features important in determining their masses. Then we discuss how the mass equation can be extended to include the bosons and quarks. Lastly we consider the physical interpretation and the implications of the mass equation.

## 2 Basic features of the particle
models.

We postulate that:

$``$ First generation massive particles are light-front solitons made of transverse harmonic oscillator excitations of U(1) gauge fields propagating in a circle at v = c.

$``$ For once around, the leptons have one node (n = 1), the bosons have two nodes (n = 2), and the quarks have three nodes (n = 3).

$``$ The gauge fields have a polarization that twists about the propagation path by $`\pi `$ radians as the front propagates from one node to the next. This is required for the wave function to consistently meet the cyclic boundary conditions.

$``$ The fermions (leptons and quarks) are spinor particles; their fields must go twice around for a complete cycle of the U(1) field. The bosons are vector particles; their fields only go once around for a complete U(1) cycle.

$``$ The twist direction may or may not reverse at a node. This distinguishes one particle from another. For example, the $`W^{}`$ has the twist pattern L-L for one cycle, while the Z has a twist pattern of L-R or R-L. So the W has a total twist of 2$`\pi `$ radians per cycle while the Z has a total twist of zero. This is the source of the mass splitting between the W (80.4 GeV) and the Z (91.2 GeV).

$``$ First generation particles have a writhe = 0. That is, they have planar propagation paths. Second generation particles have writhe = 1, and third generation particles have writhe = 2. But since the more physical “winding number” equals (writhe + 1), we use winding numbers $`\omega `$ = 1,2,3 as the generation quantum number.

These features are embodied in the model depictions of Fig. 2. At the top of the figure the one-node lepton configurations are shown for the three generations. Intense orange and light orange colors (dark gray and light gray in monochrome) represent the minus and plus parts of the electric field respectively. The central ring is an event horizon making the interior fields unobservable. Hence a distant observer sees these as negatively charged particles when time averaged over many cycles. Bright yellow (lighter gray) and light yellow (very light gray) denote the north and south parts of the magnetic field. Because the twist inversion of the geometry offsets the cyclic inversion of the field, the north half of the field is always atop the ring, and the south half is always below. Thus a distant observer sees these particles as magnetic dipoles. The central part of Fig. 2 shows the two-node charged and neutral boson configurations. It is an experimental fact that the bosons do not show generational behavior. A possible reason for this will emerge from an interpretation of the mass relationship we will develop. At the bottom of Fig. 2 the models for the first generation u and d quarks are shown, and the quantum numbers for all the quarks are given. We do not show the field configurations for the higher generation quarks. With the higher winding numbers, the quark fields become so complex that the depictions are not very useful. The quark models have three degenerate forms that accommodate the color degree of freedom, but we show only one color form.

## 3 Known and estimated particle data

The masses of the leptons plotted in Fig. 1 have been experimentally measured and are very well known . Since the quarks have never been isolated, their masses must be deduced indirectly from many measurements as interpreted by QCD theory. For the lighter quarks, only a range of values can be specified. Collected in Table 1 are these mass values along with the relevant quantum numbers for the different particles. Where a range of mass is given, we use the average value. Also, in each generation of quark, the doublet of masses appears to be split from some parent level. At this stage we do not understand the splitting interaction, so we simplify by assuming that the observed mass levels come from a parent level half way in between. The masses for these unsplit levels are given in Table 2. For the bosons, no entry is given in the winding number column. This is because the bosons do not show generational behavior.

## 4 Extending the formula to
include the bosons and quarks

For particles with more than one node, we expect a contribution to the mass from “nodal” excitations. In order to see what this might be we have manipulated Eqn. 1.1, moving everything to the left hand side except the node dependence and plotted that against the number of nodes as in Fig. 3. For the W and Z bosons the generational term has been eliminated by taking the limit as $`\omega `$ goes to infinity (This has a physical basis that will be discussed later). The linear relationship with a slope of log $`\pi `$ as seen in Fig. 3 tells us that each node represents a factor of $`\pi `$ contribution to the particle mass equation. Knowing this lets us write an extended mass equation:

$$\mathrm{ln}\left(\frac{m\pi ^2}{\pi ^n}\right)=\frac{1}{2\pi \alpha }\left(\frac{1}{2}\right)^\omega +\mathrm{ln}\left(m_W\right)$$
(2)
where the symbols are as in Eqn. 1.1 and n is the particle node number taken from Tables 1 and 2. For the bosons, we set $`\omega =\mathrm{}`$ which removes the first term on the rhs of Eqn. (4.1) for them. This equation is plotted as the line in Fig. 4, along with points for all the elementary particles using the data from Tables 1 and 2.

As expected from Fig. 1, the experimentally known masses of the leptons and bosons fit the line given by Eqn. 4.1 quite well. The currently accepted masses of the quarks , with the doublets combined as described show more scatter about the line. Yet they are reasonably well represented given the approximations made. Thus we take this mass relationship as a valid equation for the masses of the leptons, bosons and quarks. Assuming that, the natural question is: “What features and interactions of the particles could give rise to such an equation?”

## 5 Physical interpretation of the mass equation

As a start, we exponentiate Eqn. 4.1, and rearrange it into a form that is easier to interpret:

$$m=m_W\frac{1}{\pi ^2}\pi ^n\mathrm{exp}\left[\frac{1}{2\pi \alpha }\left(\frac{1}{2}\right)^\omega \right]$$
(3)
The explicit multiplication dots are given to emphasize an important aspect; the different contributions to the mass are multiplicative, not additive. This can be understood if Eqn. 5.1 is a partition function equation. In statistical mechanics the total partition function for a system having independent parts is the product of the partition functions for the parts. This also shows up in topological quantum field theories where the vacuum expectation value of the field is given by the product of Wilson loop holonomies. Indeed, in those theories the (vev) is called a partition function. From this point of view, we can think of Eqn. 5.1 as giving the partitioning of the mass-energy of the particle among the degrees of freedom open to it. The $`m_W`$ factor sets the mass scale for all the particles. The $`\pi ^2`$ factor should be considered a normalization factor. Since the W boson has two nodes (n=2), this factor cancels the next factor when the equation is applied to the W boson. The $`\pi ^n`$ factor accounts for the level of excitation of the field harmonic oscillators in terms of the node number. The last factor, the exponential, is different in several respects. It may be thought of as an “embedding factor”, giving the fraction of the intrinsic particle mass that is observable to a real observer far away. This is the factor that distinguishes the masses of the different generations according to their winding numbers. For the bosons we allow the winding number to go to infinity, and this factor becomes unity. This sets the fermions apart from the bosons in a way discussed below.

One of the surprising features of Eqn. 5.1 is its simplicity. When one considers the mass of an elementary particle, one expects contributions from many sources. For charged particles there is the energy stored in the electromagnetic field. There was a time when the entire mass of the electron was thought to reside in its electromagnetic field. We now know that this contribution to the mass is very small compared to that arising from other sources. Then there are radiative corrections and possible quantum corrections that may contribute to particle masses. They seem not to be included in Eqn. 5.1. Lastly, contributions from internal excitations should contribute to the particle mass. These are included as the $`\pi ^n`$ factor in our equation. But comparing the relative contributions, it is clear that this one is minor to that specified by the $`m_W`$ factor. Here is a clue as to why this mass equation is so simple. There is one overwhelming and dominant effect that makes all the other contributions negligible. What is this dominant effect? We believe it is the need to provide energy to embed the rotating particle in the surrounding spacetime. Equation 5.1 is a “winding number” equation in the following sense. According to our model, these particles are quantum black holes of Kerr geometry that are rotating at v = c. From studies of Kerr black holes applied to rotating stars, it is understood that near the event horizon, the whole of spacetime is set into rotation. It is like a vortex in fluid mechanics in which the surrounding fluid has a rotating velocity field surrounding the vortex core. When a stable particle is formed, it must have enough energy to set up the vortex field in the surrounding spacetime. This is the dominant energy term for these particles that Equation 5.1 describes, and why $`m_W`$ is so large. Yet somehow the fermions seem to evade the bulk of this energy cost. How? It is because they are spinor particles. The rotation of spin 1/2 particles is very special. When a spin 1/2 particle is rotated about any axis for $`2\pi `$ radians, it winds up its connections to the surroundings just like a boson. But when the spinor is rotated another $`2\pi `$ radians, it can unwind the connection to the surroundings by suitable writhing. This possibility was known to Dirac, and is well described in the book of Misner, Thorne, and Wheeler. Applied to our particle models we have this picture. For the bosons, as time passes, their fields continuously wind around, frame-dragging nearby spacetime with them. This is the origin of the infinite winding number for the bosons. For the fermion particles, as their fields rotate around, they cyclically wind the surrounding spacetime on the first 2$`\pi `$ revolution and then unwind it on the second 2$`\pi `$ revolution. The result is that the fermions have to spend very much less energy embedding themselves into their surroundings than do the bosons. One can see this effect at work in the exponential factor of Eqn. 5.1. The internal winding number $`\omega `$ specifies how many times the field goes around before the node is reached (the half cycle point) where the spinor can begin unwinding its connection. This translates into the fraction of the energy for infinite winding that the particle must have. If this picture is correct, it implies that all the massive fermions should be much less massive than any of the massive bosons. This is in fact true excepting the top quark, as can be seen by looking at Table 1. From our mass correlation, it appears that the top quark is a special case and has a higher mass than expected. That is the why the tb quark point lies above the line in both Fig. 3 and 4. One possible reason for this is that some new effect not included in our considerations (Supersymmetry mixing?) is beginning to show up at the highest winding numbers and node numbers. Perhaps the same effect keeps fourth and higher generation families from having observable particles.

## 6 Discussion

The particle mass equation 5.1 describes mainly how the various particles embed themselves into spacetime. The mass effects connected with that dominate all the other mass effects that one usually thinks about. It is surprising that the mass equation contains no adjustable parameters if one considers $`m_W`$ and $`\alpha `$ as constants of nature. Clearly, much additional work needs to be done to fully appreciate what this mass equation is telling us.


text:
# Infrared Space Observatories: How to Mitigate Zodiacal Dust Intereference

## 1 Introduction

An unprecedented interest in space-borne observation missions has arisen out of NASA’s Origins Program, a collection of missions aimed at determining and characterizing the origin and development of galaxies, stars, planets, and the chemical conditions necessary to support extraterrestrial life444http://origins.jpl.nasa.gov. Over the course of the next two decades, the Origins Program will focus on developing space-based observatories through a series of precursor missions such as the Space Infrared Telescope Facility (SIRTF); first generation missions, such as Space Technology 3 (ST3), Space Interferometry Mission (SIM) and the Next Generation Space Telescope (NGST); and second generation missions, such as the Terrestrial Planet Finder (TPF).

One of the most important constraints of such missions is the so-called zodiacal dust (zodi) or interplanetary dust. In our own solar system, dust (very tiny solid particles composed of silicates, ices, and other minerals) is ever present. This material is both a remnant of the formation of the planets and a consequence of continuing collisions among comets, asteroids, and other small bodies. The zodi has a potentially serious impact on the ability of space-borne observatories to detect and study their targets. Specifically, the zodi reveals itself as a diffuse component of the sky brightness, attributed to the scattering of sunlight in the UV, optical, and near-IR, and the thermal re-radiation of absorbed energy in the mid-IR and far-IR. At infrared wavelengths from approximately 1 $`\mu `$m, the signal from the zodi is a major contributor to the diffuse sky brightness and dominates the mid-IR (10–60 $`\mu `$m) sky in nearly all directions (except very low galactic latitudes). Consequently, the emission from the local zodi is a major noise source, considerably affecting the size and weight requirements of the mirrors. In Beichman, Woolf, & Lindensmith (1999) it is shown that for a 1 AU TPF mid-IR interferometer mission, the local zodiacal background constitutes roughly 70% of the total noise. The intensity of the zodi radiation (in terms of thermal emission) is greatly reduced when moving in a direction normal to the ecliptic plane.

The zodiacal dust problem is the principal motivation for the design of trajectories having considerable displacements normal to the ecliptic. We have examined a widely used zodi model that was measured by the Cosmic Background Explorer (COBE) and the Infrared Astronomical Satellite (IRAS) (Kelsall et al. 1998). Simulation of this model quantifies the evolution of the zodi brightness (from thermal emission only) as a function of the height above the ecliptic plane, at various Earth positions. The results are depicted in Figure 1, with $`\lambda `$ denoting Earth’s angular position relative to the Sun. It is apparent that at 0.18 AU above the ecliptic plane more than 50% of the emission is avoided, and at 0.4 AU above the ecliptic, more than 80% of the emission is avoided. While higher excursions reduce further the noise generated by the zodi emission, the energetic requirements involved may be substantial. In other words, the 0.4 AU point on the abscissa of Figure 1 represents a threshold above which the dependence of the normalized brightness on the height above the ecliptic plane is small when compared to the energy cost to get there.

A variety of orbits have been considered for other missions. The diversity of these missions and their stringent resolution and accuracy requirements necessitates the design of specialized trajectories. For example, for the NGST mission, a halo orbit around the collinear Lagrangian libration point has been considered (Lubow 2000). For the TPF and Darwin missions, orbits at 5 AU where the zodi is highly reduced are being considered (Landgraf & Jehn 2001) but cost, communication, power, and delay to data return can be problematic. For the TPF mission, both $`L_2`$ and heliocentric Earth-trailing orbits have also been examined (Beichman et al. 1999). To date, several types of trajectories associated with the collinear libration points have been reported. These trajectories are obtained by exploring either the planar circular restricted three-body problem (CR3BP) or the spatial CR3BP. The so-called Lyapunov orbits obtained in the planar case (Szebehely 1967) bifurcate into a spatial family of Lissajous orbits (Gomez, Massdemont, & Simo 1997), which generate the well-known Halo orbits as a particular case only when perturbations are ignored and specific initial conditions are selected (Howell, Barden, & Lo. 1997; Henon 1974; Goudas 1963). While libration point orbits have been widely explored — as they often offer valuable operational and scientific features — less attention has been given to other trajectories emerging from the spatial CR3BP that are more suitable for the unique constraints of mid-IR space-borne observation missions.

This paper describes work aimed at synthesizing families of out-of-the-ecliptic trajectories, which mitigate the effects of the zodi. In §2 we discuss the zodiacal dust model, in §3 we describe the optimal trajectories, and §4 follows with a summary and discussion.

## 2 The Local Zodiacal Cloud Model

The scattered light and thermal emission from the zodiacal dust cloud within our own solar system contribute significantly to the diffuse sky brightness observable at most infrared wavelengths. The zodi has been investigated extensively using observations from the IRAS satellite and the most recent detailed COBE Diffuse Infrared Background Experiment (DIRBE). A model of the zodi was used to evaluate the effectiveness of the orbits described in this paper. The parameterized physical model of the zodi that was used is a simplified version of the highly complex model developed in direct conjunction with the measurements from the COBE DIRBE satellite (Kelsall et al. 1998).

As the following equation and discussion demonstrate, the model brightness is calculated as the integral of the product of a source function and a three-dimensional dust density distribution function evaluated along the line of sight,

$$I_\nu =E_{c,\nu }B_\nu (T(R))n_c(X,Y,Z)𝑑s.$$
(1)
At infrared wavelengths of approximately 10-60$`\mu `$m, thermal reradiation of absorbed energy from the zodi is estimated to contribute 90% or more to the total sky brightness in nearly all directions (except at very low Galactic latitudes). In contrast, the scattered sunlight contribution is of significance in the UV, optical, and near-IR. Therefore, as a first approximation, it can be neglected in the case of an observing mission centered at a 12 $`\mu `$m wavelength. The simplified model thus considers the thermal emission contribution of the zodi, which is expressed in equation (1) as a blackbody $`B_\nu (T(R))`$, and its associated emissivity modification factor at 12 $`\mu `$m $`E_{c,\nu }`$ (initially a free parameter in the original DIRBE model) which is $`E_{c,12\mu \mathrm{m}}=0.958`$. The dust grain temperature $`T`$ variation with distance from the Sun is given by $`T(r)=T_0R^\delta `$, where $`\delta =0.467`$. The three dimensional dust density distribution, $`n_c(X,Y,Z)`$, is composed of several structured components: a smooth cloud, three asteroidal dust bands and a circumsolar ring at 1AU. At a wavelength of 12 $`\mu `$m, the DIRBE satellite (Kelsall et al. 1998) measured the zodi components of the smooth cloud, dust bands and circumsolar ring as 28.476, 1.938, 3.324 MJy sr-1 respectively. The simplified model adopted in this paper neglects the contribution to the zodi brightness of the dust bands and circumsolar ring.

Because of the relatively small inclination of the midplane of the smooth cloud distribution, the zodical light is the only component of the sky brightness that is not fixed on the celestial sphere. This important and unique feature, depicted in Figure 2, results in the temporal variation of zodi brightness observed in a given celestial direction by an Earth-based observer. The model calculations for a spacecraft are performed in heliocentric ecliptic coordinates $`(X,Y,Z)`$, where $`s,R_{},\lambda `$ represent respectively the height above the earth, the Earth-Sun distance (1 AU) and the heliocentric longitude of the earth.

$`X=R_{}\mathrm{cos}\lambda `$ (2)
$`Y=R_{}\mathrm{sin}\lambda `$ (3)
$`Z=s`$ (4)
The center of the smooth cloud is offset from the Sun by $`(X_0,Y_0,Z_0)`$ and the resulting translated cloud coordinates are:

$`X^{}=XX_0`$ (5)
$`Y^{}=YY_0`$ (6)
$`Z^{}=ZZ_0`$ (7)
$`R_c=\sqrt{X^2+Y^2+Z^2}`$ (8)
The vertical structure of the smooth dust cloud is determined by the height of its inclined symmetric midplane;

$$Z_c=X^{}\mathrm{sin}\mathrm{\Omega }\mathrm{sin}iY^{}\mathrm{cos}\mathrm{\Omega }\mathrm{sin}i+Z^{}\mathrm{cos}i$$
(9)
where $`i`$ and $`\mathrm{\Omega }`$ are the inclination and ascending node of the midplane.

The density of the smooth cloud is separable into radial and vertical terms;

$$n_c(X,Y,Z)=n_0R_c^\alpha f\left(\left|\frac{Z_c}{R_c}\right|\right),$$
(10)
where

$$f\left(\left|\frac{Z_c}{R_c}\right|\right)=\{\begin{array}{cc}e^{\beta \left(\frac{|Z_c/R_c|^2}{2\mu }\right)^r}\hfill & \mathrm{for}\left|\frac{Z_c}{R_c}\right|<\mu \hfill \\ e^{\beta \left(|Z_c/R_c|^2\mu /2\right)^r}\hfill & \mathrm{for}\left|\frac{Z_c}{R_c}\right|\mu \hfill \end{array}$$
(11)
Following the same procedure of the DIRBE model, integration along the line of sight was performed from the satellite to an outer radial cutoff of 5.2 AU from the Sun. A particular line of sight is defined by two angles, the projected angle in the ecliptic plane and the angle from the ecliptic. To evaluate the zodi at a particular point along the orbit a single angle from the ecliptic (usually 30 degrees) was used and the zodi brightness was averaged over 360 degrees of the projected angle, these numbers were then normalized to go from 0 to 1. Our simulations found that the normalized zodi brightness were very close for angles above the ecliptic in the range of 0 to 60 degrees. Note that the spacecraft observes away from the ecliptic plane, so it rotates 180 degrees when it is crosses the ecliptic.

## 3 Trajectories

### 3.1 The Equations of Motion

A vast amount of literature exists on the CR3BP (e.g. Gomez, Masdemont, & Simo 1998; Gomez et al. 1997; Howell 1997; Marchal 1990; Szebehely 1967 and references therein). In almost all past research, the standard rotating coordinate system has been used with the origin set at the barycenter of the large primary, $`m_1`$, and the small primary, $`m_2`$. The $`x`$-axis is positive in the direction of $`m_1`$, the $`z`$-axis is perpendicular to the plane of rotation and is positive when pointing upward, and the $`y`$-axis completes the set to yield a right-hand reference frame. Normalization is performed by setting $`m_2=\mu ,m_1=1\mu `$ where $`\mu =m_2/(m_1+m_2)`$. Thus, $`m_1`$ is located at $`(\mu ,0,0)`$ and $`m_2`$ is located at $`(\mu 1,0,0)`$. Usually, this coordinate system is sufficient to model the problem and yields a fruitful characterization of diverse families of trajectories. However, in this study we have adopted a slightly different rotating coordinate system, which we found to be particularly useful for characterization of trajectories that reduce the zodi interference. The basic notion is to choose the origin of the coordinate system at the center of the small primary, Earth in our case — rather than the barycenter — and to normalize the masses by the mass of the large primary, the Sun in our case. This coordinate system was first used by Rabe (1961) and later by Breakwell (1963) in his approximate analysis of three-dimensional trajectories.

Following the usual analysis but in this variant coordinate system, let $`𝒓`$ denote the position vector of the vehicle relative to Earth, and $`𝑹`$ that of the earth relative to the Sun. The acceleration of the spacecraft relative to earth is given by Breakwell (1963):

$$\ddot{𝒓}=\frac{\mu _E𝒓}{r^3}[\frac{\mu _S(𝑹+𝒓)}{|𝑹+𝒓|^3},\frac{\mu _S𝑹}{R^3}]$$
(12)
where $`\mu _E`$ is the gravitational constant of the earth and $`\mu _S`$ is the gravitational constant of the Sun.

The acceleration $`\ddot{𝒓}`$ is evaluated in a rotating Earth-fixed coordinate system as depicted in Figure 3. This local vertical reference frame originates at the center of the earth, with the $`x`$-axis directed radially outward along the local vertical, the $`y`$-axis lying along the direction of Earth’s motion, and the $`z`$-axis normal to the ecliptic to complete the Cartesian right-hand setup.

We shall use the following unit convention: The position of the vehicle is measured in astronomical units (AU), where the mean Earth-Sun distance, assumed constant, is $`|𝑹|=R=1`$ AU$`=1.496\times 10^8`$ km. The time unit is normalized by the earth’s mean heliocentric angular velocity, i.e. $`t=t^{}/\sqrt{R^3/\mu _S}`$ with $`t^{}`$ as the time measured in seconds. Accordingly, the velocity vector of the vehicle $`𝒗=[\dot{x},\dot{y},\dot{z}]^T`$ is normalized by $`R/\sqrt{R^3/\mu _S}`$. Also, let $`\mu =\mu _E/\mu _S`$.

The inertial frame used here, denoted by $`XYZ`$, is a heliocentric-ecliptic coordinate system, as depicted in Figure 3. By utilizing the above unit convention, the transformation from the rotating frame to the inertial frame can be found,

$$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{ccc}\mathrm{cos}t& \mathrm{sin}t& 0\\ \mathrm{sin}t& \mathrm{cos}t& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x+1\\ y\\ z\end{array}\right].$$
(13)
Also, we define a pseudo-potential function,

$$\mathrm{\Omega }(x,y,z)=\frac{1}{2}(x^2+y^2)+\frac{\mu }{r}+\frac{1}{\rho }+x,$$
(14)
where

$$r=\sqrt{x^2+y^2+z^2}$$
(15)
is the distance from Earth, and

$$\rho =\sqrt{(x+1)^2+y^2+z^2}$$
(16)
is the distance from the Sun.

To proceed, we assume that the vehicle performs a ballistic motion, i.e., no control forces are used, and that the motion of the vehicle is unperturbed. Computing the well-known expression for the acceleration in a rotating coordinate frame yields the following equations of motion:

$`\ddot{x}2\dot{y}=\mathrm{\Omega }_x`$ (17)
$`\ddot{y}+2\dot{x}=\mathrm{\Omega }_y`$ (18)
$`\ddot{z}=\mathrm{\Omega }_z.`$ (19)
where the subscript stands for partial differentiation. Substituting from equation (14), equations (17)–(19) can be equivalently written as

$`\ddot{x}x2\dot{y}={\displaystyle \frac{\mu x}{(x^2+y^2+z^2)^{3/2}}}{\displaystyle \frac{(1+x)}{[(x+1)^2+y^2+z^2]^{3/2}}}+1`$ (20)
$`\ddot{y}y+2\dot{x}={\displaystyle \frac{\mu y}{(x^2+y^2+z^2)^{3/2}}}{\displaystyle \frac{y}{[(x+1)^2+y^2+z^2]^{3/2}}}`$ (21)
$`\ddot{z}={\displaystyle \frac{\mu z}{(x^2+y^2+z^2)^{3/2}}}{\displaystyle \frac{z}{[(x+1)^2+y^2+z^2]^{3/2}}}`$ (22)
or, in state-space vector form,

$$\dot{𝐱}=𝐟(𝐱),$$
(23)
where

$$𝐱[x,\dot{x},y,\dot{y},z,\dot{z}]^T.$$
(24)
Multiplying equations (20)–(22) by $`\dot{x},\dot{y},\dot{z}`$, respectively, adding the results and performing the integration, renders the integral of motion, the well-known Jacobi constant, given by:

$$C=\mathrm{\Omega }(x,y,z)\frac{1}{2}(\dot{x}^2+\dot{y}^2+\dot{z}^2).$$
(25)
Note that for any given initial conditions, the differential equations (20)–(22) are equivariant under the following transformations:

$$(x,y,z,t)(x,y,z,t)\mathrm{or}(x,y,z,t)$$
(26)
Hence, symmetries exist only for a backward integration of the equations of motion. However, various symmetries exist for different sets of initial conditions. An obvious symmetry that occurs in the $`xz`$ plane is

$$x(t,(𝐱_0)_1)=x(t,(𝐱_0)_2),y(t,(𝐱_0)_1)=y(t,(𝐱_0)_2),z(t,(𝐱_0)_1)=z(t,(𝐱_0)_2),$$
(27)
where

$$(𝐱_0)_1=[x_0,\dot{x_0},y_0,\dot{y_0},z_0,\dot{z_0}]^T,(𝐱_0)_2=[x_0,\dot{x_0},y_0,\dot{y_0},z_0,\dot{z_0}]^T.$$
(28)
Other symmetries that may exist as functions of initial conditions are much less obvious, due to the rotation of the coordinate system. The next section describes some stability considerations, which pave the way to synthesizing the constraints for the trajectory optimization procedure.

### 3.2 Practical Stability

One of the most important aspects of designing a trajectory is to define, study and characterize its stability. Since stability is a matter of definition, one has to use an appropriate stability framework that is suitable to the specific dynamical system involved. We have found that the so-called practical stability theory renders a good means to quantify our analysis and design of the CR3BP trajectories.

For given initial conditions, a trajectory is said to be practically stable if $`r\beta =r_{max}`$ for $`t[t_1,t_2]`$, i.e., if the vehicle has not gone beyond some pre-specified distance from Earth at a given time interval. This stability definition is used to constrain the Genetic Algorithm optimization procedure (see §3.3 below). It reflects operational considerations, since at long distances from Earth communication becomes too costly.

Regions of practically stable motion for the spatial CR3BP can be semi-analytically determined by examining the curves of zero velocity, also known as Hill’s regions. These curves are generated by plotting equipotential contours of the function (14), or, equivalently, plotting the position of the vehicle for different values of the Jacobi constant evaluated along $`\dot{x}=\dot{y}=\dot{z}=0`$. The permissible regions of the vehicle’s motion are confined by different values of $`C`$. For the three-dimensional case discussed here, sufficiency conditions for practically stable motion can be determined as follows (Szebehely 1967): Let $`C(L_1,0)`$ denote the Jacobi constant for the Sun-Earth system evaluated along the zero velocity curve that intersects the Lagrangian equilibrium point $`L_1`$. If

$$C>C(L_1,0)=1.0004513$$
(29)
then the zones of possible motion of the vehicle are divided into three disconnected parts, either near the earth, near the Sun, or far away from both. The first two regions are practically stable and are known as Hill-stable regions (Marchal 1990). To distinguish between the practically stable regions and the third (not necessarily practically stable) region, we notice that the vehicle can drift away from the primaries only if the gradient of its potential function is positive. Hence, if equation (29) holds and in addition

$$\mathrm{\Omega }(x,y,z)<0$$
(30)
a practically stable motion results.

Unfortunately, in order to satisfy both equations (29) and (30), one has to choose impractical initial conditions, since, generally speaking, orbits around the earth have a small normal motion magnitude, and trajectories that start near the Sun, which potentially have large normal deflections, can impose undesirable practical constraints on the mission. Instead, we look for different values of $`C`$ that also result in practically stable trajectories (recall that the conditions above for stability are sufficient but not necessary). The basic requirement is to find initial conditions that lie within reasonable distance from Earth, but shift the vehicle to some (possibly) heliocentric orbit inclined to the ecliptic with the least possible energy. This task is fulfilled in the next section.

### 3.3 Optimal Trajectories

The discussion in the previous section stresses the complexity and counter-intuitive nature of solutions to the spatial CR3BP. We used a Genetic Algorithm (GA), specifically the Deterministic Crowding GA (Mahfoud 1995), to optimize the trajectories. This method is preferable over other optimization algorithms such as the gradient search and the simplex method because it avoids local minima (i.e. the search is performed over the entire state space) and promotes diversity of solutions.

We wish to maximize the normal displacement subject to the following

$``$ the differential equations of motion;

$``$ the practical stability constraint — the vehicle should not exceed a given distance from Earth during a given time interval;

$``$ the initial position vector should lie outside Earth’s sphere but inside some pre-determined radius;

$``$ the initial velocity vector should not exceed some pre-specified limit, derived from the overall propellant mass and capability of the launch vehicle.

See Gurfil & Kasdin (2001) for the GA parameters used for the trajectory optimizations. The trajectory search procedure is divided into two conceptual stages: characterization of optimal trajectories, where the search space is confined by rather loose bounds; and design, where the search space is narrowed in order to obtain practical results.

### 3.4 Trajectory Characterization

Eight optimization sets were carried out with different upper and lower bounds on the initial conditions and the maximum permissible distance from Earth (see Gurfil & Kasdin 2001). In order to visualize the trajectories obtained, the individual with the highest fitness in the last generation is selected as the optimal solution. The time scale selected for visualization purposes is 15 years (larger than the 5 year time scale selected in the optimization itself).

The upper left panel in Figure 4 depicts the three-dimensional trajectory that resulted from the one of the optimization sets. Note that in this case the GA generated a quasiperiodic Lissajous trajectory that spirals above the ecliptic plane. In other words, if $`\omega _z`$ denotes the frequency of the vertical motion, and $`\omega _x,\omega _y`$ denote characteristic frequencies of the radial and transverse motions, respectively, then this trajectory satisfies

$$\omega _z<\omega _x\omega _y$$
(31)
which results in a large fraction of the orbit being above the ecliptic. Trajectories satisfying equation (31) are categorized as Type I trajectories. The main deficiency of Type I trajectories is their considerable distance from Earth. Note also that the solution depicted in the upper left panel of Figure 4 is asymmetric relative to the ecliptic plane, since most of the time the vehicle remains above the ecliptic. However, due to the symmetry property (27), a mirror image of the trajectory relative to the ecliptic can easily be generated. This is true for all the trajectories considered.

The upper right panel in Figure 4 depicts the three-dimensional trajectory that resulted from another optimization set. In this case the GA also generated a quasi-periodic Lissajous trajectory. However, here the trajectory satisfies555$`\omega _x,\omega _y`$ in this case denote the frequencies of the dominant harmonics.

$$\omega _z>\omega _x\omega _y$$
(32)
which results in frequent vertical crossings of the ecliptic while maintaining reasonable distance from Earth. Trajectories satisfying equation (32) are categorized as Type II trajectories. The main deficiency of Type II trajectories is that frequent ecliptic crossing reduces the overall duration above a certain height. Note that as with Type I, the vertical crossing is asymmetric.

A third type of trajectory, shown in the lower left panel of Figure 4, was spotted when performing the characterization on yet a third optimization set. It satisfies the condition

$$\omega _z\omega _x\omega _y.$$
(33)
Thus, the trajectory is actually a quasi-periodic Lissajous trajectory as well, but it is almost closed in three dimensions, due to equation (33), as illustrated by the lower left panel of Figure 4.

A further examination shows that the other optimization sets produced Type I, II, and III trajectories as well. Generally speaking, reducing the practical stability limit reduces the frequency of the in-plane motion while keeping the frequency of the vertical motion almost unchanged. Reducing the initial velocity limit reduces the amplitude of both the vertical and the in-plane motion. The next step is to narrow the search space in order to design practical operational trajectories for space-borne observation missions. This step is carried out using GAs as well.

### 3.5 Trajectory Design: Low-Energy Optimal Trajectory

The design procedure is different from the characterization process, as the constraints of the problem are determined based upon practical engineering considerations. The major limitation of a real-life trajectory is the maximum energy available from the launcher for vehicle injection. We used twice the energy per unit mass required to inject a vehicle into an Earth departure hyperbola starting from a circular orbit of radius $`r_0`$ as our design constraint (this measure is known as C3).

To simplify the discussion, we shall describe only the first stage of the iterative design procedure. To this end, we chose to start from a 200 km parking orbit, i.e. $`(r_0)_{max}=6578`$ km. Next, an estimate of the mass should be made. We chose the benchmark value 4000 kg, which is roughly the estimated mass of the current TPF configuration (Beichman et al. 1999). Using charts of lift capabilities versus C3 for various launchers, it was found that the Atlas ARS or the Delta IV could provide a C3 of approximately 40 given this mass. This C3 is smaller than the C3 required for generation of orbits normal to the ecliptic using planetary flybys (Buglia 1973) or direct injection (Renard 1970). Consequently, assuming C3 = 40 km2/s2 and $`(r_0)_{max}=6578`$ km, we have $`(\nu _0)_{max}=12.7`$ km/s. Also, we let $`r_{max}=2`$ AU.

Using the constraints mentioned, the GA was used to generate a population of optimal trajectories. The last generation of the GA optimization comprises several solutions, but two dominant solutions stand out (due to the symmetry property). Because the mid-plane of the zodiacal cloud is inclined 2.03 degrees to the ecliptic, with the ascending node at 77.7 degrees, going above the ecliptic renders slightly better reduction in zodi noise than going below the ecliptic, for the epoch used here. Hence, the initial conditions giving a positive normal displacement were chosen.

The optimization procedure will always pick $`\nu _0`$ and C3 to lie on the constraint surface. Thus, properties of the low-energy optimal trajectories are as follows: $`\nu _0=12.7`$ km/s, C3=40 km2/s2 and accordingly $`\mathrm{\Delta }\nu =4.9`$ km/s. The trajectory starts from a 200 km parking orbit. The maximum normal deflection above the ecliptic is 0.223 AU with a maximum distance of 2 AU from Earth. Although the resulting optimal trajectory is non-Keplerian, examining it in the inertial reference frame (13) reveals that it is a heliocentric orbit with the distance from the Sun satisfying $`0.984\rho 1.124`$, assuming a mission lifetime of 5 years. This property is most important when considering solar arrays and power management for the mission, since it implies that modest-sized solar arrays can be used. Comparing the frequencies of the in-plane motion to the frequency of the normal motion classifies this trajectory as Type II.

Figure 5 depicts the normal, transverse, and radial displacements in the rotating Earth-fixed coordinate system, and the distance from Earth, for a 5 year mission for the low-energy optimal trajectory. The three-dimensional trajectory is presented as well. Because of the limited time scale, it seems as though the trajectory drifts away from Earth; however, the transverse and radial displacements are periodic having very slow frequencies relative to the frequency of the normal motion, hence a Type II trajectory results.

In order to estimate how much the proposed trajectory reduces the zodi noise, we used the zodi model described in §2 and incorporated it into the integration of the differential equations (20)–(22) using an appropriate coordinate transformation. First, note the periodic behavior resulting from the periodic normal displacement; in fact, these two variables are coherent. Second, note that for a positive normal displacement, the reduction in zodi brightness is 5% higher than for negative normal displacement. The maximum reduction in brightness is 67%. Averaging the values over time yields a mean reduction of 45%. Practically speaking, during 60% of the mission lifetime the zodi brightness is reduced by more than 50%. This allows for a considerable reduction in mirror size, permits faster data integration times and allows collection of more observations for a given mission lifetime. Figure 7 shows a plot of the trajectory in heliocentric coordinates, where the intensity of the plot represents the normalized zodi brightness. Figure 8 is the cumulative brightness distribution for this trajectory. This represents the percentage of the mission life time for which the normalized zodi brightness is below a certain value.

Using the methodology proposed in Beichman et al. (1999) for the calculation of mirror diameter, we have found that for the TPF IR interferometer the noise decrease allows a reduction of 20% in mirror diameter, which yields up to 35% reduction in payload mass. Note that this reduction is not included in the optimization. Iteration on the optimization scheme would result in a reduced energy requirement for this reduced mass and perhaps a higher normal displacement.

### 3.6 Trajectory Design: High-Energy Optimal Trajectory

The previous section described an optimal trajectory that emerges from a 200 km parking orbit. The constraints $`r_{max}2`$ AU and $`(\nu _0)_{max}=12.7`$ km/s resulted in a Type II trajectory that yielded a significant reduction in the zodi brightness. The purpose of this section is to present a different optimal trajectory which was obtained using other constraints. These constraints represent the maximum lift capability of existing launch vehicles and may even exceed them to some extent. Nevertheless, we have chosen to present this high-energy optimal trajectory because we believe that the energetic requirements for its implementation can be achieved by future launch vehicle such as the EELV or possibly via a combination of low-thrust electric propulsion and impulsive velocity changes. The principal merit of this trajectory is the outstanding reduction in zodi noise for most of the mission lifetime.

In this new trajectory, we have chosen to start from a 36,000 km Geosynchronous orbit (GEO), i.e. $`(r_0)_{max}=42,378`$ km. Anticipating that the reduction in zodi noise allows the mass to be reduced by half666This is the maximum mass reduction possible for TPF. It is obtained when the zodi noise is completely eliminated. Other noise sources, however, prevent further mirror size and mass reduction (Beichman et al. 1999)., we use a 2000 kg payload. With this mass, we chose C3=95 km2/s2, which exceeds — but not by much — the lift capability of the Atlas ARS. This dictates the constraint $`(\nu _0)_{max}`$. Also, we let $`r_{max}=3`$ AU. Again, the GA optimizer was used to synthesize an optimal trajectory.

Properties of the optimal trajectory are as follows: The trajectory starts from a 36,000 km orbit. The maximum normal deflection above the ecliptic is 0.374 AU with a maximum distance of 3 AU from Earth. As in the previous case, the resulting optimal trajectory drifts away from Earth for a given mission lifetime of 5 years. It is a non-Keplerian heliocentric orbit satisfying $`0.958\rho 2.274`$. Comparing the frequencies of the in-plane motion to the frequency of the normal motion classifies this trajectory as Type III.

The optimal trajectory is visualized in Figure 9, which exhibits some of the extraordinary features of this trajectory. Figure 9 depicts the normal, transverse and radial displacements in the rotating Earth-fixed coordinate system, the 3D trajectory, and the distance from Earth, for a 5-year mission. Note that the minimum approach to Earth is 0.2 AU and is obtained after 802 days (approximately 2.2 years). The relatively close approach (0.2 AU) means that the orbit might be able to be perturbed to come close enough to Earth for replenishment or maintenance. Thus, this Type III trajectory offers somewhat of an advantage over the previously discussed optimal Type II trajectory. This specialized feature is typical of Type III trajectories.

Similarly to the treatment in the previous section, we used the zodi model described in §2 to estimate the reduction in zodi brightness. Figure 10 describes the normalized zodi brightness as a function of time. The results are dramatic: the maximum reduction in brightness is 97%. Averaging the values over time yields a mean reduction of 75%. During 82% of the mission lifetime the zodi brightness is reduced by more than 70%. Figure 11 shows a plot of the trajectory in heliocentric coordinates, the intensity of the plot represents the normalized zodi brightness. Note that unlike in Figure 7, the zodi brightness is asymmetric; this is because one of the ecliptic crossings is farther away from the Sun than the other so the zodi brightness is lower. Figure 12 is the cumulative brightness distribution for this trajectory. This represents the percentage of the mission life time for which the normalized zodi brightness is below a certain value. The zodi reduction for this orbit allows for a considerable reduction in telescope mirror size, which not only results in a remarkable reduction in mass, but moreover, can significantly reduce the development and manufacturing costs.

Using the methodology proposed in Beichman et al. (1999) for the calculation of mirror diameter, we have found that for the TPF IR interferometer, the noise reduction due to the high-energy trajectory allows a reduction of 36% in mirror diameter, which yields up to 50% reduction in payload mass.

## 4 Summary and Discussion

The zodiacal dust is a significant problem for mid-IR missions like TPF because of its strong thermal emission. While this emission can be mitigated by a mission orbit at 5 AU from the Sun (as proposed by Landgraf & Jehn 2001), the associated design aspects renders such an orbit infeasible even though the zodi is considerably reduced. For a finite mission lifetime, the overall scientific return is dramatically reduced due to the time required to get to a 5AU orbit. Moreover, these trajectories require very large solar arrays and massive communication antennas. Also, in case of a failure, it is practically impossible to replace the damaged components or to perform other replenishment missions. Note also that going a unit distance above the ecliptic decreases the zodi intensity considerably more than going a unit distance from the Sun within the ecliptic. This was the motivation for the study of out-of-the-ecliptic trajectories for space-borne observatories. The spacecraft observes away from the ecliptic plane, so that it turns over when it crosses the ecliptic plane. We have used genetic algorithms as the parameter optimization procedure which has resulted in a fruitful probing of the complex dynamics of the restricted problem. For example, the farthest either of our trajectories goes from Earth is 2.2 AU within a 5 year mission lifetime, which means an order of magnitude smaller solar arrays can be used than for a mission at 5 AU. Also, data collection and scientific interpretation can begin right away instead of having to wait several years. A direct trip to 5 AU requires a C3 of 80 km2/s2 (Meissinger, Wertz, & Dawson 1997), equivalent to our high energy trajectory, and a travel time of two years. Galileo’s trajectory used a C3 of 13 km2/s2 (Meissinger et al.), but it took 6 years to reach Jupiter, with most of that time spent around 1 AU, and required correct planetary alignment to allow for three planetary flybys. Going to 5 AU also means long travel times through the ecliptic plane, which increases the chance of micrometeorite damage. Finally, we did not iterate on mass reduction in our optimization procedure. The large mass reduction would provide a reduction in launch costs or alternatively an orbit with a greater normal displacement above the ecliptic.

Three types of out-of-ecliptic trajectories were identified using Deterministic Crowding GAs. Type I was characterized by a high frequency in-plane motion and a slow frequency normal displacement. This type was obtained for large departures from Earth ( 5AU). Type II was obtained when the maximal Earth-departure was restricted to small values ($`<`$ 2 AU). Its main feature was low frequency in-plane motion and high frequency out-of-the-plane motion. Type III was obtained for an intermediate Earth departure distance (3 AU) and exhibited almost equal characteristic frequencies. This classification was used as the infrastructure for design of optimal trajectories. It should be stressed that Type II trajectories can also achieve large normal displacements if an appropriate C3 is chosen. (Recall that C3 is twice the energy per unit mass which is required to inject a vehicle into an Earth departure hyperbola starting from a circular orbit of radius $`r_0`$.) However, Type III trajectories offer both longer times above the ecliptic and a closer approach to Earth than the other two types of trajectories. Type II trajectories provide a smaller departure from Earth and a shorter distance from the Sun. These facts should be traded off when considering utilization of either of these trajectories for real-life space-borne observation missions. The trajectories are selected a priori to be stable; hence, stationkeeping maneuvers are reduced to correcting perturbations only, and, furthermore, escape to an undesirable region is impossible.

The design process incorporated tighter parameter constraints than the characterization procedure. In order to obtain a Type II optimal trajectory, an estimate of the payload mass was made and an appropriate C3 was chosen, based upon lift capabilities of existing launchers. We call this the low-energy trajectory. The second optimal trajectory was a high-energy trajectory. Although its required C3 exceeded the capabilities of existing launchers, it might constitute a worthwhile alternative using future launchers or a combination of impulsive and electric propulsion. The two trajectories have the properties described in Table 1.

The most important point is that both trajectories imply that the overall collecting area and hence the mass and cost can be significantly reduced. This is due to the substantial reduction in zodi which results from the out-of-the-ecliptic displacement of 0.22 and 0.37 AU for the low and high energy trajectories respectively. For comparison, in typical libration point trajectories, the normal excursion is up to an order of magnitude smaller. This high displacement also considerably reduces the probability of micrometeorite or space debris damage, which could extend the mission lifetime. The low energy trajectory has a mean distance of 1 AU from the Sun which guarantees that standard solar arrays can be used for power supply. The high-energy trajectory, on the other hand, could require larger solar arrays, though still well within current capabilities with minimal increase in the mass. The drift from Earth is 2 AU for the low energy trajectory and so communication and uplink/downlink operations can be carried out using modest-sized antennas. The high energy trajectory has a drift from Earth of 3 AU and larger antennas would be required. Although the high-energy trajectory requires more communciation and operating power, and relies on future planned launch vehicles, the expected reduction in mass, development and manufacturing costs is likely to be large and should compensate for the expected increase in launch costs.

The mass reduction figures used in this design process are rather specific for a TPF mid-IR interferometer looking for terrestrial planets. As the local zodi cloud noise is mitigated, a noise floor is reached. This minimum noise level has specific characteristics derived from the mission itself, and is determined by other noise sources that become drivers as we escape the ecliptic. If the trajectories obtained are generalized to other missions, such as filled aperture telescopes like NGST and SIRTF, even further benefits in terms of mass reduction are possible. Benefits are not as dramatic for optical missions because the zodi intensity is lower in the visible wavelength region.

The quest for optimal trajectories for space-borne observatories has by no means been exhausted by this study. While the foundation for out-of-the-ecliptic trajectories has been built, there are several additional points that need to be clarified and further investigated. These include optimal transfers from Earth to the initial orbit, use of electric propulsion, searching for optimal initial conditions in other regions of state space, including perturbations from other solar system bodies, and using lunar and planetary flybys.

This work was performed as part of the Ball Aerospace architecture studies for the Terrestrial Planet Finder mission under JPL contract 1217281. The authors gratefully acknowledge the help and support of the people at Ball and JPL as well as the other members of the science team and the Princeton Terrestrial Planet Finder group. SS is supported by the W.M. Keck Foundation.


text:
# Nematic - Isotropic Transition in Porous Media - a Monte Carlo Study

## Abstract

We propose a lattice model to simulate the influence of the porous medium on the Nematic-Isotropic phase transition of liquid crystals confined to its pores. The effects of the pore surface and the pore interconnectivity are modelled through a disorder parameter. Monte Carlo calculations based on the model leads to results that compare well with experiments on the suppression of the transition temperature with decreasing pore size, the softening of the transition from the first order to continuous when the pore size becomes small and the enthalpy of transition.

Statistical mechanics of liquid crystals in confined geometry is a subject of interest both from technological and basic science points of view . Many electro-optical devices operate based on the properties of liquid crystals confined in suitable geometries. For example, incorporating liquid crystal in gel network has led to a new type of nematic display device with a significantly improved switching property . The academic interest in this subject stems from the fact that confinement leads to several complex, unexpected and interesting phenomena that are not very well understood. The basic issue concerns the phase behaviour of the confined liquid crystals which is determined by the delicate balance established amongst the competing effects of the disorder induced by the confining surface, the finite size, the orientational order and the temperature. A porous medium with a random network of interconnected pores of random sizes and shapes provides a simple confining medium . The sizes of the pores in the medium, pore surface roughness, pore connectivity etc., influence the nature and content of the phase and its transition.

In this letter we shall focus attention on the Nematic-Isotropic (N-I) phase transition. Experiments tell us that the transition temperature, denoted by $`T_{\mathrm{NI}}`$, is lower in the porous medium as compared to the bulk sample; also $`T_{\mathrm{NI}}`$ decreases with decreasing pore size . The N-I transition is weakly first order in the bulk. In the porous medium, however, the transition softens as the pore size decreases and becomes continuous for small pores .

Several models have been proposed and investigated to understand the above experimental findings. Our interest here is to study the orientational ordering of the liquid crystals; hence lattice models with nearest neighbour interaction would suffice. There are numerous approaches that have been tried out to model the influence of the porous matrix. These include explicit construction of an independent pore with different rigid boundary conditions and director configurations , Potts spin models with the porous medium approximated by the diffusion-limited cluster-cluster aggregate , random field Ising models , random anisotropy nematic models  and models based on dilution . An important objective of all these models has been to make a quantitative prediction of the depression of $`T_{\mathrm{NI}}`$ with decreasing pore size. The models proposed so far have been successful only qualitatively: the predicted decrease of $`T_{\mathrm{NI}}`$ with decreasing pore size is atleast an order of magnitude higher, see e.g. , than what has been observed experimentally. Such a large discrepancy has been attributed to several reasons. For example, in the independent pore models, the influence of the interconnectivity of the pores is not accounted for . In the dilution models, only the excluded volume effect of the pore surface is modelled. Despite this, the inter-particle interaction seems to be excessively perturbed in these models .

In this letter, we propose a model that takes into account the effect of the disorder, induced by the pore surface as well as pore connectivity. Consider an $`L\times L\times L`$ cubic lattice with a total of $`L^3`$ lattice sites. Each lattice site holds a liquid crystal molecule interacting only with its six nearest neighbours as per the angular part of the Berne - Pechuka (BP) potential . Periodic boundary conditions are imposed. For a given value of $`L`$ we select a certain number, $`N_d`$, of lattice sites randomly and declare the liquid crystal molecules on these sites to have quenched random orientations. This disorder parameter $`N_d`$ is key to our model. It is defined in such a way that it serves as a good measure of the confining effects of a single pore as well as the connectivity of the random network of pores, wide infra.

Consider a pore that contains $`N`$ liquid crystal molecules. Of these, on the average, $`N^{2/3}`$ molecules will be on the surface and these are the ones that are directly influenced by the pore surface. In a lattice model, the substrate site perturbs one out of the six nearest neighbour interactions (bonds). However, in our model, we propose to take care of the surface effect by placing molecules with quenched random orientations inside the lattice; each such quenched molecule has six bonds. Thus it is reasonable to say that the relevant number of randomly quenched molecules required to represent the surface effects of a single pore should be $`6^1N^{2/3}`$. Next, consider the cubic lattice whose volume is $`L^3v`$, where $`v`$ is the volume of a single liquid crystal molecule. If the typical volume of a pore is $`V`$, then the number of pores in the lattice cube is $`L^3v/V`$. It immediately follows, $`N_d=L^3vV^1N^{2/3}6^1`$. This can be expressed conveniently in terms of the typical size (diameter) of a pore ($`L_p`$), the linear size of the lattice system ($`L`$) and the volume of a single liquid crystal molecule $`(v)`$ as,

$$N_d(L,L_p,v)=\frac{L^3}{L_p}\left(\frac{v}{36\pi }\right)^{1/3}.$$
(1)
It is clear that a set of quenched sites represents the effect of the wall of a single pore. Several such sets of quenched sites are distributed independently, randomly and uniformly in the lattice cube and they represent the effect of the walls of all the pores. Each active lattice site is connected to every other ones through a path consisting of nearest neighbour connecting segments. In this sense our model represents, approximately, the interconnected network of pores. Also, we can conveniently employ periodic boundary conditions which is known to have least finite size effects. In addition, wherever possible, we can employ finite size scaling to calculate the macroscopic properties in the thermodynamic limit. For a typical liquid crystal molecule of 8CB, studied extensively in different porous media, $`v=480`$ Å3. The porous media studied experimentally have typical sizes ranging from $`80`$Å to $`400`$Å . This situation can be simulated by making proper choices of $`L_p`$ and $`v`$.

Monte Carlo simulations with standard Metropolis algorithm were carried out for five pore sizes (approximately corresponding to experimental situations involving 8CB molecules) and for each pore, for five values of $`L`$. The values of pore sizes considered in our study are $`L_p=400,180,150,110,`$ and $`80`$ Å , besides the bulk liquid crystal, which corresponds to $`N_d=0`$ . For each case, we considered system sizes of $`L=10,15,17,20,`$ and $`22`$, and simulations were carried out at forty values of $`T`$, bracketing the transition temperature. Also, the simulations were repeated for three random disorder configurations and the results were averaged over these runs. For each case, the Monte Carlo simulations were carried out in three stages. In the first stage, short runs of $`20,000`$ to $`50,000`$ Monte Carlo (MC) cycles were carried out. From the data generated, we calculated the values of the specific heat $`C_V`$ at the chosen temperatures. The transition temperature was located approximately as the one at which $`C_V`$ is maximum. In the second stage, extensive simulations consisting of $`500`$ thousand MC cycles with initial hundred thousand cycles ignored for equilibration purpose, were carried out at the approximate transition temperature estimated in the first stage. From these results, we calculated the distribution of energy $`E`$. Then, employing Ferrenberg-Swendson reweighting technique , the energy distribution was analyzed to yield fresh specific heat profiles as a function of $`T`$. The location of the specific heat peak at this stage yields the next level of approximation to the transition temperature. In the third stage, large Monte Carlo runs of five hundred thousand MC cycles were carried out at the transition temperature estimated in the second stage and this process was repeated for three disorder configurations. The distribution of energy thus obtained was analyzed for each case, employing again Ferrenberg-Swendson reweighting technique to obtain profiles of specific heat $`C_V`$, ordering susceptibility $`\chi `$, and the fourth cumulant (also known as Binder cumulant) , $`V_4`$ of energy fluctuations.

$`T_{\mathrm{NI}}`$ deduced from the data on $`C_V`$ as well as from that on $`V_4`$ are plotted as function of $`L^3`$ in Fig. 1(a), for the bulk sample. The two sets of data fall on straight lines establishing finite size scaling. Extrapolating the system size to infinity, we get the transition temperature in the thermodynamic limit. Similar scaling behaviour was found for the simulations carried out for a pore size of $`400`$ Å , and is shown in Fig. 1(b). For smaller pores however, we could not get a neat finite size scaling behaviour; hence, for these cases, we have calculated the transition temperature by taking the average of $`T_{\mathrm{NI}}`$ obtained from $`C_V`$ and $`V_4`$ for the largest lattice size considered ($`L=22`$). The data on $`T_{\mathrm{NI}}`$ thus obtained and normalized to the bulk transition temperature, is plotted as a function of the inverse of the pore size ($`L_p`$), in Fig. 2, alongwith the experimentally measured data . We find that our results match remarkably well with the measured values. For comparison, we have shown the coresponding predictions of the independent pore model, in the same figure.

A set of macroscopic quantities that would reveal the nature of the phase transition through their finite size scaling patterns, includes the peak amplitude and Full Width at Half Maximum (FWHM) of the temperature profiles of $`C_V`$, $`\chi `$ and $`V_4`$. We have shown in Fig. 3 (Left) the peak amplititude of the specific heat $`C_V`$ versus $`L^3`$ for $`L_p=0`$ (bulk sample) and $`400`$Å . In the case of bulk and pore size $`L_p=400`$ Å , the peak amplitude of $`C_V`$ scales as $`L^3`$, as expected for first order transition. For $`L_p=180`$Å however we find $`L^3`$-scaling does not hold good . This failure of $`L^3`$\- scaling is an indication of the softening of the phase transition as the pore size decreases. Similar results on the orientational susceptibity are shown in Fig. 3 (Right). We find that for the bulk and pore size $`400`$Å , the susceptibility scales as $`L^3`$; for pore size $`180`$Å the scaling breaks down .

We calculated the minimum of $`V_4`$ as a function of $`L`$ for the bulk as well as for pore size $`L_p=400`$ Å . In the asymptotic ($`L\mathrm{}`$) limit, this value reaches $`0.64`$. However for the case with $`180`$Å we find that asymptotically the minimum of $`V_4`$ reaches a value of $`0.66`$, expected in a continuous transition. This gives a further indication of the cross over of the transition from the first to the continuous one with decreasing pore size, for details see .

We have analyzed the data on the FWHM of the temperature profiles of $`C_V`$, $`\chi `$ and $`V_4`$. For first order phase transition, these quantities are expected to scale as $`L^3`$. Indeed we find such a scaling behaviour for the bulk and for pore size of $`L_p=400`$ Å . The results for the pore size of $`L_p=180`$ Å , do not exhibit such a scaling giving further evidence for the cross over phenomenon. The details of these results are given in .

The energy distribution near a first order transition should show a double peak structure corresponding to a double well potential. This corresponds to the coexistance of the two phases separated by a potential barrier, representing surface energy terms arising from domain formation. Employing reweighting techniques, the energy distribution obtained as described earlier, was extrapolated to yield distributions in the neighbourhood of the transition temperature. We found that in the case of bulk as well as for pore of size $`400`$ Å , a double peak structure in the distribution of energy yielding a profile with equal peak heights. The energy distribution thus obtained, was transformed to a quantity which differs from free energy only by an additive constant. The ’free energy’ thus obtained had a double well structure for the bulk sample. Also the barrier height was found to increase with increase of $`L`$. A similar result was obtained for the pore size $`400`$Å . We then plotted the barrier height as a function of $`L^2`$ for the bulk sample as well as for $`L_P=400`$Å . The curve was found to be linear establishing $`L^2`$-scaling. This confirms that the transition is first order. We could not extract such structures (typical for first order transition) for the cases with pore sizes $`180`$Å and less. This provides a definitive evidence that the nature of the transition indeed crosses over from the first order to the continuous one, when the pore size decreases. The details are given in .

The enthalpy exchanged in the transition can be obtained as integrated area of the specific heat profile and is depicted in Fig. 4 as a function of pore size; also shown are the experimental results. The agreeement is found to be very good.

In conclusion, we have proposed a model to study the phase behaviour of liquid crystals in a porous medium. We have given a method for calculating the disorder parameter $`N_d`$ in terms of the typical pore size, the system size, and the volume of the liquid crystal molecule. The predictions of our model on the suppression of the transition temperature with decreasing pore size, on the softening of the transition as the pore size decreases and on the enthalpy of transition compare well with experimental findings.

We thank G. Bharadwajkumar and V. Satheesh for help during the early stage of the calculations, B. V. R. Tata and M. C. Valsakumar for discussions and the Board of Research in Nuclear Sciences, India (Project No. 99/37/09/BRNS) and the office of the Naval Research, USA (Grant No. N00014-97-1-0994) for support.

FIGURE CAPTIONS

Figure 1 Transition temperature deduced from Specific heat ($`\times `$) and the Binder Cumulant ( o ), plotted against $`L^3`$. (a) for the bulk (b) for $`L_p=400`$ Å . The error bars are also marked.

Figure 2 $`T_{\mathrm{NI}}(L_p)/T_{\mathrm{NI}}`$(bulk) versus $`1/L_p`$. o - Experiment; + - present simulations; $`\times `$ \- independent pore model.

Figure 3 Maximum of specific heat (Left) and of the orientational susceptibility (Right) plotted against $`L^3`$. $`\times `$ \- bulk; o - pore size 400Å .

Figure 4 Variation of relative enthalpy with pore size; o - experimental values; $`\times `$ \- our results .