text:
# On the Phase Transitions of Graph Coloring and Independent Sets

## 1 Introduction

The class of decision problems that can be solved in nondeterministic polynomial time, known as the NP class, is central to the theory of computational complexity. Informally, a decision problem is in NP if its solution can be checked to be correct in polynomial time. Obtaining that solution, however, is altogether a different matter: for some problems it can be done in polynomial time as well, while for others it is as yet unknown whether a polynomial-time procedure exists. In order to characterize such seemingly harder problems, it has proven useful to look at the problems that are “complete” for NP (the so-called NP-complete problems), that is, problems in NP for which the discovery of a polynomial-time procedure to find a solution would immediately warrant the existence of such a procedure for all the other problems in NP as well .

The NP-complete problems are thus the hardest problems in NP, and are in a sense essentially equivalent to one another in terms of how hard they are to be solved . But it has been known already for many years, both from practical experience (e.g., ) and from looking at the minutiae of the structure of NP , that some NP-complete problems are in fact harder than others, and that a given NP-complete problem may have instances that are significantly harder than other instances of the same problem. Following some initial results of about one decade ago , it is now known that, for NP-complete problems like satisfiability and its derivations  and graph coloring , sharp phase transitions with respect to some order parameter exist and are frequently correlated with the hardness of finding a solution.

We are concerned in this paper with combinatorial optimization problems. Hard optimization problems have been characterized in much the same way as their decision-form counterparts. By a minor technicality, however, they are best termed NP-hard problems (as opposed to NP-complete problems) to indicate only that they are at least as hard as the problems in NP (but not necessarily one of them) . Invariably, an optimization problem whose decision-form variant is NP-complete, is NP-hard.

Unlike decision problems, optimization problems are only now beginning to be looked at in order to explain the relative hardness of their instances, but we already have some empirical evidence of the presence of similar phase transitions . In the case of graph coloring, for example, the optimization problem asks for the graph’s chromatic number—the least number of colors needed to assign one color to each node without ever assigning the same color to neighbors in the graph (i.e., nodes that are connected by an edge).

What is known for this problem comes from considering a sequence of graphs of increasing density (number of edges per node) and what happens along this sequence at the points in which the chromatic number increases. It has been discovered that finding the chromatic number just before these points is distinctly harder than just after them. Also, sharp peaks in the size of the so-called backbone of each graph in the sequence are detected just before those points as well, thus indicating a strong positive correlation between problem hardness and backbone size. A graph’s backbone in this case is the set of node pairs that are assigned the same color by every coloring that employs a number of colors given by the graph’s chromatic number. So a graph with a large backbone presents many opportunities for an algorithm that seeks the optimum to waste time trying to assign two different colors to a node pair that belongs to the backbone.

Coloring a graph optimally is one of the problems that we treat in this paper. While the backbone size relates clearly, in an intuitive way, to why larger backbones tend to imply harder instances of the problem, we feel that it lends little combinatorial insight into the hardness of those instances, specifically into how the structure of the graph affects the hardness of coloring its nodes optimally. Our contribution in the context of this problem, presented in Section 2, is to demonstrate that other indicators exist that relate just as clearly to the appearance of phase transitions related to the hardness of graph coloring while at the same time carrying what we think is better combinatorial intuition. The indicators we use come from considering a graph’s chromatic polynomial and its set of acyclic orientations.

The latter brings us to the second problem of interest in this paper, treated in Section 3, which is the problem of finding an independent set of maximum cardinality in a graph. An independent set is a set of nodes that includes no neighbors. The cardinality of a maximum independent set of a graph is the graph’s independence number. Finding this number is also an NP-hard problem, one that shares with optimal graph coloring a clean combinatorial interpretation in terms of the graph’s acyclic orientations. This problem does not appear to have already been examined for phase transitions related to the hardness of its instances.

Throughout the paper, we use $`G(n,m)`$ (or simply $`G`$, if $`n`$ and $`m`$ can be inferred from the context) to denote an undirected graph with $`n`$ nodes and $`m`$ edges. All our empirical results are based on fixed sequences of graphs, each graph with $`n`$ nodes but increasingly more edges (up to $`N=n(n1)/2`$). For $`1sN`$, the $`s`$th graph in this sequence, denoted by $`G_s(n,s)`$ or simply by $`G_s`$, has $`n`$ nodes, $`s`$ edges, and is generated according to the random-graph model that samples uniformly from the set of all graphs having the same number of nodes and edges  (equivalently, for $`s1`$, $`G_{s+1}`$ may be regarded as being obtained from $`G_s`$ by the random addition of a new edge).

Using one single sequence as the basis of each experiment precludes the smoothing effect of taking averages over larger ensembles, known to mask the appearance of very sharp phase transitions for problems like graph coloring . We are always careful, however, to make sure that the observed phenomena are also present, qualitatively, in several other sequences.

For small values of $`s`$, a graph $`G_s`$ in the sequence $`G_1,\mathrm{},G_N`$ is likely to have isolated nodes (nodes without neighbors). Such nodes do not affect the graph’s chromatic number, and affect its independence number only trivially (every isolated node is a member of all maximum independent sets). So the sequence that is actually used in all our experiments is the sequence $`H_1,\mathrm{},H_N`$, where, for $`1sN`$, $`H_s`$ is obtained by stripping $`G_s`$ of its isolated nodes.

Before proceeding, we pause momentarily to consider this issue of isolated nodes more carefully. Let $`\nu _s`$ denote the number of isolated nodes of $`G_s`$. We have $`\nu _1=n2`$, while for $`s1`$ it is easy to see that

$$\nu _{s+1}=p_0\nu _s+p_1(\nu _s1)+p_2(\nu _s2),$$
(1)
where $`p_k`$, for $`k\{0,1,2\}`$, is the probability that the addition of the $`s+1`$st edge to $`G_s`$ incorporates $`k`$ new nodes into $`H_s`$. These three probabilities can be assessed easily as fractions of $`Ns`$ and lead to a simplification of (1) as

$$\nu _{s+1}=\left(1\frac{n1}{Ns}\right)\nu _s,$$
(2)
which is clearly solved by

$$\nu _s=n\underset{k=0}{\overset{s1}{}}\left(1\frac{n1}{Nk}\right).$$
(3)
Furthermore, for $`sN`$ we obtain

$$\nu _sn\left(1\frac{n1}{N}\right)^sne^{2s/n},$$
(4)
so $`\nu _s`$ obviously decreases rapidly with $`s`$.

We now turn to our two main sections. At the end, concluding remarks are given in Section 4.

## 2 Graph coloring

Let $`\chi (G)`$ denote the chromatic number of $`G`$. We show in Figure 1 two plots for $`n=65`$, one indicating the time needed to find $`\chi (H_s)`$ by a public-domain code  that is based on the heuristic of , the other indicating the evolution of $`\chi (H_s)`$ as $`s`$ is increased. We only show data for a certain range of $`s`$ values, since coloring instances outside this interval tend to be relatively trivial, thus requiring little time for solution. As expected, the chromatic number increases steadily as the graph gets denser (acquires more edges), and does so increasingly rapidly. Also, it is often the case that the time needed to find the chromatic number goes down significantly immediately after each increase in the chromatic number. In this section, we develop new combinatorial arguments that show that such sudden transitions are indeed to be expected.

Our study of the hardness of finding $`\chi (G)`$ starts with an investigation of how abundant optimal colorings of $`G`$ are, that is, colorings that require exactly $`\chi (G)`$ colors. In order to carry out this investigation, we resort to the chromatic polynomial of $`G`$, denoted by $`\pi (G,x)`$, which gives the number of distinct ways in which $`G`$ can be colored by at most $`x0`$ colors.

This polynomial has several interesting properties. For example, it is a degree-$`n`$ polynomial in $`x`$, the coefficient of $`x^n`$ is $`1`$, and the coefficient of $`x^{n1}`$ is $`m`$ . Also, by definition $`\chi (G)`$ is the least value of $`x`$ for which $`\pi (G,x)`$ is positive, giving the number of distinct ways in which $`G`$ can be colored optimally. So finding $`\pi (G,x)`$ is expected to be no easier than finding $`\chi (G)`$, although the following simple method can be used for relatively small graphs.

Let $`e`$ be any edge of $`G`$, denote by $`G\backslash e`$ the graph obtained from $`G`$ by removing $`e`$, and by $`G/e`$ the graph obtained from $`G`$ by contracting the end nodes of $`e`$ into one single node. We see that $`G`$ can be colored in as many distinct ways as $`G\backslash e`$ can, except for those colorings of $`G\backslash e`$ that assign to the end nodes of $`e`$ the same color. But these are precisely the colorings of $`G/e`$, so we get

$$\pi (G,x)=\pi (G\backslash e,x)\pi (G/e,x).$$
(5)
Clearly, (5) defines a simple recursion for calculating $`\pi (G,x)`$ that stops either at graphs that only contain isolated nodes or at graphs that are completely connected. If $`k`$ is the number of nodes in either case, then the former graphs admit $`x^k`$ distinct colorings and the latter $`x!/(xk)!`$ distinct colorings, so the polynomial can be calculated easily at the bases of the recursion and upward from them.

In Figure 2, three plots are given for $`n=10`$: one depicts the continual increase of the chromatic number as the number of edges is increased, while another shows the number of distinct optimal colorings for each graph $`H_s`$, that is, $`\pi (H_s,\chi (H_s))`$ (the third plot is discussed shortly in what follows). For each graph, the chromatic polynomial has been computed using public-domain code  based on the recursion of (5). This computation quickly exhausts processing and memory resources as the numbers of nodes and edges increase—thence the reason why we present data for $`n=10`$ only. So we must always bear in mind that our ability to draw general conclusions may be impaired.

As the plots of Figure 2 indicate, the number of distinct optimal colorings increases dramatically at each step in the chromatic number, and after that decreases more or less steadily until immediately before the next step. In absolute terms, then, at each increase in the chromatic number optimal colorings become strictly more abundant. If the same could be shown to hold in relative terms as well (i.e., if the fraction of optimal colorings relative to some larger universe could be shown to undergo the same transition), then we might be able to continue and investigate in more detail how this relates to the time it takes to color the graphs optimally. However, we find that it is as yet unclear how to characterize such a larger universe.

Even so, it is worth examining the matter further, because at least in part the sudden jumps in the number of optimal colorings are really the product of well understood combinatorial growth at work. For example, let us examine what happens when the chromatic number increases from $`k`$ to $`k+1`$ along the sequence of graphs. Suppose nodes $`i`$ and $`j`$ are the nodes that, when connected to each other by an edge, cause the graph’s chromatic number to increase. Clearly, each optimal coloring before the increase (i.e., with $`k`$ colors) yields at least $`k+1`$ distinct optimal colorings after the increase (i.e., with $`k+1`$ colors): viewing colors as positive integers, first choose one of $`i`$ or $`j`$ and assign to it color $`k+1`$, while the remaining nodes all retain their previous colors; then assign color $`k+1`$ to all nodes that had color $`k`$ and this color to the one of $`i`$ or $`j`$ that was selected previously; then proceed likewise until this same node has been assigned all colors ($`k+1`$ down through $`1`$).

While this may all seem like obvious combinatorial growth at the points where the chromatic number changes, we remark that such an effect may also be present in the behavior of the backbones commonly used to characterize the phase transitions of graph coloring. In the setting that we just examined, nodes $`i`$ and $`j`$ clearly constitute one of the node pairs of the backbone—or else the addition of an edge between them would not cause the chromatic number to increase, because there would be at least one optimal coloring with $`k`$ colors that would assign different colors to them. So let us consider the probability that randomly chosen nodes $`i`$ and $`j`$ constitute a backbone pair.

We do so by first conceding, just for the sake of the argument, that nodes are uniformly distributed among the colors in all optimal colorings. In this case, the probability that $`i`$ and $`j`$ have the same color in all optimal colorings (that is, that $`i`$ and $`j`$ form a backbone pair) when the chromatic number is $`k`$ is $`(k/k^2)^\rho =k^\rho `$, where $`\rho `$ is the number of optimal colorings with $`k`$ colors. But $`\rho `$ increases to at least $`(k+1)\rho `$ when the chromatic number increases to $`k+1`$, so the probability we just computed gets divided by at least $`(k+1)^{k\rho }`$. We then see that the sudden increases observed in the number of optimal colorings are probably also inherently related to what happens to backbones at the same points in the sequence of graphs. This provides a new perspective on the collapse of backbones at the coloring transitions while at the same time providing a better understanding of why it happens.

Another well understood combinatorial-growth effect is that, when the chromatic number is $`k`$, every optimal coloring is essentially equivalent to $`k!1`$ others, each corresponding to a permutation of the colors among the nodes. This brings us to the third plot of Figure 2, where the number of optimal colorings is shown normalized by the factorial of the current chromatic number. Evidently, all sudden jumps are still there, but they now possess a stronger significance, because only one optimal coloring is counted out of all colorings that are equivalent to one another by straightforward permutation of colors (that is, without implying a different partition of the node set).

But let us return to the chromatic polynomial of $`G`$. The usefulness of this polynomial goes beyond the counting of distinct colorings, as it provides the first link to yet another characterization of the hardness of optimal graph coloring, now based on the acyclic orientations of $`G`$. An orientation of $`G`$ is an assignment of directions to the edges of $`G`$; it is acyclic if no directed cycles are formed, that is, if it is impossible to reach the same node twice by following edges according to their directions exclusively.

If we let the set of the acyclic orientations of $`G`$ be denoted by $`\mathrm{\Omega }(G)`$ and consider the number of acyclic orientations of $`G`$ in terms of what happens to the graphs $`G\backslash e`$ and $`G/e`$ introduced earlier, then we have the following. Every acyclic orientation of $`G\backslash e`$ yields either one or two acyclic orientations for $`G`$ when $`e`$ is assigned a direction. The former case happens when one direction assignment for $`e`$ forms a directed cycle in $`G`$ but not the other, the latter when both assignments preserve acyclicity in $`G`$ (it cannot happen that both assignments form directed cycles because the orientation of $`G\backslash e`$ is acyclic). Similarly, every acyclic orientation of $`G/e`$ is necessarily one of those acyclic orientations of $`G\backslash e`$ from which two orientations of $`G`$ are obtained. Thus,

$$|\mathrm{\Omega }(G)|=|\mathrm{\Omega }(G\backslash e)|+|\mathrm{\Omega }(G/e)|$$
(6)
for any edge $`e`$ of $`G`$.

The recursion in (6) is strikingly similar to the one in (5), and this has been shown to give rise to the remarkable identity

$$|\mathrm{\Omega }(G)|=(1)^n\pi (G,1).$$
(7)
That is, the number of acyclic orientations of $`G`$ can be obtained from applying the chromatic polynomial of $`G`$ to the negative unit . Obtaining (7) and some of its refinements  from the relation between (5) and (6) comes as a consequence of the so-called theory of P-partitions and its order polynomials . Of interest to us is that such polynomials quantify the following very useful relationship between the acyclic orientations of $`G`$ and its colorings.

Suppose that $`G`$ can be colored by $`k`$ colors. Still viewing colors as positive integers, suppose also that we assign to the edges of $`G`$ an orientation that makes every edge point from the node with the higher color to the one with the lower. This orientation is clearly acyclic and induces no directed path in $`G`$ containing more than $`k`$ nodes. Conversely, suppose we start from an acyclic orientation of $`G`$. If $`k`$ is the number of nodes on the longest directed path in $`G`$ according to this orientation, then $`G`$ can be colored by at most $`k`$ colors, as follows. We assign color $`1`$ to the sinks (nodes whose adjacent edges are all oriented inward), then color $`2`$ to the sinks that would be formed if the original sinks were to be removed from $`G`$, then the lowest available color to the set of sinks that would appear next, and so on.

Note, in this process, that starting with a coloring yields a unique acyclic orientation. The converse is not necessarily true, however: starting with an acyclic orientation may yield more than one coloring of the nodes of $`G`$. Consider, for example, the acyclic orientation shown in Figure 3 and the two corresponding colorings shown in parentheses next to the nodes. What the process indicates, however, is that it is possible to seek optimal colorings for $`G`$ by looking for acyclic orientations of $`G`$ that are shortest in terms of how many nodes there are in a longest directed path. Letting $`\mathrm{\Omega }^{}(G)\mathrm{\Omega }(G)`$ be the set of such orientations, what we have seen is that

$$|\mathrm{\Omega }^{}(G)|\pi (G,\chi (G)),$$
(8)
so we may have a tighter characterization of how abundant optimal colorings are by looking at shortest acyclic orientations instead of the chromatic polynomial applied to $`\chi (G)`$.

The formal relationship between $`\chi (G)`$ and the acyclic orientations of $`G`$ can be stated as follows , and strengthens important earlier results . For $`\omega \mathrm{\Omega }(G)`$, let $`P_\omega `$ be the set of all directed paths in $`G`$ according to $`\omega `$. For $`pP_\omega `$, let $`|p|`$ indicate the number of nodes in $`p`$. Then

$$\chi (G)=\underset{\omega \mathrm{\Omega }(G)}{\mathrm{min}}\underset{pP_\omega }{\mathrm{max}}|p|.$$
(9)
An illustration is given in Figure 4, where two acyclic orientations are shown for the same graph, together with the corresponding partition into sinks alluded to earlier. This partition is known as the sink decomposition of the graph according to the acyclic orientation . In the figure, each such decomposition is shown with a rightmost box containing the sinks, then another box to its left containing the sinks that appear if the former sinks are eliminated, and so on. The set of nodes in each box is normally referred to as a layer of the sink decomposition. Notice that the number of layers in the sink decomposition for $`\omega `$ is precisely $`\mathrm{max}_{pP_\omega }|p|`$. In the case of Figure 4, the bottommost acyclic orientation has the smaller sink decomposition and thus corresponds to a better (in this case, optimal) assignment of colors (a different color to the nodes in each layer).

For $`n=10`$, Figure 5 depicts the relationship between the chromatic number and the number of acyclic orientations of each graph in $`H_1,\mathrm{},H_N`$. The thinner dashed plot in the figure gives the number of acyclic orientations of each graph. This number has been computed using the algorithm of , which although efficient in several aspects relevant to the analysis of enumerative algorithms, becomes prohibitive very quickly as the graph gets larger. The same enumeration process has been used to record the number of acyclic orientations that are shortest, that is, those whose sink decompositions have as many layers as the graph’s chromatic number. This small addition to the algorithm employs straightforward depth-first search , and the results are shown as the thicker solid plot of the figure.

Remarkably, an effect very similar to the one observed in Figure 2 is seen to occur now as well: the number of shortest acyclic orientations increases sharply whenever the chromatic number increases, and subsequently goes down until immediately before the next increase. These are also absolute data, but now it is obvious how to make them relative: we simply observe the percentage of all acyclic orientations that are shortest. This is shown in Figure 6, which confirms that, also in relative terms, shortest acyclic orientations become significantly more abundant right after an increase in the chromatic number, becoming increasingly rarer from there onward until the next increase.

This indicates that, in a sense, coloring a graph optimally becomes easier immediately after an increase in the chromatic number, and incrementally harder until the next transition is reached. Behind this statement is the intuitive feeling that coloring algorithms should in general fare better when optimal solutions are more abundant, be such abundance assessed as some indicator of how many optimal colorings there are or as the fraction $`|\mathrm{\Omega }^{}(G)|/|\mathrm{\Omega }(G)|`$ of shortest acyclic orientations. Of course, this intuition calls for experimental support, but notice that at least the most naïve of all random approaches—a series of Bernoulli trials inside $`\mathrm{\Omega }(G)`$ until the first member of $`\mathrm{\Omega }^{}(G)`$ is found—is certain to benefit from such abundance of shortest acyclic orientations, as clearly the expected time for its convergence is $`|\mathrm{\Omega }(G)|/|\mathrm{\Omega }^{}(G)|`$ . More serious methods with the potential to benefit from the relative abundance of optimal acyclic orientations exist , though, and a systematic effort to assess their capabilities on sequences of random graphs is under way.

The number (or fraction) of optimal acyclic orientations is, in principle, also subject to the same concerns with the disguising action of obvious combinatorial growth that we expressed earlier. To see what happens when we apply the same normalization by the factorial of the current chromatic number, we have in each of Figures 5 and 6 a plot with the results (the one in thick dashes). The pronounced jumps are still present, but this normalization is only meaningful as inherited from the relationship between optimal colorings and shortest acyclic orientations.

A better normalization may exist and we would like to digress on this possibility briefly, although several problems related to it are still open. Suppose we take an acyclic orientation and turn all its sinks into sources (nodes whose adjacent edges are all oriented outward). This necessarily yields another acyclic orientation, and the continual repetition of the process must eventually lead to a period of orientations. This attractor dynamics has several interesting properties ; in our context, the most crucial property is that the number of layers in the graph’s sink decomposition is continually nonincreasing along the process of obtaining new acyclic orientations by turning sinks into sources. So the period at the core of each attractor comprises orientations that yield sink decompositions all with the same number of layers, this number being also no larger than that resulting from any other orientation in the same attractor. This means that, in addition to the chromatic indicator originally observed when this attractor dynamics was first analyzed (the graph’s interleaved multichromatic number, cf. ), each attractor is related to finding the chromatic number as well, since finding a period whose orientations are shortest over all attractors immediately yields the chromatic number. A period, in summary, provides a means of expressing the equivalence of several acyclic orientations and may become suitable for the normalization we need if only more knowledge can be obtained on it. Unfortunately, we thus far lack this necessary additional knowledge.

To finalize, we comment on yet another interesting insight that can be gained from considering the relationship between a graph’s colorings and its acyclic orientations. Suppose that $`\omega `$ is an acyclic orientation of $`G`$, and consider the random addition of an edge to $`G`$ between two nodes not currently connected. If we look at $`\omega `$ from the perspective of its sink decomposition, then let its layers be numbered $`1`$ through $`L`$, starting at the layer that contains the sinks and onward. As illustrated in Figure 7, in some cases the addition of the new edge will preserve the sink decomposition, while in others it will not. These two possibilities are shown in the leftmost sink decomposition in the figure as two dashed directed edges. The addition of one of them preserves the sink decomposition (shown in the middle sink decomposition of the figure); the addition of the other, which connects two nodes in the same layer, forces the sink decomposition to acquire another layer, as shown in the rightmost sink decomposition of the figure.

Let us then assume that, if the new edge is added between layers, then it is oriented in the direction of the lower-numbered layer. Let also the $`\mathrm{}`$th layer have $`s_{\mathrm{}}`$ nodes, so that $`n=_{\mathrm{}=1}^Ls_{\mathrm{}}`$. The number of node pairs not yet connected by an edge is $`Nm`$; of these, the ones that have nodes in different layers amount to $`_{\mathrm{}=1}^{L1}_{\mathrm{}^{}=\mathrm{}+1}^Ls_{\mathrm{}}s_{\mathrm{}^{}}m`$. If we let $`q`$ be the probability that adding the new edge does not disturb the sink decomposition, then $`q`$ is the probability that the edge is added between layers. That is,

$$q=\frac{_{\mathrm{}=1}^{L1}_{\mathrm{}^{}=\mathrm{}+1}^Ls_{\mathrm{}}s_{\mathrm{}^{}}m}{Nm}.$$
(10)
By (10), the addition of an edge between layers (this happens with probability $`q`$) causes $`q`$ to decrease, as the only change in the formula is the concomitant subtraction of $`1`$ off both the numerator and the denominator. When the edge is added between nodes of the same layer (with probability $`1q`$), the double summation in the numerator of (10) may either remain the same or vary. If it remains the same or decreases, then $`q`$, as before, decreases. In order to verify what happens otherwise, let $`S`$ denote the double summation. Then $`q`$ is seen to vary by at least

$$\frac{S+1m1}{Nm1}\frac{Sm}{Nm}=\frac{Sm}{(Nm)^2(Nm)}0,$$
(11)
so $`q`$ increases unless $`S=m`$, in which case it was zero to begin with and remains zero.

If we now consider the number of edges that need to be randomly added to $`G`$ so that $`q`$ once again assumes the value it currently has, say $`ϵ`$, it is easy to see from (10) that this number is given by $`\delta /(1ϵ)`$, where $`\delta `$ is the increase that $`q`$ incurs along the way. The value of $`\delta `$ is very hard to quantify, but it seems reasonable to assume, at least for the sake of the argument, that it gets smaller as the graph gets denser (i.e., acquires more edges). In this case, the number of random edge additions needed for $`q`$ to return to the value $`ϵ`$ is ever smaller as the process unfolds. Overall, what we witness is a process that resembles the increase in $`\chi (G)`$ as $`G`$ becomes denser, slow at first but increasingly rapid as the graph’s density gets higher. We hope to obtain a better characterization of this resemblance as further research adds detail to the picture. If we succeed, it may be possible to obtain a generic prescription for determining all the points at which the chromatic number increases along the sequence, thus adding to what is already known .

## 3 Independent sets

Let $`\alpha (G)`$ denote the independence number of $`G`$. We have used public-domain code  based on  to find $`\alpha (H_s)`$ for each graph in $`H_1,\mathrm{},H_N`$ with $`n=75`$. What this code finds is not a maximum independent set directly, but rather a maximum clique in the graph that is complementary to the graph of interest (i.e., has an edge joining two distinct nodes if and only if the graph of interest does not). A clique is a subgraph whose nodes are all connected to one another, so the correspondence to independent sets should be clear. This is illustrated in Figure 8.

Figure 9 has two plots for $`n=75`$, one to indicate the time to find the graph’s independence number, the other to indicate how this number evolves as the graph becomes denser. We show data for a certain density interval only. To the left of what is shown the graphs have isolated nodes and the behavior is uncharacteristic, while to the right times become too small to be indicative of any particularly interesting behavior. But inside the density interval used in the figure the independence number goes down nearly steadily,<sup>1</sup><sup>1</sup>1The two increases in the independence number back to $`30`$ from $`29`$ for $`s=202`$ and $`s=225`$ can only be explained by the existence of two last isolated nodes in $`G_{201}`$, one of which gets incorporated into $`H_{202}`$, the other into $`H_{225}`$. That a node should remain isolated through $`s=224`$ for $`n=75`$ is unlikely but entirely conceivable. In fact, the probability of such an event, given approximately by $`e^{2s/n}`$ (cf. (4)), is $`0.0025`$. and at some of these downward transitions there appears to be a considerable decrease in the time for optimal solution right after the decrease in the independence number. Although this evidence is less ubiquitous than in the case of graph coloring, and for this reason certainly less compelling, this study seems to be the first one to address the issue of phase transitions related to finding a graph’s independence number and for this reason we investigate the matter further. In this section, we introduce some combinatorial arguments that appear to relate to these phenomena.

As for graph coloring, our initial step is to assess the abundance of maximum independent sets for a given graph $`G`$, that is, independent sets of size $`\alpha (G)`$. Unlike the case of graph coloring, however, it is as yet unknown how to count this number exactly. This remains true even if we settle for maximal (as opposed to maximum) independent sets, that is, independent sets that cannot be enlarged without losing the independence property, although in this case increasingly better upper bounds on the desired number have been discovered recently (cf.  and the references therein).

But at this point it helps to recall that the sequence of graphs $`G_1,\mathrm{},G_N`$ is randomly generated, each graph being drawn uniformly from the set of graphs having the same number of nodes and edges. We may therefore attempt to assess the number of maximal independent sets of a given size in each graph by computing the expected value of this number, which is in fact a random quantity. In order to accomplish this more easily, it is helpful to resort to the model of random graphs in which an edge exists between two nodes with constant probability, say $`p`$, independently of the nodes. Although this is not the model under which our graphs were generated, away from limiting situations it is safe to assume that the two models are equivalent to each other with $`p=m/N`$ for $`G(n,m)`$ .

Let $`\xi (G,k)`$ denote the expected number of maximal independent sets of size $`k`$ in $`G`$. The number of candidate sets is $`\left(\genfrac{}{}{0pt}{}{n}{k}\right)`$, and the probability that each one is an independent set is $`(1p)^{k(k1)/2}`$. If a candidate is an independent set, then the probability that it is maximal is the probability that each of the remaining $`nk`$ nodes is connected to at least one of its $`k`$ nodes, that is, $`[1(1p)^k]^{nk}`$. We then get

$$\xi (G,k)=\left(\genfrac{}{}{0pt}{}{n}{k}\right)\left(1\frac{m}{N}\right)^{\frac{k(k1)}{2}}\left[1\left(1\frac{m}{N}\right)^k\right]^{nk}.$$
(12)
Now let $`G`$ be a random graph generated in compliance with the edge density given by $`p`$, and let the value of $`\alpha (G)`$ be known. Letting $`k=\alpha (G)`$ in (12) yields the expected number of maximal independent sets in similar random graphs, whose size is, at least for one of them, the size of its maximum independent set.

For $`n=75`$, we show the evolution of this number for the sequence $`H_1,\mathrm{},H_N`$ in Figure 10, along with the evolution of the independence number. Interestingly, every transition of the independence number to a smaller value causes the expected number of maximal independent sets of size $`\alpha (H_s)`$ to increase markedly. From there onward, this number decreases until the next similar transition occurs. This is one first indication, albeit imprecise, that the abundance of such maximal sets may be related to the increased ease with which the graph’s independence number can sometimes be found immediately after a decrease in the independence number. However, as in our analysis of the chromatic polynomial in Section 2, it is not clear how to proceed and characterize the relationship more effectively. Not only this, but it is still possible, as in the case of graph coloring, that some underlying inherent equivalence among maximal independent sets exists that would indicate a way of normalizing the data shown in the figure. We still do not know how that can be achieved.

Once again, though, a deep relationship exists between a graph’s independent sets and its acyclic orientations. It is in this case both more complex and more subtle than in the case of graph coloring, so we describe it with the aid of an illustration right from the start. First consider Figure 11, where two acyclic orientations of the same graph are shown, each one alongside what is known as a chain decomposition of the graph according to it. A chain decomposition of a graph according to an acyclic orientation is a partition of the graph’s node set such that the nodes in each set of the partition are arranged by the acyclic orientation as a single chain of nodes (a directed path). In the case of Figure 11, each partition is displayed with the aid of boxes to enclose the nodes that go in each set. Furthermore, each of the chain decompositions shown in the figure employs the minimum number of chains for the corresponding acyclic orientation. The topmost acyclic orientation admits a chain decomposition comprising one single chain, while for the other no chain decomposition exists with fewer than two chains.

Understanding how chain decompositions and independent sets are related involves several technicalities that we will not discuss formally but illustrate through examples instead. Figure 12 contains two flow networks, that is, two directed graphs with two distinguished nodes, $`s`$ and $`t`$, on which we consider the maximum flow from $`s`$ to $`t`$ subject to certain edge capacities (unit for all edges that emanate from $`s`$ or arrive at $`t`$, infinite for all others). Each flow network corresponds to one of the acyclic orientations of Figure 11 and is constructed from that orientation as follows. For each node $`i`$ in $`G`$, two nodes, $`i^{}`$ and $`i^{\prime \prime }`$, are added to the flow network, along with a directed edge from $`s`$ to $`i^{}`$ and one from $`i^{\prime \prime }`$ to $`t`$. If an edge exists between nodes $`i`$ and $`j`$ and the acyclic orientation directs that edge from $`i`$ to $`j`$, then the flow network contains an edge directed from $`i^{}`$ to $`j^{\prime \prime }`$.

The maximum flow from $`s`$ to $`t`$ in each of the networks of Figure 12 is indicated by solid edges (edges that carry unit flow) and dashed edges (edges that carry zero flow). In general, the following interesting property holds for maximum flows in such networks . Of the solid edges, those whose removal disconnects $`s`$ from $`t`$ minimally (the minimum cut) induce a node cover in $`G`$, that is, a set of nodes that touches every edge. In Figure 12, an edge in the minimum cut either leads from $`s`$ to a node enclosed in a box or from such a node to $`t`$. Nodes thus represented correspond to the node sets $`\{1,2,3,4\}`$ and $`\{2,3,5\}`$ in $`G`$, each set a node cover. Note that, by definition, the complements of these sets with respect to the node set of $`G`$ are necessarily independent sets of $`G`$—in the figure, these are the sets $`\{5\}`$ and $`\{1,4\}`$, each of whose nodes corresponds to one of the chains in Figure 11.

Formally, the relationship between $`\alpha (G)`$ and the acyclic orientations of $`G`$ is as given next , and brings sharper focus to the classic result established by Dilworth’s theorem . For $`\omega \mathrm{\Omega }(G)`$, let $`D_\omega `$ be the set of all chain decompositions of $`G`$ according to $`\omega `$. For $`dD_\omega `$, let $`|d|`$ be the number of chains in $`d`$. Then

$$\alpha (G)=\underset{\omega \mathrm{\Omega }(G)}{\mathrm{max}}\underset{dD_\omega }{\mathrm{min}}|d|.$$
(13)
As one readily notices, (13) is a sort of dual of (9): it expresses $`\alpha (G)`$ as the number of chains in the minimum chain decomposition that has the most chains over $`\mathrm{\Omega }(G)`$. Henceforth, we let $`\mathrm{\Omega }^+(G)\mathrm{\Omega }(G)`$ denote the set of such widest acyclic orientations, understood as those acyclic orientations that achieve the minimum chain decomposition that is maximum over $`\mathrm{\Omega }(G)`$.

It is now instructive to return briefly to Figure 11. For each acyclic orientation $`\omega `$ displayed in the figure, the number of chains in the corresponding chain decomposition is $`\mathrm{min}_{dD_\omega }|d|`$. Of the two acyclic orientations, the bottommost achieves the maximum of this quantity over $`\mathrm{\Omega }(G)`$.

We now investigate the use of $`|\mathrm{\Omega }^+(G)|`$ as an indicator of how abundant the widest acyclic orientations of $`G`$ are and how this relates to the hardness of finding $`\alpha (G)`$. This time, our procedure to enumerate all the acyclic orientations of a graph has been coupled with a public-domain code  that implements the algorithm of  to find the maximum flow in a network (and also the minimum cut, as a by-product).

The results are shown in Figure 13, in which the number of all acyclic orientations of each graph is plotted alongside its independence number and the number of acyclic orientations that are widest (those whose chain decomposition into the fewest possible chains requires $`\alpha (H_s)`$ chains). The figure indicates that the same sharp transitions that appear in Figure 10 are also observed for the acyclic orientations: every downward transition that the independence number undergoes is accompanied by a sharp increase in the number of widest acyclic orientations.

When we assess such increases with respect to the set of all the acyclic orientations of each graph, we obtain what is shown in Figure 14. Clearly, every decrease in the independence number corresponds to a marked increase in the fraction of acyclic orientations that are widest. We hope to eventually be able to conclude that such increases are correlated with the time it takes to find a graph’s maximum independent set. As we mentioned in Section 2, obviously the unreasonable approach of random trials does benefit from an elevated $`|\mathrm{\Omega }^+(G)|/|\mathrm{\Omega }(G)|`$ ratio, but we believe this may also be the case for heuristics that exploit the role of acyclic orientations directly .

## 4 Conclusions

We have investigated the NP-hard problems of coloring the nodes of a graph optimally and of finding a maximum independent set in a graph. For each of these problems, we first displayed empirical evidence that, as the number of edges in the graph is increased by the random addition of one edge at a time, significant variations take place in the time to solve the problem optimally. Often these variations coincide with the transition of the value of the optimal solution to a new level, a higher one for graph coloring, a lower one for independent sets.

For graph coloring, we have shown that the upward transitions in the chromatic number along the sequence of increasingly denser graphs coincide with sharp increases in the abundance of distinct colorings employing the optimal number of colors. More importantly, we have shown that the same phenomenon takes place when we consider the shortest acyclic orientations of the graphs, that is, those orientations whose sink decompositions have as many layers as the graphs’ chromatic numbers. In this case, sharp increases are observed in the ratio $`|\mathrm{\Omega }^{}(H_s)|/|\mathrm{\Omega }(H_s)|`$.

Our conclusions for the maximum independent set problem are similar. The downward transitions in the independence number that occur along the sequence of graphs coincide with sharp increases in the expected number of maximal independent sets of certain sizes. As for the acyclic orientations of the graphs, we have given evidence that the same type of phenomenon takes place regarding the widest acyclic orientations of the graphs, that is, those whose minimum chain decompositions have as many chains as the graphs’ independence numbers. In this case, sharp increases take place in the ratio $`|\mathrm{\Omega }^+(H_s)|/|\mathrm{\Omega }(H_s)|`$.

We find that the two special subsets of $`\mathrm{\Omega }(G)`$, $`\mathrm{\Omega }^{}(G)`$ and $`\mathrm{\Omega }^+(G)`$, underlie an interpretation of the phase transitions observed in optimal graph coloring and optimal independent sets that is full of combinatorial meaning. This meaning is directly related to the structure of the graphs involved and that of its acyclic orientations. In addition, it provides a direct interpretation of the hardness of the problems as given by the relative abundance of the acyclic orientations that give the optima.

Except for those shown in Figures 1, 9, and 10, all our empirical results have been given for $`n=10`$ only. As we indicated earlier, the reason for this has been the severe combinatorial explosion that occurs when finding a graph’s chromatic polynomial or its set of acyclic orientations. For $`n=15,20`$, so far these computations could only be carried out through about $`40`$ edges, yielding results that fully support the conclusions we have drawn along the paper. But, instead of presenting such partial results, we opted for giving the reader data on the full evolution through the densest graphs for $`n=10`$.

We think the empirical evidence we have provided is only the beginning of a deeper investigation of how the abundance of certain acyclic orientations affects the hardness of optimal graph coloring and of finding maximum independent sets. Given the aforementioned computational difficulties, the most promising avenue for continuation seems to be the search for additional analytic properties that can explain the findings we have described. Perhaps the evolution of the probability in (10) superficially described at the end of Section 2, and of an analogous indicator in the case of independent sets, should be investigated more deeply as a first step.

### Acknowledgments

The authors acknowledge partial support from CNPq, CAPES, the PRONEX initiative of Brazil’s MCT under contract 41.96.0857.00, and a FAPERJ BBP grant.


format: text

text:
# The X-ray spectra of the flaring and quiescent states of AT Microscopii observed by XMM-Newton Based on observations obtained with XMM-Newton, an ESA scienc mission with instruments and contributions directly funded by ESA Member States and the USA (NASA)

## 1 Introduction

Hot outer atmospheres (coronae) are very common for relatively cool stars in the spectral classes F-M. Many of these coronae are characterized by temperatures up to about 20 MK and densities $`n_\mathrm{e}>10^{10}`$ cm-3. It was already known that the coronae of many stars are different from the solar corona (e.g., reviews by Pallavicini 1989 and Mewe 1991). The heating mechanism of these coronae, however, is still not well understood. A variety of possible explanations is present, such as, convection zones with MHD waves transfering energy from the photosphere into the corona, large flares heating plasma by reconnection, or a large number of continuously heating micro- and nanoflares. High-resolution spectroscopy of stellar systems available from Chandra (Brinkman et al. 2000, Canizares et al. 2000) and XMM-Newton (Brinkman et al. 2001) offers the possibility to study the coronal spectra in great detail and to determine various coronal quantities, such as densities, temperatures, abundances, emission measures, and line ratios, during the flaring states as well as during the quiescent state.

Late-type stars (M and K dwarfs) show more magnetic activity and higher coronal temperatures than the Sun, characterized by flaring (Pettersen 1989) in the optical, UV, X-ray and radio wavelengths (Pallavicini et al. 1990). To understand the mechanisms that underly the high temperatures and the variety of densities, emission measures, and abundance values observed in coronal systems, the X-ray spectra of stars of type F-M are investigated, especially focused on stars with flare activity. In very active stars, however, the “quiescent state” might be produced by a superposition of flare decays and microflares rather than by a real quiescent coronal plasma. Here we present the investigation of the active M-type binary star AT Mic (dM4.5+dM4.5) at a distance of 10.22 pc (Perryman et al 1997). Both components are known to flare frequently (Joy & Wilson, 1949). Large optical flares were observed by García-Alvarez et al. (2002), while X-ray flares were studied e.g., by HEAO-1 (Kahn et al. 1979) and by EXOSAT (Pallavicini et al. 1990). The quiescent corona was observed by EUVE (Monsignori Fossi et al. 1995).

## 2 Observations

The X-ray spectrum of AT Mic was observed by rgs and epic-mos on board XMM-Newton on 14 October 2000 during revolution 156. The observation log of the data is presented in Table 1.

The data were processed by means of the XMM-Newton SAS version 5.3.3 (June 2002). The exposure time was 28.3 ks for rgs and 25.4 ks for epic-mos1, of which 19.0 ks were free of solar flare proton “pollution” applying the criterion that the latter’s count rate should be $`<`$ 0.35 c/s. For rgs the first-order net spectra were extracted by including 95% of the cross dispersion PSF (xpsfincl = 95 in rgsproc) and the background spectra were extracted by excluding 98% (xpsfexcl = 98 in rgsproc).

The rgs spectral resolution is $`\mathrm{\Delta }\lambda `$0.07 Å (FWHM), with a maximum effective area of about 140 cm2 around 15 Å. The wavelength uncertainty is 7-8 mÅ. The total rgs spectrum runs from 5 to 37 Å. Only the range from 8 to 37 Å was used. For epic-mos1, which observed in the small-window mode, we extracted the spectrum by means of a circle centered around the source with a radius of 40 arcsec. The background was obtained by applying the same circle in a source free area at another CCD. The same SAS version was used to construct response matrix and ancillary response files. The range from 1 to 14 Å was used. epic-mos2 was operating in the timing mode. For more instrumental details on rgs and epic-mos see den Herder et al. (2001) and Turner et al. (2001), respectively. The pn data were used to produce a lightcurve only.

## 3 Analysis

### 3.1 Quiescence versus flaring state

Based on the lightcurves obtained with RGS, EPIC-MOS and pn (see Fig. 1) we have concluded that AT Mic was flaring several times during our observations. A period of about 2.9 ks related to the most prominent flare was extracted from the GTI table as flare state, leaving a pseudo-quiescent part of 14.1 ks. Three other weak flares can be seen in Fig. 1. However, they do not dominate and have therefore not been flagged as flares in the further spectral analysis. From Fig. 1 we notice that the top of the flare is not constant, but shows small variations, which seem to be periodic. However, these time intervals are too short to allow for individual spectral analysis.

The observed spectra are shown in Fig. 2 for the quiescent (top) and flaring (bottom) state. As can be seen from Fig. 2, epic-mos has a low resolution, but a high sensitivity compared to rgs. Due to the short exposure time of the flaring state the signal/noise ratio is lower for that spectrum. For both states shown in Fig. 2 it is clear that the X-ray spectrum obtained by rgs in the wavelength range from 5 to 37 Å is dominated by H- and He-like transitions of C, N, O, and Ne and by Fe xvii lines, while the low-resolution epic-mos spectrum shows the Fe K-shell transitions around 1.9 Å and the H- and He-like transitions of Ne and Si. The Fe K-shell feature (see also Table 4) as well as the continuum appear more prominent in the spectrum of the flaring state (bottom panel) than in that of the quiescent state.

#### 3.1.1 Multi-temperature fitting

The data of the quiescent state and the flaring state have been fitted with a 3-$`T`$ Collisional Ionization Equilibrium model using the spex code (Kaastra et al., 1996a) in combination with an updated version of mekal (Mewe et al. 1995, Phillips et al. 1999). The mekal database is publicly available111http://www.sron.nl/divisions/hea/spex/version1.10/line/line\_new.ps.gz  as an extended list of fluxes of more than 5400 spectral lines. Fig. 3 shows the comparison of the quiescent spectrum with the best-fit model. The fits have been performed by attaching the errors to the data (“fit weight data”), given in the columns “data”, and by attaching the errors to the model (“fit weight model”), given in the columns “model”. In the abstract the values from the “data fit” are given.

Table 2 shows the 3 temperatures, the emission measures in 1051 cm-3, and abundances together with their standard deviations within parentheses. The abundances are normalized to their corresponding solar photospheric values (Anders & Grevesse 1989), except for iron and oxygen. For Fe we have changed the value of log AFe from 7.67 (Anders & Grevesse 1989) to 7.51 (Grevesse & Sauval 1998 and 1999). For oxygen we used the new value given by Allende Prieto & Lambert (2001) (log AO = 8.69 instead of 8.93). Here log AFe and log AO is the logarithm of the Fe- and O-abundance relative to log AH=12.0.

The temperatures in Table 2 cover a range from 3 to 35 MK. It turns out that the line spectrum is mainly produced by the EM related to $`T_1`$ = 3 MK and to $`T_2`$ = 8 MK, while the hot component $`T_3>`$20 MK contributes to the continuum below 10 Å and the Fe K-shell lines in the epic-mos observation. This is especially true for the flaring state. A fourth cool component (1 MK) is possible but not well constrained. This cool component is only related to the C vi line (influencing the C abundance) and very weakly coupled to the continuum. The total emission measure of the three components is 12.2(.5)$`\times 10^{51}`$ cm-3 for the quiescent state and 19.5(.8)$`\times 10^{51}`$ cm-3 for the flare state. Based on EUVE observations, Monsignori Fossi et al. (1995) have determined an emission measure maximum of about 3$`\times 10^{51}`$ cm-3 around 10 MK (their Figure 5), detecting only the emission measure related to our second temperature component at $`T_2`$(7.9 MK).

When fitting the O VII lines a slight increase of the electron density $`n_e`$ is noticed, from the quiescent state to the flaring state. However, this increase is not statistically significant.

In the energy interval from 0.3 to 10 keV the X-ray luminoscity $`L_x`$ is 2.13(.08)$`\times 10^{29}`$erg/s for the quiescent state and 3.63(.22)$`\times 10^{29}`$erg/s for the flaring state. These values are averages of the “data” and “model” fit. Pallavicini et al. (1990) observed AT Mic with EXOSAT (ME 1-10 keV detector) during quiescence and during a flare. The quiescent X-ray luminosity was measured to be 1.8 1029 erg/s (Table 5 of Pallavicini et al. 1990) in agreement with our value. However, they observed a larger flare, but with the same temperature as ours.

The abundances have been determined by fitting to the spectrum as a whole. As a check, whether a large scatter exists between different lines of the same element, abundances were also fitted to isolated individual lines, keeping temperatures, emission measures and other abundances fixed at the values from the global fit (Table 2). In this way blends have been taken into account. Due to the weak continuum the absolute abundances cannot be coupled accurately to the continuum and therefore they are strongly anticorrelated with the emission measure. This means that if one quantity increases the other decreases and vice versa. This is clear from comparing the $`EM_{\mathrm{total}}`$ from the “data fit” with that from the “model fit” and by comparing the absolute Oxygen abundance (O/H) in these columns. The product $`EM_{\mathrm{total}}`$$`\times `$ O/H, however, is more stable. For this reason the abundance ratios, normalized to oxygen, are more robust. Oxygen was chosen for its low relative standard deviation, thanks to the strong O viii and O vii lines. For the (less active) solar corona Feldman et al. (1992) discussed the presence of a First Ionization Potential (FIP) effect. This effect implies that elements with a low value of the first ionization potential (e.g., $`<`$ 10 eV) show enhanced coronal abundances relative to photospheric values. However, the abundance ratios of the low-FIP elements Mg, Si, and Fe are relatively low in AT Mic, indicating an inverse FIP-effect (IFIP). Similar IFIP effects were observed by Brinkman et al. (2001), Güdel et al. (2001), Audard et al. (2002) and Audard et al. (2003) for other active late-type stars. This IFIP effect is illustrated in Fig. 4 (top), which shows the trend of an increasing abundance ratio with increasing FIP for the quiescent state. In the bottom panel the same quantities are shown for the flaring state. However, due to the short time interval all quantities determined in the flare spectrum have large statistical uncertainties. The poor signal-to-noise ratio of the flaring state data is also reflected in the “too good” $`\chi ^2`$/dof of the multitemperature fit. This suggests an overestimate for the errors in this low count rate spectrum. For that reason the fits have also been performed, attaching the errors to the model, resulting in reasonable $`\chi ^2`$/dof-values.

Comparing the quiescent state with the flaring state some flattening of this IFIP effect for the low FIP elements (Fe and Mg) is noticed in the flaring state. This is shown in Table 3 and Fig. 4c, in which the ratios between the A/O values of the flaring and quiescent state are given. The values used in Table 3 and Fig. 4 are averages between “data” and “model” values. Si seems not to be affected by the difference between the flaring state and the quiescent state. In the solar corona the enhancement of the low-FIP elements starts from 10 eV down. Results in Table 3 and Fig. 4c may indicate a value lower than 10 eV below which the elemental abundances are influenced. Such a shift might be related to the lower surface temperature of AT Mic.

### 3.2 DEM-modeling

To show the smooth connection between the three separated temperature components of the multi-temperature fit (Table 2) we have performed a DEM-modeling of the rgs and epic-mos spectra of AT Mic in the quiescent and flaring state, applying the regularization method and the polynomial method present in the spex-code (see Kaastra et al. (1996b)). In this DEM-modeling the abundance ratios given in Table 2 have been used in combination with an O/H value of 1.57 for the quiescent state and 1.39 for the flare state. The emission measures $`EM=n_\mathrm{e}n_\mathrm{H}V`$ (per logarithmic temperature bin) of the corona of AT Mic are shown in Fig. 5ab for both states. The emission measure is distributed from $``$1 up to about 60 MK for the flare state. The DEM distribution agrees with the three temperature components from the multi-temperature fit if we take into account the uncertainties corresponding to the standard deviations.

Table 2 and Fig. 5 show the increase of the emission measure and a shift to higher temperatures for the hot EM (40 MK) during the flare state (bold) compared to the quiescent state (weak). The emission measure below 10 MK is the same for the quiescent and flaring state.

### 3.3 Individual lines

The more prominent lines of the rgs spectrum have been measured individually. Delta functions were folded through the instrumental response matrix in order to derive integrated line fluxes. No additional width to the delta functions was needed to fit the line shapes. A “constant” level was adjusted to take into account the real continuum or the “pseudo-continuum”, created by the overlap of numerous weak lines.

The line fluxes are compared between the quiescent and flaring state (Table 4). Especially the flux difference of the Fe XXV line at 1.84 Å between the two states is clear. This line has a far higher flux in the flare state than in the quiescent state. The same is true for other highly ionized (hot) ions, such as Si XIV, Si XIII, Fe XIX, Fe XVIII, and Fe XVII. However, for these ions the effect is less significant.

### 3.4 Density measurement

The three He-like transitions, consisting of the resonance line (r) 1s2 1S0 \- 1s2p 1P1, the intercombination line (i) 1s2 1S0 \- 1s2p 3P1, and the forbidden line (f) 1s2 1S0 \- 1s2s 3S1 are density- and temperature-dependent (e.g., Gabriel and Jordan 1969). At increasing electron density the low-lying 1s2s 3S1 level will be depopulated in favour of the 1s2p 3P1 level, resulting in a decrease of the intensity of the forbidden line and an increase of the intercombination line. This makes the ratio R = f/i a valuable diagnostic tool for densities. The He-like lines of oxygen are well separated, although the intercombination line in the quiescent spectrum is affected by an innershell transition in O vi (see Fig. 6). These lines at 21.602 Å (r), 21.802 Å (i) and 22.101 Å (f) have been used to obtain temperature and density values from the G- and R-ratio using the tables by Porquet et al. (2001). For representative temperatures around 3 MK, we derive electron densities of $`3\pm 2\times 10^{10}`$ cm-3 for the quiescent state and $`8_{7.5}^{+20}\times 10^{10}`$ cm-3 for the flaring state. Within the errors these values agree with those given in Table 2.

## 4 Conclusions

The X-ray spectrum of AT Mic obtained by rgs and epic-mos is dominated by H- and He-like transitions of Ne, O, N, and C and by Fe xvii lines (rgs) and by H- and He-like transitions of Si and Mg and Fe K-shell transitions (epic-mos).

The spectra correspond to a quiescent and a flaring state as was concluded from the observed light curve. The dominant temperatures and emission measures cover a broad temperature range from 1 up to about 60 MK. The total emission measure in this temperature domain is $`12.2\times 10^{51}`$ cm-3 for the quiescent state and $`19.5\times 10^{51}`$ cm-3 for the flaring state. Both values are obtained from the “fit weight data” fits.

From the O vii lines a density of $`2\times 10^{10}`$ cm-3 is derived for the quiescent state and of $`4\times 10^{10}`$ cm-3 for the flaring state. The density of the flaring state is somewhat higher, but due to the poor statistics this is not significant.

Most line features of the rgs spectra are produced in the temperature range from 3-10 MK, dominated by quiescent emission. Lines from highly ionized atoms are enhanced in the flare state spectrum.

The abundance pattern in AT Mic shows the depletion of low-FIP elements ($`<`$10 eV) relative to high-FIP elements in the quiescent state. This indicates the presence of a possible IFIP effect in this active star, as has been shown for other active stars by Audard et al. (2002), Audard et al. (2003), Brinkman et al. (2001), Güdel et al. 2001, and Güdel et al. (2002). However, in the flare state this effect appears to be suppressed. This may imply that during the flare material from deeper layers with photosperic/chromospheric abundances is transported into the corona.

###### Acknowledgements.

The SRON National Institute for Space Research is supported financially by NWO. The PSI group acknowledges support from the Swiss National Science Foundation (grant 2100-049343). MA acknowledges support from the Swiss National Science Foundation (fellowship 81EZ-67388).


text:
# Surface plasmon toy-model of a rotating black hole.

\[

## Abstract

Recently introduced surface plasmon toy black hole model has been extended in order to emulate a rotating black hole (Kerr metric). Physical realization of this model involves a droplet of an optically active liquid on the metal surface which supports propagation of surface plasmons. Such droplets are shown to exhibit giant optical activity in the frequency range near the surface plasmon resonance of a metal-liquid interface.

PACS no.: 78.67.-n, 04.70.Bw

\]

The realization that solid-state toy models may help in an understanding of electromagnetic phenomena in curved space-time has led to considerable recent effort in developing toy models of electromagnetic  and sonic  black holes. In the case of media electrodynamics this is possible because of an analogy between the propagation of light in matter and in curved space-times: it is well known that the Maxwell equations in a general curved space-time background $`g_{ik}(x,t)`$ are equivalent to the phenomenological Maxwell equations in the presence of a matter background with nontrivial electric and magnetic permeability tensors $`ϵ_{ij}(x,t)`$ and $`\mu _{ij}(x,t)`$ . In this analogy, the event horizon corresponds to a surface of singular electric and magnetic permeabilities, so that the speed of light goes to zero, and light is ”frozen” near such a surface. In the absence of established quantum gravitation theory the toy models are helpful in understanding electromagnetic phenomena in curved space-times, such as Hawking radiation  and the Unruh effect . Unfortunately, until recently all the suggested electromagnetic black hole toy-models were very difficult to realize and study experimentally, so virtually no experimental work was done in this field. Very recently a simple surface plasmon toy black hole model has been suggested and observed in the experiment . This model is based on the surface plasmon resonance at the edge of a liquid droplet deposited on the metal surface, which supports propagation of surface plasmons. Surface plasmons are collective excitations of the conductive electrons and the electromagnetic field . They exist in ”curved three-dimensional space-times” defined by the shape of the metal-dielectric interface. Since in many experimental geometries surface plasmons are weakly coupled to the outside world (to free-space photons, phonons, single-electron excitations, etc.) it is reasonable to treat the physics of surface plasmons separately from the rest of the surface and bulk excitations, so that a field-theory of surface plasmons in a curved space-time background may be considered. For example, a nanohole in a thin metal membrane may be treated as a ”wormhole” connecting two ”flat” surface plasmon worlds located on the opposite interfaces of the membrane . On the other hand, near the plasmon resonance (which is defined by the condition that $`ϵ_m(\omega )=ϵ_d`$ at the metal-dielectric interface, where $`ϵ_m(\omega )`$ and $`ϵ_d`$ are the dielectric constants of metal and dielectric, respectively ) the surface plasmon velocity vanishes, so that the surface plasmon ”stops” on the metal surface, and the surface charge and the normal component of the electric field diverge. For a given frequency of light, the spatial boundary of the plasmon resonance (”the event horizon” of our toy model) may be defined at will using the geometry $`ϵ_d(x,y)`$ of the absorbed layer of dielectric on the metal surface. Thus, the plasmon resonance becomes a natural candidate to emulate the event horizon of a black hole. As a result, toy two-dimensional surface plasmon black holes can be easily produced and studied .

In what follows I am going to extend the surface plasmon toy black hole model in order to emulate the space-time metric of a rotating black hole (Kerr metric ). Physical realization of this model involves a droplet of an optically active liquid on the metal surface which supports propagation of surface plasmons. I will show that the droplets of such liquids on the metal surface exhibit giant optical activity in the frequency range near the surface plasmon resonance of a metal-liquid interface.

As a first step, let us consider in detail the dispersion law of a surface plasmon (SP), which propagates along the metal-dielectric interface. The SP field decays exponentially both inside the metal and the dielectric. Inside the dielectric the decay exponent is roughly equal to the SP wave vector. As a first step let us assume that both metal and dielectric completely fill the respective $`z<0`$ and $`z>0`$ half-spaces. In such a case the dispersion law can be written as

$$k^2=\frac{\omega ^2}{c^2}\frac{ϵ_dϵ_m(\omega )}{ϵ_d+ϵ_m(\omega )},$$
(1)
where we will assume that $`ϵ_m=1\omega _p^2/\omega ^2`$ according to the Drude model, and $`\omega _p`$ is the plasma frequency of the metal. This dispersion law is shown in Fig.1(b) for the cases of metal-vacuum and metal-dielectric interfaces. It starts as a ”light line” in the respective dielectric at low frequencies and approaches asymptotically $`\omega =\omega _p/(1+ϵ_d)^{1/2}`$ at very large wave vectors. The latter frequency corresponds to the so-called surface plasmon resonance. Under the surface plasmon resonance conditions both phase and group velocity of the SPs is zero, and the surface charge and the normal component of the electric field diverge. Since at every wavevector the SP dispersion law is located to the right of the ”light line”, the SPs of the plane metal-dielectric interface are decoupled from the free-space photons due to the momentum conservation law.

If a droplet of dielectric (Fig.1(a)) is placed on the metal surface, the SP dispersion law will be a function of the local thickness of the droplet. Deep inside the droplet far from its edges the SP dispersion law will look similar to the case of a metal-dielectric interface, whereas near the edges (where the dielectric is thinner) it will approach the SP dispersion law for the metal-vacuum interface. As a result, for every frequency between $`\omega _p/(1+ϵ_d)^{1/2}`$ and $`\omega _p/2^{1/2}`$ there will be a closed linear boundary inside the droplet for which the surface plasmon resonance conditions are satisfied. Let us show that such a droplet of dielectric on the metal interface behaves as a ”surface plasmon black hole” in the frequency range between $`\omega _p/(1+ϵ_d)^{1/2}`$ and $`\omega _p/2^{1/2}`$, and that the described boundary of the surface plasmon resonance behaves as an ”event horizon” of such a black hole.

Let us consider a SP within this frequency range, which is trapped near its respective ”event horizon”, and which is trying to leave a large droplet of dielectric (see Fig.2). The fact that the droplet is large means that ray optics may be used. Since the component of the SP momentum parallel to the droplet boundary has to be conserved, such a SP will be totally internally reflected by the surface plasmon resonance boundary back inside the droplet at any non-zero angle of incidence. This is a simple consequence of the fact that near the ”event horizon” the effective refractive index of the droplet for surface plasmons is infinite (according to eq.(1), both phase and group velocity of surface plasmons is zero at surface plasmon resonance). On the other hand, even if the angle of incidence is zero, it will take the SP infinite time to leave the resonance boundary. Thus, the droplet behaves as a black hole for surface plasmons, and the line near the droplet boundary where the surface plasmon resonance conditions are satisfied plays the role of the event horizon for surface plasmons.

The toy black holes described above are extremely easy to make and observe . In our experiments a small droplet of glycerin was placed on the gold film surface and further smeared over the surface using lens paper, so that a large number of glycerin microdroplets were formed on the surface (Fig.3(a)). These microdroplets were illuminated with white light through the glass prism (Fig.1(a)) in the so-called Kretschman geometry . The Kretschman geometry allows for efficient SP excitation on the gold-vacuum interface due to phase matching between the SPs and photons in the glass prism. As a result, SPs were launched into the gold film area around the droplet. Photograph taken under a microscope of one of such microdroplets is shown in Fig.3(b). The white rim of light near the edge of the droplet is clearly seen. It corresponds to the effective SP event horizon described above. Near this toy event horizon SPs are stopped or reflected back inside the droplet. In addition, a small portion of the SP field may be scattered out of the two-dimensional surface plasmon world into normal three-dimensional photons. These photons produced the image in Fig.3(b). We also conducted near-field optical measurements of the local surface plasmon field distribution around the droplet boundary Fig.3(c,d,e) using a sharp tapered optical fiber as a microscope tip. These measurements were performed similar to the measurements of surface plasmon scattering by individual surface defects described in . The droplet was illuminated with 488 nm laser light in the Kretschman geometry. The tip of the microscope was able to penetrate inside the glycerin droplet, and measure the local plasmon field distribution both inside and outside the droplet. Inside the droplet (in the right half of the images) the shear-force image (c) corresponds to the increase in viscous friction rather than the droplet topography. However, this image accurately represents the location of the droplet boundary, shown by the arrow in Fig.3(e). The sharp and narrow local maximum of the surface plasmon field just inside the droplet near its boundary is clearly visible in the near-field image Fig.3(d) and its cross-section Fig.3(e).

Unfortunately, the described toy SP black hole model does not work outside the frequency range between $`\omega _p/(1+ϵ_d)^{1/2}`$ and $`\omega _p/2^{1/2}`$. On the other hand, this is a common feature of every electromagnetic toy black hole model suggested so far. All such toy models necessarily work only within a limited frequency range. In addition, the losses in metal and dielectric described by the so far neglected imaginary parts of their dielectric functions will put a stop to the singularities of the field somewhere very near the toy event horizon (however, such ”good” metals as gold and silver have very small imaginary parts of their dielectric constants, and hence, very pronounced plasmon resonances ). Notwithstanding these limitations, the ease of making and observing such toy SP black holes makes them a very promising research object. If we forget about the language of ”black holes” and ”event horizons” for a moment, the SP optics phenomenon represented in Fig.3(b) remains a potentially very interesting effect in surface plasmon optics. Namely, this photo shows the existence of a two-dimensional SP analog of whispering gallery modes, which are well-known in the optics of light in droplets and other spherical dielectric particles. Whispering gallery modes in liquid microdroplets are known to substantially enhance nonlinear optical phenomena due to cavity quantum electrodynamic effects . One may expect even higher enhancement of nonlinear optical mixing in liquid droplets on the metal surfaces due to enhancement of surface electromagnetic field inherent to surface plasmon excitation, and in addition, due to accumulation of SP energy near the surface plasmon event horizons at the droplet boundaries. This strong enhancement of nonlinear optical effects in liquid droplets may be very useful in chemical and biological sensing applications.

Let us now try and extend this surface plasmon toy black hole model so that a rotating black hole metric (Kerr metric) can be emulated. The relationship between the $`\stackrel{}{D}`$ and $`\stackrel{}{E}`$ fields in a gravitational field with the metric $`g_{ik}`$ can be written in the form :

$$\stackrel{}{D}=\frac{\stackrel{}{E}}{h^{1/2}}+\stackrel{}{H}\times \stackrel{}{g},$$
(2)
where $`h=g_{00}`$ and $`g_\alpha =g_{0\alpha }/g_{00}`$. For an electromagnetic wave this relationship is very similar to the relationship between the $`\stackrel{}{D}`$ and $`\stackrel{}{E}`$ fields in the optically active medium, if we identify $`1/h^{1/2}`$ as $`ϵ`$ of the medium, and recall how optical activity is introduced in media electrodynamics.

There are different ways of introducing optical activity (gyration) tensor in the macroscopic Maxwell equations. It can be introduced in a symmetric form, which is sometimes called Condon relations :

$$\stackrel{}{D}=ϵ\stackrel{}{E}+\gamma \frac{\stackrel{}{B}}{t}$$
(3)
$$\stackrel{}{H}=\mu ^1\stackrel{}{B}+\gamma \frac{\stackrel{}{E}}{t}$$
(4)
Or it can be introduced only in an equation for $`\stackrel{}{D}`$ (see ). In our consideration I will follow Landau and Lifshitz , and for simplicity use only the following equation valid in isotropic or cubic-symmetry materials:

$$\stackrel{}{D}=ϵ\stackrel{}{E}+i\stackrel{}{E}\times \stackrel{}{g},$$
(5)
where $`\stackrel{}{g}`$ is called the gyration vector. If the medium exhibits magneto- optical effect and does not exhibit natural optical activity $`\stackrel{}{g}`$ is proportional to the magnetic field $`\stackrel{}{H}`$:

$$\stackrel{}{g}=f\stackrel{}{H},$$
(6)
where the constant $`f`$ may be either positive or negative. For metals in the Drude model at $`\omega >>eH/mc`$

$$f(\omega )=\frac{4\pi Ne^3}{cm^2\omega ^3}=\frac{e\omega _p^2}{mc\omega ^3},$$
(7)
where $`\omega _p`$ is the plasma frequency and $`m`$ is the electron mass . In any case, the local difference in refractive index for the left $`n_{}`$ and right $`n_+`$ circular polarizations of light is proportional to the local value of the gyration vector $`\stackrel{}{g}`$.

Thus, from the relationships between the $`\stackrel{}{D}`$ and $`\stackrel{}{E}`$ fields in a gravitational field and in media electrodynamics we conclude that space-time behaves as a chiral optical medium if $`\stackrel{}{g}0`$ (we must take into account different conventions for the use of imaginary numbers implemented in  and , and hence equations (2) and (5)).

A space-time region with $`\stackrel{}{g}0`$ can be described locally as a rotating coordinate frame with an angular velocity

$$\stackrel{}{\mathrm{\Omega }}=\frac{ch^{1/2}}{2}\stackrel{}{}\times \stackrel{}{g}$$
(8)
Similar vector $`\stackrel{}{\mathrm{\Omega }}`$ field can be defined for any chiral medium regardless of the nature of its optical activity (natural or magnetic field induced). Thus, electrodynamics of such medium may be understood as if there is a distribution of local angular rotation $`\stackrel{}{\mathrm{\Omega }}`$ field inside the medium.

Now it is clear that in order to emulate a rotating black hole we must use a droplet of liquid which exhibits optical activity. Let us consider a medium with $`\stackrel{}{g}=\stackrel{}{g}(r)`$ directed along the $`\varphi `$-coordinate (Fig.4(b)). Such a distribution may be produced, for example, in a twisted layer of liquid crystal or in a liquid droplet exhibiting magneto-optical effect if a droplet is formed around a small cylindrical wire and a current is passed through the wire. Such a geometry will result in different refractive indices $`n_+`$ and $`n_{}`$ of the liquid as seen by the surface plasmons, which propagate in the clockwise (right) and counterclockwise (left) directions near the droplet boundary. As a result, these plasmons will have slightly different dispersion laws (Fig.4(a)). This difference produces two main effects: the position of the effective event horizon for the left and right plasmons will be different inside the liquid droplet (Fig.4(b)); and there will be a narrow frequency range (located near the surface plasmon resonance of the metal-liquid interface) where the event horizon will exist for only one (left or right) kind of surface plasmons (Fig.4(a)).

The first effect closely resembles the properties of the Kerr metric, which characterizes a rotating black hole. The Kerr metric has two characteristic surfaces around a rotating black hole. There exists a spherical event horizon similar to the case of a non-rotating black hole, and there is a so called ergosphere  just outside the black hole event horizon, which corresponds to the area of space where every particle must rotate in the direction of the black hole rotation. There are no particle states, which rotate in the opposite direction inside the ergoshere. Similar to the case of a real rotating black hole, our toy surface plasmon black hole has a toy-ergosphere: between the two event horizons shown in Fig.4(b) there are no surface plasmon states which rotate in the direction opposite to the direction of the effective $`\stackrel{}{\mathrm{\Omega }}`$ field of the toy black hole. This area of space is quite unusual from the point of view of surface plasmons. We may say that even the surface plasmon vacuum rotates in this area, similar to the description of the electromagnetic vacuum inside the band gap of a photonic crystal as being ”colored” vacuum. Further study of the linear and nonlinear optical properties of such rotating vacuum may be quite interesting from the fundamental optics point of view.

The second effect may be quite important from the point of view of practical applications. In the frequency range where the effective event horizon exists only for one (left or right) kind of plasmons, the array of droplets will behave as a medium which exhibits giant planar optical activity. Such media are being actively developed at the moment .

In conclusion, recently introduced surface plasmon toy black hole model has been extended in order to emulate a rotating black hole (Kerr metric). Physical realization of this model involves a droplet of an optically active liquid on the metal surface which supports propagation of surface plasmons. Such droplets are shown to exhibit giant optical activity in the frequency range near the surface plasmon resonance of a metal-liquid interface.

This work has been supported in part by the NSF grants ECS-0210438 and ECS-0304046.

Figure captions.

Fig.1 (a) Experimental geometry of surface plasmon toy black hole observation. (b) Surface plasmon dispersion law for the cases of metal-vacuum and metal-dielectric interfaces

Fig.2 A surface plasmon trapped inside a droplet near the effective ”event horizon”: The projection of surface plasmon momentum parallel to the droplet edge must be conserved. Due to effectively infinite refractive index near the droplet edge surface plasmons experience total internal reflection at any angle of incidence.

Fig.3 Far-field (a,b) and near-field (c,d) images of a toy surface plasmon black hole: (a) Droplet of glycerin on a gold film surface (illuminated from the top). The droplet diameter is approximately 15 micrometers. (b) The same droplet illuminated with white light in the Kretschman geometry, which provides efficient coupling of light to surface plasmons on the gold-vacuum interface (Fig.1(a)). The white rim around the droplet boundary corresponds to the effective surface plasmon ”event horizon”. (c) and (d) show $`10\times 10\mu m^2`$ topographical and near-field optical images of a similar droplet boundary (droplet is located in the right half of the images) illuminated with 488 nm laser light. Cross-sections of both images are shown in (e). Position of a droplet boundary is indicated by the arrow.

Fig.4 (a) Surface plasmon dispersion law for the metal-vacuum and metal-dielectric interfaces in the case of optically active dielectric. (b) Difference in refractive indices for the left and right surface plasmons results in different locations of the effective event horizons for these plasmons. The local directions of the gyration filed are shown by the arrows.


text:
# 1 We plot the DE analysis for different synthetic series generated with the GP model and three different exponents of the Pareto law, i.e. 𝛾={3.5,2.25,1.25}. The values of 𝛿 obtained by fitting the asymtotic behavior are also reported with the corresponding error. These values increase as 𝛾 diminishes passing respectively from 0.52±0.02 (for 𝛾=3.5) to 0.78±0.02 (for 𝛾=2.25) and finally to 0.99±0.02 (for 𝛾=1.25).

Reply to the Comment “Intercluster correlation in Seismicity” by Helmstetter and Sornette

The comment by Helmstetter and Sornette  addresses an important point of our paper  that needs to be clarified. Any scaling measure, including also the diffusion entropy (DE) method, when applied to the earthquake time series, yields anomalous scaling for two reasons, long-range time correlation (intermittency) and the Pareto’s statistics. In fact, the DE random walker position is determined by the number of events (earthquakes) in a given time interval, which, in turns, is related to both the number of clusters and the cluster size. The distribution of the time distances between clusters, $`\psi (\tau ^{[m]})`$ in ref., determines the number of clusters, while the Pareto’s law

$$p(n)n^\gamma $$
(1)
gives the probability to get a cluster of size $`n`$. The resulting asymptotic scaling detected by the DE is the scaling of the most anomalous of the two processes, i.e. the process with the lowest power law exponent.

The authors of the comment  discuss the Generalized Poisson (GP) model with an exponent $`\gamma =2.25`$ , which generates $`\delta =0.8`$. In fact, according to the Continous Time Random Walk theory , the scaling is $`\delta =1/(\gamma 1)`$ if $`2<\gamma <3`$, $`\delta =0.5`$ for $`\gamma >3`$ and $`\delta =1`$ for $`\gamma <2`$.

We have repeated the calculations with $`\gamma =2.25`$ and the results, shown in fig.1, confirm the theoretical prediction. In the figure we plot also two other cases in order to illustrate how $`\delta `$ changes as a function of $`\gamma `$. In Ref. we explored the range $`\gamma =2.53`$.

The experimental data yield a value $`\delta =0.94`$ confirming that $`\psi (\tau ^{[m]})`$ has a power law tail with an exponent $`\mu =2.06`$. Moreover, the fact that the scaling observed in fig. 2 of our previous paper  is not influenced by the threshold of magnitude considered, is an additional evidence that what we are actually observing is the real scaling determined by the clusters distances.

The LR model  is a minimal generalization of the GP model (only $`\psi (\tau ^{[m]})`$ changes). Therefore it is very interesting that a more sophisticated model, like the ETAS, which reproduces also other properties of seismicity, yields the same $`\delta `$ observed in the LR and in the experimental data. However we would like to point out that we define $`\tau ^{[m]}`$ as the distance between two consecutive clusters, while in  this symbol refers to distance between two consecutive earthquakes with magnitude greater than $`M=5`$. So there is no contraddiction between the correlations found in  and our results.

Mirko S. Mega1, Paolo Allegrini2, Paolo Grigolini1,3,4, Vito Latora5, Luigi Palatella1, Andrea Rapisarda5 and Sergio Vinciguerra5,6

1Dipartimento di Fisica dell’Università di Pisa and INFM, via Buonarroti 2, 56127 Pisa, Italy

2Istituto di Linguistica Computazionale del CNR, Area della Ricerca di Pisa, Via G. Moruzzi 1, 56124, Pisa, Italy

3Center for Nonlinear Science, University of North Texas, P.O. Box 311427, Denton, Texas 76203-1427

4 Istituto di Biofisica del CNR, Area della Ricerca di Pisa, Via Alfieri 1, San Cataldo, 56010, Ghezzano-Pisa, Italy

5 Dipartimento di Fisica e Astronomia, Università di Catania, and INFN sezione di Catania, Via S. Sofia 64, 95123 Catania, Italy

6 Osservatorio Vesuviano - INGV, Via Diocleziano 328, 80124 Napoli, Italy

PACS numbers: 91.30.Dk,05.45.Tp,05.40.Fb