text:
# Tiling space and slabs with acute tetrahedra

## 1 Problem definition

Triangulations of two and three-dimensional domains find numerious applications in scientific computing, computer graphics, solid modeling and medical imaging. Most of these applications impose a quality constraint on the elements of the triangulation. Among the most popular quality criteria for elements  are the aspect ratio (circumradius over inradius), the minimum dihedral angle, and the radius-edge ratio (circumradius over shortest edge). However, many other quality criteria have been considered, including maximum dihedral angle. Bern et al. for instance, studied nonobtuse triangulations , where domains are meshed with simplices having no obtuse angles. In this paper, we consider a slightly stronger quality constraint: all the dihedral angles in the mesh are forced to be acute (strictly less than $`90^{}`$). Although acuteness seems only slightly stronger than nonobtuseness, this problem turns out to be considerably harder than the nonobtuse triangulation problem, as we observe below in Section 3.

###### Definition.

An angle is acute if it is strictly less than a right angle ($`\frac{\pi }{2}=90^{}`$). A simplex is acute if all its (interior) dihedral angles are acute. A triangulation is acute if all of its simplices are acute.

###### Problem 1.

Given a domain $`\mathrm{\Omega }`$, compute an acute triangulation of $`\mathrm{\Omega }`$.

There has been extensive work on the two-dimensional version of this problem, for the special cases where the domain $`\mathrm{\Omega }`$ is a triangle, square, quadrilateral, or a finite point set . We review those results in Section 3. In three-dimensional space, however, almost nothing has been known about acute triangulations before now. To the best of our knowledge, even the following relaxed form of the problem, where the the input domain is the entire space, had not been addressed in the literature.

###### Problem 2.

Is it possible to tile three-dimensional Euclidean space using acute tetrahedra?

We present an affirmative answer to this question, by several different constructions. The two-dimensional analog of this problem has a trivial positive answer: congruent copies of any single triangle will tile the plane. However, this idea does not extend to 3D, as the regular tetrahedron (for instance) cannot tile space. All tetrahedra known to tile space have right angles, as further discussed in Section 3.

We started this research on acute triangulations because of a method developed for space-time meshing which required an acute base mesh. This and our other motivations are discussed in Section 2. Section 3 surveys previous research in acute triangulations. Section 4 investigates what acuteness means for a three-dimensional simplex and gives a comparison of acute and Delaunay triangulations. Constructions tiling three-dimensional space, and hence solving Problem 2, are given in Section 5. The paper concludes in Section 6 with a quality assesment of these constructions.

## 2 Motivation

We were originally motivated to study acute triangulations by the space-time meshing algorithm of Üngör and Sheffer . This *tent-pitcher* algorithm was designed to discretize space-time domains into meshes that obey a certain cone constraint, which requires all faces in the mesh to have smaller slopes than the cones that define the domain of influence imposed by the numerical (engineering) problem. (For instance, we might require simply that all faces make at most a $`45^{}`$ angle with the horizontal.) Because there is then a well-defined direction of information flow across element boundaries, such meshes enable the use of very efficient element-by-element methods (including space-time discontinuous Galerkin methods) to solve a wide variety of numerical problems, for instance in elastodynamics. The tent-pitching algorithm starts with a space mesh of the two- or three-dimensional input domain and constructs the space-time mesh using an advancing front approach. The algorithm is known to generate a valid space-time mesh if the initial space mesh is an acute triangulation , but may fail if there is an obtuse angle or even a right angle.

Later, Erickson *et al.*  proposed an improved version of the tent-pitching algorithm. By removing the acute angle requirement, the new space-time algorithm works over arbitrary spatial domains. However, there is a loss of efficiency (more elements are required) whenever there is a nonacute angle.

Thus the study of Problem 1 is motivated by current space-time meshing algorithms. But even a solution to Problem 2 is useful, since it leads to a better understanding of the acute triangulation problem for more general input domains, and it also finds some direct applications in mesh generation.

Spacial tilings of high quality have been used for designing meshing algorithms: Fuchs , Field and Smith  and Naylor  built meshes by overlaying standard tilings onto the given polyhedral domain. They used tilings known at the time, such as Sommerville constructions (Figure 5) and subdivided cubes (Figure 7), which we discuss in Section 3.2. Their approach has three steps, illustrated in Figure 1:

* Overlay the chosen tiling with the given domain. The main challenge in this step is finding the right scaling, location and orientation for the tiling so that it matches the domain boundary as closely as possible.
* Adjust the points to get a better fit. For this purpose, one of the standard smoothing techniques  can be used. Alternatively, Fuchs  suggested minimizing a function which penalizes configurations that produce irregular vertices.
* Construct the mesh by computing the conforming Delaunay triangulation of the adjusted point set and the domain boundary.

Fuchs  reports good performance of his experiments when he used the second Sommerville construction (Figure 5(b)) as the space tiling. (This tiling is the Delaunay triangulation of the body-centered cubic lattice.) The dihedral angles of his mesh in Figure 1(c) range between $`7.6^{}`$ and $`168.2^{}`$. However, most of the angles (here and also in his meshes of similar geometric domains) cluster around $`60^{}`$ and $`90^{}`$, which are exactly the dihedral angles of the BCC tetrahedron in the input tiling. Some of the constructions we propose in Section 5 are considerably better in terms of dihedral angles and also other quality measures. Our new constructions can find immediate use to improve the results of this previous research  on tiling-based meshing.

Bossavit suggests  that acute triangulations may also be useful in computational electromagnetics.

## 3 Background

### 3.1 Acute and nonobtuse triangulations

There has been considerable research  on the nonobtuse triangulation problem, which imposes a slightly weaker constraint than the acute triangulation problem. Angles in a nonobtuse triangulation are less than or equal to $`90^{}`$. Bern et al.  showed that any $`d`$-dimensional point set can be triangulated with $`O(n^{d/2})`$ simplices, none of which has any obtuse dihedral angles. However, they also proved that no similar result is possible if all angles are required to be at most $`90^{}ϵ`$. This indicates that the acute triangulation problem is much more challenging than nonobtuse triangulation. To appreciate this difference, consider the two problems for a square domain in two dimensions. A single diagonal cuts a square into two nonobtuse triangles, as in Figure 2(a). Finding an acute triangulation, however, can be a challenging recreational math problem.

Lindgren  showed that at least eight triangles, as in Figure 2(b), are needed. Later, Cassidy and Lord  showed that for any $`n10`$ (but not for $`n=9`$) there is an acute triangulation with exactly $`n`$ triangles. Figure 2(c) shows the solution with ten triangles. We can use the maximum angle in a triangulation as a quality measure. The triangulations in Figure 2(b,c) can be realized with maximum angles about $`85^{}`$ and $`80.3^{}`$, respectively. Eppstein  improved this angle to $`72^{}`$ using fourteen acute triangles, as shown in Figure 2(d). Using Euler’s formula, Eppstein also showed that any acute triangulation of a square must have an interior vertex of valence five, implying that $`72^{}`$ is the best possible. It is unknown whether there is a triangulation achieving this with fewer than fourteen triangles.

The acute triangulation problem has been studied for other simple polygons as well. Gardner  asked the question for triangles. Manheimer proved that seven acute triangles are necessary and sufficient to subdivide a nonobtuse triangle . Recently, Maehara  showed that an arbitrary quadrilateral can be tiled by $`10`$ (but perhaps not by any fewer) acute triangles. Gerver  considered the problem of finding triangulations with a stricter upper bound (between $`72^{}`$ and $`60^{}`$) on their angles, and gave necessary conditions for a polygonal domain to have such a triangulation. If we restrict ourselves to two-dimensional point sets, a solution to the acute triangulation problem is given by Bern et al. . Their approach starts with a quadtree, and replaces the squares by tiles with protrusions and indentations. Figure 3 shows sample tiles together with an acute triangulation resulting from their algorithm.

Krotov and Krizek  studied refinement methods to subdivide a nonobtuse tetrahedral partition into another finer one. Unfortunately, they called the resulting triangulations acute type instead of nonobtuse even though $`90^{}`$ dihedral angles were ubiquitous in them. Another related work is by Hangan, Itoh and Zamfirescu  who studied acute surface triangulations of certain special shapes such as a cube, sphere and icosahedron.

### 3.2 Acute and nonobtuse tilings

Aristotle claimed that regular tetrahedra could meet five-to-an-edge to tile space, and this claim was repeated over the centuries (see ). This of course is false, because the dihedral angle of a regular tetrahedron is not $`72^{}`$ but $`\mathrm{arccos}\frac{1}{3}70.53^{}`$. Figure 4 shows the small gap left when five tetrahedra are placed around an edge.

There are, however, tetrahedral shapes which can tile space. Sommerville  found four such tetrahedra, shown in Figure 5.

Four decades later, Davies  and Baumgartner  independently rediscovered three of the Sommerville tetrahedra; Baumgartner also found a new example. Goldberg  surveyed the list of all known space-tiling tetrahedra, and found three infinite families, including the one shown in Figure 6(a).

The construction of this family is based on a tiling of the plane by equilateral triangles of side-length $`e`$. The infinite prism over each triangle is filled with tetrahedra whose sides are $`3a,b,b,b,c,c`$, where $`b^2=a^2+e^2`$ and $`c^2=4a^2+e^2`$, as shown in Figure 6. Since the ratio $`a/e`$ is arbitrary, there is a continuous family of tetrahedral space-fillers of this type. Goldberg’s two other families can be derived simply by cutting these tetrahedra into two conguent pieces, either by the plane through $`C`$, $`D`$ and the midpoint of $`\overline{AB}`$, or by the plane through $`A`$, $`B`$ and the midpoint of $`\overline{CD}`$. Whether the list is complete or not is still an open problem . None of the known space-tiling tetrahedra is acute (although several are nonobtuse). In fact, the tilings all contain edges of valence four. (In Goldberg’s family, these are the edges of length $`c`$.)

Since it seems likely that there is no tiling of space by congruent acute tetrahedra, we will now consider tilings with several shapes of tetrahedra. There are now many more ways to fill space, for instance by subdividing the cube into five or six tetrahedra as in Figure 7. These tilings also, of course, have $`90^{}`$ dihedral angles, and so are nonobtuse but not acute.

There are many results (like minmax and maxmin angle results) known about optimality of Delaunay triangulations in the plane. But these do not extend to three dimensions, and little is yet known about optimum triangulations in space. Thus it is not surprising that the constructions of acute triangulations do not easily extend from two to three dimensions. It is perhaps remarkable that acute triangulations of space can be constructed.

## 4 Acute Tetrahedra

An acute tetrahedron does not necessarily have high quality in terms of either aspect ratio or radius-edge ratio. Low-quality tetrahedra have been classified into nine types , and three of these (the spire, splinter and wedge) can have all their dihedral angles acute. However, fortuitously, the tetrahedra in our constructions are mostly quite close to regular, and are high-quality for use in mesh generation and numerical simulations.

### 4.1 Acuteness test

By definition, a tetrahedron is acute if each of its six dihedral angles is less than $`90^{}`$.

###### Lemma 1.

Consider an edge $`ab`$ of a tetrahedron $`abcd`$, and let $`\mathrm{\Pi }`$ denote projection to a plane normal to $`ab`$. The dihedral angle along $`ab`$ is acute if and only if $`\mathrm{\Pi }(a)=\mathrm{\Pi }(b)`$ lies strictly outside the circle with diameter $`\mathrm{\Pi }(c)\mathrm{\Pi }(d)`$.

###### Proof.

The dihedral angle along $`ab`$ is by definition the angle $`\mathrm{}\mathrm{\Pi }(c)\mathrm{\Pi }(a)\mathrm{\Pi }(d)`$; the lemma follows from standard plane geometry (Thales’ theorem, see Figure 8(c)). ∎

This lemma can be applied to each of the edges of a tetrahedron. We now examine some alternate criteria for acuteness.

###### Lemma 2.

A tetrahedron is acute if and only if the orthogonal projection of each vertex onto the plane of the opposite facet lies strictly inside that facet. An acute tetrahedron has acute facets, but this is not a sufficient condition.

###### Proof.

Suppose the projection of a vertex $`d`$ is not inside the opposite triangle $`\mathrm{}abc`$. Then, the dihedral angle along any edge that separates the region it projects to from $`\mathrm{}abc`$ must be nonacute, as in Figure 8(a). (If the projection is on a triangle edge, then the corresponding dihedral angle is exactly $`90^{}`$.) Conversely, if the dihedral angle along edge $`ab`$ is nonacute, then $`d`$ projects outside $`\mathrm{}abc`$, as in Figure 8(b).

We prove the second statement in contrapositive form, while noting that nonacute sliver tetrahedra can have acute face angles. Suppose tetrahedron $`abcd`$ has a nonacute face angle $`\mathrm{}bac`$; we will show the tetrahedron is nonacute. If the projection of $`d`$ onto $`\mathrm{}abc`$ is not in the interior we are done by the first part of the lemma. Otherwise, we claim the dihedral angle along $`ad`$ is larger than $`\mathrm{}bac`$ and thus is nonacute. To check the claim, remember the spherical dual law of cosines (see ):

$$\mathrm{cos}d^{}=\mathrm{cos}b^{}\mathrm{cos}c^{}+\mathrm{sin}b^{}\mathrm{sin}c^{}\mathrm{cos}\mathrm{}bac$$

where $`b^{}`$, $`c^{}`$ and $`d^{}`$ are the dihedral angles along the edges $`ab`$, $`ac`$ and $`ad`$, respectively. (See Figure 8(d).) Assuming $`b^{},c^{}<\frac{\pi }{2}`$, this gives $`\mathrm{cos}d^{}<\mathrm{cos}\mathrm{}bac`$ as desired. ∎

### 4.2 Acuteness of Delaunay Triangulations

Given a set of vertices, the Delaunay triangulation is optimal in many ways. However, a Delaunay triangulation in any dimension can have obtuse angles. In this section, we investigate the converse, whether an acute triangulation is necessarily the Delaunay triangulation for its vertices. The answer is positive in the plane, but negative in three-space.

###### Lemma 3.

Any acute two-dimensional triangulation $`𝒯`$ is Delaunay.

###### Proof.

Since $`𝒯`$ is acute, the diametral circle of each edge is empty of other vertices. By definition, this means the edge is in the Gabriel graph of the vertex set, which is a subgraph of the Delaunay triangulation. But since the edges of $`𝒯`$ form a triangulation, it must be the entire Delaunay triangulation. See also Figure 9. ∎

###### Corollary.

If an acute triangulation of a two-dimensional vertex set exists then it is unique. ∎

###### Lemma 4.

There is an acute triangulation $`𝒯`$ in three dimensions which is not Delaunay.

###### Proof.

Consider a “cube corner” tetrahedron, and glue it to a copy of itself across the equilateral face. Then move the two corner vertices away from each other a tiny amount to make the tetrahedra acute. In coordinates, take a suitable small $`ϵ>0`$, and let $`a=(ϵ,ϵ,ϵ)`$, $`b=(1,0,0)`$, $`c=(0,1,0)`$, $`d=(0,0,1)`$, and $`e=(2/3+ϵ,2/3+ϵ,2/3+ϵ)`$. The two tetrahedra $`abcd`$ and $`bcde`$ are acute, but the Delaunay triangulation of these five points consists of three tetrahedra: $`abce`$, $`acde`$, and $`abde`$, as in Figure 10, because $`e`$ is inside the circumsphere of $`abcd`$. Note that the three Delaunay tetrahedra are obtuse, having $`120^{}`$ dihedral angles along edge $`ae`$. The acute triangulation we started with is obtained by performing a 3-to-2 flip on the Delaunay triangulation. ∎

## 5 Constructions for Acute Tilings

### 5.1 TCP triangulations

Our first set of acute triangulations basically come from the crystallography literature. Chemists studying alloys of two transition metals have often found that since the two types of atoms are similar (but slightly different) in size, the Delaunay triangulation of their positions is built of nearly regular tetrahedra. These TCP (tetrahedrally close packed) structures were first described by Frank and Kasper  and have been studied extensively by the Shoemakers  among others.

A combinatorial definition of the TCP class was given by Sullivan : A triangulation is called TCP if every edge has valence 5 or 6, and no triangle has two 6-valent edges. This definition includes all the chemically known TCP structures, but also allows some new structures  not yet seen in nature.

It is not hard to check that the definition allows exactly four types of vertex star in a TCP triangulation. Dually, the voronoi cell around any vertex has one of the four combinatorial types shown in Figure 11: these are the polyhedra with pentagonal and hexagonal faces but no adjacent hexagons. (It is interesting that these dual structures are seen in some other crystal structures: in some zeolites, silicon dioxide outlines the voronoi edges, while in clathrates, water cages along the voronoi skeleton trap large gas molecules.)

All known TCP structures can be viewed as convex combinations of the three basic ones (called A15, Z and C15) shown in Figure 12.

There are many ways to understand these structures . To construct A15, we can start with a BCC lattice. Its Delaunay triangulation is one of the Sommerville tilings discussed above; since the edges have even valence the tetrahedra can be colored alternately black and white. If we take the BCC lattice together with the circumcenters of all black tetrahedra, we have the vertices of A15: their Delaunay tetrahedra are all now nearly regular. Similarly, the C15 structure arises from the diamond lattice by adding selected circumcenters, and the Z structure can be obtained similarly starting with hexagonal prisms.

The C15 structure (also known as the cubic Friauf–Laves phase) is shown in Figure 13, where the red spheres are centered on a diamond lattice (FCC together with a certain translate) and the blue spheres are at selected circumcenters.

In any triangulation of space, the average dihedral angle multiplied by the average edge valence is exactly $`360^{}`$. If a tiling could be made of regular tetrahedra, the average edge valence would thus be $`n_0:=360^{}/\mathrm{arccos}(\frac{1}{3})5.1043`$. But by symmetry, the regular tetrahedron is a critical point for average dihedral angle, so any tiling made of nearly regular tetrahedra should have average valence quite close to $`n_0`$. Indeed, all known TCP structures have average valence between $`5\frac{1}{10}`$ and $`5\frac{1}{9}`$, the values for C15 and A15.

Sullivan  has formalized a construction suggested by Frank and Kaspar for mixing the basic TCP structures. Start with any tiling of the plane by copies of an equilateral triangle and a square, like one of the four shown in Figure 14.

Suppose the side length of the square and triangle is $`4`$. Mark black and white dots on the tilings as shown. (The dots are at edge midpoints, at triangle centers, and at distance $`1`$ from the sides of the sqaures.) Then the vertices of the corresponding TCP structure are at heights $`4k1`$ above the black dots, at heights $`4k+1`$ above the white dots, and at heights $`4k`$ and $`4k+2`$ above the vertices of the square/triangle tiling. (Here $`k`$ ranges over all integers.)

Again, in each case the TCP triangulation is simply the Delaunay triangulation of this periodic point set. See Figure 16. The triangulations constructed in this way are all combinations of the A15 and Z structures. A variant of this construction  builds combinations of the Z and C15 structures, again starting from an arbitrary square/triangle tiling.

Especially in the mixed structures like $`\sigma `$ and H, the particular geometry we have described here may differ slightly from that found in the actual crystals with the same combinatorics. Presumably, these slight adjustments do not affect the shapes of the tetrahedra very much. The quality figures we present below are measured using the exact geometry we have just described.

Sullivan’s original interest in these structures was for the mathematical study of foam geometry. The Kelvin problem asks for the most efficient partition of space into unit-volume cells, that is, for the partition with least surface area. Lord Kelvin’s suggested solution was a slightly relaxed form of the BCC voronoi cells (truncated octahedra). But in 1994, Weaire and Phelan  discovered that a relaxed form of the voronoi cells for the TCP structure A15 is more efficient than the Kelvin foam .

It is perhaps not surprising that TCP structures are related to foams: Plateau’s rules for singularities in soap films minimizing their surface area imply that a foam is combinatorially dual to some triangulation, preferably one with nearly regular tetrahedra. It is an interesting question whether any triangulation meeting the combinatorial definition of TCP can be built with tetrahedra close to regular, but certainly for the known TCP structures this seems always to be the case. Thus our acute triangulations arising from this construction have high quality by almost all measures.

### 5.2 Icosahedral Construction of the Z Structure

An alternate construction for the TCP Z structure is inspired by the work of Field . His tilings involved right-angled tetrahedra, but by selectively adjusting the point set, we obtain a tiling with only acute tetrahedra. A regular icosahedron can be subdivided into 20 acute (and nearly regular) tetrahedra simply by coning to the center point. We place icosahedra in a hexagonal lattice in the plane, each in the same orientation, touching edge to edge, as in Figure 15.

This layer then gets repeated vertically, with each icosahedron sharing a horizontal face with the ones just above and below it. Our point set is then the vertices and centers of all the icosahedra. Its Delaunay triangulation, shown in Figure 16(d), is combinatorially the TCP Z structure, but with slightly different geometry than that constructed before. The horizontal faces (seen head-on as equilateral triangles in Figure 15) are shared by two icosahedra. Each other face separates an icosahedron from one of the four types of Delaunay tetrahedra that fill the gaps. There are two type of gaps. The deeper gaps are defined by the points $`a,b,c,d,e,f,g,`$ and $`g^{}`$ (opposite to $`g`$) and are filled with two type of tetrahedra, e.g., $`bdfg`$ and $`abfg`$. The shallower gaps are defined by the points $`d,f,h,i,j,k,e,`$ and $`e^{}`$ (opposite of $`e`$) and are filled with two type of tetrahedra, e.g., $`ee^{}df`$ and $`ee^{}fh`$.

### 5.3 An Acute Triangulation of a Single Slab

The acute triangulations we have described so far, though periodic, do not have any planar boundaries within them. Here we describe an acute triangulation of a slab (which can of course be repeated to fill all of space). We view this as partial progress towards the problem of triangulating an arbitrary domain, although it seems much harder to find an acute triangulation of a cube or even an infinite square prism. We triangulate the bottom half of the slab in the following eight steps. Let $`h`$ be the height of the slab and $`\gamma =h/14.2`$.

1. Start with a grid of equilateral triangles of side length $`6\gamma `$ on the base plane, as in Figure 17(a).
2. Place a near-regular tetrahedron (with height $`4\gamma `$) over each triangle, as in Figure 17(b).
3. Add a tetrahedron in the gap between each pair of adjacent tetrahedra, as in Figure 17(c). The resulting surface has deep hexagonal dimples at the original vertices in the base plane.
4. Add six tetrahedra in each dimple, each with one vertex on the starting plane, two vertices at height $`4\gamma `$, and one new vertex at height $`4.6\gamma `$ over the starting vertex, as in Figure 17(d). Now we have a surface with shallow hexagonal bumps.
5. Place a vertex at height at $`7.1\gamma `$ over the midpoint of the edge between each pair of adjacent bumps. Each such vertex and edge form a vertical triangle; let this separate two new tetrahedra whose fourth vertices are the the two nearby bump vertices, as in Figure 17(e). The surface is now covered by tall diamond-shaped bumps.
6. Place a tetrahedron between each adjacent pair of bumps, as in Figure 17(f). We now have an alternating grid of medium-depth six-sided holes (over each of the shallow hexagonal bumps) and deep tetrahedral holes (over the points where three of the shallow hexagonal bumps meet).
7. Fill each tetrahedral hole, to form a surface alternating between six-sided holes and flat triangles, as in Figure 17(g).
8. Place six tetrahedra into each medium-height hexagonal hole to turn it into a medium-height hexagonal bump, as in Figure 17(h). In order to make the bumps equal to the holes, the height of the new vertices is chosen as $`[2(7.14.6)+4.6]\gamma =9.6\gamma `$.

To complete the triangulation of the slab, we now repeat the first seven steps in reverse order.

Any of the constructions given in this section serves to prove our main result:

###### Theorem 5.

It is possible to tile three-dimensional Euclidean space with acute tetrahedra. ∎

## 6 Evaluation of the constructions

Our constructions use tetrahedra of quite good quality (as summarized in Table 1) and so they are quite suitable for mesh generation.

There are still some challenges left to make use of these tilings in real-life meshing techniques. A strategy is required to fit the tilings into a planar projection of the spatial domain. Malkevitch studies this problem in . He describes the conditions for a polygon to be tiled by squares and equilateral triangles. Also, even though one of our constructions fits between two parallel planes in a slab, all of them have dimples (cavities on the surface) in most directrions, making them not very suitable for meshing domains with flat surfaces.

Open problems related to this work include the following

1. Given a point set in 3D that has an acute triangulation, how can we compute it?
2. Are these constructions the best possible? For instance, which tiling of space with tetrahedra minimizes the maximum dihedral angle?
3. Is it possible to tile the space with congruent copies of some single acute tetrahedron?
4. Is it possible to subdivide a cube (or even an acute tetrahedron) into acute tetrahedra?
5. Is it possible to extend the planar acute triangulation algorithm of  into three dimensions, using constructions similar to ours?

### Acknowledgments

Many thanks to Jeff Erickson for discussions and asking how to subdivide a cube into acute tetrahedra, Marshall Bern, Olaf Delgado, Herbert Edelsbrunner, and Sariel Har-Peled for discussions, and Damrong Guoy for his help on visualization.

This work was partially supported by NSF ITR grant DMR-0121695. David Eppstein was supported by NSF grant CCR-9912338. John Sullivan was also partially supported by NSF grant DMS-0071520. Alper Üngör was also partially supported by a UIUC Computational Science and Engineering Fellowship.


format: text

text:
# Timing noise, glitches and the braking index of PSR B0540-69

## 1 Introduction

PSR B0540$``$69 was discovered in the soft X-rays by Seward et al. (1984) in a field of the Large Magellanic Cloud observed with the Einstein Observatory. Pulsations at optical frequencies were soon detected by Middleditch & Pennypacker (1985) with a mean pulsed magnitude of 22.5. In the radio band PSR B0540$``$69 is a quite faint source and pulsed signals were first observed only in 1989-90 (Manchester et al. 1993).

PSR B0540$``$69 is one of the youngest rotation powered pulsars. It has a period of about 50 ms and a large period derivative of $`4.79\times 10^{13}`$ s s-1, comparable to that of the Crab pulsar, and a spin down age of about 1500 years. The pulse shape, at X and optical wavelengths, is broad and almost sinusoidal. Several estimates of the braking index $`n`$ have been reported in the literature (see Boyd et al. 1995 for a compilation of older results) ranging from $`2.01\pm 0.02`$ (Manchester & Peterson 1989; Nagase et al. 1990) to $`2.74\pm 0.10`$ (Ögelman & Hasinger 1990). Deeter et al. (1998) analyzed an extended set of GINGA observations and found a braking index of 2.08$`\pm `$0.02. This result was substantially confirmed by Mineo et al. (1999), who combining a BeppoSAX frequency measurement with earlier ASCA results derived a value of 2.10$`\pm `$0.1. A recent analysis of all ASCA pointed observations gave a braking index of 2.10 (Hirayama et al. 2002). A glitch in the timing history of PSR B0540$``$69, the only one reported up to now for this young Crab-like pulsar, has been recently pointed out by Zhang et al. (2001). These authors based their analysis on a collection of RXTE observations, spanning 1.2 years. The approximate epoch of the glitch, according to their estimates, occurred on MJD 51325 $`\pm `$ 45 and the change in frequency and its first derivative were $`\mathrm{\Delta }\nu /\nu =(1.90\pm 0.04)\times 10^9`$ and $`\mathrm{\Delta }\dot{\nu }/\dot{\nu }=(8.5\pm 0.07)\times 10^5`$, respectively. Zhang et al. (2001) measured also the braking index after the glitch which resulted equal to 1.81 $`\pm `$ 0.07, significantly lower than the values reported from other analyses.

In this paper we present an exhaustive timing analysis of PSR B0540$``$69 performed on RXTE data set spanning aboUt 5 years. In particular, we extend the set of RXTE observations used by Zhang et al. (2001), adding more observations for a total time interval of about 3 years before and 2 years after the epoch of the glitch episode claimed by these authors. Moreover, to further increase the time interval data set, public ASCA and BeppoSAX observations are also included in the analysis extending the length up to about 8 years.

## 2 Observation and Data Reduction

The RXTE observations considered in the present paper were performed between March 8, 1996 and March 14, 2001. We used only data obtained with the PCA (Jahoda et al. 1996) accumulated in “Good Xenon” telemetry mode, time-tagged with a 1$`\mu `$s accuracy with respect to the spacecraft clock, which is maintained to the UTC better than 100 $`\mu `$s. The pulsar position inside the instrument Field of View (FoV) is different in the various pointings and ranges between 0 and 25 arcminutes far from the centre. Following Zhang et al. (2001) we considered only events with pulse-height channels between 5 and 50, corresponding to the 2-18 keV energy interval. Furthermore, we verified that the selection of all PCA detector layers, instead of those from the top layer only, increased significantly the S/N ratio of the pulsation and adopted this choice for all the observations.

PSR B0540-69 was observed by ASCA (Tanaka et al. 1994) 14 times between 1993-06-13 and 1999-11-03. Only data from GIS (Ohashi et al. 1996, Makishima et al, 1996) were used in our analysis. The pulsar position inside the FoV ranges between 0 and 22 arcminutes off-axis. The Narrow Field Instruments (NFIs) onboard BeppoSAX (Boella et al. 1997a) observed PSR B0540-69 five times between 1996-10-26 and 2000-02-19. We considered only data from MECS (Boella et al. 1996) because they are in the same energy range used for the analysis of the other two satellites. The position of the source inside the instrument FoV lies between 0 and 17 arcminutes. Data from both imaging instruments (GIS and MECS) were extracted from regions centered at the source position and the shapes and sizes of these regions were optimized for each observation, depending upon the different off-axis locations, to achieve the highest S/N ratio. The log of all these observations is given in Table A1, available in electronic form at CDS (http://cdsweb.u-strasbg.fr/).

## 3 Timing analysis and pulse profile

### 3.1 Frequencies and derivatives evaluation

The UTC arrival times of all selected events were first converted to the Solar System Barycentre using the (J2000) pulsar position $`\alpha =5^\mathrm{h}40^\mathrm{m}10\stackrel{s}{.}95`$ and $`\delta =69\mathrm{°}19\mathrm{}55\stackrel{}{.}1`$ (Caraveo et al. 1992) and the JPL2000 ephemeris (DE200; Standish 1982). For each observation we searched the pulsed frequency $`\nu `$ using the folding technique in a range centered at the values computed by means of the pulsar ephemeris reported by Deeter et al. (1999). The central time of each observation was chosen as reference epoch and the corresponding frequency was estimated by fitting the $`\chi ^2`$ peak with a gaussian profile. Both these data are listed in Table A1 (columns 2 and 7, respectively). Frequency errors at 1 $`\sigma `$ level were computed with two different methods. The first evaluation was performed computing the frequency interval corresponding to a unit decrement with respect to the maximum in the $`\chi ^2`$ curve ($`err=\nu (\chi _{max}^2)\nu (\chi _{max}^21)`$). In the second method we produced a template profile folding the longest RXTE observation (rows 33$`+`$34 in Table A1) with its best frequency (see Fig.1). The zero phase was taken at the centre of the peak. Then, we compared the folded profiles of other pointings with an accurate analytical model of the template changing the frequency within the searching range and computed the relative $`\chi ^2`$ values. For a correct comparison we subtracted the off-pulse levels and normalized the total pulsed counts. The frequency uncertainties were then set equal to the frequency interval corresponding to an increase of the $`\chi ^2`$ by a unit with respect to the minimum in the curve ($`err=\nu (\chi _{min}^2+1)\nu (\chi _{min}^2)`$). No significant difference was found between the error estimates obtained with the two methods, we then considered as proper statistical uncertainties of the pulsar frequencies those computed with the first one (column 8). All errors, thereafter, refer to one standard deviation. The times of arrival (TOA), referred to the peak centre, were also determined for all RXTE observations by cross-correlating the folded profiles with the analytical model of the template. The resulting phase offset was added to the central time of the observation. Errors on TOAs were assumed equal to the statistical uncertainties of the peak centre derived from a best fit of the pulse profile. TOA and errors are reported in Table A1 (columns 5 and 6.)

A first estimate of the frequency derivatives was obtained by the best fit of the frequency history listed in Table A1 with the second-degree polynomial:

$$\nu (t)=\nu _0+\dot{\nu }_0(tt_0)+\frac{1}{2}\ddot{\nu }_0(tt_0)^2.$$
(1)
This procedure was first applied separately to the two ASCA+BeppoSAX and RXTE frequency sets and the resulting parameters’ values were found in a good agreement within their uncertainties. We then evaluated the ephemeris for the whole set of data (RXTE+ASCA+BeppoSAX); a summary of these results is given in Table 1. For easy comparison, the same reference epoch ($`t_0`$) was taken, in all the fittings, equal to MJD 50372.5481748585. Fig. 2 shows the residuals obtained fitting the entire frequency set with Eq.(1): the number of frequency values in excess of 2 standard deviations is about 6% of the total, confirming the right evaluation of the frequency errors.

### 3.2 Pulse phase analysis

The precision of the timing parameters can be enhanced by maintaining a pulse coherence over the entire time interval covered by the observations. This analysis was performed only on the RXTE data. The BeppoSAX timing, in fact, does not maintain the indispensable accuracy to the UTC, fundamental to apply a coherent phase analysis. The ASCA observations are rather sparse and the systematics in the absolute time assignment (Hirayama et al. 1996) makes more complicated the proper phase recycle correction.

The arrival times of every event for RXTE observations were folded by using the values of $`\nu _o`$, $`\dot{\nu }_o`$, $`\ddot{\nu }_o`$ reported in Table 1 (line 3). Phase shifts of the pulse profiles, expected because of the ephemeris accuracy, were computed for all observations by a cross-correlation with the analytical model of the template (Fig. 1). Phase errors were taken equal to the TOA uncertainties multiplied by the corresponding frequencies. Throughout this paper phases are measured in cycle units. The resulting values of the phase shifts, corrected for the presence of recycles, are shown in Fig. 3 as a function of the observation epoch.

The large phase variations, of the order of several cycles in a few hundred days, indicate that input ephemerides of Table 1 are not the best ones and that it is possible to improve them. The corrections, in absence of frequency irregularities, can be obtained by fitting the phase lags by the simple third-degree spin-down model:

$$\mathrm{\Delta }\varphi (t)=\mathrm{\Delta }\nu _o(tt_0)+\frac{1}{2}\mathrm{\Delta }\dot{\nu }_o(tt_0)^2+\frac{1}{6}\mathrm{\Delta }\ddot{\nu }_o(tt_0)^3,$$
(2)
where $`\mathrm{\Delta }\varphi (t)`$ is the measured phase difference at time $`t`$, $`t_0`$ is the reference epoch and $`\mathrm{\Delta }\nu _o`$, $`\mathrm{\Delta }\dot{\nu }_o`$, $`\mathrm{\Delta }\ddot{\nu }_o`$ are the ephemeris corrections for the frequency, first and second derivatives, respectively. We found that this equation fails to describe the entire set of phase data showing very large systematic deviations (panel $`a`$ of Fig. 3). Such deviations indicate the presence of a timing noise characteristic of several other pulsars (Arzoumanian et al. 1994, Lyne 1999). We tried also to correct the phases by fitting higher order polynomials, up to the sixth degree. The best fit for this model is shown in the panel $`b`$ of Fig. 3: the residuals are of course largely reduced with respect to the other fit and the general behaviour is much better described, but smaller amplitude systematic deviations remain apparent.

To quantify the timing noise we can calculate the $`\mathrm{\Delta }_8`$ parameter defined by Arzoumanian et al. (1994) as:

$$\mathrm{\Delta }_8=\mathrm{log}(10^{24}|\ddot{\nu }|/6\nu ).$$
(3)
For PSR 0540$``$69 instead of the measured $`\ddot{\nu }`$, whose large value is related to the regular spin down rather than the timing noise, we used the correction that takes into account the residuals shown in Fig. 3a. A reliable estimate of this correction is given by the difference between the third degree coefficient of the best fit parameters of the sixth and third degree polynomials. This difference was found equal to $`2.27\times 10^{22}`$ Hz s-2 and gives $`\mathrm{\Delta }_8`$=0.28, higher than all the values of the sample of Arzoumanian et al. (1994). Note that PSR B0540-69 lies well in the upper part of the plot $`\mathrm{\Delta }_8`$ $`vs`$ $`\dot{P}`$ given by these authors.

For comparing our results with those of Zhang et al. (2001) we divided our data into two subsets, before and after the time of the glitch, assumed at MJD 51325 – the interval of possible epochs given by Zhang et al. (2001), spans 45 days – and applied a cubic polynomial fit to each subset. We obtained acceptable fits, characterized by nearly zero residuals; resulting best fit values for the two intervals are reported in Table 2, the ephemerides are given in Table 3. Propagating the parameters of Table 3 to the epoch of the glitch we obtained a marginal detection for the frequency change $`\mathrm{\Delta }\nu /\nu `$ = (1.8$`\pm `$1.0) $`\times `$ $`10^9`$, whereas more significant differences were found for the first and second derivative: $`\mathrm{\Delta }\dot{\nu }/\dot{\nu }`$ = (1.69$`\pm `$0.01) $`\times `$ $`10^4`$ and $`\mathrm{\Delta }\ddot{\nu }/\ddot{\nu }`$ = (2.043$`\pm `$0.001) $`\times `$ $`10^3`$. The significance of these parameters’ differences can be largely affected by the timing noise systematics. To verify this hypothesis we performed a similar analysis on other data subsets selected by changing the separation epochs and values of $`\mathrm{\Delta }\nu /\nu `$ of the same order and similar significance were found. This last finding addressed us to interpret the marginal detection of frequency jump at the glitch epoch claimed by Zhang et al. (2001) as a not genuine result but due to the pulse phase analysis in presence of strong timing noise.

### 3.3 The braking index

As discussed in the Introduction, a relevant difference between our previous results (Mineo et al. 1999) and those of Zhang et al. (2001) is in the value of the braking index: we essentially confirmed the estimates of Deeter et al. (1998), Hirayama et al. (2002), while Zhang et al. (2001) gave the value of 1.81, about 14% smaller, but significantly different when considering the associated uncertainties. Using the new ephemerides given in Table 3, we computed for the two intervals the values of $`n`$ and obtained 2.1272$`\pm `$0.0003 and 2.122$`\pm `$0.001 for the first and second interval, respectively.

A change in the timing parameters implies a variation of the braking index:

$$\frac{\mathrm{\Delta }n}{n}=\frac{\mathrm{\Delta }\nu }{\nu }2\frac{\mathrm{\Delta }\dot{\nu }}{\dot{\nu }}+\frac{\mathrm{\Delta }\ddot{\nu }}{\ddot{\nu }}.$$
$`(4)`$
From the above results it is clear that the parameter responsible for the much lower value of $`n`$ measured by Zhang et al. (2001) is the second frequency derivative, for which these authors reported the value of (3.23$`\pm `$0.12)$`\times 10^{21}`$ s-2, just 15% smaller than our result.

This discrepancy is likely due to the presence of such high timing noise, because the data set considered by Zhang et al. (2001) spans a rather narrow time interval with respect to that considered by us. To further investigate this point we computed the braking index in several independent time intervals, evaluating the best corrections to the pulsar ephemerides in each of them. The resulting values are shown in Fig. 4: they ranges from 1.97 to 2.47. Again the largest variation is due to the estimates of second derivative of the pulsar frequency, which differs of $``$20%, between the various intervals. We conclude, that reliable estimates of the braking index can be obtained only considering the longest time intervals in which a good fit of phases (or TOA) with Eq.(2) can be obtained, likely spanning several years. The use of shorter intervals can introduce a bias due to the timing noise.

## 4 Discussion

The only way to study the timing noise of PSR B0540-69 is the use of a dense set of X-ray observations, like that of RXTE, because this young pulsar is in the Large Magellanic Cloud and its flux is too weak to be monitored in the radio band. Our analysis on a large database of X-ray observations of PSR B0540-69, covering more than 5 years, provided a good evidence for a relevant timing noise affecting the phase of the pulsed signal. In particular, we showed that the best fit of a third degree polynomial, including up to the second frequency derivative, gives phase residuals up to a few cycles and that residuals as large as 0.4 remain even when a sixth degree polynomial is used. Assuming that the difference in the second derivative obtained from these best fits is a measure of its fluctuations, we evaluated the $`\mathrm{\Delta }_8`$ parameter (Arzoumanian et al.1994) which was found equal to 0.28, confirming the high level of timing noise. Taking into account that even larger variations of this derivative are also found when polynomial best fits are performed over shorter time intervals, this result could be considered a lower limit. Furthermore, it supports the finding that timing noise is stronger in young pulsars with a high $`\dot{P}`$.

A consequence of this high timing noise is that the glitch claimed by Zhang et al. (2001) cannot be confirmed. The frequency difference of this event given by these authors is very small, about 0.04 $`\mu `$Hz. The glitch cannot be detected directly as a sudden frequency jump because it is quite less than the uncertainties in the frequency measurements which are typically of the order of a few $`\mu `$Hz, and in the best conditions of a few tenth of $`\mu `$Hz. Zhang et al. (2001) also excluded that this effect can be due to timing noise, but their conclusion is affected by the use of a shorter time interval that does not allows an accurate analysis of the timing noise.

At variance, our results show that frequency differences of the same order of that given by Zhang et al. (2001) are usually found when different selections of time intervals are considered and they do not depend upon a well defined episode.

From the timing noise analysis we were also able to show how much the first and second derivative of the pulsar frequency are stable in time. We found that the former can have fluctuations of amplitude of about 10-4, while for the latter fluctuations can be much higher, and in some cases the estimates can differ of $``$ 10–20%, depending upon the length of the time interval taken into account. Such large variations affect the evaluation of the braking index, particularly when it is found by the fitting of Eq.(2) to the pulse phases. We showed that when the longest possible intervals are considered, $`n`$ turns out to be very close to 2.12, in agreement with several other previous estimates, while values as that given by Zhang et al. (2001) are obtained in shorter intervals. We verified this interpretation by fitting Eq.(2) to the same subset of RXTE observations used by these authors (more specifically, the fit was performed on RXTE observations from row 82 to row 114 of Table A1) and derived from the best fit ephemeris $`n`$ = 1.854 $`\pm `$ 0.003. However, when these parameters are used to extrapolate the phase shift to the entire subset data after the epoch of the glitch claimed by Zhang et al. (2001) they produce a systematic deviation of the residuals which increases with the elapsed time, as shown in plot of Fig. 5.

Finally, we note that the agreement of our estimate of braking index with those derived from ASCA (Hirayama et al. 2002) and BeppoSAX (Mineo et al. 1999) data can be easily understood on the basis of the rather long time intervals covered by the observations. They are too sparse to provide a good information on the timing noise but can give a right evaluation of the mean second derivative.

###### Acknowledgements.

The authors are very grateful to the referee, R.N.Manchester, for the relevant comments and suggestions that greatly improved the scientific content of the paper. This work has been partially supported by the Italian Space Agency (ASI).

## Appendix

In this appendix we report a table containing the central epoch of the observation (column 2), the total duration (column 3), exposure time (column 4), TOA and error (column 5 and 6; only for RXTE observations) and the frequency with the relative error (column 6,7) for each observation analyzed.


text:
# Focusing of electromagnetic waves by periodic arrays of dielectric cylinders

## Abstract

By numerical simulations, we show that properly arranged two dimensional periodic arrays, formed by dielectric cylinders embedded in parallel in a uniform medium, can indeed act as an optical lens to focus electromagnetic waves, in accordance with the recent conjecture in the literature. The numerical simulations are based on an exact multiple scattering technique. The results suggest that the E-polarized waves are easier to be focused than the H-polarized waves. The robustness of the focusing against disorders is also studied. Comparison with the corresponding cases for acoustic waves is also discussed.

Photonic crystals (PCs) review ; rev2  are made of periodically modulated dielectric materials, and most sonic crystals sigal  (SC) are made up of materials with periodic variation of material compositions. Photonic and sonic crystals have been studied both intensively and extensively. When passing through photonic crystals, the propagation of electromagnetic (EM) waves can be significantly affected by the photonic crystals in the same way as that electrons are controlled by the crystals. This interesting phenomenon has stimulated a variety of possible applications of PCs in controlling photons. In particular, considerable efforts have been devoted to finding photonic crystals that can completely block propagation of electromagnetic waves in all directions within a certain range of frequencies, termed as photonic bandgap. It has been suggested that PCs may be useful for various applications such as antennae brown , optical filters chen , lasers evans , prisms lin1 , high-Q resonant cavities lin2 , wave-guides kraus1 , mirrors kraus2 , left-handed materialsNotomi ; LHM1 ; LHM2 , and second harmonic generations mart . These applications mostly rely on the existence of photonic bandgaps, and a majority of them is not concerned with the linear dispersion region well below the first gap. In other words, most of earlier studies were focused on the formation of bandgaps and the inhibited propagation of waves.

Recently, the interest in the low frequency region, where the dispersion relation is linear, has just started. Since the wavelength in this region is very large compared to the lattice constant, the wave sees the medium as if it were homogeneous, in analogy with wave propagation in normal media. Consequently, a possible new application of PCs has been suggested by a number of authorsHalevi1 ; Halevi2 . These authors suggested that PCs could also be employed as custom-made optical components in the linear regime below the first bandgapHalevi2 . However, no physical realization of optical lenses has been made so far. Along the same line of thought, it was suggested that SCs may also be used to build acoustic lenses to converge the acoustic waves. A necessary condition to be satisfied for constructing an acoustic lens is that the acoustic impedance contrast between the SC and the air should not be large; otherwise acoustic waves will be mostly reflected. The recent experimentcervera  and the corresponding numerical simulationBikash  on acoustic waves propagation through a lenticularly shaped SC have confirmed that acoustic lenses by SCs are indeed possible. Encouraged by these findings, in this paper we would like to further explore the possibility of using PCs as an optical lens to focus electromagnetic waves, following the line of the simulation of acoustic lensesBikash .

In this paper we carry out numerical simulations on the focusing of EM waves by PCs. We wish to theoretically realize the particular predictions made in Halevi1 ; Halevi2 . Since the multiple scattering technique has been successfully applied earlier ye  to reproduce some experimental results on acoustic propagation and scattering in SCs and this technique can be fully adopted to EM waves, we will use this technique to study the focusing effect of EM waves by PCs in detail. To the best of our knowledge, there has been no earlier attempt in using the multiple scattering theory to investigate the focusing phenomenon of PCs.

The system considered here is similar to what has been presented in Halevi1 ; Halevi2 . Assume that $`N`$ uniform dielectric cylinders of radius $`a`$ are placed in parallel in a uniform medium, perpendicular to the $`xy`$ plane. The arrangement can be either random or regular. The scattering and propagation of EM waves can be solved by using the exact formulation of TwerskyTwersky . While the details can be found in ye , here we only brief the main procedures. A unit pulsating line source transmitting monochromatic waves is placed at a certain position. The scattered wave from each cylinder is a response to the total incident wave, which is composed of the direct contribution from the source and the multiply scattered waves from each of the other cylinders. The response function of a single cylinder is readily obtained in the form of the partial waves by invoking the usual boundary conditions across the cylinder surface. The total wave ($`E`$ or $`H`$ for the E- or H- polarization respectively) at any space point is the sum of the direct wave ($`E_0`$ or $`H_0`$) from the transmitting source and the scattered wave from all the cylinders. The normalized field is defined as $`TE/E_0\text{or}H/H_0`$; thus the trivial geometrical spreading effect is eliminated.

When the cylinders are placed regularly, the phenomenon of band structures prevails, and can be evaluated by the standard plane-wave expansion method. Fig. 1 shows the band structures for both E- and H- polarized EM waves when propagating through an array of square lattice of dielectric cylinders with radius $`a=0.38`$ cm, placed in the air. The dielectric constant for the cylinders is 10, which is smaller than that in Halevi1 . The fractional area occupied by the cylinders for a unit area , i. e. the filling factor, is 0.13. In the simulation, the frequency is made non-dimensional by scaling as $`ka`$. Here is shown that a complete bandgap appears for the E-polarized wave. Following Halevi1 ; Halevi2 , we also calculate the phase speed from the band structure for the first band as a function of the filling factor, and the results are shown in Fig. 1 (c). It is the linear region of the first band that we will consider in the following. The results indicate that the phase speed for the E-wave is more significantly reduced by the PC.

A PC made optical lens is illustrated in Fig. 2, which is in line with the prediction in Halevi1 . The source is placed at a distance far enough from the lens so that the incident waves can be regarded nearly as plane waves. In this way, the focusing point of the lens can be inferred from the transmitted field on the other side of the array.

In the rest computation, the following parameters and arrangements are used: (1) Square lattices of cylinders and the propagation along  direction are considered; (2) the dielectric constant for the cylinders is 10; (3) the filling fraction is 0.13; (4) the radius of cylinders is 0.38 cm; (5) the frequency of the waves is taken as $`ka=0.152`$, well within the linear region; (6) the lens thickness, i. e. the max range along the $`x`$-axis, is 10, and the height, i. e. the max span along the $`y`$-axis, is 20; (7) the distance between the source and the center of the lens is 100. All the lengths are scaled non-dimensionally by the lattice constant.

Fig. 3 shows the two-dimensional spatial distributions of the transmitted intensities. Here we see that although there are some focusing effects for both E- and H- polarized waves, the focusing effect is mostly prominent for the E-polarized wave in particular. This agrees with the expectation from the phase speed estimate in Fig. 1 (c): the phase speed for the E- polarized wave is more reduced, rendering a bigger contrast to the outside medium, thereby yielding a bigger refractive index. The focusing point of the E-polarized waves is at about $`x=9`$, and the waves are better focused along the $`y`$-axis. From the results in Fig. 3, we see that both E- and H- polarized waves have certain spreading along the x-axis, i. e. the propagation direction. This is very similar to what has been observed in the acoustic lensescervera ; Bikash . In any event, the fact that the focusing features are in certain qualitative agreement with the earlier predictionHalevi1 ; Halevi2  is encouraging. Another note should be made here. In Halevi2 , Halevi et al. conjectured an elliptically shaped lens using a 2D PC (See Fig. 1 in the paper). According to our earlier simulation on the acoustic lenses, a lens of such a shape is less efficient than the diamond-like lens illustrated by Fig. 2.

For comparison, we have also considered the E-polarized EM wave transmission in the  direction through a slab of rectangular array of dielectric cylinders. The size of the slab is $`10\times 20`$: the length along the $`x`$-axis is 10 and the height along the $`y`$-axis is 20. The transmission results are shown in Fig. 4. Interestingly, there are some focusing effects. For example, there is a focused field centered around $`x=26,y=0`$ and there are other two focusing centers at $`x=11`$ and $`y=5`$ and $`5`$ respectively. This is quite different from the corresponding case with the acoustic arraysBikash . This seems not in the expectation. The reason follows. As the source is quite far away from the slab, the incident wave on the slab could be thought as a plane wave. Then it is expected that the transmitted wave should not be focused in the space when there is no such focusing effect as shown in LHM2 . The results in Fig. 4 together with that in Fig. 3 imply that although the focusing is quite a general feature of a lattice arrangement of dielectric cylinders, the shape of the lens as shown in Fig. 2 may be essential for a unique focusing.

We have also examined the robustness of the focusing against disorders. Here we consider the positional disorder of the cylinders. Since the wavelength is larger than the lattice constant in the present cases, one might intuitively conclude that the positional disorder has no effects. The results for the E-polarized wave are shown in Fig. 5. Here the shape of the array is the same as that in Fig. 2, except that the cylinders are placed in a complete randomness. Comparing to Fig. 3, it is shown that such a disorder does not destroy the focusing, in accordance with the intuition. That the disorder has little effect on the focusing of EM waves differs from the situation with the acoustic systems. In Bikash , it was shown that the disorders can completely destroy the focusing phenomenon in the acoustic lenses. All these results tend to support the homogenization which has been carried out in Halevi1 .

In summary, here we report the results of EM waves transmission through lenticular structures made of dielectric cylinders. Complying with the previous conjecture, the EM wave focusing effects are indeed observable by such structures.

###### Acknowledgements.

This work received support from National Science Council of Republic of China and the Department of Physics at the National Central University.


text:
# Lattice constant in diluted magnetic semiconductors (Ga,Mn)As

## Acknowledgment

This work has been done within the project AVOZ1-010-914 of the ASCR. The financial support provided by the Grant Agency of the ASCR (Grant No. A1010214), by the Grant Agency of the Czech Republic (202/00/0122), and by RTN project No. HPRN-CT-2000-00143 the EC acknowledged.