text:
# Moving black holes via singularity excision

## I Introduction

Recent advances in numerical relativity are expanding the dynamical range covered by non-linear, three-dimensional simulations of general relativistic systems. However, we are still far from having evolutions that successfully simulate a binary black-hole system starting from the in-spiral of the holes, continuing through their merger and ending in the ring-down of the resulting single black hole. The obstacles defining this challenge are multiple and not all understood or even known. For example, it has become increasingly evident that the form of the Einstein equations implemented in numerical codes is one of the important aspects determining the behavior of the simulations. In principle, there are an infinite number of ways to recast the Einstein equations as a Cauchy problem. The search for the formulation or family of formulations that will lead to long-lived evolutions has, in some situations, become a tour-de-force. In spite of the many “3+1 flavors” of the Einstein equations so far introduced, there has not been a formulation exhibiting a clear superiority.

Besides the form of the evolution equations, there are other essential ingredients that have been demonstrated to yield improvements in the duration of the simulations. Gauge or coordinate conditions, methods to handle the singularities, non-reflective boundary conditions and discretization methods are some examples of these elements Luis\_review ; Thomas\_review .

Another major focus in numerical relativity has been producing long-lasting, three-dimensional evolutions of the simplest case of a space-time containing black-hole singularities, namely that of a single, excised black hole in which the location of the black hole remains fixed on the computational grid Pitt\_stable . It is safe to say that this milestone has been achieved Alcubierre1 ; KST ; Yo\_stable . One should keep in mind, however, that some of the elements used to obtain simulations lasting forever cannot be extended or applied to binary black-hole orbits and coalescences. Specifically, the elements in question, as we will later discuss in detail, are those that hinge on keeping fixed the coordinate location of the black-hole singularities and thus the excision boundary.

For sufficiently large separations, even in the non-linear regime, it should be possible to construct an approximate co-rotating coordinate system in which the coordinate drift of the black holes is avoided via an appropriate set of gauge or coordinate conditions IBBH . However, once a given co-rotating frame can no longer keep the location of the black holes fixed without introducing extreme grid stretching, construction of a new co-rotating frame will be needed. It is then likely that in this new co-rotating frame, the coordinate location of the black holes, and therefore the excision regions, will be different. Consequently, there is a natural need to develop an excision algorithm that is able to accommodate motion of the black holes. This implies handling, among other things, grid points that emerge from inside of the excised regions and become part of the computational domain.

There are examples of simulations with moving excision already in the literature. Using a characteristic formulation, long-lasting evolutions of a wobbling black hole have been obtained wobble . Other examples of simulations involving motion of the excision boundary include boosted, single black holes moving  and grazing collisions Brandt . Recently, Yo, Baumgarte and Shapiro Yo\_moving  presented a simple scheme for a moving excision boundary when solving the scalar field equation in the background of a boosted black hole.

The objective of this paper is to enhance the excision algorithm used in long-term stable simulations of single non-moving black holes to handle the case of moving black holes. Our approach can be viewed as an extension of the work in Ref. Yo\_moving . The paper is structured as follows: the formulation of the Einstein equations and the method used in our code, called Maya, to solve these equations are summarized in § II. The crux of this paper, namely the implementation of excision, is given in § III. Tests of the excision method for a moving black hole are presented in § IV. Finally, we discuss the ramifications of our results and future developments in § V.

## II The Maya Code

When viewed as a Cauchy problem, the Einstein equations provide a mechanism to construct the time-history of the geometry of three-dimensional space-like hypersurfaces. The geometry of each hypersurface is characterized by the intrinsic metric $`g_{ij}`$ and extrinsic curvature $`K_{ij}`$. The space-time foliation is glued together with the help of the lapse function $`\alpha `$ and the shift vector $`\beta ^i`$. In principle, the problem reduces to finding $`g_{ij}`$ and $`K_{ij}`$ for all the hypersurfaces by solving the Einstein equations written explicitly in terms of $`g_{ij}`$ and $`K_{ij}`$, namely the Einstein equations in a 3+1 form. An example of these equations are the ADM equations ADM ; York . However, evolutions obtained with the ADM equations have had very limited success. Nakamura and Shibata SN , and later Baumgarte and Shapiro BS , introduced a formulation (BSSN) of the Einstein equations that has clearly yielded remarkable improvements. Instead of $`g_{ij}`$ and $`K_{ij}`$ as primary variables, the BSSN formulation introduces new variables $`\mathrm{\Phi },\widehat{g}_{ij},K,\widehat{A}_{ij}`$ and $`\widehat{\mathrm{\Gamma }}^i`$. The relationships between the BSSN and ADM variables are:

$`\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{1}{6}}\mathrm{ln}g^{1/2}`$ (1)
$`\widehat{g}_{ij}`$ $`=`$ $`e^{4\mathrm{\Phi }}g_{ij}`$ (2)
$`K`$ $`=`$ $`g^{ij}K_{ij}`$ (3)
$`\widehat{A}_{ij}`$ $`=`$ $`e^{4\mathrm{\Phi }}A_{ij}`$ (4)
$`\widehat{\mathrm{\Gamma }}^i`$ $``$ $`\widehat{g}^{jk}\widehat{\mathrm{\Gamma }}_{jk}^i=_j\widehat{g}^{ij},`$ (5)
where $`A_{ij}=K_{ij}g_{ij}K/3`$. Above, Eqs. (1) and (2) imply that the conformal metric $`\widehat{g}_{ij}`$ has unit determinant, which in turn yields the second equality in Eq. (5). In addition to the standard constraints, Hamiltonian and momentum constraints, the BSSN equations require satisfying the conditions $`\sqrt{\widehat{g}}=1`$ and $`\widehat{A}_i^i=0`$ if one chooses to evolve all the tensor components of $`\widehat{g}_{ij}`$ and $`\widehat{A}_{ij}`$. Explicitly enforcing $`\sqrt{\widehat{g}}=1`$ is not necessary, but insuring that $`\widehat{A}_{ij}`$ remains trace-free is crucial in extending the life of numerical evolutions. In terms of the new variables (1-5), the Einstein evolution equations take the following form:

$`_t\mathrm{\Phi }`$ $`=`$ $`_\beta \mathrm{\Phi }{\displaystyle \frac{1}{6}}\alpha K`$ (6)
$`_t\widehat{g}_{ij}`$ $`=`$ $`_\beta \widehat{g}_{ij}2\alpha \widehat{A}_{ij}`$ (7)
$`_tK`$ $`=`$ $`_\beta K_i^i\alpha +\alpha (\widehat{A}_{ij}\widehat{A}^{ij}+K^2/3)`$ (8)
$`_t\widehat{A}_{ij}`$ $`=`$ $`_\beta \widehat{A}_{ij}+e^{4\mathrm{\Phi }}\left(_i_j\alpha +\alpha R_{ij}\right)^{TF}`$ (9)
$`+`$ $`\alpha (K\widehat{A}_{ij}2\widehat{A}_{il}\widehat{A}_j^l)`$
$`_t\widehat{\mathrm{\Gamma }}^i`$ $`=`$ $`_\beta \widehat{\mathrm{\Gamma }}^i2\widehat{A}^{ij}_j\alpha `$ (10)
$`+`$ $`2\alpha \widehat{\mathrm{\Gamma }}_{jk}^i\widehat{A}^{jk}+12\alpha \widehat{A}^{ij}_j\mathrm{\Phi }{\displaystyle \frac{4}{3}}\alpha \widehat{g}^{ij}_jK`$
$`+`$ $`\left[\left(\chi +{\displaystyle \frac{2}{3}}\right)\widehat{\gamma }^{kl}\widehat{\mathrm{\Gamma }}_{kl}^i\chi \widehat{\mathrm{\Gamma }}^i\right]_j\beta ^j.`$
The term in Eq. (9) with superscript $`TF`$ denotes the trace-free part of the tensor terms between brackets. Eq. (10) has been modified from the original BSSN form ($`\chi `$-terms) following Ref. Yo\_stable . This modification is crucial to achieve stability in single, non-moving black-hole simulations. Finally, the first terms in the r.h.s. of Eqs. (6-10) involving $`_\beta `$ denote Lie derivatives of the scalar $`K`$ and tensor densities $`\widehat{g}_{ij}`$ and $`\widehat{A}_{ij}`$, as well as their extension to quantities related to tensor densities, such as $`\mathrm{\Phi }`$ and $`\widehat{\mathrm{\Gamma }}^i`$.

We have developed a numerical code, called Maya, to solve Eqs. (6-10). The code views this system of equations as having the following structure:

$$_tu=\beta ^i_iu+T(u,\beta )+S(u,u,^2u)\rho $$
(11)
where $`u=\{\mathrm{\Phi },\widehat{g}_{ij},K,\widehat{A}_{ij},\widehat{\mathrm{\Gamma }}^i\}`$. Notice that the terms $`_\beta u`$ have been explicitly separated into $`\beta ^i_iu+T(u,\beta )`$ terms. In the interior of the computational domain, that is, excluding the outer and excision boundaries, spatial derivatives appearing in $`T`$ and $`S`$ are approximated using second-order, centered finite differencing. The terms of the form $`\beta ^i_iu`$ are commonly called “advection” terms. They are approximate via a second order or higher upwind differencing; for details, see Ref. Kelly .

At the outer boundary, the traditional approach for BSSN-based evolutions has been to assume that

$$u=u_o+\frac{w(tr)}{r},$$
(12)
where $`u_o`$ represents an analytic solution. For the case of single black holes, this analytic solution is the exact solution to the Einstein equations. If an exact solution is not available, the flat space-time solution is used provided the outer boundary is sufficiently far away from the holes. Eq. (12) is then implemented in a differential form, namely

$$_tu=\frac{x^i}{r}_i(uu_o)+\frac{(uu_o)}{r}.$$
(13)
The r.h.s. of Eq. (13) is discretized as in the advection term in Eqs. (6-10). It was found Alcubierre1 ; Yo\_stable  that in order to achieve stability for single, non-moving black-hole evolutions, condition (13) should not be used for the connection $`\widehat{\mathrm{\Gamma }}^i`$. For this field, one needs instead to impose $`_t\widehat{\mathrm{\Gamma }}^i=0`$.

Notice that in deriving Eq. (13), it was assumed that the analytic solution is time-independent. This will not be the case for the moving black hole in our tests. For simplicity, we have chosen instead to directly use the Dirichlet condition $`u=u_o`$. This will turn out to be one of the two factors limiting the duration of the simulations.

At the excision boundary, the discretization of spatial derivatives also requires special treatment. This constitutes the focus of the paper and is addressed in detail in the following section.

Once all spatial operators in Eq. (11) are discretized, $`u`$ is updated using a second-order, iterative Crank-Nicholson scheme Teukolsky . Parallelization, IO and parameter manipulations in the Maya code are handled via the Cactus toolkit developed by the Albert Einstein Institute in Golm, Germany Cactus .

Finally, several studies have shown Kelly ; Alcubierre1 ; nct ; Yo\_stable  that stable evolutions of single black holes are not possible if both the lapse function and shift vector are analytic functions from the exact solutions. Instead, one needs “driver” gauge conditions such as the $`1+\mathrm{log}`$ condition for $`\alpha `$ or the $`\mathrm{\Gamma }`$-driver for $`\beta ^i`$; for a review of gauge conditions, see Ref. AEI\_review . Unfortunately, these driver-gauge conditions have been designed having in mind time-independent solutions. These gauge conditions will not be applicable for the excision tests of a moving black hole under our consideration. We are then limited to using the exact solution to specify $`\alpha `$ and $`\beta ^i`$. Here resides the second factor limiting the duration of our simulations.

## III Black Hole Excision

From causality considerations, one argues that it is completely unnecessary to evolve the system of equations inside the event horizon of a black hole – no physically meaningful information can possibly propagate out to affect the space-time outside the horizon. Doing so has the desirable effect of removing the black-hole singularity from the computational domain. If one chooses to ignore, or excise, the interior of a black hole, the immediate problem one must face is that of finding its event horizon. The event horizon, however, is a global object requiring a complete knowledge of the space-time, namely the solution to the problem itself. Unruh unruh  suggested using instead the apparent horizon. Apparent horizons wald  are defined in terms of quantities local-in-time. That is, to find apparent horizons, one only needs the spatial metric $`g_{ij}`$ and extrinsic curvature $`K_{ij}`$ of the space-like hypersurfaces in the 3+1 foliation. Thus, following the history of apparent horizons naturally adapts to the evolution of the system since only information at a given instant of time is required. The importance of apparent horizons in connection with excision lies on the fact that an apparent horizon will always be coincident or contained within the event horizon. Therefore, as long as the region that is ignored or excised is bounded by or contained within the apparent horizon, this region will also be interior to the event horizon.

Once the interior of an apparent horizon has been removed from the calculation, the next step to address is whether boundary conditions are required at the excision boundary. If all the fields involved in the calculation have at the excision boundary outgoing characteristics in the direction of the singularity, away from the computational domain, there is no need for imposing boundary conditions. In this case, the fields at the excision boundary are evolved in the same form as the fields in the interior of the computational domain.

The problem in general relativity is that the only characteristics that one can blindly assume to be outgoing at the excision boundary are those of physical modes since only those modes are causally constrained. For other modes, such as gauge modes, the characteristics depend on the particular structure of the evolution equations as well as the set of gauge conditions imposed. Manifestly hyperbolic formulations of the Einstein equations facilitate determining the characteristics of all the fields, thus allowing the identification of those fields requiring conditions at the excision boundary. Here is one of the main reasons for the popularity of explicitly hyperbolic formulations of the Einstein equations hyper\_review . This would suggest that numerical simulations involving excision are only possible using a hyperbolic formulation. However, studies by several groups using the BSSN formulation, a formulation that is not explicitly hyperbolic, have shown BS ; AEI\_BSSN ; BSSN\_hyper  that this is not the case. BSSN-based codes, assuming all fields having outgoing characteristics at the excision boundary, have produced results similar to those based on hyperbolic formulations.

In summary, we work under the assumption that there is no need for boundary conditions at the excision boundary. The only task is to design a discretization of the evolution equations that is appropriate at the excision boundary where centered finite differencing is no longer feasible.

A possible approach to handle spatial derivatives at the excision boundary is to modify all finite-difference stencils at this boundary so as to avoid using grid-points from within the excised region (e.g. one-sided differences), while still keeping the same order of accuracy. This typically involves getting lost in a veritable jungle of logic, especially when the excision region is that of the “Lego sphere”, a representation of a sphere using only the points on a fixed Cartesian grid. Each type of corner must be treated with its own finite difference stencil. This approach was followed in early simulations moving ; Brandt . To circumvent the complications of this stencil-modification method, Alcubierre and Brügmann Alcubierre1  introduced a very simple excision method which is ideally suited for stationary space-times such as those involved with single black-hole evolutions. In Alcubierre1 , simulations of single black holes that did not assume any symmetries had limited duration. It has been recently demonstrated Yo\_stable , however, that the factor restricting the duration of these simulations was not the excision method used, but a term in the evolution equation of the connection $`\widehat{\mathrm{\Gamma }}^i`$.

Anticipating that for more general cases, e.g. binary black holes, the excision boundary and the r.h.s. in Eqs. (6-10) will have strong temporal dependence, we present a method for handling the excision boundary that generalizes to non-stationary space-times. Our approach extends the method used by Yo, Baumgarte and Shapiro in their study of the dynamics of a scalar field in the background of a single, boosted black hole Yo\_moving . We have carried out tests to show that our method, when applied to the case of a non-moving black hole, reproduces the observed long-term stability observed in Ref. Yo\_stable .

Once an apparent horizon is found, the first step is to choose the shape and size of the excision region. A brief comment is warranted regarding the choice of the shape of the excision region. Historically, a sphere was used. The rationale for this is that the apparent horizon will tend to be spherical or ellipsoidal in all but extremely dynamical space-times. As mentioned before, on the Cartesian grids typically used in three-dimensional numerics, one achieves only a crude approximation of the sphere, colloquially known as a “Lego sphere”. Alcubierre and Brügmann simplified matters Alcubierre1  by using a cube as an excision shape. We have experimented with a variety of excision shapes, all possessing at the very least octahedral symmetry (because of the frequent use of octant, quadrant, and equatorial symmetry boundary conditions in performing three-dimensional simulations), and in the end, have settled on the sphere, the cube, and the cuboctahedron. These shapes are shown in Fig. 1.

To facilitate the excision method, we adopt the common practice of carrying an extra grid function, called the “mask”, to indicate the state of a grid-point in the computational domain. The typical states of a grid-point are: excised, excision boundary, interior, and outer boundary. Fig. 2 depicts a schematic representation of the mask variable. The large, dark circle is the location of the apparent horizon. The excision region is a “Lego sphere” represented by the thick line. The filled circles are the excision boundary points, while the empty circles are the excised points. We also use the mask function to label the outer boundary points, but this is irrelevant in the present discussion on excision.

Given $`u^n`$ at the interior, excision boundary and outer boundary points, the r.h.s. $`\rho `$ of the evolution equations (11) for the interior grid-points are constructed and used to update $`u^nu^{n+1}`$ at these points. Before any subsequent updating takes place, including intermediate iterated Crank-Nicholson steps, updating of $`u^n`$ at the excision boundary is required.

Our approach for updating data at the excision boundary is designed to avoid altering the finite difference stencils used for the interior points. This is crucial for preserving a simple code. It consists of obtaining the data at the excision boundary via polynomial extrapolation. In order to simplify matters even further, we perform one-dimensional extrapolations along normal directions to the excision boundary. Because in general the excision boundary is only a discretized approximation to the geometric shape desired, the normal to the excision boundary is not aligned with the grid-points in the computational domain. We select the direction of extrapolation to be as close as possible to the direction of the normal for which all data involved on the extrapolation lie on grid-points.

We consider two approaches to construct data at the excision boundary. One approach involves the extrapolation of $`\rho `$ from the interior points to the excision boundary (rhs-extrapolation method). The other approach is based on extrapolation of the solution $`u^{n+1}`$ (sol-extrapolation method). Fig. 3 shows an illustration of both methods. Excised points are labeled with an X, the excision boundary points with a box and interior points with a solid dot. As mentioned before, for each interior point we can use centered finite differencing to calculate $`\rho `$. That is, $`\rho _2`$ through $`\rho _6`$ in Fig. 3 are known, but $`\rho _1`$ is unknown. The value of $`\rho _1`$ is, in principle, required to obtain $`u_1^{n+1}`$. The rhs-extrapolation method finds the value $`\rho _1`$ using the values of $`\rho _2`$ through $`\rho _5`$, depending on the order of the extrapolation.

Given that for the interior points, the truncation error in the discretizations to obtain $`\rho `$ is second-order, the rhs-extrapolation must be such that the overall truncation error is preserved. This implies rhs-extrapolations of at least second-order. That is,

$$\rho _1=2\rho _2\rho _3+O(h^2),$$
(14)
with $`h`$ the grid spacing. It is important to notice that the simple excision method first used by Alcubierre1  and later by Yo\_stable  consists of first-order extrapolation, namely $`\rho _1=\rho _2`$. This choice yields truncation errors of $`O(h)`$. Our numerical experiments with a moving black hole using zero-order extrapolations quickly became unstable, in agreement with similar results in Yo\_moving . A possible reason why zero-order extrapolation works for non-moving single black holes is that, from Eq. (11), the $`\rho _1=\rho _2`$ condition is equivalent in the continuum to

$$0=_r\rho =_r_tu=_t_ru,$$
(15)
which is consistent with the time-independent nature of the solution.

The sol-extrapolation method uses the values $`u_2^{n+1}`$ through $`u_5^{n+1}`$ to obtain $`u_1^{n+1}`$, as depicted by the solid black arrow. In order to preserve second-order accuracy of the solution in this case, third-order extrapolations are needed. This is because the extrapolated value $`u_1^{n+1}`$ will be used in computing $`\rho _2`$ for the updating $`u_2^{n+1}u_2^{n+2}`$.

As the coordinate location of a black hole changes, so does the position of the excision region. We adjust the numerical evolution such that the excision region does not move by more than a grid-point. The consequences of this change of location are that points previously labeled excision boundary are now interior and points previously labeled excised are now excision boundary. Fig. 4 depicts this process. The movement is indicated by the shift in the apparent horizon (dark circles), which causes the excision region to shift in order to remain centered on the apparent horizon. The dotted circles are the original apparent horizon and excision region. The filled circles are the new excision boundary points that were previously excised points. While nothing needs to be done with interior points, we need to populate the previously excised points with values. This is done following the sol-extrapolation method previously described.

## IV Dancing Black Hole

In this section, we present numerical experiments that demonstrate the ability of the Maya code to handle moving excision boundaries. Because the algorithms and data structures presented here are general, they can provide a preliminary treatment for non-trivial time-dependent systems such as black-hole binaries.

In order to test our excision algorithm, we introduce the simplest time-dependency into the system that does not change the form of the metric. We consider a single, non-rotating black hole in ingoing Eddington-Finkelstein coordinates and apply a coordinate transformation of the form:

$$x^ix^i+d^i,$$
(16)
where

$`d^i`$ $`=`$ $`[A\mathrm{sin}\left({\displaystyle \frac{Vt}{A}}\right),0,0]\text{or}`$ (17)
$`d^i`$ $`=`$ $`[A\mathrm{cos}\left({\displaystyle \frac{Vt}{A}}\right),A\mathrm{sin}\left({\displaystyle \frac{Vt}{A}}\right),0].`$ (18)
Notice that this is a coordinate transformation of the spatial coordinates only. The above transformation induces the following change on the shift vector:

$$\beta ^i\beta ^i+b^i,$$
(19)
where

$`b^i`$ $`=`$ $`[V\mathrm{cos}\left({\displaystyle \frac{Vt}{A}}\right),0,0]\text{or}`$ (20)
$`b^i`$ $`=`$ $`[V\mathrm{sin}\left({\displaystyle \frac{Vt}{A}}\right),V\mathrm{cos}\left({\displaystyle \frac{Vt}{A}}\right),0].`$ (21)
The functional form of $`\alpha ,g_{ij}`$ and $`K_{ij}`$ remains the same. The only effect of this coordinate transformation is to change the coordinate location of the black hole. For the case (17), the black hole bounces back and forth with an amplitude $`A`$ and maximum coordinate velocity $`V`$ (bouncing-hole). Similarly, for the case (18), the black hole moves in a circle of radius $`A`$ and coordinate velocity $`V`$ (circling-hole).

The size of the excision region depends on the shape chosen. In the case of excising a spherical region, the excision radius was set to $`1.5M`$ with $`M`$ the mass of the black hole and $`2M`$ its radius, while the cubical excision region is inscribed to fit inside a sphere of radius $`1.5M`$. That is, we are excising a larger region of the computational domain when we use a sphere. We have empirically found that sol-extrapolation performs a little better than rhs-extrapolation for the tests under consideration. We use third-order ($`O(h^3)`$) sol-extrapolation for the excision boundary and fourth-order for repopulating. The simulations were carried out with grid-spacings $`h=0.2M`$ and imposing equatorial symmetry about the $`z`$-axis. The computational domain has a total length of $`20M`$ in the $`x`$ and $`y`$ directions and $`7M`$ in the $`z`$ direction. We used values of $`A=1.0M`$ and $`V=0.5`$ for the simulations presented here.

The driver gauge conditions commonly used to produce long-term stability in single, non-moving black-hole evolutions are the $`1+\mathrm{log}`$ slicing and $`\mathrm{\Gamma }`$-driver; that is,

$`_t\alpha `$ $`=`$ $`_i\beta ^i\alpha K`$ (22)
$`_t\beta ^i`$ $`=`$ $`\lambda _t\widehat{\mathrm{\Gamma }}^i,`$ (23)
respectively. Unfortunately, these gauge condition do not work for our moving black-hole tests. The black-hole solution that one obtains under the coordinate transformation (16) are not compatible with these conditions. New or modified driver gauge conditions would have to be constructed. These new driver conditions would likely be dependent on the particular nature of our moving black-hole solution, thus not generalizable to astrophysically interesting cases such as black-hole binaries. Because the focus of the work in this paper is on excision, we have chosen to set $`\alpha `$ and $`\beta ^i`$ to the analytic values given by the exact solution. It is well known that using analytic lapse and shift gauge conditions yields evolutions of single, non-moving black hole lasting $`<100M`$. As we will see next, the limits in duration of our moving black-hole simulations are consistent with the stationary cases.

As mentioned before, a similar problem arises with the boundary condition (13). This condition assumes that the analytic solution is time-independent. For the case of time-dependent solutions, the condition becomes

$$_tu=_tu_o\frac{x^i}{r}_i(uu_o)+\frac{(uu_o)}{r}.$$
(24)
The modified condition (24) should in principle work for our time-dependent case. However, one needs to remember that the outer boundary condition for the connection $`\widehat{\mathrm{\Gamma }}^i`$ that yields long-term stability for non-moving holes was not (13), but $`_t\widehat{\mathrm{\Gamma }}^i=0`$. We were unable to obtain a condition applicable to our time-dependent case that would lead to significant improvements in the duration of the evolutions. Therefore, for simplicity, we have chosen instead to directly use for all fields the analytic solution as the boundary condition.

Fig. 5 is a $`xt`$-plot of the variable $`\mathrm{\Phi }`$ for the bouncing-hole simulation. The sinusoidal flat region is a $`x`$-cross section of the world-tube of the excision boundary. Fig. 6 is also a $`xt`$-plot but in this case for the normalized Hamiltonian constraint. The reflections at the outer boundary are clear from this figure. These effects become stronger with time and eventually are one of the main causes of the simulation terminating. The other contributing factors are the analytic gauge conditions. The circling case is depicted in Fig. 7. The figure consists of stacked time slices of the variable $`\mathrm{\Phi }`$. The time slices are separated by $`15M`$ increments in time. For clarity, the domain displayed only represents $`\pm 5M`$ in the $`x`$ and $`y`$ dimensions. The “orbit” of the black hole is $`4\pi `$, so each slice represents an evolution of approximately one and a quarter orbits from the slice preceding it. Animations of both bouncing and circling holes can be found at movies .

For comparison of evolutions, we monitor the Hamiltonian constraint, the area of the apparent horizon, and a “mass” function $`M_\mathrm{\Phi }`$ obtained from

$$M_\mathrm{\Phi }=\frac{r}{2}e^{12\mathrm{\Phi }}1.$$
(25)
For the exact solution, $`M_\mathrm{\Phi }`$ is equal to the mass of the black hole. In all the following figures (Figs. 9-11), the top panels depict results using a spherical excision region and the bottom cubical. We also show results from a stationary black-hole case with a solid line. This is our reference case as we do not expect the moving cases to supersede the stationary case. The dotted lines represent the bouncing-hole and the dashed the circling-hole.

Fig. 9 contains plots of the L2-norm of the Hamiltonian constraint versus time for each case. It is evident from this figure that the simulation is unstable even in the case of a stationary black hole. The instability is likely due to the gauge and boundary conditions. When we turned on the driver conditions (22) and (23), as well as the boundary condition (13), long-term stable runs for a stationary hole were obtained, as one can see from Fig. 8, similar to those in Yo\_stable . Fig. 10 shows the area of the apparent horizon versus time; the apparent horizon tracker is described in schnetter . The apparent horizon tracker fails to give reasonable results at approximately $`t=90M`$ because the coordinate distortions are such that the horizon intersects the excision region at one or more points. Finally, Fig. 11 shows plots of the L2-norm of $`M_\mathrm{\Phi }`$ versus time. The reason why the L2-norm of $`M_\mathrm{\Phi }`$ does not show the jagged behavior is because the L2-norm of $`M_\mathrm{\Phi }`$ was done with data that did not include excision boundary data. The data were taken from within a shell centered at the location of the black hole a distance $`0.5M`$ from the excision boundary. The jaggedness in the Hamiltonian constraint plots is due to the artificial jump of values because the center of the excision region is forced to be at a grid point, thus producing discontinuous motion of the excision region.

## V Conclusions

Ref Alcubierre1  described a simple excision algorithm that extended the lifetimes of time-independent, single black-hole evolutions. This work was expanded in Ref. Yo\_moving , showing that higher order extrapolation is necessary for stability in studying the dynamics of a scalar field in the background of a boosted Kerr-Schild background. We have generalized and applied the excision algorithm presented in Yo\_moving  to the case of a single black hole moving its coordinate location throughout the computational domain. We showed that this algorithm is generalizable to non-stationary black-hole space-times. We also find that, despite the fact that a cubical excision region is all that is necessary for long-lived stationary black-hole runs, a spherical excision boundary is viable. A spherical shape may be more desirable in dynamic simulations because it allows the excision region to fit more easily within the apparent horizon. It is important to emphasize that the algorithm in this work for excising the black hole as it moves through the computational domain behaves in a similar way to excising a stationary black hole. This suggest that as gauge and boundary conditions suitable for time-dependent solutions are developed, the excision method discussed in this work can be directly implemented without requiring further or with minimal modifications.

## VI Acknowledgments

We thank Miguel Alcubierre, Bernd Brügmann and Jorge Pullin for helpful discussions. Parallel and I/O infrastructure provided by Cactus. We acknowledge the support of the Center for Gravitational Wave Physics funded by the National Science Foundation under Cooperative Agreement PHY-0114375. Work partially supported by NSF grants PHY-9800973 to Penn State and PHY-9800737 PHY-9900672 to Cornell University.


format: text

text:
# Three New Long Period X-ray Pulsars Discovered in the Small Magellanic Cloud

## 1 Introduction

Our program of routinely evaluating the pulsar content of archival observations of the Small Magellanic Cloud has previously revealed five new X-ray pulsars (Macomb et al. 1999, Finger et al. 2001, Lamb et al. 2002a, Lamb et al. 2002b). To date, there are at least 31 known X-ray pulsars in the Small Magellanic Cloud, 25 of which were discovered in the last six years (Yokogawa 2002, Lamb et al. 2002a, Laycock et al. 2002, Corbet et al. 2002, Lamb et al. 2002b). This explosion in the number of known pulsars is due mainly to four X-ray missions, ASCA, ROSAT, RXTE, and Chandra, each of which has devoted time to long, sensitive observations of the SMC.

The advantages to studying the compact object content of the Magellanic Clouds are well known: a relatively small angular size, at a reasonably close and known distance, located at a high galactic latitude with little obscuration by the interstellar medium of the Galaxy. In addition, there is evidence for recent star formation in the SMC over the last few 10’s of millions of years which creates an environment in which X-ray pulsars are expected to be plentiful (Haberl & Sasaki 2000, Maragoudaki et al. 2001).

Our analysis of a recent Chandra X-ray Observatory observation of the SMC reveals evidence for four X-ray pulsars. One, CXOU J011043.1-721134, is possibly an anomalous X-ray pulsar (Lamb et al. 2002b). This paper describes the discovery of X-ray pulsations from the other three sources: CXOU J005750.3-720756 (152.10 s period), CXOU J010102.7-720658 (304.49 s period) and CXOU J005736.2-721934 (564.83 s period). After a description of the discoveries, we address the characteristics of each source, followed by a discussion of the general significance of these detections.

## 2 Discovery of Pulsations

The discovery data were obtained from CHANDRA obsid 1881, a 100ks ACIS-I observation which began 2001 May 15. The ACIS-I read-out time for this observation was set at a coarse 3.241 seconds, giving a Nyquist frequency of 0.154 Hz. Processing of the ACIS field reveals 140 point sources detected at the 3 sigma level using a standard sliding box detection algorithm. We extracted photons for each point source in the field using an ellipse of varying size, defined by the source detection software, which encompassed $`90\%`$ of the source photons.

The photon arrival times were barycentered and a fast Fourier transform performed for each individual source. With the average fft power normalized to one, we set a power threshold of 20 to flag potential signals. For a 100 ksec observation binned at 3.241 seconds, there are about 31000 independent trials. This gives a probability of around $`1\%`$ that any peak power beyond 20 arises by chance for all of the 140 sources analyzed.

Concentrating on frequencies above 0.001 Hz (a 1000 second period), we find 21 of the 140 sources with powers beyond 20. Of these, all but four sources had powers at two set frequencies that appeared in multiple sources. The first frequency was near 0.001 Hz and the second at 0.001415 Hz, periods of 1000 and 706 seconds, both of which are attributable to known dithering effects (as described at http://cxc.harvard.edu/ciao/why/dither.html). Of the four remaining sources, one is CXOU J011043.1-721134 which is a 5.44 s pulsar (Lamb et al. 2002b), leaving us with three new candidate pulsars. Table 1 lists the Chandra centroid for each source along with the proposed counterpart.

The Fourier power distributions for the three remaining sources are shown in Figure 1. These pulsar timing distributions are based upon a slightly different photon set than for the general search. These refined photon selections encompassing the energy range from 0.6 - 7.5 keV are close to maximizing the Fourier power for all three sources. In addition, we chose position regions corresponding to an 80% encircled energy to further optimize pulsar timing parameters. These selections resulted in 5429 photons for the 152 s pulsar, 325 photons for the 304 second pulsar, and 3141 photons for the 565 second pulsar.

The first source, centered on a frequency of 0.0065747 Hz (152.10 seconds) has a peak power of 90 while the second is at a frequency of 0.0032842 Hz (304.49 seconds) with peak power 30.4. While this frequency is close to being twice the frequency of the 152 s pulsar, they are not precisely doubled. In addition, the fact these sources are in different ACIS chips and that neither of these frequencies show up in any other field source no matter the source strength indicates that they are in fact distinct signals. The third source is detected at a power of 57 at a frequency of 0.00177045 Hz (564.83).

In the case of the 152 and 565 second pulsars, a third harmonic of the fundamental frequency is clearly seen, and there is also a strong second harmonic for the 152 second pulsar. No other sources in the field of view have significant Fourier power at any of these frequencies. The implication is that these three sources are long-period X-ray pulsars. As their individual characteristics make clear, they are both likely to be examples of high-mass X-ray binary systems with Be star companions.

There is little evidence of variability from any of the three new sources over the 100 ksec Chandra observation. While a $`\chi ^2`$ test of variability for the pulsars binned over single pulse periods for the 152 and 565 second pulsars and 5 pulse periods for the 304 second pulsar shows little gross variability, there is modest variability for the 565 second pulsar (3% chance of being consistent with a flat rate). There is no evidence for any of the pulsars dramatically changing their flux during the observation, i.e. an onset or end of an emission episode. Consequently, we have taken all spectral and timing information for these three pulsars from the full observation.

## 3 Individual Source Characteristics

### 3.1 The 152 second pulsar: CXOU J005750.3-720756

A search of the ROSAT observation catalog shows many ROSAT HRI or PSPC observations overlapping the position of CXOU J005750.3-720756. In addition, the ROSAT results archive lists a source, RX J0057.8-7207, a source with an 11 arcsecond error radius that is within 6 arcsecond of the Chandra position centroid. No previous reports of X-ray pulsations from this source have been reported, however it was suggested as a Be binary pulsar candidate by Haberl & Sasaki (2000) based upon positional coincidence with an H-alpha region cataloged by Meyssonnier & Azzopardi (1993). Given the small positional difference between the CHANDRA position (absolute position good to approximately 0.6”) and the centroid of the ROSAT/PSPC source (localization 5.1”, Haberl & Sasaki 2000), it is likely that these two sources are the same. Unfortunately, a check of the three strongest ROSAT detections of this source does not reveal evidence for pulsed emission.

This is not surprising, however, when one considers the energy spectrum of this source. The spectrum of CXOU J005750.3-720756 has a hard photon index of $`1.0\pm 0.1`$. The complete power law fit parameters are listed in Table 2, along with the estimated pulse period and period error. The error on the pulse period is calculated from the Rayleigh power, based upon the frequency points corresponding to a drop in power equal to the square root of the peak power. The source flux is the highest of these three sources. Given this very hard spectrum, it is likely that a soft X-ray telescope with reduced effective area such as the ROSAT/PSPC and ROSAT/HRI are less effective at detecting the periodicity.

Figure 2 shows the folded Chandra lightcurve for the 152.1 second pulsar. The pulse profile is complicated, showing a sharp dip with a duration of 10% of the pulse period. The pulsed fraction is 63$`\pm 3`$%. A search of the full dataset over frequency derivative (binning at 1.0x10-13 Hz/s) and frequency provides a small but statistically insignificant improvement in detected power. This is the case for all three pulsars, so all quoted pulse periods and frequencies are uncorrected for frequency derivative.

### 3.2 The 304 second pulsar: CXOU J010102.7-720658

The weakest of the three sources by an order of magnitude, CXOU J010102.7-720658 has a less well-defined pulse profile as shown by Figure 2. The source has a pulsed fraction of 90$`\pm `$8%, a fortuitously large value since if the pulsed fraction were more like the other two, we would not have detected the periodicity. This source, with only 470 total photons considered for spectroscopy, has a spectral shape which is more difficult to quantify. Table 2 lists the power-law fit parameters. The spectrum is slightly softer than the 152 second pulsar, but the power-law spectral shape is not the best fit. A blackbody fit gives a temperature of 1.0 keV, consistent with it being a hard spectrum.

As with the 152 second pulsar, this previously known X-ray emitter was postulated to be an X-ray binary pulsar by Haberl & Sasaki and by Kahabka & Pietsch (1996) based upon correlations of the ROSAT X-ray source, RX J0101.0-7206, with an H-alpha emission source. The putative ROSAT counterpart is only 0.96 arcseconds from the CHANDRA detection, well within the combined source location errors for the two detectors indicating they are the same source. The most remarkable thing about this source is the extremely weak flux. Assuming a distance of 57 kpc for the SMC, the X-ray luminosity is only 1.1x1034 ergs/s. This is a very small discovery flux, making it the weakest accreting SMC source for which X-ray pulsations have been discovered.

### 3.3 The 565 second pulsar: CXOU J005736.2-721934

The pulse profile of CXOU J005736.2-721934 is shown in Figure 2. The pulse profile is quite complicated, showing three distinct peaks, a characteristic indicated by the strong presence of the first three harmonics in the Fourier transform spectrum. This source is also formally the most variable of the three with a $`\chi ^2`$ test for variability based upon binning the flux over a single pulse period shows that the flux is consistent with being constant at the 3% level unlike the other two sources, which are fully consistent with constant flux. This is not, however, conclusive evidence for pulse to pulse variability.

Unlike the previous two sources, there has been no previous suggestion that this is an accreting X-ray pulsar. This source is only 18 arcseconds away, however, from source 19 of the ASCA SMC source catalog, well within the half-power radius of the ASCA XRT (Yokogawa et al. 2000). Its pulsed fraction is 48$`\pm `$5%, and it is has the softest spectrum of the three sources. There are no reported ROSAT detections of this source. A search of stellar databases does show that there is an H-alpha emission line source located within 0.8 arcsecond of the CHANDRA position (Meyssonnier & Azzopardi - Star 1020 with a J2000 location of 00 57 36.0 -72 19 34. This lends support to the interpretation of this source as yet another SMC Be-binary.

## 4 DISCUSSION

There are two independent arguments that these pulsars have high mass companions. The first is their apparent association with H-alpha sources which are themselves markers of regions containing massive stars (Kennicutt 1983).

We have assessed the possibility that the apparent association of these 3 pulsars with known H-alpha sources is due to an accidental overlap in position. A cross-check of the Chandra positions with the H-alpha catalogs of Meyssonnier & Azzopardi (1996) and Murphy & Bessell (2000) shows that of the 140 Chandra sources in our field, only 4 were within 5 arcseconds of an H-alpha source. The three closest coincidences were for these 3 pulsars - CXOU J005736.2-721934 (565 s) at a distance of 1.64”, CXOU J010102.7-720658 (304 s) at a distance of 1.14”, and CXOU J005750.3-720756 (152 s) at a distance of 0.67”. The binomial probability of these three particular sources having the smallest 3 distances to an H-alpha source is $`10^6`$. Meyssonnier & Azzopardi do not quote precise error radii for their sources, but indicate that overall “seeing” for the observations was on the order of 1 arcsec. In addition, the CHANDRA positions are expected to be good to within about 0.6 arcsec. Thus, the distances for all three pulsars from the putative Halpha source are within a sum of the estimated Chandra and Meyssonnier & Azzopardi error radii ($``$1.6”).

The second argument which supports the identification of these pulsars as high mass systems is given by the empirical correlation noted by Corbet (1986) which, for X-ray pulsars with spin periods in excess of 100 s, implies a companion which is either a Supergiant or a Be star.

Thus, like virtually all other SMC pulsars (31 total), these three are associated with high mass binaries. In contrast, from the Milky Way with approximately 100 times the mass of the SMC only about 40 such objects have been detected. This overabundance has been discussed in some detail by Schmidtke et al. 1999 and Yokogawa et al. 2000. This overabundance points to a rather recent star formation activity $``$ 107 yr ago.

Further evidence for this episode of star burst activity comes from examination of the spatial distribution of the SMC X-ray binaries. Maragoudaki et al. (2001) have recently analyzed the star content of the SMC using UV photometry over a large (6 by 6) field. The distribution of the SMC X-ray binaries follows very closely (with 3 exceptions) the their distribution of the youngest stars, those having ages less than 8$`\times `$106 yr (figure 9, Maragoudaki et al.) This distribution differs considerably from the general distribution of 517 SMC X-ray sources catalogued by Haberl et al. (2000). For example, 29% of the X-ray sources (149 of 517) fall in the south-eastern region of the SMC (RA greater than 01h and Dec. less than -72.5). On the other hand only 10% (3 of 31) of the X-ray pulsars are located in this region. This is another reflection of the fact that the population of the SMC pulsars are prominently (if not, exclusively) high mass and therefore generally younger then other types of X-ray sources. Further detailed comparisons between the star maps of Maragoudaki et al. and the distribution of the SMC X-ray binaries may provide useful constraints on evolutionary scenarios of high mass X-ray binaries.

Despite the relatively poor timing resolution (3.241 s) of this particular observation, the single 100 ksec exposure reveals four new X-ray pulsars, increasing the number of known SMC pulsars by over 10%. A longer exposure with the Chandra instruments at much higher time resolution would probably provide even more discoveries and allow important assessments of the underlying distributions of SMC pulsars.

This work made use of software and data provided by the High-Energy Astrophysics Archival Research Center (HEASARC) located at Goddard Space Flight Center. Additional use was made of the SIMBAD catalogue at CDS, Strasbourg France, and NASA’s Astrophysics Data System (ADS) .


text:
# Spectroscopic analysis of newborn massive stars in SMC N81

## 1. Introduction

N81 belongs to the class of the High Excitation Blobs (HEB) first introduced by Heydari-Malayeri & Testor (1982) (see Heydari-Malayeri et al., these proceedings). It displays cavities, shocks, ionisation ridges, turbulent structures typical of massive star forming environments. Several stars just emerging from their parental cloud are grouped within the 2” constituting the core region. Thus, N81 represent a unique opportunity to study both the earliest phases of the evolution of massive stars and the metallicity dependence of their wind properties.

## 2. Results

A qualitative analysis of the STIS spectra (Fig. 1) reveals that the N81 stars are O dwarfs (presence of NV($`\lambda 1240`$), OV($`\lambda 1371`$), CIV($`\lambda 1550`$) and HeII($`\lambda 1640`$); absence of strong SiIV($`\lambda 1400`$)) probably lying near the ZAMS as indicated by their low luminosity and the weakness of their wind features (see Heydari-Malayeri et al 2002) : they likely belong to the Vz class (Walborn & Parker 1992).

The preliminary quantitative analysis of one of the two brightest stars of N81 using the non-LTE spherically expanding line blanketed code CMFGEN (Hillier & Miller 1998) gives : $`T_{eff}=40000`$ K, $`\mathrm{log}g=4.1`$, $`\mathrm{log}(L/L_{})=5.15`$, $`\mathrm{log}\dot{M}=9.0`$, $`v_{\mathrm{}}=2000`$ km/s. Fig. 1 shows the wind momentum - luminosity relation (WLR) for O stars : the N81 star has a wind $``$ 2 orders of magnitude weaker than normal dwarfs. This may indicate either a breakdown of the WLR at low luminosity, or a steeper slope at low metallicity.

## 3. Conclusion

Work is in progress to improve our determination of the parameters of these N81 O stars. If confirmed, such weak winds for O stars are puzzling because they are not consistent with the prediction of radiation driven winds theory (e.g. Vink et al. 2000, Hoffmann et al, Puls, these proceedings).

## References

Crowther, P.A. et al., 2002, ApJ, submitted

Herrero, A., Puls, J., Villamariz, M. R., 2000, A&A, 354, 193

Heydari-Malayeri, M., Testor, G., 1982, A&A, 111, L11

Heydari-Malayeri, M. et al., 2002, A&A, 381, 951

Hillier, D. J., Miller, D. L., 1998, ApJ, 496, 407

Hillier, D.J. et al., 2002, ApJ, submitted

Kudritzki, R.P., Puls, J., 2000, ARA&A, 38, 613

Lamers, H.J.G.L.M., Casinelli, J.P., 1996, ASP conference serie, 98, 162

Puls, J. et al., 1996, A&A, 305, 171

Vink, J., de Koter, A., Lamers, H.J.G.L.M., 2000, A&A, 362, 295

Walborn, N. R., Parker, J. W., 1992, ApJ, 399, L87


text:
# 1 Introduction.

## 1 Introduction.

The present article is a continuation of . The notation of  will be used. Now, we will solve numerically six separate monocriterial optimization problems O1-O6. We assumed the data from real street canyon of the Krasiński Avenue in Cracow . We solve the set of nonlinear partial differential equations E1-E8 with given boundary B0-B8, and the initial conditions C0-C8, and sources D0-D2, by the finite difference method, using the C language programme written by the author. We solve this set in the cuboid $`\mathrm{\Omega }`$, starting from initial conditions, and we iterate it over the time period using the direct finite difference method taking into account the boundary B0-B8 conditions, and the sources D0-D2, and initial conditions C0-C8, at each time step. The functionals in F1-F6 are iterated with the same steps as the equations E1-E8 are iterated. The first time derivative is approximated by the first differential quotient using forward two-point first difference in the direction of time coordinate, whereas the spatial first derivatives are approximated by the first differential quotients using central three-point first differences in the directions of spatial coordinates. In general, the numerical results are in good agreement with the measured data from .

## 2 General Traffic Management Inferences.

On the basis of the performed simulations, we deduce the traffic management inferences. In general, some vehicular traffic and hydrodynamical parameters influence the solutions of optimization problems O1-O6, but not all of them.

R1. The direction of velocity of air mixture is important. The optimal pollutant concentrations for both the single canyon and the canyon with its two substitute nearest neighbour canyons, O3, O6, are the lowest ones if the velocity components of the boundary and the initial value problems are equal to zero; further, they are grater if the velocity has only nonzero vertical component $`v_z`$; then, they are greater if velocity has only nonzero x-component $`v_x`$; further, they are again greater if the velocity has two nonzero components $`(0,v_y,v_z)`$; next, they are again greater if the velocity has three nonzero components $`(v_x,v_y,v_z)`$; finally, they are the greatest if the velocity has only nonzero y-component $`v_y`$. We infer that in case R1.6, the optimal green time $`g_{2,\mathrm{C}}^{}`$ is different from other cases. It means that velocity of the mixture influences optimal vectors of control.

R2. The optimal values for cases O1, O4, for cases O2, O5, and for cases O3, O6, decrease with “uniformization” of the vehicles, when we pass from nonuniform vehicles $`\mathrm{UNIFORM}=0`$ to uniform ones $`\mathrm{UNIFORM}=1`$. It is a result of reduction of the number of vehicles moving in the canyon. For $`\mathrm{UNIFORM}=1`$ the values of maximum free flow speed, jam, saturation, threshold, green, and red densities reach the minimal values .

R3. The long vehicular queues decrease the total travel times O1, O4, and they increase both the optimal emissions O2, O5, and optimal concentrations of pollutants O3, O6 . The decrement of total travel times O1, O4, with long vehicular queues is a result of clustering of the vehicles.

R4. The constant of temperature scale $`\mathrm{TH}`$ does not differentiate the values of optimal concentrations of pollutants O3, O6 in the temperature range near standard temperature and pressure STP conditions. However, it diminishes them even hundredfold for very high temperatures .

R5. The functional form of initial and boundary conditions affects the optima. If they are constant, then optimal concentrations O3, O6, are two times higher than in the case when they are changing exponentially in space in three dimensions. The index $`\mathrm{FORM}=\mathrm{CONSTANT}`$ corresponds to constant initial and boundary conditions, whereas $`\mathrm{FORM}=\mathrm{EXPONENTIAL}`$ corresponds to the spatially exponentially changing initial and boundary conditions. For the latter case, we chose the following functional form of initial conditions: the temperature is equal to $`T_0(x,y,z,t)=\mathrm{exp}(\frac{x}{a}\frac{y}{b}\frac{z}{c})\mathrm{TH}`$, the only non-zero coordinates of velocities of mixture are the x-coordinates and they are equal to $`v_{x,0}(x,y,z,t)=\mathrm{exp}(\frac{x}{a}\frac{y}{b}\frac{z}{c})\mathrm{VX}`$, the density is equal to $`\rho _0(x,y,z,t)=\mathrm{exp}(\frac{x}{a}\frac{y}{b}\frac{z}{c})\rho _{\mathrm{STP}}`$, the concentrations are equal to $`c_{i,0}(x,y,z,t)=\mathrm{exp}(\frac{x}{a}\frac{y}{b}\frac{z}{c})c_{i,\mathrm{STP}}`$, and the pressure is equal to $`p_0(x,y,z,t)=\mathrm{exp}(\frac{x}{a}\frac{y}{b}\frac{z}{c})p_{\mathrm{STP}}`$. The boundary conditions B0-B8 are similar.

R6. The presence of vehicles on both the left and right lanes is important. The optimal total travel times and emissions are halved in absence of vehicles on the left or right lanes when compared to situation when they circulate on both the left and right lanes .

R7. The values of saturation, arrival, or jam vehicular density, and of vehicular free flow velocities also affect the optima O1-O6 .

R8. The assumption of energy conservation equation, of the thermodiffusion effect, of the chemical potential and of the Grand Canonical ensemble, as well as of influence of gravity on intrinsic energy and on the chemical potential, drastically changes the optimal concentrations O3, O6, towards the measured ones .

R9. The value of time of simulation and of discretization in time considerably affects the optima (compare ). The values of optimal solutions O1-O6, increase from tenfold to hundredfold. Also the optimal vectors (5-tuples) of control for O1-O6, change their values. It is a result of cumulative effect of length of the period of simulation $`T_S`$ on integral functionals F1-F6.

We use the following notation: $`\mathrm{yL}=0/\mathrm{yR}=0`$ means that the lengths of all queues on the left/right lanes at the beginning were equal to zero, $`\mathrm{VX}=0`$ is meant for zero initial and boundary velocity of mixture, while $`\mathrm{VX}=1`$ is put for the velocity of $`1[\mathrm{m}\mathrm{s}^1]`$. The same holds for $`\mathrm{VY}`$ and $`\mathrm{VZ}`$. If $`\mathrm{VX}0`$, then the left lanes are leeward and the right ones windward. If there are no vehicles on the left or right lanes, then $`\mathrm{LON}=0`$ or $`\mathrm{RON}=0`$, respectively. $`\mathrm{UNIFORM}=0`$ stands for non-homogeneous (different) values of maximum free flow speed $`w_{vt,f}`$, jam $`k_{vt,\mathrm{jam}}`$, saturation $`k_{vt,\mathrm{sat}}`$, threshold $`k_{vt,\mathrm{threshold}}`$, green $`k_{vt,\mathrm{GREEN}}`$, and red $`k_{vt,\mathrm{RED}}`$, vehicular densities for $`VT=4`$ types of vehicles: passenger cars, 8-ton, 12-ton, and 16-ton trucks, where threshold vehicular densities $`k_{vt,\mathrm{threshold}}`$ are the numerical values for which we assume in numerical simulations that there is a vehicle on the given lane. We also assume that these values are the same for both the left and right lanes. $`\delta _x,\delta _y,\delta _z,\delta _t`$, are the steps in $`x,y,z,t,`$ directions, respectively. $`\delta _{C_1},\delta _{g_1},\delta _{C_2},\delta _{g_2},\delta _F`$, are the steps in $`C_1,g_1,C_2,g_2,F`$, directions, respectively. In the boundary and initial conditions B0a, B0b, B0c, C0, we assume that the temperature is equal to a given constant $`\mathrm{TH}`$.

## 3 Discretization scheme.

In order to solve numerically the set of nonlinear partial differential equations E1-E8, one uses the apparatus of finite differences. Firstly, we discretize the points of domain $`\mathrm{\Sigma }`$ in the standard way:

$$(x,y,z,t)(x_i,y_j,z_k,t_l),i=0,\mathrm{},N_x,j=0,\mathrm{},N_y,k=0,\mathrm{},N_z,l=0,\mathrm{},N_t,$$
(3.1)
$$x_i=x_0+\frac{i1}{N_x}(x_{N_x}x_0),$$
(3.2)
$$y_i=y_0+\frac{j1}{N_y}(y_{N_y}y_0),$$
(3.3)
$$z_k=z_0+\frac{k1}{N_z}(z_{N_z}z_0),$$
(3.4)
$$t_l=t_0+\frac{l1}{N_t}(t_{N_t}t_0).$$
(3.5)
Secondly, the value of given function $`f`$ at point $`(x,y,z,t)`$ of its domain $`\mathrm{\Sigma }`$ is approximated by:

$$f(x,y,z,t)f(x_i,y_i,z_k,t_l)f_{i,j,k,l}.$$
(3.6)
Thirdly, the first finite differences in spatial and temporal directions:

$`\mathrm{\Delta }_xf_{i,j,k,l}={\displaystyle \frac{1}{2}}(f_{i+1,j,k,l}f_{i1,j,k,l}),`$ (3.7)
$`\mathrm{\Delta }_yf_{i,j,k,l}={\displaystyle \frac{1}{2}}(f_{i,j+1,k,l}f_{i,j1,k,l}),`$ (3.8)
$`\mathrm{\Delta }_zf_{i,j,k,l}={\displaystyle \frac{1}{2}}(f_{i,j,k+1,l}f_{i,j,k1,l}),`$ (3.9)
$`\mathrm{\Delta }_tf_{i,j,k,l}=f_{i,j,k,l+1}f_{i,j,k,l},`$ (3.10)
allow to approximate the first partial derivatives:

$`{\displaystyle \frac{f(x,y,z,t)}{x}}{\displaystyle \frac{\mathrm{\Delta }_xf_{i,j,k,l}}{\mathrm{\Delta }x}},`$ (3.11)
$`{\displaystyle \frac{f(x,y,z,t)}{y}}{\displaystyle \frac{\mathrm{\Delta }_yf_{i,j,k,l}}{\mathrm{\Delta }y}},`$ (3.12)
$`{\displaystyle \frac{f(x,y,z,t)}{z}}{\displaystyle \frac{\mathrm{\Delta }_zf_{i,j,k,l}}{\mathrm{\Delta }z}},`$ (3.13)
$`{\displaystyle \frac{f(x,y,z,t)}{t}}{\displaystyle \frac{\mathrm{\Delta }_tf_{i,j,k,l}}{\mathrm{\Delta }t}}.`$ (3.14)
Since the vector and scalar partial differential equations E1-E8 are of the following form:

$$\frac{f^\alpha (x,y,z,t)}{t}=F^\alpha (x,y,z,t,f^\alpha (x,y,z,t),\frac{f^\alpha (x,y,z,t)}{x_m},\frac{^2f^\alpha (x,y,z,t)}{x_mx_n},f^\beta (x,y,z,t)),$$
(3.15)
hence we approximate them as follows:

$$\frac{\mathrm{\Delta }_tf_{i,j,k,l}^\alpha }{\mathrm{\Delta }t}=F^\alpha (x_i,y_j,z_k,t_l,f_{i,j,k,l}^\alpha ,\frac{\mathrm{\Delta }_{x_m}f_{i,j,k,l}^\alpha }{\mathrm{\Delta }x_m},\frac{\mathrm{\Delta }_{x_m}\mathrm{\Delta }_{x_n}f_{i,j,k,l}^\alpha }{\mathrm{\Delta }x_m\mathrm{\Delta }x_n},f_{i,j,k,l}^\beta ).$$
(3.16)
Hence, we can rewrite it in following forward scheme:

$$f_{i,j,k,l+1}^\alpha =f_{i,j,k,l}^\alpha +\mathrm{\Delta }tF^\alpha (x_i,y_j,z_k,t_l,f_{i,j,k,l}^\alpha ,\frac{\mathrm{\Delta }_{x_m}f_{i,j,k,l}^\alpha }{\mathrm{\Delta }x_m},\frac{\mathrm{\Delta }_{x_m}\mathrm{\Delta }_{x_n}f_{i,j,k,l}^\alpha }{\mathrm{\Delta }x_m\mathrm{\Delta }x_n},f_{i,j,k,l}^\beta ).$$
(3.17)
We give now the examples of dicretized equations. The discretization of continuity equation E2 reads:

$`\rho _{i,j,k,l+1}=\rho _{i,j,k,l}{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}}v_{i,j,k,l}^x(\rho _{i+1,j,k,l}\rho _{i1,j,k,l})`$ (3.18)
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }y}}v_{i,j,k,l}^y(\rho _{i,j+1,k,l}\rho _{i,j1,k,l})`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }z}}v_{i,j,k,l}^z(\rho _{i,j,k+1,l}\rho _{i,j,k1,l})`$
$`+\mathrm{\Delta }t{\displaystyle \frac{T_S}{\rho _0}}S_{i,j,k,l}`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}}\rho _{i,j,k,l}(\rho _{i+1,j,k,l}\rho _{i1,j,k,l})`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }y}}\rho _{i,j,k,l}(\rho _{i,j+1,k,l}\rho _{i,j1,k,l})`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }z}}\rho _{i,j,k,l}(\rho _{i,j,k+1,l}\rho _{i,j,k1,l}).`$
For the x-component of Navier-Stokes equation E1 we have:

$`v_{i,j,k,l+1}^x=v_{i,j,k,l}^x{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}}v_{i,j,k,l}^x(v_{i+1,j,k,l}^xv_{i1,j,k,l}^x)`$ (3.19)
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }y}}v_{i,j,k,l}^y(v_{i,j+1,k,l}^xv_{i,j1,k,l}^x)`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }z}}v_{i,j,k,l}^z(v_{i,j,k+1,l}^xv_{i,j,k1,l}^x)`$
$`\mathrm{\Delta }t{\displaystyle \frac{T_S}{\rho _0}}{\displaystyle \frac{v_{i,j,k,l}^x}{\rho _{i,j,k,l}}}S_{i,j,k,l}`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}}BT_0{\displaystyle \frac{T_S^2}{a^2}}(T_{i+1,j,k,l}T_{i1,j,k,l})`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}}BT_0{\displaystyle \frac{T_S^2}{a^2}}{\displaystyle \frac{T_{i,j,k,l}}{\rho _{i,j,k,l}}}(\rho _{i+1,j,k,l}\rho _{i1,j,k,l})`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{(\mathrm{\Delta }x)^2}}{\displaystyle \frac{\eta T_S}{\rho _0a^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i+1,j,k,l}^x2v_{i,j,k,l}^x+v_{i1,j,k,l}^x)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{(\mathrm{\Delta }y)^2}}{\displaystyle \frac{\eta T_S}{\rho _0b^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i,j+1,k,l}^x2v_{i,j,k,l}^x+v_{i,j1,k,l}^x)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{(\mathrm{\Delta }z)^2}}{\displaystyle \frac{\eta T_S}{\rho _0c^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i,j,k+1,l}^x2v_{i,j,k,l}^x+v_{i,j,k1,l}^x)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{(\mathrm{\Delta }x)^2}}(\xi +{\displaystyle \frac{\eta }{3}}){\displaystyle \frac{\eta T_S}{\rho _0a^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i+1,j,k,l}^x2v_{i,j,k,l}^x+v_{i1,j,k,l}^x)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{4\mathrm{\Delta }x\mathrm{\Delta }y}}(\xi +{\displaystyle \frac{\eta }{3}}){\displaystyle \frac{\eta T_S}{\rho _0b^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i+1,j+1,k,l}^xv_{i+1,j1,k,l}^xv_{i1,j+1,k,l}^x+v_{i1,j1,k,l}^x)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{4\mathrm{\Delta }x\mathrm{\Delta }z}}(\xi +{\displaystyle \frac{\eta }{3}}){\displaystyle \frac{\eta T_S}{\rho _0c^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i+1,j,k+1,l}^xv_{i+1,j,k1,l}^xv_{i1,j,k+1,l}^x+v_{i1,j,k1,l}^x)`$
$`+\mathrm{\Delta }t{\displaystyle \frac{g_0T_S^2}{a}}g^x.`$
Similarly for the y-component:

$`v_{i,j,k,l+1}^y=v_{i,j,k,l}^y{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}}v_{i,j,k,l}^x(v_{i+1,j,k,l}^yv_{i1,j,k,l}^y)`$ (3.20)
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }y}}v_{i,j,k,l}^y(v_{i,j+1,k,l}^yv_{i,j1,k,l}^y)`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }z}}v_{i,j,k,l}^z(v_{i,j,k+1,l}^zv_{i,j,k1,l}^z)`$
$`\mathrm{\Delta }t{\displaystyle \frac{T_S}{\rho _0}}{\displaystyle \frac{v_{i,j,k,l}^y}{\rho _{i,j,k,l}}}S_{i,j,k,l}`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }y}}BT_0{\displaystyle \frac{T_S^2}{b^2}}(T_{i,j+1,k,l}T_{i,j1,k,l})`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }y}}BT_0{\displaystyle \frac{T_S^2}{b^2}}{\displaystyle \frac{T_{i,j,k,l}}{\rho _{i,j,k,l}}}(\rho _{i,j+1,k,l}\rho _{i,j1,k,l})`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{(\mathrm{\Delta }x)^2}}{\displaystyle \frac{\eta T_S}{\rho _0a^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i+1,j,k,l}^y2v_{i,j,k,l}^y+v_{i1,j,k,l}^y)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{(\mathrm{\Delta }y)^2}}{\displaystyle \frac{\eta T_S}{\rho _0b^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i,j+1,k,l}^y2v_{i,j,k,l}^y+v_{i,j1,k,l}^y)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{(\mathrm{\Delta }z)^2}}{\displaystyle \frac{\eta T_S}{\rho _0c^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i,j,k+1,l}^y2v_{i,j,k,l}^y+v_{i,j,k1,l}^y)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{4\mathrm{\Delta }x\mathrm{\Delta }y}}(\xi +{\displaystyle \frac{\eta }{3}}){\displaystyle \frac{\eta T_S}{\rho _0b^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i+1,j+1,k,l}^xv_{i+1,j1,k,l}^xv_{i1,j+1,k,l}^x+v_{i1,j1,k,l}^x)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }y^2}}(\xi +{\displaystyle \frac{\eta }{3}}){\displaystyle \frac{\eta T_S}{\rho _0b^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i+1,j+1,k,l}^yv_{i+1,j1,k,l}^yv_{i1,j+1,k,l}^y+v_{i1,j1,k,l}^y)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{4\mathrm{\Delta }y\mathrm{\Delta }z}}(\xi +{\displaystyle \frac{\eta }{3}}){\displaystyle \frac{\eta T_S}{\rho _0b^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i,j+1,k+1,l}^xv_{i,j+1,k1,l}^xv_{i,j1,k+1,l}^x+v_{i,j1,k1,l}^x)`$
$`+\mathrm{\Delta }t{\displaystyle \frac{g_0T_S^2}{b}}g^y.`$
Finally, for the z-component:

$`v_{i,j,k,l+1}^z=v_{i,j,k,l}^z{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}}v_{i,j,k,l}^x(v_{i+1,j,k,l}^zv_{i1,j,k,l}^z)`$ (3.21)
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }y}}v_{i,j,k,l}^y(v_{i,j+1,k,l}^zv_{i,j1,k,l}^z)`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }z}}v_{i,j,k,l}^z(v_{i,j,k+1,l}^zv_{i,j,k1,l}^z)`$
$`\mathrm{\Delta }t{\displaystyle \frac{T_S}{\rho _0}}{\displaystyle \frac{v_{i,j,k,l}^z}{\rho _{i,j,k,l}}}S_{i,j,k,l}`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }z}}BT_0{\displaystyle \frac{T_S^2}{c^2}}(T_{i,j,k+1,l}T_{i,j,k1,l})`$
$`{\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }z}}BT_0{\displaystyle \frac{T_S^2}{c^2}}{\displaystyle \frac{T_{i,j,k,l}}{\rho _{i,j,k,l}}}(\rho _{i,j,k+1,l}\rho _{i,j,k1,l})`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{(\mathrm{\Delta }x)^2}}{\displaystyle \frac{\eta T_S}{\rho _0a^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i+1,j,k,l}^z2v_{i,j,k,l}^z+v_{i1,j,k,l}^z)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{(\mathrm{\Delta }y)^2}}{\displaystyle \frac{\eta T_S}{\rho _0b^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i,j+1,k,l}^z2v_{i,j,k,l}^z+v_{i,j1,k,l}^z)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{(\mathrm{\Delta }z)^2}}{\displaystyle \frac{\eta T_S}{\rho _0c^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i,j,k+1,l}^z2v_{i,j,k,l}^z+v_{i,j,k1,l}^z)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{4\mathrm{\Delta }x\mathrm{\Delta }z}}(\xi +{\displaystyle \frac{\eta }{3}}){\displaystyle \frac{\eta T_S}{\rho _0c^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i+1,j,k+1,l}^xv_{i+1,j,k1,l}^xv_{i1,j,k+1,l}^x+v_{i1,j,k1,l}^x)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{4\mathrm{\Delta }y\mathrm{\Delta }z}}(\xi +{\displaystyle \frac{\eta }{3}}){\displaystyle \frac{\eta T_S}{\rho _0b^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i,j+1,k+1,l}^xv_{i,j,k1,l}^xv_{i,j,k+1,l}^x+v_{i,j1,k1,l}^x)`$
$`+{\displaystyle \frac{\mathrm{\Delta }t}{(\mathrm{\Delta }z)^2}}(\xi +{\displaystyle \frac{\eta }{3}}){\displaystyle \frac{\eta T_S}{\rho _0c^2}}{\displaystyle \frac{1}{\rho _{i,j,k,l}}}(v_{i,j,k+1,l}^z2v_{i,j,k,l}^z+v_{i,j,k1,l}^z)`$
$`+\mathrm{\Delta }t{\displaystyle \frac{g_0T_S^2}{c}}g^z.`$
## 4 Method of comparison.

We briefly describe the method of comparison of the numerical results with the experimental ones. The comparison of simulations with the measured data, if authorized, is aimed at showing the correctness of the model (conservation of the rank of simulated parameters), and of presenting the correct direction of deviations (the deviations always tend to secure the direction). We transform both measurements and optima $`J_\mathrm{C}^{}`$ O3, in order to compare them. The measured pollutant concentrations are multiplied by the volume of the canyon, and we obtain pollutant masses (we assume that the measured concentrations are the same in the entire canyon and in time). The pollutant masses are added up, assuming $`\mathrm{HC}`$ mass to be equal to zero ($`\mathrm{HC}`$ were not measured in Ref. ). The calculated optima $`J_\mathrm{C}^{}`$ O3 are divided by time of simulation $`T_S`$, giving approximation of pollutant masses in every moment of simulation (we assume that $`J_\mathrm{C}^{}`$ O3 are homogeneous in time, in order to make possible the comparison with measurements). Thence, we have two comparable magnitudes. Also, only insufficient vehicular flow data are given in Ref. . Neither the jam, saturation vehicular densities, nor the maximum vehicular speeds are present. However, the initial and boundary conditions of traffic flow B5-B8, C5-C8 are not constant functions because they are parametrized by measured traffic parameters. The data are not instantaneous, but they are measured in intervals of hours. Moreover, the measured data are not accompanied by measurements of wind velocity, temperature, density. If the measurements had been performed instantly using more specialized equipment (e. g., lidar), and if they had been stored in computer facility, then, it would be possible to make more serious comparisons. The simulations were performed for the time period $`\mathrm{T}=60[\mathrm{s}]`$, whereas the measurement were averaged over the time period $`\mathrm{T}=1800[\mathrm{s}]`$. The integral functionals F3, F6, are cumulative with respect to the time period.

## 5 Computer simulations.

From thousands of performed optimizations O1-O6 we selected some representable optimizations. The C language programme of more than 14000 lines written by the author was being run on workstations of Hewlett-Packard, Sun, Silicon Graphics, under UNIX operating system. It was being run also on Pentium personal computers under Linux operating system. We used C language UNIX/Linux compilers cc. The time of simulations varied from a couple of hours to several days or even weeks, depending on the discretization parameters. The discretization steps were tested in order to obtain finite solutions to the optimization problems in reasonable time of simulation. They were also studied from the point of view of sensitiveness of optima appearing on them.

## 6 Mathematical questions.

The problem of uniqueness of solutions of optimization problems O1-O6 is complex. If there are no vehicles allowed in the canyon, then all the six functionals O1-O6 are constant functions of parameters, the optima exist, and the number of optimal solutions is continuum. The values of functionals are then zero for total travel time and emissions, and are positive for concentrations. These results are analytical. For other cases there are no such analytical results for both the uniqueness and existence known to the author. The numerical solutions of optimization problems are approximately global within the error induced by discretization of the physical domain $`\mathrm{\Sigma }`$ and of the domain of control parameters $`U^{\mathrm{adm}}`$. The method of optimization was a full search.

## 7 Conclusions.

The proecological traffic control idea and advanced model of the street canyon have been developed. It has been found that the proposed model represents the main features of complex air pollution phenomena.

## 8 Acknowledgements.

It is my pleasure to thank Professor Andrzej Adamski for formulating the problem, Professors Stanisław Białas, Wojciech Mitkowski, and Gwidon Szefer, for many constructive criticisms. I also thank Professor Włodzimierz Wójcik for his giving me access to his computer facilities.