text:
# Stress transmission through three-dimensional granular crystals with stacking faults

## I Introduction

Static packings of granular particles support stress in a nontrivial manner . The forces on any given particle from all neighboring particles and from gravity must balance vectorally. Stress on a random granular pack is supported by a disordered network of these interparticle contacts in local force balance with each other . It has been suggested that elastic theory is capable of accounting for these inhomogenaities . However, the indeterminacy when the particles are infinitely hard, so that particle distortions are not allowed, has led to a search for alternative descriptions of the mechanical properties of granular materials in the hard-sphere limit .

The proposed alternatives have had varying degrees of success. Coppersmith et al. proposed a diffusive q-model for the transmission of force through granular packings , and a hyperbolic model was created by Cates et al. . Experimentally it has been difficult to rule out any of these descriptions conclusively, but recently progress has been made in studying the conditions under which these models apply.

Much understanding and additional surprises have been gained from experiments studying the response to a localized external force, similar to a Green’s function. DaSilva and Rajchenbach  studied a two dimensional packing of rectangular photoelastic bricks. When applying an external force at the top of the packing they measured the response forces between neighboring bricks within the pack. The results showed that at a given depth the maximum interparticle force was located directly below the position of the applied force. Deeper in the packing (farther from the applied force) they found the width of this peak of maximum force to broaden as the square root of the depth. Similar results were found by Moukarzel et al. in recent displacement response experiments . Such findings support diffusive models. However, a cross-over to a linear widening of the response peak cannot be ruled out for deeper packings. Reydellet and Clément looked at the response forces at the bottom surface of an amorphous three-dimensional packing of spheres . Again they found the maximum force situated directly beneath the applied force, but now the width of the peak grew linearly with depth from the surface. This linear growth lends credence to the elastic theories. Similar elastic-type behavior was found in two-dimensional packings by Geng et al.. In addition those researchers also found a strong dependence on the spatial ordering of the particles .

The importance of spatial ordering was further shown in three dimensional granular crystals in Ref. . Within face-centered-cubic (FCC) crystals, the force was found to be supported along straight lines of contacts between beads, resulting in three areas of large force at the bottom surface in response to a locally applied external load at the top surface. In contrast, an hexagonal-close-packed (HCP) arrangement resulted in a ring of large force at the bottom surface. In both cases there was a local minimum of force directly beneath the point of application. Such patterns are reminiscent of the hyperbolic models, which predict straight-line propagation of forces. Bouchaud et al. have extended these hyperbolic models to include the splitting of force chains at defects within a granular packing . Such splittings are predicted to lead to the elastic-like behavior that has been observed in amorphous packings.

It is, therefore, of interest to explore the regime between perfectly-ordered crystalline bead packs and completely amorphous arrangements. As a first step in this direction, we present an experimental study of the effect on the force propagation of a small number of stacking fault defects within otherwise perfectly ordered FCC crystals.

## II Experimental Methods

Employing the triangular acrylic cell described by Blair et al.  and Mueggenburg et al.  and a ’by hand’ construction technique, we produced large FCC crystals of approximately $`20,000`$ soda lime glass spheres of diameter $`3.06\pm 0.04`$mm with an extremely low number of defects (estimated to be fewer than ten beads significantly out of place in the entire pack) oriented as horizontal planes of triangular order. Furthermore, we carefully controlled the stacking of the planes in order to create true three-dimensional crystals and to introduce stacking faults at specific locations.

In this manner, we constructed five different granular FCC crystals, each having one or two defects in the stacking order. We then applied an impulsive force to a small region, approximately two beads in diameter, centered on the top of the crystal . Using a carbon paper technique, we characterized the normal forces on the top and bottom surfaces of the crystal in response to this force . A piece of carbon paper was sandwiched between the crystal and a piece of white paper on both the top and bottom surfaces of the crystal. The applied force caused individual beads to press into the carbon paper and leave marks on the white paper. These marks were then digitized and image analysis software was used to calculate the area and intensity of each mark .

We conducted twenty experimental runs on each crystal. For each run the resultant carbon marks were aligned with the expected triangular lattice. The areas and intensities were averaged over experimental runs and, in order to improve statistics, were averaged over symmetries of the crystal (one reflection and two rotations) when such averaging did not qualitatively alter the intensitiy patterns.

## III Results

A pure FCC crystal consists of layers of triangular order stacked in such a way that every third layer lies on top of the first. We annotate this stacking as $`\underset{¯}{ABCABCABCABCABC\mathrm{}}`$ A stacking fault corresponds to a disruption of this order and may be considered to be a boundary between two pure FCC crystals with different stacking orders.

We constructed a $`19`$ layer crystal with a stacking fault centered at the tenth layer. This stacking is represented as $`\underset{¯}{ABCABCABC}\overline{\underset{¯}{A}}\overline{CBACBACBA}`$ The underline represents one pure FCC crystal, while the overline corresponds to a second pure FCC crystal with a different orientation. At the boundary between the two crystals there exists an $`\overline{\underset{¯}{A}}`$ layer which correctly matches both stacking orders and can thus be considered to be in both crystals. Note that, by necessity, the three layers centered around the stacking fault $`\underset{¯}{C}\overline{\underset{¯}{A}}\overline{C}`$ follow an hexgonal-close-packed stacking order. When applying a quick impulse to a small area at the center of the top of this crystal, the bottom surface displayed six regions of large force at the vertices of a regular hexagon as shown in figure 1. The regions of large force are approximately the same size as the area over which the impulse was applied on the top surface. We observed another region of larger than background force inside the hexagon, but closer to one face. At other symmetric places of the pattern, areas of larger than background force are present, but with much lower intensity. It is not clear if this area of heightened force is a real feature of the crystal with stacking fault or if this is merely an extreme fluctuation. Averaging over crystal symmetries would accentuate such fluctuations, and thus we also show the intensity pattern without this averaging.

Additionally, we explored possible assemblies of crystals with two stacking faults. Figure 2 represents a $`13`$ layer FCC crystal with stacking faults at the seventh and eleventh layers. The stacking order (from top to bottom) was $`\underset{¯}{ABCABC}\overline{\underset{¯}{A}}\overline{CBA}\overline{\underset{¯}{C}}\underset{¯}{AB}`$ as shown schematically in figure 2(a). We found that this crystal produced a twelve-peaked intensity plot (Figure 2b). The plot shows two pairs of peaks of nearly equal height aligned along the faces of an equilateral triangle. Each pair of peaks is approximately four beads wide and the separation between pairs is approximately three beads.

Another crystal with two stacking faults consisted of eleven layers with faults centered at the fourth and eighth layers ($`\underset{¯}{ABC}\overline{\underset{¯}{A}}\overline{CBA}\overline{\underset{¯}{C}}\underset{¯}{ABC}`$). For this crystal, we found three tall peaks flanked by peaks half the height of the taller ones arranged along the faces of an equilateral triangle (Figure 3). The width of each of the taller peaks was slightly larger than that of the smaller flanking peaks which correspond to the width of the area over which the force was applied at the top surface.

The next crystal we built contained thirteen layers with faults at the fifth and ninth layers. The stacking order was $`\underset{¯}{ABCA}\overline{\underset{¯}{B}}\overline{ACB}\overline{\underset{¯}{A}}\underset{¯}{BCAB}`$. We found an intensity pattern (Figure 4) similar to that of the previous crystal with two faults. This crystal produced three tall peaks flanked by peaks about one third the height of the taller peaks arranged along the faces of an equilateral triangle. The taller peaks were approximately three beads wide and the shorter peaks were approximately two beads wide. In contrast to figure 3, which corresponds to a crystal with fewer layers above and below the stacking faults, the flanking peaks are now farther away from the taller peaks.

Another possible type of stacking fault involves the omission of one plane from a large FCC crystal. This fault maintains the orientation of the crystal above and below the fault and thus might be considered a more basic defect than the grain boundary described above. However, this fault could be described as a special case of two grain-boundary faults. Figure 5 shows the results for a crystal with layer stacking $`\underset{¯}{ABCABC}\overline{\underset{¯}{A}}\overline{\underset{¯}{C}}\underset{¯}{ABCABC}`$. Although one might suppose that this crystal exhibited the least amount of stacking disorder, we found that the force intensity pattern consisted of six approximately equal-height peaks, one at each of the vertices of an equilateral triangle and one at each of the midpoints of the sides (Figure 5). The peaks at the vertices of the triangular pattern are four beads in width while those on the sides are three beads wide. We describe this pattern below as a special case of the nine or twelve peaks shown in the crystal packings with two faults.

## IV Discussion

As in the experiments performed on pure FCC crystals by Mueggenburg et al. , we can reconstruct the force chain paths through the crystal based on the force intensity plots from the experiments. More importantly, the simple force balance model used to explain the force patterns in that study, also explains the present results.

Forces will travel through pure FCC sections of the crystal, along the spines of a tetrahedron. Thus, to understand an elementary stacking fault, all we need to determine is what happens at the interface between two FCC crystals with different orientations. Again, the force balance model can be applied. In order to maintain force balance at the interface, the force chain must split. Note that locally this is the same geometry as an HCP crystal. Below the stacking fault each of the new forces will travel through the lower portion of the crystal in straight line paths. Figure 6(a) illustrates this force splitting at a grain boundary.

We found a regular hexagonal pattern of force in our first crystal because the fault was placed in the middle of the packing. In general, we would expect a six-peaked pattern in which two peaks reside along each side of a triangle. As the defect is moved higher or lower in the crystal, these pairs of peaks would move towards or away from each other.

It is simple to extend this idea to crystals containing two stacking faults. The crystal can be thought of as three pure FCC crystals with top and bottom sections having the same stacking order and the middle section having a different stacking order. At each interface the force splitting process is repeated resulting in twelve peaks in the intensity pattern (3 split into 6 and then split into 12). Figure 6(b) illustrates this multiple splitting process. However, most of our results do not demonstrate the expected twelve-peaked pattern. This can be explained by force recombination.

For certain combinations of stacking faults, the forces will converge upon the same lattice point on the bottom of the crystal. In other cases, two forces could travel along parallel lines through the crystal and arrive at neighboring lattice points on the bottom surface. In the former scenario, the force intensity at the convergence point should be twice the intensity of a lattice point with only one line of force. In the latter scenario, the forces should produce a wider than normal peak (the exact width would depend upon the width of the applied impulse and the separation of the force chains within the packing). These two effects account for our observed force patterns from crystals with two faults.

In the second and third crystals with two stacking faults (Figures 3 and 4), the taller central peaks are places where forces have recombined and the smaller flanking peaks correspond to single force chains. The widths of the peaks corroborate this explanation in both cases. The fourth example, where one plane of the lattice has been removed, can be categorized as a crystal with two stacking faults in which the twelve peaks have been reduced to six due to this recombination. The corner peaks in the intensity plot (Figure 5) are caused by parallel forces terminating at adjacent lattice points and the side peaks are caused by two forces terminating at the same lattice point as shown in figure 6(c). The widths of the peaks support this explanation. This is the extreme case of the trend from figure 3 to figure 4 where the flanking peaks moved towards the corners of the triangular pattern as the spacing between defects relative to the overall crystal height decreased. In the case, where the two defects occur at successive layers, the flanking peaks are positioned very near the corners and thus overlap with each other.

## V conclusions

In this study, we have established how stacking faults affect the transmission of stress through an FCC crystalline granular assembly. We have observed that the force chains split into two equal components at the boundary between different FCC crystal orientations. However, we note that the forces are transmitted as in a pure FCC crystal before and after splitting. Our findings reinforce those in reference  and the force balance model proposed therein.

We have identified the elementary stacking fault defect as the boundary between two pure FCC crystals of different stacking orders. All further stacking faults can be described as repetitions of this basic defect. Alternately, we can express the stacking faults as sections of HCP crystals within the overall FCC crystal. Although we find the grain-boundary description the most lucid, this other interpretation is geometrically equivalent.

As the number of stacking faults increases, the intensity pattern looks increasingly like that of an HCP crystal, and indeed the local stacking order around a fault is HCP. This leads to the conclusion that the force pattern of an HCP crystal can be viewed as the extreme limit of an FCC crystal with a stacking fault at every layer.

## VI Acknowledgments

This work was supported by NSF-CTS 0090490 and by the NSF MRSEC Program under DMR-0213745. MJS acknowledges support by the University of Chicago MRSEC Summer 2002 REU program.


format: text

text:
# Iron Line Diagnostics for the GRS 1915+105 Black Hole

## 1 INTRODUCTION

Martocchia et al. (2002, Paper I ) and Martocchia et al. (2003, Paper II ) reported the discovery of intense iron K$`\alpha `$ fluorescent emission in two BeppoSAX observations of the microquasar GRS 1915+105, which took place respectively on April 19-20, 1998, and April 21-22, 2000. During these observations the Fe line is broad and asymmetric, best fitted with a relativistic disc model (see e.g. , ,  and Paper I for references on the relativistic line model).

These are the only two cases, to our knowledge, in which the iron K$`\alpha `$ line is so prominent in this source. The fact that the line is not always such strong in GRS 1915+105 observations should not surprise, since finding the line out of the thermal emission high-energy tail can be awkward in galactic Black Hole (BH) candidates, because of the high disk temperature. Nevertheless, strong, broad, possibly relativistic iron K$`\alpha `$ lines have recently been detected in various spectra of galactic BH candidates, e.g. Cyg X-1, XTE J1650-500, XTE J1550-564, V4641 Sgr, GRO J1655-40, XTE J2012+381, GX 339-4, XTE J1748-288, XTE J1908+094. They may help estimating the rotation parameter of the BHs hosted at the center of these sources by means of the iron K$`\alpha `$ diagnostics.

## 2 THE SOURCE

GRS 1915+105 is a superluminal jet source, well-known since its discovery by the WATCH experiment on Granat . The source distance is about 10 kpc, and its inclination – based on observations of its jets – 70 degrees. The optical counterpart was found only recently, via infrared observations, yielding evidence that GRS 1915+105 belongs to the class of low-mass X-ray binaries; the same observations allowed to determine the mass of the central compact object, which has been constrained to $`M_\mathrm{c}=14\pm 4M_{}`$, i.e. well above the standard neutron star mass limit (see  and references therein). GRS 1915+105 is thus believed to host a BH with a gravitational radius $`r_g=\frac{GM}{c^2}21`$ km.

However, the issue of the central BH spin remains open. In the case of the two BeppoSAX observations here discussed, we could profit of the Fe line diagnostics, which allows to determine the accretion disk innermost stable orbit $`r_{\mathrm{ms}}`$, a known function of $`J/M`$ in Kerr metric, customarily put equal to the inner boundary of the optically thick portion of the disk itself.

GRS 1915+105 is a notoriously variable source, with the inner accretion disk being blown-up and becoming optically thin for a relevant fraction of the time. This is one more reason for the line not to be always visible in BeppoSAX observations. We had to select appropriate time intervals, out of both observations, in which the variability is less dramatic.

## 3 DATA ANALYSIS

Spectral fits have been performed with the XSPEC software package. We assumed an optically thick, neutral, and geometrically thin Keplerian accretion disk with solar elemental abundances. For both data sets, results of a preliminary fit showed that the dependence of the disc emissivity on the radius is rather flat: we assumed a phenomenological dependence $`r^2`$. “Cold” absorption was also included: its column density parameter $`n_\mathrm{H}`$ always stabilizes at about $`5.4`$ and $`5.45\times 10^{22}`$ cm-1, in the two observations respectively; we thus choose these as fiducial values.

### 3.1 First BeppoSAX observation: April 19-20, 1998 (Paper I)

This observation started at 11:33:19 UT and ended on the next day at 20:14:52 UT. Out of the entire $`\mathrm{\Delta }t120`$ ks, six time intervals of $`2000`$ s each were sorted out where the source variability is less pronounced. These are six orbits in which the source is in a ”quiescent” state. Data from the LECS, MECS and PDS instruments have been used in the intervals 0.1–4, 1.7–10 and 15–150 keV, respectively. The basic ingredient to model the continuum is a powerlaw with cutoff. As a model for the locally emitted fluorescent and reflection features we initially considered numerically computed spectra, including the Fe line and the underlying reflected continuum derived from Monte Carlo computations  and processed through our additive routine kerrspec. However, we verified that the model is quite insensitive to the Compton-reflected emission: to save computational time we adopted for it a simple, non-relativistic model (pexrav in XSPEC), and used kerrspec only to fit the line profile.

We applied these ingredients on MECS and PDS data only. The iron line parameters do not significantly change with more refined assumptions for the continuum, for instance accounting for LECS data and including a multicolor blackbody component (diskpn). Indeed, the disc thermal emission is not clearly recognizable. Furthermore, the bestfit disk temperature comes out unrealistically high ($`kT3.3`$ keV), while the BH mass results far too small from a physical point of view. These facts indicate that the standard interpretation of the continuum is not fully appropriate here. Here are the bestfit parameters: $`F_{210\mathrm{keV}}1.8\times 10^8`$ erg cm-2 s-1; $`\mathrm{\Gamma }2.6÷2.7`$; $`E_{\mathrm{cut}}25`$ keV. The reflection parameter $`R`$ is not well constrained and can vary from $`0.15`$ up to $`1.2\times 2\pi `$. The innermost line-emitting radius $`r_{\mathrm{in}}`$ stays always above $`6r_\mathrm{g}`$ (from $`7.8r_\mathrm{g}`$ in the fifth interval up to $`79r_\mathrm{g}`$ in the second one); the Fe line equivalent width (EW) is much larger than 100 eV in all cases, up to 300 eV in one interval.

Even though the $`\chi ^2`$ statistics is not very good – mainly due to the difficulty in modeling the continuum components – the distorted profile (cp. Fig. 1) is always much better fitted with a relativistic line than with a gaussian. We therefore infer that the relativistic solution for the iron line, besides being the most natural physical explanation for the very broad line observed, does not depend from the choice of the continuum model.

### 3.2 Second BeppoSAX observation: April 21-22, 2000 (Paper II)

In this case we analyzed separately two selected – “lower” and “higher flux” – spectral phases, lasting $`10`$ ks each, out of the entire (19:56:55–11:16:55 UT) observation. Only data from the MECS and PDS instruments were available, and were used in the same energy intervals as above. The flux is $`F_{210\mathrm{keV}}0.7`$ and $`1.1\times 10^8`$ erg cm-2 s-1, in the two states respectively. Again, we started our data analysis using a powerlaw with cutoff, disk multi-blackbody emission and Compton reflection “hump”. Contrary to the 1998 observation, thermal emission from the disk is not negligible this time. Broad and skewed residuals are apparent around 6.4 keV. The Fe feature is unphysically broad when modelled with a bare Gaussian ($`\sigma 2`$ keV, $`\chi _{\mathrm{red}}^2>2`$), suggesting relativistic broadening. We thus tried with a fully relativistic model wabs * ( diskline + diskpn + refsch ) where refsch accounts for relativistic distortions in the reflected continuum. Since on this way the best-fit innermost line-emitting orbit always results to be equal to the inner hard limit – $`6r_\mathrm{g}`$, the last stable orbit in Schwarzschild metric – it has been straightforward to test the spinning BH assumption with these data. We employed the Kerr line model kerrspec (cp. Paper I) to replace diskline. Best-fit parameters are shown in the Table. The iron line rest energy is $`E_06.4`$ keV. The value of the farthest emitting radius $`r_{\mathrm{out}}`$ is not well constrained by the iron line models, and always tends to the outer hard limit set by the model itself.

## 4 DISCUSSION

The Fe K$`\alpha `$ fluorescent emission feature detected in both 1998 and 2000 BeppoSAX observations of the microquasar GRS 1915+105 is strong, broad and asymmetric, best-fitted with a relativistic disc model. The line is emitted from neutral or low ionized iron. In the case of the 1998 observation we found evidence of emission from a region of the disk which is near, but still outside the Schwarschild innermost stable orbit; on the other hand, in the 2000 spectrum emission from $`r<6r_\mathrm{g}`$ (the last stable orbit in Schwarzschild metric) is compatible with the data, even if the results of our data analysis do not allow a firm conclusion on this regard. Emission from inside $`6r_\mathrm{g}`$ would indicate a Kerr spacetime, i.e. that a rotating BH is hosted at the center of the system. Of course, the disc parameters are somewhat dependent on the adopted model for the continuum, and must be taken with caution. Since the BH spin stays constant between the two observations, differences in the innermost emitting radii must reflect changes in the accretion flow.

Most of the observed broad, intense Fe lines in galactic BH candidates have beed detected during “very high”, “intermediate” or “quiescent” spectral states. A plausible explanation for this has still to be found: to produce the line, the accretion disc must efficiently reflect primary photons, even at the innermost radii. This is unlikely in “hard” states, since such states correspond to absent or optically-thin disks, but could be the case of “soft/high” states.

In the GRS 1915+105 1998 observation we registered very strong flux and an extremely steep powerlaw, with a cutoff at rather low energies. In the 2000 observation we have a similar powerlaw continuum, but thermal emission is clearly seen, and “cold” reflection is better constrained. Both observations may correspond to “very high” states, with Comptonization tails originating from the disk blackbody seed photons. BH mass estimates inferred using the thermal luminosity are physically unplausible, also due to the fact that current thermal emission models are not adequate to describe radiation transfer very near to the event horizon.

More detailed discussions on the data analysis and results are reported in Paper I and Paper II. Various issues still remain open, with regard to the accretion disk geometrical and physical states. Data with enhanced spectral resolution and sensitivity are strongly needed, and will be hopefully achieved by the next observational campaigns, first of all the ones scheduled for the satellite XMM-Newton.


text:
# 1 Introduction

## 1 Introduction

This communication is inspired by the ongoing spate of papers published on the inappropriately designated left-handed materials which are macroscopically homogeneous and display negative phase velocities in relation to the time–averaged Poynting vector, but are not chiral . Nominally, such a material is deemed to possess a relative permittivity scalar $`ϵ_r=ϵ_r^{}+iϵ_r^{\prime \prime }`$ and a relative permeability scalar $`\mu _r=\mu _r^{}+i\mu _r^{\prime \prime }`$, both dependent on the angular frequency $`\omega `$, such that both $`ϵ_r^{}<0`$ and $`\mu _r^{}<0`$ in some spectral regime.222The condition for the phase velocity and the time–averaged Poynting vector to be oppositely directed is $`\left(|ϵ_r|ϵ_r^{}\right)\left(|\mu _r|\mu _r^{}\right)>ϵ_r^{\prime \prime }\mu _r^{\prime \prime }`$, which permits — more generally — $`ϵ_r^{}`$ and/or $`\mu _r^{}`$ to be negative . An $`\mathrm{exp}(i\omega t)`$ time–dependence having been assumed here, $`ϵ_r^{\prime \prime }>0`$ and $`\mu _r^{\prime \prime }>0`$ at all $`\omega >0`$ for all passive materials. Originally conceived by Veselago , these materials have been artificially realized quite recently . Their fascinating electromagnetic properties can have technological implications of massive proportions , but those implications remain speculative at this time. Although these materials have been variously named , we prefer the name negative phase–velocity (NPV) materials as the least ambiguous of all extant names.

Using transmission–line analysis and lumped parameters, Alu and Engheta  have recently suggested a new route to realizing NPV materials: Take a thin layer of an epsilon–negative (EN) material: it has a negative real permittivity scalar but a positive real permeability scalar. Stick it to a thin layer of a mu–negative (MN) material, which has a negative real permeability scalar and a positive real permittivity scalar. Provided the two layers are sufficiently thin, the paired EN–MN layers could function effectively as a NPV material. The clear attraction of this scheme is that EN and MN layers are easier to manufacture, very likely, than the NPV materials fabricated thus far .

Our objective here is to examine the suggested scheme using continuum field theory, and to establish a restricted equivalence of an EN–MN bilayer to a NPV material. Implications for parallel–plate waveguides  and perfect lenses  are deduced therefrom.

A note about notation: Vectors are in boldface, column vectors are boldface and enclosed in square brackets, while matrixes are denoted by Gothic letters; $`ϵ_0`$ and $`\mu _0`$ are the free–space permittivity and permeability, respectively; and $`k_0=\omega (ϵ_0\mu _0)^{1/2}`$ is the free–space wavenumber. A cartesian coordinate system is used, with $`𝐮_x`$, $`𝐮_y`$ and $`𝐮_z`$ as the unit vectors.

## 2 Bilayer Theory in Brief

Consider the layers $`0<z<d_a`$ and $`d_a<z<d_a+d_b`$. Their constitutive relations are as follows:

$$\begin{array}{c}𝐃(𝐫)=ϵ_0ϵ_a𝐄(𝐫)\hfill \\ 𝐁(𝐫)=\mu _0\mu _a𝐇(𝐫)\hfill \end{array}\},0<z<d_a,$$
(1)
$$\begin{array}{c}𝐃(𝐫)=ϵ_0ϵ_b𝐄(𝐫)\hfill \\ 𝐁(𝐫)=\mu _0\mu _b𝐇(𝐫)\hfill \end{array}\},d_a<z<d_a+d_b.$$
(2)
The constitutive parameters present in the foregoing equations are complex–valued with positive imaginary parts (as befits any passive medium). The two half–spaces on either sides of the bilayer are vacuous.

Without loss of generality, the electromagnetic field phasors everywhere can be written as

$$\begin{array}{c}𝐄(𝐫)=\stackrel{~}{𝐞}(z)e^{i\kappa x}\hfill \\ 𝐇(𝐫)=\stackrel{~}{𝐡}(z)e^{i\kappa x}\hfill \end{array}\},\mathrm{}<z<\mathrm{},$$
(3)
where the transverse wavenumber $`\kappa [0,\mathrm{})`$. The fields inside the bilayer must follow the 4$`\times `$4 matrix ordinary differential equation

$$\frac{d}{dz}\left[𝐟(z)\right]=i𝒫(z)\left[𝐟(z)\right],0<z<d_a+d_b.$$
(4)
In this equation, $`\left[𝐟(z)\right]=\mathrm{col}[\stackrel{~}{e}_x(z),\stackrel{~}{e}_y(z),\stackrel{~}{h}_x(z),\stackrel{~}{h}_y(z)]`$ is a column vector, while the 4$`\times `$4 matrix function $`𝒫(z)`$ is piecewise uniform as

$$𝒫(z)=\{\begin{array}{c}𝒫_a,0<z<d_a\hfill \\ 𝒫_b,d_a<z<d_a+d_b\hfill \end{array},$$
(5)
where

$$𝒫_{a,b}=\left[\begin{array}{cccc}0& 0& 0& \frac{\kappa ^2}{\omega ϵ_0ϵ_{a,b}}+\omega \mu _0\mu _{a,b}\\ 0& 0& \omega \mu _0\mu _{a,b}& 0\\ 0& \omega ϵ_0ϵ_{a,b}+\frac{\kappa ^2}{\omega \mu _0\mu _{a,b}}& 0& 0\\ \omega ϵ_0ϵ_{a,b}& 0& 0& 0\end{array}\right].$$
(6)
The only nonzero elements of the matrixes $`𝒫_{a,b}`$ appear on their antidiagonals, of which the $`(2,3)`$ and the $`(3,2)`$ elements are relevant to $`s`$–polarized fields, and the $`(1,4)`$ and the $`(4,1)`$ elements to the $`p`$–polarized fields.

The solution of (4) is straightforward, because the matrix $`𝒫(z)`$ is piecewise uniform . Thus, the algebraic relation

$$\left[𝐟(d_b+d_a)\right]=e^{i𝒫_bd_b}e^{i𝒫_ad_a}\left[𝐟(0)\right]$$
(7)
is sufficient to solve both reflection/transmission problems as well as guided–wave propagation problems. The two matrix exponentials on the right side of (7) cannot be interchanged — unless the matrixes $`𝒫_a`$ and $`𝒫_b`$ also commute, which is possible with dissimilar materials only in quite special circumstances .

## 3 Analysis

Matrixes $`𝒫_{a,b}`$ have $`\pm \sqrt{k_0^2ϵ_{a,b}\mu _{a,b}\kappa ^2}=\pm \alpha _{a,b}`$ as their eigenvalues. Provided that $`|\alpha _{a,b}|d_{a,b}1`$ (i.e., both layers are electrically thin), the approximations

$$e^{i𝒫_{a,b}d_{a,b}}+i𝒫_{a,b}d_{a,b}$$
(8)
can be made, with $``$ as the 4$`\times `$4 identity matrix. Then

$$e^{i𝒫_bd_b}e^{i𝒫_ad_a}+i𝒫_ad_a+i𝒫_bd_be^{i𝒫_ad_a}e^{i𝒫_bd_b},$$
(9)
and the two layers in the bilayer can be interchanged without significant effect .

Let us now consider a single layer of relative permittivity $`ϵ_{eq}`$, relative permeability $`\mu _{eq}`$ and thickness $`d_{eq}`$. Quantities $`𝒫_{eq}`$ and $`\alpha _{eq}`$ can be defined in analogy to $`𝒫_a`$ and $`\alpha _a`$. Two thickness ratios are defined as

$$p_{a,b}=\frac{d_{a,b}}{d_{eq}}0,$$
(10)
in order to compare the single layer with the previously described bilayer. There is no hidden restriction on the non–negative real numbers $`p_a`$ and $`p_b`$.

Provided that $`|\alpha _{eq}|d_{eq}1`$ (i.e., the single layer is electrically thin as well), the approximation

$$e^{i𝒫_{eq}d_{eq}}+i𝒫_{eq}d_{eq}$$
(11)
can be made. Equations (9) and (11) permit us to establish the following equivalences between a bilayer and a single layer:

* s–polarization: The only nonzero field components are $`E_y`$, $`H_x`$ and $`H_z`$. Therefore, the equality of the $`(2,3)`$ elements of $`p_a𝒫_a+p_b𝒫_b`$ and $`𝒫_{eq}`$, and likewise of the $`(3,2)`$ elements, has to be guaranteed for equivalence; thus, the equations

$`p_a\mu _a+p_b\mu _b=\mu _{eq},`$ (12)
$`p_aϵ_a+p_bϵ_b\left({\displaystyle \frac{\kappa }{k_0}}\right)^2\left({\displaystyle \frac{p_a}{\mu _a}}+{\displaystyle \frac{p_b}{\mu _b}}\right)=ϵ_{eq}\left({\displaystyle \frac{\kappa }{k_0}}\right)^2{\displaystyle \frac{1}{\mu _{eq}}}`$ (13)
have to solved for $`ϵ_{eq}`$ and $`\mu _{eq}`$. We conclude therefrom that, for a given value of $`\kappa `$ and subject to the thickness restrictions $`|\alpha _{a,b,eq}|d_{a,b,eq}\stackrel{<}{}0.1`$, the bilayer and a single layer are equivalent with respect to the transformation of the $`x`$– and $`y`$– components of the fields from one exterior face to the other exterior face if

$`ϵ_{eq}`$ $`=`$ $`p_aϵ_a+p_bϵ_b`$ (14)
$`\left({\displaystyle \frac{\kappa }{k_0}}\right)^2\left[{\displaystyle \frac{p_ap_b(\mu _a\mu _b)^2}{\mu _a\mu _b(p_a\mu _a+p_b\mu _b)}}+{\displaystyle \frac{(p_a+p_b+1)(p_a+p_b1)}{p_a\mu _a+p_b\mu _b}}\right],`$
$`\mu _{eq}`$ $`=`$ $`p_a\mu _a+p_b\mu _b.`$ (15)
* p–polarization: The only nonzero field components being $`H_y`$, $`E_x`$ and $`E_z`$, the equality of the $`(1,4)`$ elements of $`p_a𝒫_a+p_b𝒫_b`$ and $`𝒫_{eq}`$ suffices, along with the equality of the $`(4,1)`$ elements of the two matrixes. For a given value of $`\kappa `$ and subject to the thickness restrictions $`|\alpha _{a,b,eq}|d_{a,b,eq}\stackrel{<}{}0.1`$, the bilayer and a single layer are equivalent if

$`\mu _{eq}`$ $`=`$ $`p_a\mu _a+p_b\mu _b`$ (16)
$`\left({\displaystyle \frac{\kappa }{k_0}}\right)^2\left[{\displaystyle \frac{p_ap_b(ϵ_aϵ_b)^2}{ϵ_aϵ_b(p_aϵ_a+p_bϵ_b)}}+{\displaystyle \frac{(p_a+p_b+1)(p_a+p_b1)}{p_aϵ_a+p_bϵ_b}}\right],`$
$`ϵ_{eq}`$ $`=`$ $`p_aϵ_a+p_bϵ_b.`$ (17)

Clearly, the constitutive parameters of the equivalent layer are functions of both $`p_a`$ and $`p_b`$; and we must point out that the sum of these two ratios need not equal unity. Furthermore, except for normal incidence (i.e., $`\kappa =0`$), the constitutive parameters of the equivalent layer depend on the incident linear polarization state. Finally, the constitutive parameters of the equivalent layer change with the transverse wavenumber $`\kappa `$.

The foregoing equations can be manipulated to yield negative values of both $`ϵ_{eq}^{}`$ and $`\mu _{eq}^{}`$ for either

* an EN–MN bilayer $`\left\{ϵ_a^{}<0,ϵ_b^{}>0,\mu _a^{}>0,\mu _b^{}<0\right\}`$ or
* a MN–EN bilayer $`\left\{ϵ_a^{}>0,ϵ_b^{}<0,\mu _a^{}<0,\mu _b^{}>0\right\}`$.

An EN–MN bilayer is equivalent to a NPV layer for both polarization states when $`\kappa =0`$, provided the condition

$$\frac{|\mu _b^{}|}{\mu _a^{}}>\frac{p_a}{p_b}>\frac{ϵ_b^{}}{|ϵ_a^{}|}$$
(18)
holds true. The inequality (18) is applicable for a MN–EN bilayer also, if the subscripts $`a`$ and $`b`$ are interchanged therein, i.e.,

$$\frac{|\mu _a^{}|}{\mu _b^{}}>\frac{p_b}{p_a}>\frac{ϵ_a^{}}{|ϵ_b^{}|}.$$
(19)
A further specialization of $`p_a+p_b=1`$ leads to the inequality

$$\frac{|\mu _b^{}|}{\mu _a^{}+|\mu _b^{}|}>p_a>\frac{ϵ_b^{}}{|ϵ_a^{}|+ϵ_b^{}}$$
(20)
for EN–MN bilayers, and

$$\frac{|\mu _a^{}|}{\mu _b^{}+|\mu _a^{}|}>p_b>\frac{ϵ_a^{}}{|ϵ_b^{}|+ϵ_a^{}}$$
(21)
for MN–EN bilayers.

The inequalities (18)–(21) should be adequate for both $`s`$– and $`p`$–polarization states when $`\kappa /k_01`$. In general, however, a given EN–MN (or MN–EN) bilayer is equivalent to a different NPV material for a different linear polarization state and/or transverse wavenumber. Thus, the equivalence between an EN–MN (or a MN–EN) bilayer and a NPV layer is restricted.

The restricted equivalence has an interesting implication for perfect lenses . A perfect lens of thickness $`d>0`$ is defined by the fulfillment of the condition $`\left[𝐟(d)\right]=\left[𝐟(0)\right]`$ for all $`\omega `$ and $`\kappa `$. Because of dispersion and dissipation, at best, this condition is fulfilled approximately. Let us imagine that the condition is fulfilled by some NPV constitutive parameters for some $`\omega `$ and all $`|\kappa |\widehat{\kappa }`$. Then, the implementation of the acceptably imperfect lens as a cascade of thin EN-MN (or MN–EN) bilayers would require that the successive bilayers have different constitutive parameters and that the entry as well as the exit faces be curved, and even those steps may not suffice. In contrast, $`\kappa `$ is fixed for any single–mode parallel–plate waveguide, and so is the range of operating frequencies; and the emulation of a NPV material via EN–MN (or MN–EN) bilayers may not be onerous.


text:
# A new approach to relevancy in Internet searching - the “Vox Populi Algorithm”

In this paper we will derive a new algorithm for Internet searching. The main idea of this algorithm is to extend the existing algorithms by a component, which reflects the interests of the users more than existing methods. The “Vox Populi Algorithm” (VPA)  creates a feedback from the users to the content of the search index. The information derived from the users query analysis is used to modify the existing crawling algorithms. The VPA controls the distribution of the resources of the crawler. Finally, we also discuss methods of suppressing unwanted content (spam). This is necessary in order to enable an efficient performance of the VPA.

The retrieval of relevant information from data sources with a very complex structure has become a challenging task since the number of documents in the Internet has reached a level of about multi billions of documents. Only a small part of them is visible in search engines. The problem of organizing and structuring these data into catalogues or searchable databases is of theoretical and significant practical (commercial) interest.
Let us define the basic components for the mathematical description of the interests of the users, the relevancy of the search results and the crawling process. The users of search engines express their needs for information through the queries which they address to a searchable database (index) $`I`$. Each of the $`k`$ queries consists of one or more keywords $`q`$ addressed to this index. It will be presented as:

$$\stackrel{}{q}_k=(q_1,\mathrm{},q_n)_k$$
(1)
$`n`$ is the length of the query $`k`$. The number of keywords per average query is $`n2`$ (status in 2003). The users are searching for documents $`d_j`$ (HTML pages, tables, text processing documents, pictures, multimedia files, …) containing information. These documents are grouped (organized) in domains $`D_k`$ presenting sets of documents under a common editorial responsibility and address (URL):

$$D_k=\stackrel{(n_k)}{}d_j^{(k)}n_k\text{ = number of documents in }D_k$$
(2)
The number of domains is about 6.4 million in Germany  and the number of documents per domain $`n_k`$ is in the interval $`10^{0\mathrm{}8}`$.
Each document $`d`$ contains searchable information, today limited to text information. Content, which is hidden for the today s search technology in non indexable formats (bitmaps, scripts etc.) will be neglected here and in the following. A document is characterized by the content of keywords $`q`$ and the position of the keyword in certain format elements $`e_i`$ (metatags, headers, tables, link text etc.):

$$d^{(k)}=f(q_1,q_2,\mathrm{},e_1,e_2,\mathrm{})$$
(3)
During the crawling and indexing process, the image of the document $`\widehat{d}`$ in the searchable index $`I`$ contains a reduced set of information - the keywords and their position in the format elements $`e`$ of the document. When a query is addressed to the index $`I`$ a ranking algorithm generates a set of documents (links) which is ordered by the relevancy of the found documents. In order to describe the document ranking process which generates the set of results on each query, one has to introduce the density $`\rho `$ of keywords within the documents:

$$\rho _i^j=\frac{n_{q_i}}{n_{e_j}}$$
(4)
where $`n_{q_i}`$ is the number of the occurrences of the keyword $`q_i`$ in the format element $`e_j`$ and $`n_{e_j}`$ is the total number of words in this format element.
Today there exist two basic types of ranking algorithms - the dynamic and the static ranking algorithms. The dynamic rank of a document depends on two factors only - the keywords $`q`$ of the query and the information content of the documents. Expressed in a ”thumb rule”: the higher the keyword density in the document the higher is the dynamic rank of this document. The relevancy function $`R_d`$, defining the dynamic rank of a document, can be written as:

$$R_d(q_1)\underset{k=1}{\overset{N}{}}\mu _k\rho ^k(q_1)\text{N - number of format elements}e$$
(5)
for a single keyword query. The coefficients $`\mu _k`$ are free parameters, defining the importance or weight of each format element. For example, the occurrence of a keyword in an URL is usually much more important than in the text itself $`\mu _{URL}>\mu _{text}`$. Queries with multiple keywords can be written as superpositions of single keyword queries:

$$R_d^n(q_1,q_2,\mathrm{},q_n)=R^1(q_1)R^1(q_2)\mathrm{}R^1(q_n)$$
(6)
Usually these functions become modified for different purposes, such as suppression of unwanted information (spam). Other modifications can take into account the freshness of the document, the type of the format or other technical parameter.
The practical work on search engines has shown that using only a document related, dynamical ranking algorithm is insufficient. In order to also include the importance or the popularity of a domain (popularity among the webmasters not necessarily among Internet users), a new type of algorithms was invented - the static ranking . The static rank $`R_s`$ of a document $`d_i`$ is related to the importance of the corresponding domain, where it is located. The idea of the static rank of a domain $`D`$ can be expressed symbolically in the following form:

$$R_s(D)\underset{j=1}{\overset{N_j}{}}R_s^j$$
(7)
where the $`R_s^j`$ is the static rank of the sites linking to the domain $`D`$. $`N_j`$ is the total amount of external links to a Domain. In  a more detailed definition of the page rank formula is given:

$$R_s(D)=(1d)+d\underset{j=1}{\overset{N_j}{}}R_s^jM_j^1$$
(8)
where $`d`$ is a free parameter (usually in the region $`d0.85`$ ) and $`M_j`$ is the total number of outgoing links of the referring site. A detailed discussion of the page rank algorithm used by Google is also found in  and .
The resulting rank of a document is a function of the the dynamic rank (5) and the static rank (7). There is no unique or even optimal way of constructing this function. A reasonable way is to choose the resulting relevancy $`R_{ds}`$ as a product of the dynamic and static rank:

$$R_{ds}=R_d(q)R_s(d_i)$$
(9)
Analyzing (9) a usual approach would be using $`R_s(D_i)`$ instead of $`R_s(d_i)`$. In practice the static rank of a document depends not only on the static rank of the domain $`D`$ containing $`d_i`$, but also on the position in the domain (link topology of the domain). At present this kind of search algorithms is in use in every major internet search engine.
The algorithms described above do indeed meet the needs of the users. This approach is reasonable from an academic point of view and it has produced remarkable results in the past. Today it has become more difficult to make use of the link topology - very often the links are not set according to the content relevancy, but for other (economic) reasons. To the extent that the search engines have become the most important information retrieval tool, they have also become a target of spamming (site owners try to fake the search engines, virtually presenting more important content than there really is). An effective method of detecting a certain type of spam is described in the appendix. Applying filter mechanisms and modifying the parameters of the dynamic and the static relevancy algorithms, one can “fine tune” the quality of the Internet search engines.
The two methods described above explicitly do not take into account the most important factor, the interest of the users searching for information. The dynamic and the static relevancy of a document are influenced by the content of the site and by the “citation” by other sites. There is no methodical component, that reflects the voice of the searching people. This will be done by the “Vox Populi Algorithm” (people‘s voice).
The main idea of the VPA is to use the information that is extractable from the user query analysis to enhance the quality of the search. This can be done in two different ways, by modifying either the ranking or the crawling algorithm. In this paper the focus is not on the ranking, but on the crawling algorithm. The crawling algorithm defines which domain and how much of the content will be included into the search index. Sites which are not included cannot be found by the best ranking algorithm. At present there is only a small fraction ($`<10\%`$) of the Internet sites indexed by the search engines. The much bigger part of the Internet (“Deep Web”) is not visible in any of the search engines.
The source of information is the analysis of the queries $`\stackrel{}{q}`$, reflecting the users interests and needs. The query set $`Q`$ may contain all single and multiple keyword queries of the users (1). Based on these queries a multidimensional tensor $`\mathrm{\Omega }`$ can be defined, containing the information of the multiple keyword correlations with the dimension $`N_{max}`$.

$$dim[\mathrm{\Omega }(Q)]=N_{max}$$
(10)
$`N_{max}`$ is the maximum length of a query - theoretically it can be infinite. Practically the amount of queries having $`>6`$ keywords is $`<1\%`$, while the average query consists of about $`N=2`$ keywords. In order to simplify the further calculations one can reduce the dimension of (10) in the following way:

$$\mathrm{\Omega }^{N_{max}}(Q)\mathrm{\Omega }^{N=2}(Q)\mathrm{\Omega }$$
(11)
In this reduction algorithm, the queries with more than two keywords are replaced by two keyword queries, containing all possible paired combinations. For example, a three keyword query is equivalent to 3 two keyword queries and so on.
The matrix $`\mathrm{\Omega }`$ is a correlation matrix of all keywords of the query set $`Q`$, which is analyzed. $`\mathrm{\Omega }`$ is a positive and symmetric matrix 111The analysis of the order of the keywords shows a statistical asymmetry for the order of keywords $`N(1,2)N(2,1)`$. Users interested in the explicit order of the keywords can use the option called “Exact Phrase”, which is available on any modern search engine. Therefore it is reasonable to assume that the order of the keywords is not important for the users when they make simple queries (more than 90% of all queries are of this type). We will use here the approximation $`(1.2)=(2.1)`$. One can calculate the eigenvectors and eigenvalues of $`\mathrm{\Omega }`$, transforming it into the diagonal form:

$$K^1\mathrm{\Omega }K=\mathrm{\Omega }^{diag}$$
(12)
The details of the diagonalization procedure are well known, see  or any other standard textbook on mathematics. It is now important to understand the practical meaning of the matrices $`K`$ and $`\mathrm{\Omega }^{diag}`$. The matrix $`K`$ consists of eigenvectors which are keyword combinations:

$$K=\left(\begin{array}{c}\stackrel{}{e}^1\\ \stackrel{}{e}^2\\ \stackrel{}{e}^3\\ \mathrm{}\end{array}\right)$$
(13)
where each eigenvector has the coordinates

$$\stackrel{}{e}^j=(c_1q_1,c_2q_2,\mathrm{})^j$$
(14)
similar to the definition (1) the $`q_i`$ are the keywords and the coefficients $`c_i^j`$ are positive numbers, giving each keyword some ”weight” compared to the other ones (How frequent do the users ask for this keyword?). The coefficients determine the relative importance of a keyword within an eigenvector. A typical eigenvector (or better “eigenquery”) has the form (based on the data , Aug. 2003).

$$\stackrel{}{e}^j=(\mathrm{`}\mathrm{`}mp3\mathrm{"},\mathrm{\hspace{0.33em}0.73}\mathrm{`}\mathrm{`}downloads\mathrm{"},\mathrm{\hspace{0.33em}0.43}\mathrm{`}\mathrm{`}free\mathrm{"},\mathrm{})$$
(15)
This query shows how the average user is asking, when he is searching for mp3 downloads at no cost. The reduced $`(N=3)`$ keyword matrix of the example above has the form :

$$\mathrm{\Omega }=\begin{array}{cccc}& mp3& download& free\\ mp3\hfill & 37.2\%& 8.8\%& 2.7\%\\ download\hfill & 8.8\%& 19.2\%& 3.6\%\\ free\hfill & 2.7\%& 3.6\%& 13.4\%\end{array}$$
(16)
The difference between the typical keyword search at present and our approach is that the words here have different weights, determining their relative importance for the users.
Another important information about the significance of keyword combinations is contained in the matrix $`\mathrm{\Omega }^{diag}`$.

$$\mathrm{\Omega }^{diag}=\left(\begin{array}{ccccc}\lambda _1& 0& 0& \mathrm{}& 0\\ 0& \lambda _2& 0& \mathrm{}& 0\\ 0& 0& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& \lambda _{N1}& 0\\ 0& 0& \mathrm{}& 0& \lambda _N\end{array}\right)$$
(17)
Each eigenvalue $`\lambda _i`$ corresponds to an eigenvector in (14). The eigenvalue can be interpreted as the importance of the corresponding eigenvector - it defines the importance of an eigenquery for the users.
Finally, we have developed the tools for defining how a search engine can use the information of the users to determine, which content should be enhanced or reduced in the index. Based on the described algorithm it is possible to define which content is the “most wanted” content and which sites deliver this type of content:

$$(c_1q_1+c_2q_2+\mathrm{})\text{search engine}\text{list of ranked domains}$$

Crawling the Internet, each domain is given certain resources by the search engine, such as CPU time and memory in the index (alternatively also the number of crawled documents or other parameters, depending on the settings of the search engine).
The practical realization of the VPA as an extension of an existing Internet search could be performed using the following procedure:

1. Generate a ranking of domains, addressing the eigenqueries (14) to the existing (old) search index, the priority of those domains is defined by the size of eigenvalues. (17).
2. Modify the existing resource ranking list with respect to these eigenvalues.
3. Use the new determined ranking of the domains for crawling the Internet according to the modified resource distribution.
4. Repeat the cycle.

In order to determine which sites best fit the eigenqueries, it is useful to calculate a dynamic rank for a whole domain, not just for a single document. A simple method would be to summarize the total score of all documents in one domain:

$$𝐑_𝐃(e_i)\underset{k=1}{\overset{N_D}{}}R_d^k(e_i)$$
(18)
Let us assume, that the amount of resources (CPU time, number of documents, data volume etc.) given to each domain, when crawling it, can be expressed in a function $`M`$, with

$$M=M(D_k,R_s,\mathrm{})$$
(19)
In order to apply the VPA one can modify (19) in the following way:

$$M\widehat{M}=MR_{VPA}$$
(20)
The function $`R_{VPA}`$ defines the VPA correction with regard to the old crawling algorithm. The function $`R_{VPA}`$ can be presented in different ways. The basic requirement for the function is that it is monotone concerning the parameters $`\lambda _i`$, which define quantitatively how relevant a query is for the users. Following Occam‘s principle of simplicity (Pluralitas non est ponenda sine neccesitate - Entities should not be multiplied unnecessarily) this function should use only a minimum set of free parameters, which will allow the adoption (or “fine tuning”) the algorithm to the local requirements:

$$R_{VPA}(D_k)=\left(1+\alpha \lambda _k^\beta \right)\alpha ,\beta >0$$
(21)
The parameter $`\alpha `$ and $`\beta `$ can be chosen freely. In the limit, the new algorithm generates the existing results in (20).

$$\underset{\lambda 0}{lim}\widehat{M}=M$$
(22)
In this paper we have shown how the analysis of queries can be used to enhance the relevant and “most wanted” content in a search index. In this way the relevancy, experienced by the users of the search should grow - the users will find more of what they are interested in. The existing system of the relevancy ranking of documents or domains can remain unchanged. The algorithm will not replace existing crawling and ranking algorithms, but the VPA will extend them by a qualitatively new component.

Appendix
The static rank algorithm has also become the target of spamming (for example, “Google bombing” ). This means that webmasters are creating clusters of domains, which consist of very similar sites, referring to a single domain or a document. This kind of spam cluster can consist of many domains, which do not contain any valuable content at all. Because of this the static rank consequently is becoming more and more a measure of the marketing budget or the cleverness of the webmaster of a domain, rather than a measure of “real” reputation or content quality. As a result of this development, the importance of the static rank as a tool for determining the quality or the relevancy of a site is decreasing.
We want to propose an algorithm which identifies this kind of spamming. The basic idea of the static rank is reasonable - the more important sites refer (link) to a site, the more important is the site. There is a way to discriminate between “natural grown” link clusters and “artificial” ones (spam).
In order to find a quantitative method which can discriminate between these two types of link clusters, one can introduce the function which describes the statistical distribution of the relevancy $`R_s^j`$ of the links, pointing to the document $`d_i`$:

$$\varphi (R_s^j)=e^{\frac{(R_s^jR_0)^2}{\sigma ^2}}$$
(23)
here $`R_0`$ is the average static rank of all sites, linking to the center of this cluster $`d_i`$. The parameter $`\sigma `$ defines the width of the distribution.
The above mentioned types of clusters can be discriminated using the distribution $`\varphi `$ \- natural grown clusters contain links from an inhomogeneous set of sites, for example, the links to a site of a well known university will come from very small (amateur) sites of students, employees and alumnies (with a low page rank), via semi professional institutional sites (spin offs, research partners, …) up to sites of other high ranked universities or institutes. The artificial link cluster consists of automatically generated sites, each of them usually optimized for different keywords, but having approximately the same static rank. As a result of this it is possible to introduce a “cut off” criteria based on formula (23). A cluster is most likely spam, if the condition

$$\sigma _{spam}<\sigma _{critical}$$
(24)
is fulfilled. Here $`\sigma _{critical}`$ is an empirical parameter, which can be determined from the analysis of known natural and artificial clusters (or from the software generating the sites of the spam cluster). Estimates have shown that one can expect a result like $`\sigma _{natural}>>\sigma _{artificial}`$. A short test example can demonstrate this: the distribution of the page ranks of sites linking to the homepage of Steven Hawking  analyzed based on formula (23) have a width of $`\sigma ^2=1.1`$ , while the sites belonging to a typical spam cluster have a page rank distribution with $`\sigma ^2=0.5\mathrm{}0.7`$ 222The data of this example are based on the page rank indicator of Google .. The parameter $`\sigma `$ can be used for separating between these two type of link clusters. The data of this example are based on the indications of the page rank indicator of Google‘s toolbar .


text:
# Route to nonlocality and observation of accessible solitons

## Abstract

We develop a general theory of spatial solitons in a liquid crystalline medium exhibiting a nonlinearity with an arbitrary degree of effective nonlocality. The model accounts the observability of accessible solitons and establishes an important link with parametric solitons.

In various areas of applied nonlinear science, nonlocality plays a relevant role and radically affects the underlying physics. Some striking evidences are found in plasma physics,Litvak and Sergeev (1978); Pecseli and Rasmussen (1980); Turitsyn (1985) or in Bose-Einstein condensates (BEC), Parola et al. (1998); Perez-Garcia et al. (2000) where, contrary to the prediction of purely local nonlinear models, nonlocality may give rise to, or prevent, the collapse of a (plasma- or matter-)wave. In nonlinear-optics, particularly when dealing with self-localization and solitary waves, nonlocality is often associated to time-domain phenomena through a retarded response (see e.g. Vanin et al. (1994); Akhmediev et al. (1998)); spatially nonlocal effects have been associated to photorefractive Duree et al. (1993); Mamaev et al. (1997); Abe and Ogura (1998); DelRe et al. (2001) and thermal or diffusive responses Dreischuh et al. (1999); Henninot et al. (2002). To assess the role of nonlocality, theoretical studies tend to distinguish between highly and weakly nonlocal behaviors Snyder and Mitchell (1997); Mitchell and Snyder (1999); Krolikowski and Bang (2000), by comparing the spatial extent of the material response (the so-called kernel-function) and the optical beam waist. Specific kernel functions, however, strongly depend on the physical system and, as in the case of BEC, Perez-Garcia et al. (2000) they are hard to determine and apply to experimental results. On the other hand, they are at the basis of the theory of spatial optical solitons (SOS) in highly nonlocal media.
SOS have become the subject of intense theoretical and experimental investigations, both on the grounds of their packet nature and in view of applications, particularly in the exploitation of their wave-guiding character. Boardman and Sukhorukov (2001) Diverse material properties have been studied in conjunction with SOS existence and properties, including various mechanisms able to counteract diffraction in one or both transverse dimensions. Stegeman et al. (2000); Trillo and Torruellas (2001) Spatial solitons due to a local nonlinearity have been known since the original work of Chiao, Garmire and Townes with reference to Kerr media. Chiao et al. (1964) In 1997 Snyder and Mitchell Snyder and Mitchell (1997) investigated SOS in a highly nonlocal system, i. e., a medium exhibiting a power - rather than intensity - dependent refractive index. They introduced the term accessible solitons for those spatial solitary waves, owing to the simplicity of the theory and transverse profiles obeying the 2-dimensional equation of a quantum harmonic oscillator (which gives a gaussian profile like the so-called “gaussons” Oficjalski and Bialynicki-Birula (1978)). Shen pointed out to this connection as an intriguing one between distinct fields of modern physics and, provided the required large correlation lengths could be made available, to the demonstration of accessible solitons as a challenge well worth undertaking. Shen (1997)
In this Letter we introduce a model able to describe optical spatial solitons and the smooth transition from the purely local (in the limit of a Kerr nonlinearity) to the entirely nonlocal case. While pursuing a general theory, however, we chose to address a specific and available nonlocal system, i. e., nematic liquid crystals (NLC) in a planar cell. NLC have been proven to exhibit a substantial nonlocal nonlinearity of molecular origin Khoo (1995) and to support (2+1)-dimensional spatial solitons, Peccianti et al. (2000) even in the case of spatially incoherent excitations. Peccianti and Assanto (2002) After deriving the ruling equations and defining a suitable nonlocal parameter able to span the soliton-family from pure-Kerr or Townes-like (T-) to higly nonlocal or accessible (A-) solitons, we will outline the rather remarkable connection between our model and the equations describing quadratic two-color solitons (or simultons) in parametric media. Karamzin and Sukhorukov (1976)
Let us consider the simple geometry sketched in figure 1: a planar glass cell containing an undoped NLC with a preset orientation of its molecular director. The aligned liquid crystal, anchored at the bounding interfaces, behaves as a positive uniaxial with $`n_{//}`$=$`n_{\widehat{z}}`$ and $`n_{}`$=$`n_{\widehat{x}}`$. In the presence of an external quasi-static (electric or magnetic) field or special anchoring at the interfaces, the refractive index $`n(\theta )`$ in the (x, z) principal plane can exhibit a distribution along $`x`$, with $`n_{}<n(\theta (x))<n_{//}`$. Tabirian et al. (1986) For a light beam linearly polarized along $`x`$ and propagating along z, with transverse size well below the cell thickness L, neglecting vectorial effects and adopting the paraxial approximation, the evolution of the optical envelope $`A`$ is described by the Foch-Leontovich equation:

$$2ik\frac{A}{Z}+_{XY}^2A+k_0^2n_a^2[sin(\theta )^2sin(\theta _0)^2]A=0$$
(1)
where $`k_0`$ is the vacuum wavenumber, $`n_a^2`$=$`n_{//}^2`$ \- $`n_{}^2`$ the optical anisotropy, $`k^2=k_0^2(n_{}^2+n_a^2sin(\theta _0^2)^2)`$, $`\theta `$ the tilt angle of the NLC director, and $`\theta _0`$ the tilt in the absence of a light beam. When an external electric field applied along $`x`$ and an optical excitation as in equation (1) are present, $`\theta `$ is subject to re-orientation according to Khoo (1995)

$$K\frac{^2\theta }{Z^2}+K_{XY}^2\theta +\frac{\mathrm{\Delta }ϵ_{RF}E^2}{2}sin(2\theta )+\frac{ϵ_0n_a^2|A|^2}{4}sin(2\theta )=0$$
(2)
with $`K`$ the relevant elastic constant taken equal for splay, bend and twist, $`\mathrm{\Delta }ϵ_{RF}`$ the low-frequency anisotropy, and $`E`$ the rms value of the quasi-static field.
In the absence of a light wave, therefore, the orientation angle $`\widehat{\theta }`$ is determined exclusively by $`E`$ and, due to simmetry, only depends on X:Khoo (1995)

$$K\frac{d^2\widehat{\theta }}{dX^2}+\frac{\mathrm{\Delta }ϵ_{RF}E^2}{2}sin(2\widehat{\theta })=0\text{.}$$
(3)
In our case the boundary conditions correspond to the planar alignment: $`\theta (X=L/2)=\theta (X=L/2)=0`$. In the general case, the angle distribution can be written as:

$$\theta (X,Y,Z)=\widehat{\theta }(X)+\frac{\widehat{\theta }(X)}{\theta _0}\mathrm{\Psi }(X,Y,Z)\text{.}$$
(4)
Taking the cell much larger than the beam waist, we can use (3) and (4) in (2) and neglect the derivative of $`\widehat{\theta }`$. For $`\widehat{\theta }\theta _0`$, at the first order in $`\mathrm{\Psi }`$ we obtain:

$$\begin{array}{c}2ik\frac{A}{Z}+_{XY}^2A+k_0^2n_a^2\mathrm{\Psi }A=0\hfill \\ K_{XYZ}^2\mathrm{\Psi }\frac{2\mathrm{\Delta }ϵ_{RF}E^2}{\pi }\mathrm{\Psi }+\frac{ϵ_0n_a^2}{4}|A|^2=0\text{,}\hfill \end{array}$$
(5)
having chosen $`\theta _0=\pi /4`$ in order to maximize the nonlinear response. Tabirian et al. (1986) We write (5) in a dimensionless form by setting $`A=(A_c/\alpha )a(R/R_c\sqrt{\alpha },Z/\alpha Z_c)exp(iZ/\alpha Z_c)`$, $`\mathrm{\Psi }=(\mathrm{\Psi }_c/\alpha )\psi (R/R_c\sqrt{\alpha },Z/\alpha Z_c)`$ with $`A_c^2=8\mathrm{\Delta }ϵ_{RF}^2E^4/\pi ^2ϵ_0k_0^2n_a^4K`$, $`Z_c=2kR_c^2`$,$`R_c^2=\pi K/2\mathrm{\Delta }ϵ_{RF}E^2`$, $`\mathrm{\Psi }_c=2\mathrm{\Delta }ϵ_{RF}E^2/\pi k_0^2n_a^2K`$, and $`ϵ=\mathrm{\Delta }ϵ_{RF}E^2/2\pi k^2K`$. $`\alpha `$ is a free parameter, to be used hereafter to trace the family of spatial solitary waves, and $`(x,y,z)=(X/\sqrt{\alpha }R_c,Y/\sqrt{\alpha }R_c,Z/\alpha Z_c)`$ are normalized coordinates. The resulting system is:

$$\begin{array}{c}i\frac{a}{z}+^2aa+a\psi =0\hfill \\ \frac{ϵ}{\alpha }\frac{^2\psi }{z^2}+^2\psi \alpha \psi +\frac{1}{2}|a|^2=0\text{.}\hfill \end{array}$$
(6)
To underline the physical meaning of the parameter $`\alpha `$, let us consider the case $`ϵ=0`$ and formally write the solution of the re-orientation equation as $`\psi =(1^2/\alpha )^1|a|^2/(2\alpha )`$. For large $`\alpha `$, we obtain for the field envelope:

$$i\frac{a}{z}+^2a+a(1+\frac{^2}{\alpha })\frac{|a|^2}{2\alpha }=0\text{.}$$
(7)
Equation (7) rules collapse-free weakly-nonlocal media,Turitsyn (1985), while the neglecting of terms such as $`O(1/\alpha ^2)`$ describes light propagation in Kerr media which, on the contrary, are subject to catastrophic self-focusing. Reducing $`\alpha `$ increases the degree of nonlocality in the interaction between medium and optical beam. In the following we will employ $`\alpha `$ as an arbitrary parameter spanning the whole family of SOS; we expect large $`\alpha `$ values to be associated to T-solitons, whereas small $`\alpha `$ will address A-solitons. Note that, according to Turitsyn (1985), when $`ϵ=0`$ the fundamental solitary wave solutions of Eqs.(6) are stable because they realize an absolute minimum of the Hamiltonian.

Solitary solutions of (6) are defined by $`_z=0`$. Without loss of generality, taking $`a`$ real-valued we have:

$$\begin{array}{c}^2aa+a\psi =0\hfill \\ ^2\psi \alpha \psi +\frac{a^2}{2}=0\text{.}\hfill \end{array}$$
(8)
Noteworthy enough, system (8) is identical to what determines the profile of parametric spatial solitary waves in $`\chi ^{(2)}`$-media. Trillo and Torruellas (2001); Karamzin and Sukhorukov (1976) This makes an unexpected connection between self-trapped beams in two distinct physical systems encompassing rather diverse nonlinear optical responses: the ultrafast electronic nonlinearity of quadratic crystals and the slow molecular re-orientation of liquid crystals. The similarity is only limited to the profiles of the solitary waves, while their dynamic properties, such as stability, are quite different. As shown below in a practical case, to a first approximation $`ϵ`$ is negligible, thus the angle “adiabatically-follows” the light beam; in general this is not true for the harmonic-field of a parametric soliton. It is well known, in fact, that for small $`\alpha `$ the $`\psi `$-field (the second-harmonic for $`\chi ^{(2)}`$ crystals, and the re-orientation for NLC) is much wider than the $`a`$-field (see Torruellas et al. in Trillo and Torruellas (2001)). The opposite holds true for $`\alpha \mathrm{}`$: in the $`\chi ^{(2)}`$-SOS literature, this is the Kerr-limit, its dynamics resembling that of $`\chi ^{(3)}`$ materials. In NLC, for optical beams much wider than the re-orientation profile, T-solitons approximate well the solution (in the framework of the validity of the large cell approximation). The angle-profile perturbation is localized close to the optical beam axis, as typically pointed out when addressing the Kerr-like response of NLC. Tabirian et al. (1986); Khoo (1995) Here we are rather interested in the opposite limit.

In figure 2 we graph the numerical solutions of eqs. (6) as obtained by a relaxation procedure. The ratio $`\rho `$ between SOS beam- and angle- waists (standard deviation) is plotted versus $`\alpha `$. When $`\alpha `$ approaches zero, i. e., for a response-length extending well beyond the optical waist, the $`\psi `$-field is much wider than the $`a`$-field. This is the high-nonlocality regime. In the figure we indicate the two opposite limits, namely A- and T-solitons. In the A-limit, however, it is worth proving that eqs. (6) reduce to the model in Snyder and Mitchell (1997). The authors in Snyder and Mitchell (1997) postulated a highly nonlocal medium in which the index of refraction could be expressed as $`n^2=n_0^2\alpha _0^2(𝒫)R^2`$, with $`𝒫`$ the optical power. Here we show that the NLC re-orientational nonlinearity does indeed exhibit such feature, as speculated rather skeptically in Shen (1997).

To study the highly non-local regime ($`\alpha 0`$), first we solve eqs. (8), in the regions around $`r=0`$ ($`r=\sqrt{x^2+y^2}`$) and $`r\mathrm{}`$, and then match the resulting expressions. For $`r0`$ we introduce the expanded variables Nayfeh (1993) $`\zeta =r/\alpha ^{1/2}`$ and $`\xi =r/\alpha ^{1/4}`$. Next, we express the field as $`a=\frac{\widehat{a}}{\sqrt{\alpha }}f(\xi )`$ with $`f(0)=1`$, while $`\psi =\psi ^{(0)}(\zeta )+(1/\sqrt{\alpha })\psi ^{(1)}(\zeta )`$. At order $`O(1/\alpha ^{1/2})`$ we obtain for $`\psi `$:

$$\psi _{\zeta \zeta }^{(0)}+(1/\zeta )\psi _\zeta ^{(0)}+(1/2)\widehat{a}^2f^2(\alpha ^{1/4}\zeta )=0$$
(9)
For $`\alpha 0`$ with $`\zeta `$ fixed, the solution reads $`\psi ^{(0)}=\psi _0\widehat{a}^2r^2/8\alpha =\psi _0\widehat{a}^2\xi ^2/8\sqrt{\alpha }`$. At order $`O(1/\alpha ^{3/2})`$ we have $`\psi _\zeta ^{(1)}=0`$; at order $`O(1/\alpha ^{1/2})`$, the field equation reduces to the harmonic oscillator equation in $`f`$:

$$f_{\xi \xi }+(1/\xi )f_\xi +\psi ^{(1)}f(1/8)\widehat{a}^2\xi ^2f=0$$
(10)
corresponding to the result in Snyder and Mitchell (1997), and $`\psi _0=1`$ at order $`O(1)`$. After some algebra, the perturbative approach provides the SOS profile near the origin ($`\psi ^{(1)}=2/\xi _0^2`$):

$$\begin{array}{c}a=\frac{2\sqrt{2}}{\sqrt{\alpha }\xi _0^2}e^{\frac{r^2}{2\sqrt{\alpha }\xi _0^2}}\hfill \\ \psi =1+\frac{2}{\sqrt{\alpha }\xi _0^2}\frac{r^2}{\alpha \xi _0^4}\hfill \end{array}$$
(11)
with $`\xi _0`$ an arbitrary parameter. The appearance of another free parameter reflects that, when $`\alpha =0`$, eqs. (6) are invariant with respect to the transformation $`(x,y)(x/\mu ,y/\mu )`$,$`a\mu a`$,$`\psi 1\mu +\mu \psi `$, with $`\mu `$ arbitrary. The physical meaning of (11) can be elucidated by expressing them in terms of the normalized beam power $`P=a^2𝑑x𝑑y`$: $`a=(\sqrt{2}P/4\pi )exp(r^2P/16\pi )`$, and $`\psi =1+(P/4\pi )(P/8\pi )r^2`$. The latter results can be also obtained by using a multiple scales expansion in $`\sqrt{P}r`$, $`Pr`$, $`\mathrm{}`$. For the field intensity in the original variables($`R=\sqrt{X^2+Y^2}`$), we have:

$$=\frac{|A|^2}{2\eta }=\frac{𝒫^2}{\pi R_c^2𝒫_c}exp(\frac{R^2}{R_c^2}\frac{𝒫}{𝒫_c})$$
(12)
In (12) $`Z_0`$ and $`\eta =Z_0/n`$ are vacuum and medium impedances, respectively, $`c`$ the speed of light in vacuum, $`n=k/k_0=\sqrt{n_{}+n_a^2/2}`$ when $`\theta _0=\pi /4`$, and $`𝒫_c`$ a reference power which depends on material properties and cell-polarization: $`𝒫_c=16cn\mathrm{\Delta }ϵ_{RF}E^2/k_0^2n_a^4`$. The soliton profile is gaussian and completely determined by $`𝒫`$. Rendering explicit the relation between $`𝒫`$ and the (intensity) waist $`𝒲`$, its existence curve is:

$$\frac{𝒫}{𝒫_c}=\frac{R_c^2}{𝒲^2}\text{.}$$
(13)
Let us now investigate the region of large $`r`$. According to eqs. (6) the asymptotic behavior of $`\psi `$ as $`r\mathrm{}`$ is given by the modified Bessel function: $`\psi GK_0(\sqrt{\alpha }r)`$, governing the angle decay far from the beam axis ($`G`$ is a constant to be determined). The angle in (11) and its derivative must match the expression at infinity, providing the turning point $`r_T`$ between the two regions: in the presence ($`r0`$) and in the absence ($`r\mathrm{}`$) of the optical excitation, respectively. We end up with

$$1+\frac{P}{4\pi }r_T^2(\frac{P}{8\pi })^2=(\frac{P}{8\pi })^2\frac{2r_TK_0(\sqrt{\alpha }r)}{\sqrt{\alpha }K_1(\sqrt{\alpha }r)}\text{,}$$
(14)
which can be further simplified by taking into account that, for large arguments, the Bessel functions $`K_0`$ and $`K_1`$ have the same asymptotic behavior. The approximated solution is $`r_T=\sqrt{\alpha }(8\pi /P)`$; being $`_T=\sqrt{\alpha }R_cr_T`$ it reads: $`_T/_c=𝒫_c/𝒫`$. When $`𝒫>>𝒫_c`$ the profile of the angle is dominated by the modified Bessel function, and decays with typical length $`R_c`$. Therefore, it seems natural to label “highly nonlocal” a regime when the angle profile is much larger than the beam waist, such that $`_c>>𝒲`$. When the power is much greater than $`𝒫_c`$ or, equivalently, when the soliton waist is much smaller that the extent $`_c`$ of the elastic response, we are in the A-soliton regime. The power dependent perturbation of the refractive index is:

$$\mathrm{\Delta }n(R,𝒫)=n_c(\frac{2𝒫}{𝒫_c}\frac{R^2}{R_c^2}\frac{𝒫^2}{𝒫_c^2})\text{,}$$
(15)
with $`n_c=\mathrm{\Delta }ϵ_{RF}E^2/\pi nk_0^2K`$. Moreover, our approach enables to go beyond the harmonic oscillator. By introducing the new transverse scale $`\sigma =r/\alpha ^{3/8}`$ we can solve the resulting equation for $`\psi ^{(0)}`$ at order $`O(\alpha ^{1/4})`$: $`\psi _{\sigma \sigma }^{(0)}+(1/\sigma )\psi ^{(0)}4\sigma ^2/\xi _0^6=0`$. After some algebra we obtain the new approximation in normalized units: $`\psi =1+(P/4\pi )(P/8\pi )r^2+(P/8\pi )^3r^4/4`$. Then, the corresponding power-dependent refractive index perturbation reads:

$$\mathrm{\Delta }n(R,𝒫)=n_c[\frac{2𝒫}{𝒫_c}\frac{R^2}{R_c^2}(\frac{𝒫}{𝒫_c})^2+\frac{R^4}{4R_c^4}(\frac{𝒫}{𝒫_c})^3]\text{.}$$
(16)
Thus higher-order approximations imply anarmonicity of the nonlocal potential, with a refractive index still depending on power. As the latter increases, higher powers of the ratio $`𝒫/𝒫_c`$ must be taken into account, similar to local media with regards to the powers of the intensity.

Finally, comparing our theory with an actual experimental geometry, such as employed in Peccianti et al. (2000), typical parameters for a $`514nm`$ wavelength and the E7 NLC (in SI units) are: $`n_a=1`$, $`K=10^{11}`$, $`E=1.3\times 10^4`$,$`L=75\times 10^6`$, $`\mathrm{\Delta }ϵ_{RF}=20ϵ_0`$. Correspondingly, $`𝒫_c=2\times 10^6W`$, $`R_c=22\mu m`$, and $`ϵ=10^6`$. Since $`𝒫`$ inside the cell was of the order of $`0.1mW`$, we may state that in experiments the highly nonlocal regime ($`𝒫>>𝒫_c`$) is being addressed. The soliton waist, after (13), is of the order of $`3\mu m`$ and in agreement with the reported results. Note that $`\mathrm{\Delta }n(0,𝒫)5\times 10^4`$, and the angle perturbation $`\mathrm{\Psi }`$ of the order of $`10^3`$ radians, thus justifying the adopted model.

In conclusions, for the first time to our knowledge, we have presented a self-consistent analytical theory of two-dimensional spatial solitary waves in nonlocal media. Our model has an intriguing unifying character, as it embraces several physical systems in which light self-trapping has recently been investigated. We believe that most of the observed SOS’ in nematic liquid crystals are indeed accessible solitons, in-as-much-as NLC are highly nonlocal. This shines new light on self-localization in liquid crystals. Furthermore, we have presented the first derivation of a power-dependent constitutive relation for a real physical system, never reported elsewhere. We confide that our results will stimulate new experiments towards a deeper understanding of self-trapping in (highly) nonlocal nonlinear media and the development of novel all-optical devices. Peccianti et al. (2002)

###### Acknowledgements.

C. C. aknowledges the Tronchetti Provera Foundation for a generous grant. This work was supported by INFM (Advanced Research Project “SPASONELIC”).


text:
# Unexpected non-Wigner behavior in level-spacing distributions of next-nearest-neighbor coupled XXZ spin chains

## I Introduction

Random matrix theories have been successfully applied to the energy spectra of not only chaotic systems but also quantum spin systems. Montam ; Hsu ; Poil ; Pals ; 8vertex ; Meyer ; Meyer2 ; Angles ; AM  If a given Hamiltonian is integrable by the Bethe ansatz, the level-spacing distribution should be described by the Poisson distribution:

$$P_{\mathrm{Poi}}(s)=\mathrm{exp}(s).$$
(1)
If it is not integrable, the level-spacing distribution should be given by the Wigner distribution:

$$P_{\mathrm{Wig}}(s)=\frac{\pi s}{2}\mathrm{exp}\left(\frac{\pi s^2}{4}\right).$$
(2)
These behaviors of the level-spacing distribution have been observed in one-dimensional (1D),Hsu ; Poil ; Angles  2D,Pals ; Montam ; Meyer2  and 3DMeyer2  spin systems. They are also confirmed for strongly correlated systems,Alt  and applied to the recent study of quantum dots. Held

The numerical observations Hsu ; Poil ; Angles ; Pals ; Montam ; Meyer2 ; Alt ; Held  are important. In fact, for the quantum spin systems, there has been no theoretical or analytical derivation of the suggested behaviors of the level-spacing distribution. Random matrix theories do not necessarily extend to them. However, it seems that no counterexample has been shown explicitly for quantum spin systems. The suggested behaviors have been numerically confirmed for quite a large number of quantum spin systems and statistical lattice systems. Hsu ; Poil ; Angles ; Pals ; Montam ; Meyer2 ; AM  Moreover, the behaviors are often considered as empirical rules, by which we can practically determine whether a given lattice system is integrable or not AM . Thus, it could be nontrivial if one would find such a quantum spin system that does not show the standard behaviors of the level-spacing distribution.

Recently, it has been found that an extraordinary symmetry appears for special cases of the integrable XXZ spin chain: the XXZ Hamiltonian commutes with the $`sl_2`$ loop algebra at some particular values of the XXZ coupling associated with roots of unity. DFM  The loop algebra is an infinite-dimensional Lie algebra, which holds only at the particular values of the XXZ coupling and not for generic values. It is suggested numerically FM  that the standard Bethe ansatz does not hold for the special cases. Furthermore, the dimensions of degenerate eigenspaces of the loop algebra are given by some exponential functions of the system size, Deguchi  and they can be extremely large. It should therefore be nontrivial how the large degeneracies are resolved by nonintegrable perturbative terms in the spin Hamiltonian. Thus, the $`sl_2`$ loop algebra symmetry may motivate us to reconsider not only the Bethe ansatz solvability of the integrable XXZ chain but also the standard statistical behaviors of energy levels for various XXZ chains close to the integrable one. We note that some level crossings of the $`sl_2`$ loop algebra are shown in Fig. 1 of Ref. KD, .

In this paper, we discuss level-spacing distributions for the spin $`\frac{1}{2}`$ XXZ chains on finite sites under periodic boundary conditions. We mainly discuss nonintegrable cases, although the present study has been originally motivated by the $`sl_2`$ loop algebra symmetry of the integrable XXZ chain. We first confirm that integrable spin chains show the Poisson distribution for generic cases. Here we exclude the special cases related to roots of unity. Then, we consider the next-nearest-neighbor (NNN) coupled Heisenberg spin chain (or XXX spin chain), which is nonintegrable. We now note that the XXX chain with NNN couplings has the Wigner distribution, as shown in Ref. Poil, . For the XXZ spin chains with NNN couplings, however, we find that the level-spacing distributions are not given by the standard Wigner distribution. The observation should be nontrivial since the systems are nonintegrable. Finally, we discuss possible reasons why the non-Wigner distribution is obtained for the NNN coupled XXZ spin chains.

## II Numerical procedure

The Hamiltonian matrices can be separated into some sectors; in each sector, the eigenstates have the same quantum numbers. This is because the eigenvalues with different symmetries are uncorrelated. The XXZ chains have three trivial symmetries, namely, reflections, translations, and spin rotations around the $`z`$ axis. Therefore, to desymmetrize the Hamiltonians, we consider three quantum numbers: parities, momenta, and total $`S^z`$. We calculate the eigenvalues of the largest sectors. The largest sector of each Hamiltonian has 440 eigenvalues for the lattice size $`L=16`$. To calculate the eigenvalues, we use standard numerical methods, which are contained in the LAPACK library.

To find universal statistical properties of the Hamiltonians, one has to deal with unfolded eigenvalues instead of raw eigenvalues. The unfolded eigenvalues are renormalized values, whose local density of states is equal to unity everywhere in the spectrum. In this paper, the unfolded eigenvalues $`x_i`$ are obtained from the raw eigenvalues $`E_i`$ in the following method. Let us define the integrated density of states as

$$n(E)=\underset{i=1}{\overset{N}{}}\theta (EE_i).$$
(3)
Here $`\theta (E)`$ is the step function and $`N`$ is the number of the eigenvalues. We choose some points of coordinates: $`(E_i,n(E_i))`$ for $`i=1,21,41,\mathrm{},N`$. The average of integrated density of states $`n(E)`$ is approximated by the spline interpolation through the chosen points. The unfolded eigenvalues are defined as

$$x_i=n(E_i).$$
(4)
The level-spacing distributions are given by the probability function $`P(s)`$, where $`s=x_{i+1}x_i`$.

## III Next-nearest-neighbor coupled XXZ chain

Let us introduce the spin $`\frac{1}{2}`$ XXZ chain with NNN couplings on $`L`$ sites by

$``$ $`=`$ $`J_1{\displaystyle \underset{l=1}{\overset{L}{}}}\left(S_l^xS_{l+1}^x+S_l^yS_{l+1}^y+\mathrm{\Delta }_1S_l^zS_{l+1}^z\right)`$
$`+`$ $`J_2{\displaystyle \underset{l=1}{\overset{L}{}}}\left[\alpha \left(S_l^xS_{l+2}^x+S_l^yS_{l+2}^y\right)+\mathrm{\Delta }_2S_l^zS_{l+2}^z\right],`$
where $`S^a=(1/2)\sigma ^a`$ and $`(\sigma ^x,\sigma ^y,\sigma ^z)`$ are the Pauli matrices; periodic boundary conditions are imposed. For simplicity, we put $`J_1=1`$ hereafter in the paper. The Hamiltonian is nonintegrable for $`J_20`$, while it is integrable for $`J_2=0`$.

Let us confirm the Poissonian behavior for the generic case of the integrable XXZ spin chain. When $`J_2=0`$, the level-spacing distribution $`P(s)`$ mostly shows the Poisson distribution as shown in Fig. 1. We confirmed numerically the standard result for some generic values of the XXZ coupling $`\mathrm{\Delta }_1`$ ($`0\mathrm{\Delta }_1<1`$). We note that we exclude the special values of the XXZ coupling that are given by $`\mathrm{\Delta }_1=\mathrm{cos}(m\pi /N)`$ for some integers $`m`$ and $`N`$, where $`\mathrm{\Delta }_1`$ is related to a root of unity, $`q`$, through the relation $`\mathrm{\Delta }_1=(q+1/q)/2`$. DFM

In nonintegrable cases ($`J_20`$), however, $`P(s)`$ shows a non-Wigner distribution against expectations. The non-Wigner behavior of the NNN coupled XXZ spin chains is generic. We have calculated $`P(s)`$ for various combinations of $`J_2`$ and $`\mathrm{\Delta }_2`$ when $`\alpha =1`$ and for various values of $`\mathrm{\Delta }_2`$ when $`\alpha =0`$. Here, $`\mathrm{\Delta }_1=1`$ for simplicity. When $`\alpha =1`$ and $`\mathrm{\Delta }_2=0.98`$ , $`P(s)`$ looks similar to the Poisson distribution rather than the Wigner distribution for any value of $`J_2`$, as shown in Figs. 2(a) and 2(c). The Hamiltonian is close to but not exactly the same as the Heisenberg chain, when $`\mathrm{\Delta }_2=0.98`$. We find such a Poisson-like distribution when $`\mathrm{\Delta }_2`$ is very close to $`1.0`$ such as more than about $`0.95`$. For the NNN coupled Heisenberg chain ($`\mathrm{\Delta }_2=1`$), however, it was shown that $`P(s)`$ has the Wigner distribution. Hsu ; Poil

When $`\mathrm{\Delta }_2`$ decreases, $`P(s)`$ becomes a non-Poisson distribution as shown in Figs. 2(b) and 2(d). However, it is not a Wigner-like distribution, either. The level-spacing distribution $`P(s)`$ strongly depends on $`\mathrm{\Delta }_2`$. On the other hand, it seems that $`P(s)`$ is almost independent of $`J_2`$.

When only one parameter is nonzero, namely, for the case of either $`\alpha =0`$ or $`\mathrm{\Delta }_2=0`$, the nonzero parameter does not change $`P(s)`$ very much. In that case, $`P(s)`$ is given by neither a Poisson-like nor Wigner-like distribution.

For all the investigated nonintegrable cases except for the case where $`\mathrm{\Delta }_2`$ is close to $`1.0`$ and $`J_20`$, $`P(s)`$ can be approximated by the arithmetic average of the Poisson and Wigner distributions:

$$P_{\mathrm{av}}(s)=\left[P_{\mathrm{Poi}}(s)+P_{\mathrm{Wig}}(s)\right]/2,$$
(6)
as shown in Fig. 3.

Let us remark on the special case of the integrable XXZ chain related to a root of unity. At the special values of $`\mathrm{\Delta }_1`$, we find that $`P(s)`$ shows a novel peak at $`s=0`$ as shown in Fig. 4. For $`L=16`$, large degeneracies remain for $`\mathrm{\Delta }_1=\mathrm{cos}(\pi /i)`$, $`i=2,3,4`$, even after the desymmetrization procedure mentioned in Sec. II is performed. We recall that for generic values of $`\mathrm{\Delta }_1`$, no degeneracy remains after the desymmetrization procedure is completed and the Poisson distribution is obtained as shown in Fig. 1. The appearance of the peak in $`P(s)`$ at $`s=0`$ should be consistent with the $`sl_2`$ loop algebra symmetry, which holds for the transfer matrices of the XXZ and XYZ spin chains only at the special values of $`\mathrm{\Delta }_1`$. DFM ; FM ; Deguchi  The result could be related to the observation in Ref. 8vertex,  that the eigenvalue spacing distribution $`P(s)`$ has a peak at $`s=0`$ for the transfer matrix of the eight-vertex model under the free-fermion conditions. Details should be discussed in forthcoming papers.

## IV Discussion

Let us discuss some possible reasons why the nointegrable models give the non-Wigner distributions. We may consider two reasons: extra symmetries or finite-size effects.

The non-Wigner distributions, such as shown in Figs. 2 and 3, suggest that the Hamiltonian should have some extra symmetries. It must be a possible standard interpretation according to Refs. Pals,  and Angles, . However, for the NNN coupled XXZ chains, it is not clear whether there indeed exist some extra symmetries other than reflections, translations, and spin rotations around the $`z`$ axis.

We may also consider some finite-size effects, since the distribution $`P(s)`$ is always obtained for some finite systems. All the calculations in this paper are performed on 16-site chains. Recall that the largest sector of the NNN coupled XXZ Hamiltonian has 440 eigenvalues. The number is not small. If the level-spacing distribution calculated on hundreds-site or thousands-site chains could give a Wigner-like distribution, then the finite-size effects should be nontrivial.

###### Acknowledgements.

The present study was partially supported by the Grant-in-Aid for Encouragement of Young Scientists (A): No. 14702012.