format: text

text:
# On the problem of computing the well-founded semantics11footnote 1A preliminary version of this paper appeared in the Proceedings of Computational Logic – CL 2000, Lecture Notes in Artificial Intelligence, 1861, Springer Verlag, 2000.

## 1 Introduction

The well-founded semantics was introduced in \[Van Gelder et al., 1991\] to provide 3-valued interpretations to logic programs with negation. Since its introduction, the well-founded semantics has become one of the most widely studied and most commonly accepted approaches to negation in logic programming \[Alferes et al., 1995, Fitting, 1991, Chen et al., 1995, Chen & Warren, 1996, Zukowski et al., 1997, Brass & Dix, 1998\]. It was implemented in several top-down reasoning systems, most prominent of which is XSB \[Rao et al., 1997\].

The well-founded semantics is closely related to the stable-model semantics \[Gelfond & Lifschitz, 1988\], another major approach to logic programs with negation. The well-founded semantics approximates the stable-model semantics \[Van Gelder et al., 1991, Fitting, 2001\]. Moreover, computing the well-founded model of propositional programs is polynomial \[Van Gelder, 1989\] while computing stable models is NP-hard \[Marek & Truszczyński, 1991\]. Consequently, evaluating the well-founded semantics can be used as an effective preprocessing technique in algorithms to compute stable models \[Subrahmanian et al., 1995\]. In addition, as demonstrated by smodels \[Niemelä & Simons, 1996\], at present the most advanced and most efficient system to compute stable models of DATALOG<sup>¬</sup> programs, the well-founded semantics can be used as a powerful lookahead mechanism.

Despite the importance of the well-founded semantics, the question of how fast it can be computed has not attracted significant attention. Van Gelder \[Van Gelder, 1989\] described the so called alternating-fixpoint algorithm. Van Gelder’s algorithm runs in time $`O(|\text{At}(P)|\times size(P))`$, where $`\text{At}(P)`$ is the set of atoms occurring in a logic program $`P`$, $`|\text{At}(P)|`$ denotes the cardinality of $`\text{At}(P)`$, and $`size(P)`$ is the size of $`P`$ (the total number of atom occurrences in $`P`$). Improving on this algorithm turned out to be difficult. The first progress was obtained in \[Berman et al., 1995\]. The algorithm described there, when restricted to programs whose rules contain at most two positive occurrences of atoms in their bodies, runs in time $`O(|\text{At}(P)|^{4/3}|P|^{2/3}+size(P))`$, where $`|P|`$ is the number of rules in $`P`$. For programs whose rules have no more than one positive atom in the body a better estimate of $`O(|\text{At}(P)|^{3/2}|P|^{1/2}+size(P))`$ was obtained. For some classes of programs this is an asymptotically better estimate than the $`O(|\text{At}(P)|\times size(P))`$ estimate that holds for the algorithm by Van Gelder.

A different approach to computing the well-founded model was proposed in \[Zukowski et al., 1997, Brass et al., 2001\]. It is based on the notion of a program transformation \[Brass & Dix, 1998\]. The authors describe there several transformations that can be implemented in linear time and that simplify a program while (essentially) preserving the well-founded semantics. These transformations are: the positive reduction, success, negative reduction, and failure (PSNF transformations, for short). They allow one to compute in linear time the Kripke-Kleene semantics \[Fitting, 1985\] of the program. To compute the well-founded semantics one also needs to detect the so-called positive loops. The complexity of this task dominates the asymptotic complexity of the well-founded semantics computation. No improved algorithms for the positive-loop detection are offered in \[Brass et al., 2001\] so the worst-case asymptotic complexity of the algorithm presented there remains the same as that of the alternating-fixpoint method. However, due to the use of PSNF transformations, that simplify the program, the algorithm based on program transformations may in practice run faster. In contrast to the approach studied in \[Brass et al., 2001\], we focus here on the positive-loop detection task.

The alternating-fixpoint algorithm works by successively improving lower approximations $`T`$ and $`F`$ to the sets of atoms that are true and false (under the well-founded semantics), respectively. The algorithm starts with $`T=\mathrm{}`$. Using this estimate, it computes the first estimate for $`F`$. Next, using this estimate, in turn, it computes a better estimate for $`T`$. The algorithm continues until further improvements are not possible. It returns the final sets $`T`$ and $`F`$ as the well-founded semantics. A dual version of the alternating-fixpoint algorithm, starting with $`F=\mathrm{}`$ and then alternatingly computing approximations to $`T`$ and $`F`$, is also possible. The most time-consuming part of this algorithm is in computing estimates to the set of atoms that are false (in this part, in particular, positive loops are detected). In the Van Gelder algorithm, the best possible approximation (given the current estimate for $`T`$) is always computed by using a bottom-up approach.

In this paper we focus on the problem of detecting positive loops and computing new false atoms. We restrict our attention to the class of programs that have at most one positive atom in the body. We denote this class of programs by $`LP_1`$. We show that for programs from $`LP_1`$, false atoms can be computed by means of a top-down approach by finding atoms that do not have a proof. Moreover, we show that it is not necessary to find all atoms that can be established to be false at a given stage. Finding a proper subset (as long as it is not empty) is also sufficient and results in a correct algorithm. We apply these techniques to design a version of an alternating-fixpoint algorithm computing the well-founded semantics of programs from the class $`LP_1`$. We demonstrate that the resulting algorithm is asymptotically better than the original alternating-fixpoint algorithm by Van Gelder. Specifically, we show that our algorithm runs in time $`O(|\text{At}(P)|^2+size(P))`$. Thus, for programs with $`size(P)|\text{At}(P)|^2`$, our algorithm runs in linear time and is asymptotically optimal! It is also easy to see that when $`|P|>|\text{At}(P)|`$, the asymptotic estimate of the running time of our algorithm is better than that of algorithms by Van Gelder \[Van Gelder, 1989\] and Berman et al. \[Berman et al., 1995\].

As mentioned above, our approach is restricted to the class $`LP_1`$. Applicability of our method can, however, be slightly extended. Let us denote by $`LP_1^+`$ the class of these logic programs that, after simplifying by means of PSNF transformations (or, equivalently, with respect to the Kripke-Kleene semantics) fall into the class $`LP_1`$. Since PSNF transformations (the Kripke-Kleene semantics) can be computed in linear time, the asymptotic estimate of the running time of our method extends to all programs in the class $`LP_1^+`$.

The paper is organized as follows. In the next section we provide a brief review of the key notions and terminology. In Section 3 we describe several modifications to the original Van Gelder algorithm, we show their correctness and estimate their running time. The ultimate effect of our considerations there is a general template for an algorithm to compute the well-founded semantics. Any algorithm computing some (not necessarily all) atoms that can be established as false given a current estimate to the well-founded can be used with it. One such algorithm, for programs from the class $`LP_1`$, is described and analyzed in Section 4. It constitutes the main contribution of the paper and yields a new, currently asymptotically most efficient algorithm for computing the well-founded semantics for programs in $`LP_1`$. The last section contains conclusions.

## 2 Preliminaries

We start by reviewing basic concepts and notation related to logic programs and the well-founded semantics, as well as some simple auxiliary results. In the paper we consider the propositional case only.

Let $`P`$ be a normal logic program. By $`\text{At}(P)`$ we denote the set of atoms occurring in $`P`$. Let $`M\text{At}(P)`$ (throughout the paper we often drop a reference to $`P`$ from our notation, whenever there is no danger of ambiguity). By $`P_M`$ we denote the program obtained from $`P`$ by removing all rules whose bodies contain negated literals of the form $`\text{not}(a)`$, where $`aM`$. Further, by $`P^h`$ we denote the program obtained from $`P`$ by removing from the bodies of its rules all negative literals. Clearly, the program $`(P_M)^h`$ coincides with the Gelfond-Lifschitz reduct of $`P`$ with respect to $`M`$ (throughout the paper, we write $`P_M^h`$ for $`(P_M)^h`$, to simplify notation). The Gelfond-Lifschitz operator on the algebra of all subsets of At, $`GL`$ (following our convention, we omit the reference to $`P`$ from the notation), is defined by

$$GL(M)=LM(P_M^h),$$

where $`LM(Q)`$ stands for a least model of a Horn program $`Q`$.

We now present characterizations of the well-founded semantics. We phrase them in the language of operators and their fixpoints. All operators considered here are defined on the algebra of subsets of $`\text{At}(P)`$. We denote a least fixpoint (if it exists) of an operator $`O`$ by $`\text{lfp}(O)`$.

It is well known that $`GL`$ is antimonotone. Consequently, $`GL^2=GLGL`$ is monotone and has a least fixpoint. The set of atoms that are true with respect to the well-founded semantics of a program $`P`$, denoted by $`T_{wfs}`$, is precisely the least fixpoint of the operator $`GL^2`$, that is, $`T_{wfs}=\text{lfp}(GL^2)`$ \[Van Gelder, 1989, Fitting, 2001\]. The set of atoms that are false with respect to the well-founded semantics of a program $`P`$, denoted by $`F_{wfs}`$, is given by $`\overline{GL(T_{wfs})}`$ (throughout the paper, $`\overline{X}`$ denotes the complement of a set $`X`$ with respect to $`\text{At}(P)`$).

One can define a dual operator to $`GL^2`$ by

$$A(M)=\overline{GL(GL(\overline{M}))}.$$

It is easy to see that $`A`$ is monotone and that its least fixpoint is $`F_{wfs}`$. Thus, $`F_{wfs}=\text{lfp}(A)`$ and $`T_{wfs}=GL(\overline{F_{wfs}})`$.

We close this section by discussing ways to compute $`GL(M)`$ for a given finite propositional logic program $`P`$ and a set of atoms $`M\text{At}(P)`$. A straightforward approach is to compute the Gelfond-Lifschitz reduct $`P_M^h`$ and then to compute its least model. The resulting algorithm is asymptotically optimal as it runs in time linear in the size of the program. However, in this paper we will use a different approach, more appropriate for the computation of the well-founded semantics. Let $`P`$ be a logic program with negation. We define $`\text{At}^{}(P)=\{\text{not}(a):a\text{At}(P)\}`$. For every set $`M\text{At}(P)\text{At}^{}(P)`$, we define $`true(M)=M\text{At}(P)`$. If we interpret literals of $`\text{At}^{}(P)`$ as new atoms, then for every set $`M\text{At}(P)`$, the program $`P\text{not}(M)`$ can be viewed as a Horn program. Thus, it has a least model. It is easy to see that

$$GL_P(M)=true(LM(P\text{not}(\overline{M}))).$$

Here, $`P`$ appearing at the left-hand side of the equation stands for the original logic program, while $`P`$ appearing at the right-hand side of the equation stands for the same program but interpreted as a Horn program. Thus, using the algorithm of Dowling and Gallier \[Dowling & Gallier, 1984\], the Gelfond-Lifschitz reduct can be computed in time $`O(size(P)+|M|)=O(size(P))`$ (since $`M\text{At}(P)`$, $`|M|=O(size(P))`$).

## 3 Algorithms

The departure point for our discussion of algorithms to compute the well-founded semantics is the alternating-fixpoint algorithm of Van Gelder \[Van Gelder, 1989\]. Using the terminology introduced in the previous section it can be formulated as follows.

|  |  |  |  |
| --- | --- | --- | --- |
| Algorithm 1 (Van Gelder) |
|  | $`F:=\mathrm{}`$; |
|  | repeat |
|  |  | $`T:=true(LM(P\text{not}(F))`$; (\* or equivalently: $`T:=GL(\overline{F})`$; \*) |
|  |  | $`F:=\overline{LM(P_T^h)}`$; (\* or equivalently: $`\overline{GL(T)}`$; \*) |
|  | until no change in $`F`$; |
|  | return $`T`$ and $`F`$. |

Let $`F^{}`$ and $`F^{\prime \prime }`$ be the values of the set $`F`$ just before and just after an iteration of the repeat loop in Algorithm 1. Clearly,

$$F^{\prime \prime }=\overline{GL(GL(\overline{F^{}}))}=A(F^{}).$$

Thus, after iteration $`i`$ of the repeat loop, $`F=A^i(\mathrm{})`$. Consequently, it follows from our earlier remarks that when Algorithm 1 terminates, the set $`F`$ that is returned satisfies $`F=F_{wfs}`$. Since there is no change in $`F`$ in the last iteration, when the algorithm terminates, we have $`T=T_{wfs}`$. That is, Algorithm 1 is correct.

We will now modify Algorithm 1. The basis for Algorithm 1 is the operator $`A`$. This operator is not progressive. That is, $`M`$ is not necessarily a subset of $`A(M)`$. We will now introduce a related progressive operator, say $`B`$, and show that it can be used to replace $`A`$. Let $`P`$ be a logic program and let $`T`$ and $`F`$ be two subsets of $`\text{At}(P)`$. By $`P_{F,T}`$ we denote the program obtained from $`P`$ by removing

1. all rules whose heads are in $`F`$
2. all rules whose bodies contain a positive occurrence of an atom from $`F`$
3. all rules whose bodies contain a negated literal of the form $`\text{not}(a)`$, where $`aT`$.

Clearly, $`P_{F,T}P_T`$.

We define an operator $`B(F)`$ as follows:

$$B(F)=\overline{LM(P_{F,T}^h)},$$

where $`T=GL(\overline{F})`$ and $`P_{F,T}^h`$ abbreviates $`(P_{F,T})^h`$. The following result gathers key properties of the operator $`B`$.

###### Theorem 3.1

Let $`P`$ be a normal logic program. Then:

1. $`B`$ is monotone
2. For every $`F\text{At}(P)`$, $`A(F)B(F)`$
3. For every $`FF_{wfs}`$, $`B(F)F_{wfs}`$
4. $`\text{lfp}(B)=F_{wfs}`$
5. For every $`F\text{At}(P)`$, $`B(F)=F(\overline{F}LM(P_{F,T}^h))`$, where $`T=GL(\overline{F})`$.

Proof: (1) Assume that $`F_1F_2`$. Set $`T_i=GL(\overline{F_i})`$, $`i=1,2`$. Clearly, $`\overline{F_2}\overline{F_1}`$ and, by antimonotonicity of $`GL`$, $`T_1T_2`$. By the definition of $`P_{F,T}`$, $`P_{F_2,T_2}P_{F_1,T_1}`$. Consequently, $`LM(P_{F_2,T_2}^h)LM(P_{F_1,T_1}^h)`$ and, so, $`B(F_1)B(F_2)`$.
(2) Let $`T=GL(\overline{F})`$. Clearly, $`P_{F,T}P_T`$. Thus, $`A(F)=\overline{LM(P_T^h)}\overline{LM(P_{F,T}^h)}=B(F)`$.
(3) We have, $`LM(P_{T_{wfs}}^h)=\overline{F_{wfs}}`$. It follows that removing from $`P_{T_{wfs}}^h`$ rules with heads in $`F_{wfs}`$ and those that contain an atom from $`F_{wfs}`$ in their bodies does not change the least model. That is,

$$LM(P_{F_{wfs},T_{wfs}}^h)=LM(P_{T_{wfs}}^h).$$

Since, $`T_{wfs}=GL(\overline{F_{wfs}})`$, $`B(F_{wfs})=\overline{LM(P_{F_{wfs},T_{wfs}}^h)}`$. Let $`FF_{wfs}`$. Then, by (1), $`B(F)B(F_{wfs})`$. Thus, we have

$$B(F)B(F_{wfs})=\overline{LM(P_{F_{wfs},T_{wfs}}^h)}=\overline{LM(P_{T_{wfs}}^h)}=F_{wfs}.$$

(4) The least fixpoint of $`B`$ is given by $`\text{lfp}(B)=B^i(\mathrm{})`$. By (3), $`\text{lfp}(B)F_{wfs}`$. On the other hand, by (1) and (2), $`A^i(\mathrm{})B^i(\mathrm{})`$. Thus, $`F_{wfs}=\text{lfp}(A)\text{lfp}(B)`$. It follows that $`\text{lfp}(B)=F_{wfs}`$.

(5) Let $`T=GL(\overline{F})`$. Since $`P_{F,T}`$ has no rules with head in $`F`$, $`LM(P_{F,T}^h)\overline{F}`$ and, consequently, $`FB(F)`$. Thus, the assertion follows. $`\mathrm{}`$

Theorem 3.1 allows us to prove the correctness of the following modification of Algorithm 1.

|  |  |  |  |
| --- | --- | --- | --- |
| Algorithm 2 |
|  | $`F:=\mathrm{}`$; |
|  | repeat |
|  |  | $`T:=true(LM(P\text{not}(F))`$; |
|  |  | $`\mathrm{\Delta }F:=\overline{F}LM(P_{F,T}^h)`$; |
|  |  | $`F:=F\mathrm{\Delta }F`$; |
|  | until no change in $`F`$; |
|  | return $`T`$ and $`F`$. |

By Theorem 3.1, each iteration of the repeat loop computes $`B(F)`$ as the new value for the set $`F`$. More formally, the set $`F`$ just after iteration $`i`$, satisfies $`F=B^i(\mathrm{})`$. Thus, when the algorithm terminates, the set $`F`$ that is returned is the least fixpoint of $`B`$. Consequently, by Theorem 3.1(4), Algorithm 2 is correct.

We will now modify Algorithm 2 to obtain a general template for an alternating-fixpoint algorithm to compute the well-founded semantics. The key idea is to observe that it is enough to compute a subset of $`\mathrm{\Delta }F`$ in each iteration and the algorithm remains correct.

Let us assume that for some operator $`\mathrm{\Delta }_w`$ defined for pairs $`(F,Q)`$, where $`F\text{At}(P)`$ and $`Q`$ is a Horn program such that $`\text{At}(Q)\overline{F}`$ (the complement is, as always, evaluated with respect to $`\text{At}(P)`$), we have:

$`\mathrm{\Delta }_w(F,Q)\overline{F}LM(Q)`$

$`\mathrm{\Delta }_w(F,Q)=\mathrm{}`$ if and only if $`\overline{F}LM(Q)=\mathrm{}`$.

Let $`F\text{At}(P)`$. By the definition of $`P_{F,T}`$, $`\text{At}(P_{F,T}^h)\overline{F}`$. Thus, we define $`B_w(F)=F\mathrm{\Delta }_w(F,P_{F,T}^h)`$, where $`T=true(LM(P\text{not}(F)))`$. It is clear that for every $`F\text{At}(P)`$, $`FB_w(F)B(F)`$, the latter inclusion follows from Theorem 3.1(5) and (W1). Consequently, for every $`i`$,

$$B_w^i(\mathrm{})B^i(\mathrm{}).$$

It follows that $`B_w^i(\mathrm{})\text{lfp}(B)=F_{wfs}`$. It also follows that there is the first $`i`$ such that $`B_w^i(\mathrm{})=B_w^{i+1}(\mathrm{})`$. Let us denote this set $`B_w^i(\mathrm{})`$ by $`F_0`$. Then $`F_0F_{wfs}`$. In the same time, by condition (W2), $`B(F_0)=F_0`$. Since $`F_{wfs}`$ is the least fixpoint of $`B`$, $`F_{wfs}F_0`$. It follows that a modification of Algorithm 2 in which line

$$\mathrm{\Delta }F:=\overline{F}LM(P_{F,T}^h);$$

is replaced by

$$\mathrm{\Delta }F:=\mathrm{\Delta }_w(F,P_{F,T}^h);$$

correctly computes the well-founded semantics of a program $`P`$. Thus, we obtain the following algorithm for computing the well-founded semantics.

|  |  |  |  |
| --- | --- | --- | --- |
| Algorithm 3 |
|  | $`F:=\mathrm{}`$; |
|  | repeat |
|  |  | $`T:=true(LM(P\text{not}(F))`$; |
|  |  | $`\mathrm{\Delta }F:=\mathrm{\Delta }_w(F,P_{F,T}^h)`$; |
|  |  | $`F:=F\mathrm{\Delta }F`$; |
|  | until no change in $`F`$; |
|  | return $`T`$ and $`F`$. |

We will now refine Algorithm 3. Specifically, we will show that the sets $`T`$ and $`F`$ can be computed incrementally.

Let $`R`$ be a Horn program. We define the residual program of $`R`$, $`res(R)`$, to be the Horn program obtained from $`R`$ by removing all rules of $`R`$ with the head in $`LM(R)`$ and by removing from the bodies of the remaining rules those elements that are in $`LM(R)`$. We have the following technical result.

###### Lemma 3.2

Let $`R`$ be a Horn program and let $`M`$ be a set of atoms such that $`Mhead(R)=\mathrm{}`$. Then $`LM(RM)=LM(R)LM(res(R)M)`$. $`\mathrm{}`$

Lemma 3.2 implies that (we treat here negated literals as new atoms and $`P`$ as Horn program over the extended alphabet)

$$LM(P\text{not}(F\mathrm{\Delta }F))=LM(P\text{not}(F))LM(res(P\text{not}(F))\text{not}(\mathrm{\Delta }F)).$$

Thus, if the set $`F`$ is expanded by new elements from $`\mathrm{\Delta }F`$, then the new set $`T`$ can be computed by increasing the old set $`T`$ by $`\mathrm{\Delta }T=true(LM(res(P\text{not}(F))\text{not}(\mathrm{\Delta }F)))`$. Important thing to note is that the increment $`\mathrm{\Delta }T`$ can be computed on the basis of the residual program and the increment $`\mathrm{\Delta }F`$. Similarly, we have

$$P_{F\mathrm{\Delta }F,T\mathrm{\Delta }T}=(P_{F,T})_{\mathrm{\Delta }F,\mathrm{\Delta }T}.$$

Thus, computing $`P_{F,T}`$ can also be done incrementally on the basis of the program considered in the previous iteration by taking into account most recently computed increments $`\mathrm{\Delta }F`$ and $`\mathrm{\Delta }T`$.

This discussion implies that Algorithm 3 can be equivalently restated as follows:

|  |  |  |  |
| --- | --- | --- | --- |
| Algorithm 3 |
| 1 | $`T:=F:=\mathrm{\Delta }T:=\mathrm{\Delta }F:=\mathrm{}`$; |
| 2 | $`R:=P`$; (\*$`R`$ will be treated as a Horn program \*) |
| 3 | $`Q:=P`$; |
| 4 | repeat |
| 5 |  | $`\mathrm{\Delta }T:=true(LM(R\text{not}(\mathrm{\Delta }F))`$; |
| 6 |  | $`R:=res(R\text{not}(\mathrm{\Delta }F))`$; |
| 7 |  | $`T:=T\mathrm{\Delta }T`$; |
| 8 |  | $`Q:=Q_{\mathrm{\Delta }F,\mathrm{\Delta }T}`$; |
| 9 |  | $`\mathrm{\Delta }F:=\mathrm{\Delta }_w(F,Q^h)`$; |
| 10 |  | $`F:=F\mathrm{\Delta }F`$; |
| 11 | until no change in $`F`$; |
| 12 | return $`T`$ and $`F`$. |

We will now estimate the running time of Algorithm 3. Clearly line 1 requires constant time. Setting up appropriate data structures for programs $`R`$ and $`Q`$ (lines 2 and 3) takes $`O(size(P))`$ steps. In each iteration, $`\mathrm{\Delta }T`$ is computed and the current program $`R`$ is replaced by the program $`\text{res}(R\text{not}(\mathrm{\Delta }F))`$ (lines 5 and 6). By modifying the algorithm from \[Dowling & Gallier, 1984\] and assuming that $`R`$ is already stored in the memory (it is avaliable either as the result of the initialization in the case of the first iteration or as a result of the computation in the previous iteration), both tasks can be accomplished in $`O(size(R^o)+|\mathrm{\Delta }F|size(R^n))`$ steps. Here $`R^o`$ denotes the old version of $`R`$ and $`R^n`$ denotes the new version of $`R`$. Consequently, the total time needed for lines 5 and 6 over all iterations is given by $`O(size(P)+|\text{At}(P)|size(R^t))=O(size(P))`$ (where $`R^t`$ is the program $`R`$, when the algorithm terminates). The time needed for all lines 7 is proportional to the number of iterations and is $`O(|\text{At}(P)|)=O(size(P))`$.

Given a logic program $`Q`$ and sets of atoms $`\mathrm{\Delta }T`$ and $`\mathrm{\Delta }F`$, it takes $`O(size(Q)size(Q_{\mathrm{\Delta }F,\mathrm{\Delta }T})+|\mathrm{\Delta }T|+|\mathrm{\Delta }F|)`$ steps to compute the program $`Q_{\mathrm{\Delta }F,\mathrm{\Delta }T}`$ in line 8. We assume here that $`Q`$ is already in the memory as a result of the initialization in the case of the first iteration, or as the result of the computation in the previous iteration, otherwise. It follows that the total time over all iterations needed to execute line 8 is $`O(size(P)+|\text{At}(P)|)=O(size(P))`$.

Thus, we obtain that the running time of Algorithm 3 is given by $`O(size(P)+m)`$, where $`m`$ is the total time needed to compute $`\mathrm{\Delta }_w(F,Q^h)`$ over all iterations of the algorithm.

In the standard (Van Gelder’s) implementation of Algorithm 3, we compute the whole set $`\overline{F}LM(Q^h)`$ as $`\mathrm{\Delta }_w(F,Q^h)`$. In addition, computation is performed in a bottom-up fashion. That is, we first compute the least model of $`Q^h`$ and then its complement with respect to $`\overline{F}`$. Such approach requires $`O(size(Q^h))=O(size(P))`$ steps per iteration to execute line 9 and leads to $`O(|\text{At}(P)|\times size(P))`$ running-time estimate for the alternating-fixpoint algorithm.

## 4 Procedure $`\mathrm{\Delta }_w`$

In this section we will focus on the class of programs, $`LP_1`$, that is, programs whose rules have no more than one positive atom in their bodies. We assume that we have a procedure false that, given a Horn program $`QLP_1`$, returns a subset of the set $`\text{At}(Q)LM(Q)`$. We also assume that false returns the empty set if and only if $`\text{At}(Q)=LM(Q)`$. For every pair $`(F,Q)`$, where $`F\text{At}(P)`$ and $`Q`$ is a Horn program such that $`\text{At}(Q)\overline{F}`$, we define

$$\mathrm{\Delta }_w(F,Q)=\text{false}(Q).$$

It is easy to see that this operator $`\mathrm{\Delta }_w(F,Q)`$ satisfies conditions (W1) and (W2). Consequently, it can be used in Algorithm 3. Clearly, the procedure $`\mathrm{\Delta }_w`$ and its computational properties are determined by the procedure false. In the remainder of the paper, we will describe a particular implementation of the procedure false and estimate its running time. We will use this estimate to obtain a bound on the running time of the resulting version of Algorithm 3.

A straightforward way to compute the least model of $`Q`$ and so, to find $`At(Q)LM(Q)`$, is ”bottom-up”. That is, we start with atoms which are heads of rules with the empty bodies and use the rules of $`Q`$ to compute all atoms in $`LM(Q)`$ by iterating the van Emden-Kowalski operator. An efficient implementation of the process is provided by the Dowling-Gallier algorithm \[Dowling & Gallier, 1984\].

The approach we follow here in the procedure false is ”top-down” and gives us, in general, only a part of the set $`At(Q)LM(Q)`$. More precisely, for an atom $`a`$ we proceed “backwards” attempting to construct a proof or to demonstrate that no proof exists. In the process, we either go back to an atom that is the head of a rule with empty body or we show that no proof exists. In the former case, $`aLM(Q)`$. In the latter one, none of the atoms considered while searching for a proof of $`a`$ are in $`LM(Q)`$ (because $`QLP_1`$ and each rule has at most one antecedent). The problem is that we may find an atom $`a`$ that does not have a proof only after we look at all other atoms first. Thus, in the worst case, finding one new false atom may require time that is proportional to the size of $`Q`$.

To improve the time performance, we look for proofs simultaneously for all atoms and grow the proofs “backwards” in a carefully controlled way. Namely, we never let one search to get too much ahead of the other searches. This controlled way of looking for proofs is the key idea of our approach and leads to a better performance. We will now provide an informal description of the procedure false followed later by a formal specification and an example.

In the procedure, we make use of a new atom, say $`s`$, different from all atoms occurring in $`Q`$. Further, we denote by $`head(r)`$ the atom in the head of a rule $`rQ`$ and by $`tail(r)`$ the atom which is either the unique positive atom in the body of $`r`$, if such an atom exists, or $`s`$ otherwise. We call an atom $`aAt(Q)`$ accessible if there are rules $`r_1,\mathrm{},r_k`$ in $`Q`$ such that $`tail(r_{i+1})=head(r_i)`$, for $`i=1,\mathrm{},k1`$, $`tail(r_1)=s`$ and $`head(r_k)=a`$. Clearly, the least model $`LM(Q)`$ of $`Q`$ is precisely the set of all accessible atoms.

In each step of the algorithm, the set of atoms from $`\text{At}(Q)`$ is partitioned into potentially false sets or pf-sets, for short. We say that a set $`vAt(Q)`$ is a pf-set if for each pair of distinct atoms $`a,bv`$ there are rules $`r_1,\mathrm{},r_k`$ in $`Q`$ such that $`tail(r_{i+1})=head(r_i)v`$, for $`i=1,\mathrm{},k1`$, $`tail(r_1)=b`$ and $`head(r_k)=a`$. It is clear that if $`v`$ is a pf-set then either all its elements are accessible (belong to the least model of $`Q`$) or none of them does (they are all false). Clearly, singleton sets consisting of individual atoms in $`At(Q)`$ are pf-sets. In the algorithm, with each pf-set we maintain its cardinality.

Current information about the state of all top-down searches and about the dependencies among atoms, that were discovered so far, is maintained in a directed graph $`G`$. The vertex set of this graph, say $`S`$, consists of $`\{s\}`$ and of a family of pf-sets forming a partition of the set $`\text{At}(Q)`$. The edges of $`G`$ are specified by a partial function $`\text{pred}:SS`$. We write $`\text{pred}(v)=𝐮ndefined`$ if pred is undefined for $`v`$. Thus, the set of edges of $`G`$ is given by $`\{(\text{pred}(v),v):\text{pred}(v)𝐮ndefined\}`$. Since pred is a partial function, it is easy to see that the connected components of the graph $`G`$ are unicyclic graphs or trees rooted in those vertices $`v`$ for which $`\text{pred}(v)`$ is undefined. Throughout the algorithm we always have $`\text{pred}(\{s\})=𝐮ndefined`$. Thus, the connected component of $`G`$ containing $`\{s\}`$ is always a tree and $`\{s\}`$ is its root.

If $`w`$ and $`v`$ are two different pf-sets, the existence of the edge $`(w,v)`$ in $`G`$ means that we have already discovered a rule in the original program whose head is in $`v`$ and whose tail is in $`w`$. Thus, if vertices in $`w`$ are accessible, then so are the vertices in $`v`$. A pf-set that is the root of a tree forming a component of $`G`$ is called an active pf-set. If $`v`$ is an active pf-set then no rule $`r`$ with $`head(r)v`$ and $`tail(r)v`$ has been detected so far. Thus, $`v`$ is a candidate for a set of atoms which does not intersect the least model of $`Q`$. Let us note that even though $`\{s\}`$ is a root of a tree in $`G`$ it is never active as it is not a pf-set in the first place.

We let active pf-sets grow by gluing them with other pf-sets. However, we allow to grow only these active pf-sets whose cardinalities are the least. In each iteration of the algorithm the value of the variable $`size`$ is a lower bound for the cardinalities of active pf-sets. To grow an active pf-set $`v`$, we look for rules with heads in $`v`$ and with tails in pf-sets other than $`v`$ (not necessarily active) or in $`\{s\}`$. The dependencies between pf-sets discovered in this way are represented as new directed edges in $`G`$. Pf-sets that appear in the same cycle are glued together (in the procedure cycle). Since $`\{s\}`$ is not an active pf-set, it never becomes an element of a cycle in $`G`$.

If, when attempting to grow a pf-set $`v`$ we discover a rule with head in $`v`$ and with the tail in a vertex of the tree of $`G`$ rooted in $`\{s\}`$, then $`v`$ is from now on ignored (all its vertices belong to the least model of $`Q`$). Indeed, $`v`$ gets connected to a tree of $`G`$ rooted in $`\{s\}`$. Consequently, it cannot become a member of a cycle in $`G`$ in the future and is never again considered by the procedure cycle.

The main loop (lines 6-23) of the algorithm false below starts by incrementing $`size`$ followed by a call to the procedure $`cycle(S,\text{pred},size,L)`$. This procedure scans the graph $`G`$ and identifies all its cycles. It then modifies $`G`$ by considering each cycle and by gluing its pf-sets into a single pf-set. To this end, it modifies the vertex set $`S`$ of $`G`$ and the function pred defining the edges of $`G`$. Each such new pf-set becomes the root of its tree in $`G`$ and so, it becomes active. The procedure cycle computes the cardinality of each new active pf-set. Finally, it creates a list $`L`$ so that it consists of active pf-sets of cardinality $`size`$. If no such set is found ($`L`$ is empty), we move on to the next iteration of the main loop and increment $`size`$ by 1. We give a more detailed description of the procedure $`cycle`$ later in the paper when we analyze the time complexity of our method.

For each active pf-set $`vL`$ we consider the tail of each rule with head in $`v`$ (lines 9-22). If there is a rule $`r`$ with $`head(r)v`$ and $`tail(r)v`$ then it is detected (line 15). The value $`\text{pred}(v)`$ is set to this element in $`S`$ that contains $`tail(r)`$ (it may be that this set is $`\{s\}`$). We also set the variable $`success`$ to true (line 16). The pf-set $`v`$ stops to be active. We move on to the next active pf-set on $`L`$.

If such a rule $`r`$ does not exist then $`success=𝐟alse`$ and $`v`$ is a set of cardinality $`size`$ consisting of atoms which are not in the least model of $`Q`$. This set is returned by the procedure false (line 21). Hence, for an active pf-set considered in the loop 6-23, either we find a pf-set $`\text{pred}(v)S\{v\}`$ (and we have to consider the next pf-set on $`L`$) or $`v`$ is returned as a set of atoms which are not in the least model of $`Q`$ (and the procedure false terminates). Thus, the procedure false is completed if either a nonempty set $`v`$ of atoms which are not in the least model of $`Q`$ is found or, after some passes of the loop 6-23, the graph $`G`$ has no active pf-sets. In the latter case $`G`$ is a tree with the root in $`\{s\}`$. Thus, $`At(Q)=LM(Q)`$ and $`v=\mathrm{}`$ is returned (line 24).

In the procedure false, as formally described below, an input program $`Q`$ is represented by lists $`IN(a)`$, $`aAt(Q)`$, of all atoms $`b`$ such that $`b`$ is the body of some rule with the head $`a`$. If there is a rule with the head $`a`$ and empty body, we insert $`s`$ into the list $`IN(a)`$.

We also use an operation $`next`$ on lists and elements. Let $`l`$ be a list and $`w`$ be an element, either belonging to $`l`$ or having a special value undefined. Then

$$next(w,l)=\{\begin{array}{cc}\mathrm{t}henextelementafterw\mathrm{i}nl\hfill & ifwl\hfill \\ \mathrm{t}hefirstelementinl\hfill & ifw\mathrm{i}s𝐮ndefined.\hfill \end{array}$$

The value undefined should not be mixed with nil which indicates the end of a list.

Finally, we use a procedure $`\text{findset}(w,S)`$ which, for an atom $`w`$ and a collection $`S`$ of disjoint sets, one of which contains $`w`$, finds the name of the set in $`S`$ containing $`w`$ (it follows from our assumptions that such a set is unique). Elements of $`S`$ are maintained as linked lists. Each element on such a list has a pointer to the head of the list. The head serves as the identifier for the list. When the procedure $`\text{findset}(w,S)`$ is called, it returns the head of the list to which $`w`$ belongs.

|  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 1 | procedure $`\text{false}(Q)`$; |
| 2 |  | $`S:=\{\{x\}:xAt(Q)\}\{\{s\}\}`$; |
| 3 |  | for $`vS`$ do $`\text{pred}(v):=𝐮ndefined`$; |
| 4 |  | for $`xAt(Q)`$ do {$`w(x):=𝐮ndefined`$; $`cardinality(x):=1`$}; |
| 5 |  | $`size:=0`$; |
| 6 |  | while $`size<|At(Q)|`$ do |
| 7 |  |  | $`\{size:=size+1`$; |
| 8 |  |  | $`cycle(S,\text{pred},size,L)`$; |
| 9 |  |  | for all $`vL`$ do |
| 10 |  |  |  | $`\{success:=𝐟alse`$; |
| 11 |  |  |  | $`u:=next(u,v)`$; |
| 12 |  |  |  | while $`u`$ nil and not $`success`$ do |
| 13 |  |  |  |  | $`w(u):=next(w(u),IN(u))`$; |
| 14 |  |  |  |  | while $`w(u)`$ nil and not $`success`$ do |
| 15 |  |  |  |  |  | $`\{`$if $`\text{findset}(w(u),S)v`$ |
| 16 |  |  |  |  |  |  | then $`\{success:=𝐭rue`$; $`\text{pred}(v):=\text{findset}(w(u),S)\}`$ |
| 17 |  |  |  |  |  |  | else $`w(u):=next(w(u),IN(u))`$ |
| 18 |  |  |  |  | end while (14)$`\}`$; |
| 19 |  |  |  |  | if not $`success`$ then $`u:=next(u,v)`$ |
| 20 |  |  |  | end while (12)$`\}`$; |
| 21 |  |  |  | if not $`success`$ then return $`v`$ (\* the procedure terminates \*) |
| 22 |  |  | end for (9)$`\}`$ |
| 23 |  | end while (6)$`\}`$; |
| 24 |  | return $`v=\mathrm{}`$ |
| 25 | end false; |

We will now illustrate the operation of the algorithm. Let us consider the following Horn logic program $`Q`$:

$$\begin{array}{cccccc}a\hfill & ba\hfill & ac\hfill & ca\hfill & ae\hfill & de\hfill \\ fd\hfill & ef\hfill & df\hfill & eg\hfill & gj\hfill & jg\hfill \\ ij\hfill & jh\hfill & kj\hfill & kh\hfill & hk\hfill & \end{array}$$

This program is represented as a graph, $`G^Q`$, in Fig. 1. The vertices of this graph correspond to the atoms of the program. In addition, $`G^Q`$ has an auxiliary vertex $`s\text{At}(Q)`$. An edge $`(x,y)`$, where $`x,y\text{At}(Q)`$, represents the clause $`yx`$ from $`Q`$. An edge $`(s,y)`$, where $`y\text{At}(Q)`$, represents the clause $`y`$. When illustrating the algorithm, we assume that atoms from $`\text{At}(Q)`$ (atoms $`a,\mathrm{},k`$ in our example) appear on the lists $`IN(x)`$, $`x\text{At}(Q)`$, in the alphabetical order. We also assume that whenever $`s`$ belongs to a list $`IN(x)`$, it appears as the first atom on the list.

In the algorithm false, the current state of knowledge about the possibility of proving an atom from $`Q`$ is represented by the graph $`G`$. Initially, $`G`$ consists of isolated vertices. Indeed, line 3 of the algorithm sets $`\text{pred}(x)`$ to undefined, for every vertex $`x`$ of $`S`$ (see Fig. 2 (left)). All of the vertices of $`G`$, except for $`\{s\}`$ are active pf-sets. The procedure $`cycle`$ (line 8), called with $`size=1`$, puts all of them on the list $`L`$.

The algorithm considers next (line 9) all elements on the list $`L`$, that is, all vertices of $`G`$ that are active pf-sets and have cardinality equal to size. During the first iteration of the loop 6-23, $`L`$ consists of all vertices of $`G`$, except for $`\{s\}`$ (that is, singleton sets $`\{x\}`$, where $`x\text{At}(Q)=V(G^Q)\{s\}`$). For each vertex $`v`$ of $`G`$ on $`L`$, the algorithm looks for a back rule for $`v`$, that is, a rule in $`Q`$ with the head in $`v`$ and the tail in a pf-set other than $`v`$ or in $`\{s\}`$. In our graphical representation of $`Q`$ by means of the graph $`G^Q`$, a back rule for $`v`$ corresponds to an edge (referred to as a back edge) in $`G^Q`$ with the head in $`v`$ and the tail in a vertex of $`G`$ other than $`v`$ (possibly in $`\{s\}`$). To find a back rule (edge) for $`v`$, all atoms $`u`$ of $`Q`$ (equivalently, all vertices $`u`$ of $`G^Q`$) that belong to $`v`$ are considered (the loop 12-20). For each such atom $`u`$, the algorithm searches for the first atom on the list $`IN(u)`$ that does not belong to $`v`$. Let us recall that $`IN(u)`$ is the list of atoms that are the tails of rules with the head $`u`$ or, in the terms of the graph $`G^Q`$, that are the tails of edges with the head $`u`$. If such an atom is found, together with $`u`$ it determines a back rule (edge) $`r`$ for $`v`$. The algorithm sets $`\text{pred}(v)`$ to be equal to the pf-set containing the tail of $`r`$ (line 16). That is, an edge from $`\text{pred}(v)`$ to $`v`$ is added to $`G`$. The algorithm moves then on to the next element of the list $`L`$.

In our example, in the first iteration of the loop 6-23, a back rule is found for every element on $`L`$, that is, for every vertex of $`G`$ other than $`\{s\}`$. For instance, for the vertex $`\{d\}`$, the algorithm considers atoms on the list $`IN(d)=(e,f)`$ (let us recall that atoms on lists $`IN(x)`$ are arranged alphabetically with the exception of the special atom $`s`$ which, if present on a list, is always its first element). The first atom on the list, $`e`$ does not belong to $`\{d\}`$. Thus, it defines, together with $`d`$ a back rule for $`\{d\}`$, $`de`$. The resulting graph $`G`$ is shown in Fig. 2 on the right.

Let us note that when scanning the list $`IN(d)`$ in subsequent iterations the algorithm resumes the scan with the first atom that has not been looked at yet (cf. the definition of the operation next). Thus, the next time $`d`$ is considered as an element of an active pf-set for which a back rule is searched for, the scan of $`IN(d)`$ will start with $`f`$. The same holds true for all lists $`IN(x)`$, $`x\text{At}(Q)`$. Consequently, each atom on each of these lists is considered just once. Such an approach still guarantees that finding back rules works correctly (that is, that they are found by the algorithm whenever they exist). Indeed, when an atom on a list $`IN(x)`$ is considered, it either defines a back rule with the head $`x`$ (and, thus, cannot define any new back rule with the head $`x`$ in the future) or it is in the same active pf-set as $`x`$ (and, thus, it neither defines a back rule now nor it will define it in the future, as it will remain in the same pf-set as $`x`$ till the algorithm terminates).

The second iteration of the loop 6-23 starts with the procedure cycle contracting each cycle in the graph $`G`$ to a single vertex. The resulting graph is shown in Fig. 3 on the left. The procedure cycle then creates a new list $`L`$. It consists of all active pf-sets of cardinality 2. In our case, $`L`$ contains $`\{g,j\}`$ and $`\{h,k\}`$ ($`\{d,e,f\}`$ is also active but has cardinality 3).

Continuing with the second iteration, the algorithm next considers each vertex on $`L`$ (the loop 9-22) and looks for back rules. In this iteration, a back rule is found for each of the nodes on $`L`$ and the modified graph $`G`$ is given in Fig. 3 on the right.

In the third iteration, the procedure cycle contracts the only cycle in $`G`$ to a single active pf-set of cardinality 4 (Figure 4, left side). It also creates a new list $`L`$. This time it consists of active pf-sets of cardinality 3. There is just one such set - $`\{d,e,f\}`$. Subsequently, the algorithm false looks for a back rule for $`\{d,e,f\}`$. It starts by considering edges ending in $`d`$ (line 11; we assume that $`v`$ is represented by the list $`(d,e,f)`$). It scans the list $`IN(d)`$ starting at the first atom that has not been inspected so far, that is, $`f`$. However, since $`f`$ belongs to the same pf-set as $`d`$, $`f`$ does not specify a back rule. Since there are no more atoms on the list $`IN(d)`$, we move on to the next iteration of the loop 12-20 and consider atom $`e`$. We have $`IN(e)=(f,g)`$. Since $`f`$ was already considered (and yielded a back rule for $`\{e\}`$) in the first iteration, we consider $`g`$. Since $`g\{d,e,f\}`$, it defines a back rule for $`\{d,e,f\}`$, $`dg`$.

The resulting graph $`G`$ is shown in Figure 4 (on the right). It has no cycles. So, the only thing done by the procedure cycle in the iteration 4 is that it puts on $`L`$ active pf-sets of cardinality 4. There is just one such set in $`G`$, $`\{g,h,j,k\}`$. The algorithm false looks for a back edge for $`\{g,h,j,k\}`$ and does not find any. The variable success remains false. The algorithm returns $`\{g,h,j,k\}`$ and terminates (line 21). Let us note that this set is a proper subset of the set $`\text{At}(Q)LM(Q)`$.

The following theorem formally establishes two key properties of the procedure false.

###### Theorem 4.1

1. The procedure false returns a set $`v`$ such that $`vAt(Q)LM(Q)`$.
2. false returns the empty set if and only if $`At(Q)LM(Q)=\mathrm{}`$.

Proof: (1) The statement is trivially true if false returns the empty set. Thus assume that the returned set $`v\mathrm{}`$. It means that the value of the variable $`success`$ is false after all passes of the loop 12-20 for some active pf-set $`v`$ in the list $`L`$. Thus every rule in $`Q`$ with the head in $`v`$ has been considered.

Suppose there is a rule $`r`$ in $`Q`$ with $`head(r)=uv`$ and $`tail(r)=bv`$. This rule was considered by the procedure false when $`u=head(r)`$ was a member of some active pf-set, say $`y`$. Since larger pf-sets are obtained by gluing smaller ones, $`yv`$. While $`r`$ was being considered, the value of $`w(u)`$ in the loop 14-18 was $`b`$ and the value of $`v`$ was $`y`$. Consequently, $`\text{findset}(b,S)y`$ in line 15 because $`yv`$ and $`bv`$ so $`by`$. Hence the value of $`success`$ was set to true and $`\text{pred}(y)`$ was defined to be, say, $`z=\text{findset}(b,S)`$ in line 16. The pf-set $`y`$ stopped to be active. Recall that $`v`$ is active when the procedure stops. Hence $`y`$ had to be glued with other pf-sets to obtain $`v`$. This is, however, impossible because if $`y`$ were glued with some other pf-sets to form a larger pf-set $`x`$ then $`\text{pred}(y)=zx`$. Notice that $`bzxv`$. We have got a contradiction with $`bv`$.

Hence, there are no rules $`r`$ in $`Q`$ with $`head(r)v`$ and $`tail(r)v`$. Thus no atom in $`v`$ is accessible so $`vAt(Q)LM(Q)`$.
(2) Suppose false returns the empty set and consider the last pass of the loop 6-23, for $`size=|At(Q)|`$. If the list $`L`$ is empty then no vertex of $`G`$ is an active pf-set. Hence, $`G`$ is a tree with the root $`\{s\}`$. Thus all atoms in $`At(Q)`$ are accessible and consequently $`LM(Q)=At(Q)`$.

If the list $`L`$ is nonempty then it contains one pf-set $`v=At(Q)`$. The empty set is returned by the procedure false so the value of the variable $`success`$ in line 16 is true for $`v=At(Q)`$. It means that for some rule $`r`$ in $`Q`$ with $`head(r)=u`$, $`w(u)=tail(r)v=At(Q)`$ so $`w(u)=s`$. Hence, $`u`$ is accessible and, consequently, all atoms in $`At(Q)`$ are accessible. That is, we have $`At(Q)LM(Q)=\mathrm{}`$.

The converse of the implication proved above follows immediately from the first part of the theorem. $`\mathrm{}`$

We shall now consider the procedure $`cycle`$ a little bit more carefully. The procedure can be informally written in the following form.

|  |  |  |  |
| --- | --- | --- | --- |
| procedure $`cycle(S,\text{pred},size,L)`$ |
| 1. | Initialize $`L`$ to empty. |
| 2. | Find all cycles $`C_1,C_2,\mathrm{},C_p`$ in the graph $`G`$. Put $`C=\{C_1,C_2,\mathrm{},C_p\}`$. |
| 3. | For every cycle $`C=\{v_1,\mathrm{},v_q\}`$, $`CC`$, do (i)-(iv). |
|  | (i) | set $`v_C:=v_1\mathrm{}v_q`$; |
|  | (ii) | compute $`cardinality(v_C)`$ (sum up the cardinalities of all vertices in $`C`$); |
|  | (iii) | update the set $`S`$: set $`S:=(S\{v_1,\mathrm{},v_q\})\{v_C\}`$; (\* $`v_C`$ becomes |
|  |  | an active pf-set \*) |
|  | (iv) | update the function pred: for every $`i=1,\mathrm{},q`$, if $`\text{pred}(z)=v_i`$ (for |
|  |  | some $`zS`$) then $`\text{pred}(z):=v_C`$; |
| 4. | For every vertex of $`G`$ that is an active pf-set, if $`cardinality(v)=size`$, insert $`v`$ |
|  | into the list $`L`$. |

Since $`G`$ is a directed graph whose connected components are either unicyclic graphs or trees, step 2 of the procedure $`cycle`$ can be implemented in $`O(|S|)`$ time. Since pf-sets are represented as linked lists, with each node on the list pointing to the head of the list, step (i) can be implemented to take $`O(|v_C|)`$ steps. The time needed for step (ii) is, clearly, $`O(|C|)`$. Each execution of step (iii) takes also $`O(|C|)`$. Finally, the running time of each execution of step (iv) is $`O(m_C)`$, where $`m_C`$ is the size of the connected component of the graph $`G`$ containing $`C`$. Thus, an iteration of the loop 3 for a cycle $`CC`$ takes $`O(|C|+m_C+|v_C|)`$. Clearly, $`|C|m_C`$. Moreover, $`_{CC}m_C|S|1|\text{At}(Q)|`$ and $`_{CC}|v_C||\text{At}(Q)|`$ (they are all disjoint subsets of $`\text{At}(Q)`$). Thus, the total time needed for the loop 3 is $`O(|\text{At}(Q)|)`$. It is easy to see that the time needed for the loop 4 is also $`O(|\text{At}(Q)|)`$. Consequently, the running time of the procedure $`cycle`$ is $`O(|\text{At}(Q)|)`$.

We are now in a position to estimate the running time of the procedure false.

###### Lemma 4.2

If the procedure $`\text{false}(Q)`$ returns a nonempty set $`v`$, then the running time of false is $`O(|v|\times |\text{At}(Q)|)`$. If $`\text{false}(Q)`$ returns the empty set then its running time is $`O(|\text{At}(Q)|^2)`$.

Proof: Let $`|\text{At}(Q)|=n`$ and $`|v|=k`$. As we have already observed the procedure $`cycle`$ runs in time $`O(n)`$. It is not hard to see that, since we represent all sets occurring in the procedure false as linked lists, with each node on a list pointing to the head of the list, the operations: findset and $`next`$ require a constant time.

First assume that the output $`v`$ of the procedure false is nonempty. Let us estimate the number of passes of the while and for loops in the procedure. Clearly, the loop 6-23 is executed $`k`$ times. Hence the total running time of all calls of the procedure $`cycle`$ is $`O(kn)`$. The number of passes of the loop 9-22 is not larger than $`|L_1|+|L_2|+\mathrm{}+|L_k|`$, where $`L_i`$ denotes the list $`L`$ in an iteration $`i`$ of the loop. Since $`L_i`$ is a list of disjoint pf-sets of cardinality $`i`$, $`|L_i|n`$, for each $`i=1,2,\mathrm{},k`$. Hence the number of passes of the loop 9-22 can be very roughly estimated by $`kn`$. The loop 12-20 is executed at most

$$\underset{i=1}{\overset{k}{}}\underset{vL_i}{}|v|kn$$

times. This inequality follows from the fact that the sets $`v`$ in the lists $`L_i`$ are disjoint subsets of atoms so $`_{vL_i}|v|n`$. The estimation of the number of passes of the loop 14-18 is a little bit more complicated. First notice that in each execution of the loop we check a rule of the program $`Q`$ and rules are checked only one time. The rules $`r`$ checked in the loop have either both the head and the tail in some pf-set $`vS`$ or $`head(r)v`$ and $`tail(r)`$ is in some other pf-set $`uS`$. In the latter case $`\text{pred}(v)`$ is defined in line 16. The number of executions of line 16 is not larger than the number of passes of the loop 9-22 so it is bounded by $`kn`$. When the procedure returns the output, the pf-sets have cardinalities not larger than $`k`$. Hence the number of rules with both the head and the tail in the same pf-set that has been checked before the procedure stops is not larger than

$$\underset{uS}{}|u|(|u|1)(k1)\underset{uS}{}|u|(k1)n.$$

Thus the number of passes of the loop 14-18 in the whole procedure false is less than $`2kn`$. It follows that if the output $`v`$ of false is nonempty then the running time of false is $`O(|v|\times |\text{At}(Q)|)`$.

Now consider the case when the procedure false returns the empty set. Clearly the number of passes of the loop 6-23 is $`n`$ so it takes $`O(n^2)`$ time for all executions of the procedure $`cycle`$. Since the rules are checked in the loop 14-18 only one time, the number of passes of this loop is not larger than the number $`m`$ of rules in $`Q`$. Obviously $`mn^2`$ so the running time of false in this case is $`O(|\text{At}(Q)|^2)`$. $`\mathrm{}`$

By Lemma 4.2 and considerations in Section 3 we get an estimation of the running time of Algorithm 3.

###### Theorem 4.3

If $`P`$ is a program whose rules have at most one positive atom in the body then Algorithm 3 can be implemented so that its running time is $`O(|At(P)|^2+size(P))`$. $`\mathrm{}`$

## 5 Conclusions

The method for computing the well-founded semantics described in this paper is a refinement of the basic alternating-fixpoint algorithm. The key idea is to use a top-down search when identifying atoms that are false. Our method is designed to work with programs whose rules have at most one positive atom in their bodies (class $`LP_1`$). Its running time is $`O(|\text{At}(P)|^2+size(P))`$ (where $`P`$ is an input program). Thus, our algorithm is an improvement over other known methods to compute the well-founded semantics for programs in the class $`LP_1`$. Our algorithm runs in linear time for the class of programs $`PLP_1`$ for which $`size(P)|\text{At}(P)|^2`$. However, it is not a linear-time algorithm in general. It is an open question whether a linear-time algorithm for computing the well-founded semantics for programs in the class $`LP_1`$ exists.

Our results extend to the class $`LP_1^+`$. However, the extension is straightforward and the class $`LP_1^+`$ is still rather narrow. Moreover, it is not specified syntactically (it is described by means of the Kripke-Kleene semantics). The question arises whether our top-down approach to positive-loop detection can be generalized to any class of programs significantly extending the class $`LP_1`$ and possessing a simple syntactic description.

Finally, let us note that the general problem of computing the well-founded semantics still remains a challenge. No significant improvement over the alternating-fixpoint algorithm of Van Gelder has been obtained for the class of arbitrary finite propositional logic programs.

## Acknowledgments

This research was supported by the NSF grants CDA-9502645 and IRI-9619233.


text:
# Phenomenological theory of magnetic anisotropy in thin films and multilayers

## Abstract

Within a phenomenological approach the density of induced uniaxial anisotropy in nanostructures $`K(𝐫)`$ is treated as a physical field additional to the magnetization $`𝐌(𝐫)`$. The equilibrium distributions of $`K(𝐫)`$ are formed under the influence of the surfaces and internal interactions. The functions $`K(𝐫)`$ and $`𝐌(𝐫)`$ are calculated for a magnetic layer between two nonmagnetic spacers. It is shown that the transition from the phase with magnetization in plane to the perpendicular phase occurs continuously via an intermediate phase which is inhomogeneous across the layer-thickness. The theory explains experimentally observed thickness dependences of the effective anisotropy and gives a method to obtain values of the characteristic parameters for the theory from experimental data. It is shown that the separation of the induced anisotropy into volume and surface contribution is valid only for thick films.

Magnetic anisotropy is one of the key properties exploited in ferromagnetic nanostructures. Reduced dimensionality may greatly alter the effective magnetic anisotropy of such structures. For layered systems, numerous investigations demonstrate that reduced thickness and complex interactions on surfaces and interfaces induce uniaxial magnetic anisotropy which, by orders of magnitude, exceeds the values of the intrinsic anisotropy in the corresponding bulk materials . During the last decade a large variety of magnetic thin films and multilayers have been synthesized and investigated experimentally. They have already found applications in modern magnetoelectronics such as spin-valves or magnetoresistive heads .

Recent first principle numerical calculations have proved the reliability of the quantum-mechanical microscopic theory in its applications to magnetic nanostructures . They are offering considerable insight into the phenomena of induced anisotropy. However, these ab initio calculations are still unable to give a complete description of the magnetization structures and processes in real layered systems. Hence, our understanding of and control over the magnetic anisotropy of nanostructures is rather incomplete. Up to now, the analysis of experimental data on anisotropy effects in magnetic layered systems is mostly based on the effective volume anisotropy method introduced by Néel about fifty years ago . Within this approach the value of average induced uniaxial anisotropy in a magnetic layer of thickness $`d`$ is given by an empirical ansatz

$$K_{\text{eff}}=K_V+\frac{K_S}{d}$$
(1)
where $`K_V`$ is the volume anisotropy (per unit volume), and $`K_S`$ is the surface (interface) contribution (per unit area). For many layered systems with reduced thickness, the function $`K(1/d`$) becomes strongly nonlinear which is at variance with eq. (1) . Moreover, the general validity of the Néel theory has been questioned in a number of publications . The separation of the uniaxial anisotropy into volume and surface contribution and the reduction of surface anisotropy into an effective volume contribution should be considered as an excessive simplification. Furthermore, equation (1) implies that the anisotropy is constant within a magnetic layer and, thus, stabilizes a homogeneous distribution of the magnetization. It was shown, however, that the competition between volume and surface anisotropies may induce inhomogeneous states in magnetic films and multilayer structures . The surface induced anisotropy can also gradually relax into the depth of the layer . On the whole, the heuristic Néel approach as well as similar models have succeeded to give qualitative explanations of some general features of the induced anisotropy. However, they do not account for many other observed effects and they fail to give a quantitative description of the magnetization structures in magnetic nanostructures.

In this paper we propose a simple phenomenological theory which gives a consistent method to calculate the anisotropy and the equilibrium magnetization in thin films and multilayers. Let us consider a magnetic nanostructure embedded into nonmagnetic media. It may be a thin magnetic film confined by the interface with the substrate and the surface, or a magnetic layer between nonmagnetic spacers in a multilayer, or any other magnetic inclusion in a nonmagnetic matrix (such as nanowires, magnetic dots, or magnetic clusters). The magnetic energy of such a nanostructure can be written as an integral over its volume

$$W_m=\left[A\underset{i}{}\left(\frac{𝐦}{x_i}\right)^2+K(𝐫)(𝐦𝐧)^2𝐇𝐌\frac{1}{2}𝐇_d𝐌\right]𝑑𝐫.$$
(2)
Here $`𝐦=𝐌/M_0`$ is the normalized value of the magnetization vector $`𝐌`$ ($`M_0=\left|𝐌\right|`$), $`𝐧`$ is a unity vector perpendicular to the surface of the layer, $`A`$ is the exchange stiffness constant, $`𝐇`$ is an external magnetic field, and $`𝐇_d`$ is a demagnetizing field. The second term in (2) is the energy density of the induced uniaxial anisotropy. To describe the distribution of $`K`$($`𝐫`$) within the nanostructure we introduce an interaction functional of Landau-Ginzburg type

$$W_A=\left[\alpha \underset{i}{}\left(\frac{K}{x_i}\right)^2+f(K)\right]𝑑𝐫.$$
(3)
The first term in (2) represents the stiffness energy and $`f`$($`K`$) is the energy of a homogeneous induced anisotropy. In many cases the function $`f`$($`K`$) may be written as $`f(K)=aK^2+2bK`$. The stiffness parameter $`\alpha >0`$ together with $`a>0`$ characterize the resistance of the system against the uniaxial anisotropy imposed by the boundaries. The linear term in $`f`$($`K`$) reflects the tendency for symmetry breaking and the ensuing rise of homogeneous uniaxial anisotropy. Minimization of $`f`$($`K`$) yields the value of this anisotropy $`K_b=b/a`$ . The scalar function $`K(𝐫)`$ may be treated as a physical field additional to the magnetization field $`𝐌(𝐫)`$. Both fields are coupled by the anisotropy energy in (2). In most cases, this coupling energy should be negligible compared to the strong surface and volume interactions inducing the uniaxial anisotropy. Thus, the equilibrium distribution $`K`$($`𝐫`$) can be calculated independently from the magnetization field by minimization of the functional (3). The anisotropy on the confining surfaces is supposed to have values $`K(𝐫)_S=K_0`$ fixing the boundary conditions for the variation problem. Then, the equilibrium distributions of the magnetization are determined by variation of the functional (2) with this definite function $`K`$($`𝐫`$). Here, the boundary conditions $`K_0`$ are kept constant. But, in general, $`K_0`$ may vary within the surfaces, and all phenomenological constants in (3) may depend on characteristic sizes of the nanostructure. The values $`K_0`$ (volume anisotropy density on the surfaces) should not be confused with the “surface anisotropy” $`K_S`$ in (1) which is an anisotropy energy per area. The relation between these quantities will be discussed below.

As an application of the theory we calculate the induced anisotropy and the magnetization in a magnetic layer of thickness $`d`$ sandwiched between two identical nonmagnetic spacers. The layer is supposed to be infinite in $`x`$ and $`y`$ directions and bounded by parallel planar surfaces $`z=\pm d/2`$. On the surfaces the anisotropy is $`K(d/2)=K(d/2)=K_0`$. The vector $`𝐧`$ is directed along the $`z`$-axis and $`K_0>0`$ (perpendicular, surface induced anisotropy). Within the magnetic plate $`K`$($`𝐫`$) varies along $`z`$-direction and is homogeneous in ($`x,y`$) plane. Variation of the functional (3) yields the following solution for $`K`$($`z`$):

$$K(z)=K_b+\frac{\left(K_0K_b\right)\mathrm{cosh}\left(z/\lambda _a\right)}{\mathrm{cosh}\left(d/(2\lambda _a)\right)}$$
(4)
where $`\lambda _a=\sqrt{\alpha /a}`$ is a characteristic length. To derive the equilibrium functions $`𝐌(𝐫`$) one has to minimize the functional for the magnetic energy (2) with $`K`$($`z`$) given by (4). We consider structures that are homogeneous in the plate plane (i.e. monodomain structures). Further, we ignore the variation of $`M_0`$ observed in some multilayers . Under these assumptions the magnetization in the layer is described by the polar angle $`\theta (z)`$ between the vector $`𝐌`$ and $`z`$-axis. For such a one-dimensional problem the magnetostatic equations for $`𝐇_d`$ have a rigorous solution, and the corresponding stray field energy is expressed by an anisotropy-like term $`w_d=𝐇_d𝐌/2=2\pi M_0^2\mathrm{cos}^2(\theta )`$ . This energy contribution is commonly named “shape anisotropy”. Thus, the functional (2) with $`K(z)`$ given by (4) may be expressed in the following reduced form

$$W=2\pi M_0^2_{d/2}^{d/2}\left\{\lambda _{ex}^2\left(\frac{d\theta }{dz}\right)^2\frac{\stackrel{~}{K}(z)}{2\pi M_0^2}\mathrm{cos}^2\theta h\mathrm{cos}(\theta \psi )\right\}𝑑z$$
(5)
where $`h=H/(2\pi M_0)`$ is a reduced value of the magnetic field and $`\psi `$ is the angle between the vector $`𝐇`$ and $`z`$-axis, $`\lambda _{ex}=\sqrt{A/\left(2\pi M_0^2\right)}`$ is the exchange length, and the volume density of the total anisotropy $`\stackrel{~}{K}(z)`$ is

$$\stackrel{~}{K}(z)=2\pi M_0^2+\frac{K_0\mathrm{cosh}\left(z/\lambda _a\right)}{\mathrm{cosh}\left(d/(2\lambda _a)\right)}.$$
(6)
For simplicity we suppose that the homogeneous part of the induced anisotropy $`K_b`$ (see (4)) is included into the coefficients $`K_0`$ and $`2\pi M_0^2`$.

The energy (5) has the functional form of an uniaxial ferromagnet with anisotropy $`\stackrel{~}{K}(z)`$ (6) varying along the layer thickness. The distribution of $`\stackrel{~}{K}(z)`$ within the layer strongly depends on the layer thickness (Fig. 1). The length $`\lambda _a`$ characterizes the extension of the surface anisotropy into the layer depth and may be named “penetration depth of the induced anisotropy”. When $`d\lambda _a`$ the induced anisotropy is localized near surfaces and rapidly relaxes into the depth of the layer. For thinner layers, where $`\lambda _a`$ and $`d`$ have the same order of magnitude, the induced anisotropy attains considerable values within all the layer (Fig. 1). The possible magnetic phases and the regions of their stability are determined by optimization of the energy (5) with free boundary conditions that read for this case $`\left(d\theta /dz\right)_{z=\pm d/2}=0.`$ A complete analysis of the functional (5) will be done elsewhere. Here we only discuss some general properties of the solutions at zero field as representative for the problem. In this case the layer thickness $`d`$ and the parameters $`q=K_0/\left(2\pi M_0^2\right)`$ and $`\Lambda =\lambda _{ex}/\lambda _a`$ span the phase space of the solutions. A typical ($`q`$, $`d`$) phase diagram for fixed values of $`\Lambda `$ is shown in Fig. 2. It consists of two regions where homogeneous phases exist with either perpendicular ($`\theta =0`$) or parallel  ($`\theta =\pi /2`$) magnetization. These are separated by an inhomogeneous (twisted) phase where the angle $`\theta `$ gradually increases from a definite value $`\theta _10`$ on the surfaces to a largest value in the center $`\theta _1<\theta _2\pi /2`$.

The perpendicular phase is stable for small values of the thickness and when $`q>1`$. The parallel phase exists only in rather thick layers and when the parameter $`q`$ is smaller than a critical value $`q_c(\Lambda )`$. In the region $`q>q_c`$ the perpendicular phase transforms continuously into the twisted phase at the critical line $`d_1(q)`$. For $`1<q<q_c`$ the region of the existence of the twisted phase is bound by two lines of second order transitions into the homogeneous phases $`d_1(q)`$ and $`d_2(q)`$. Thus, according to our model the transition between the parallel and the perpendicular phases always occurs continuously via the transformation into an intermediary (twisted) phase. Fig. 3 shows the evolution of the magnetic structures during such a transition.

By integrating $`\stackrel{~}{K}(z)`$ in (6) over $`z`$ one obtains for the average value of the uniaxial anisotropy density

$$K_{\text{eff}}=\frac{1}{d}_{d/2}^{d/2}\stackrel{~}{K}(z)𝑑z=2\pi M_0^2+\frac{K_0\mathrm{tanh}\left(d/(2\lambda _a)\right)}{d/(2\lambda _a)}.$$
(7)
Values of the anisotropy $`K_{\text{eff}}`$ have been measured in many layered systems . Usually, experimental values of $`K_{\text{eff}}`$ are plotted as product $`K_{\text{eff}}d`$ versus $`d`$. In Fig. 4 the thickness dependence of the function $`\Phi (d)=K_{\text{eff}}d`$ corresponding to (7) is shown for different values $`q`$.

In the limit of large thickness ($`d\lambda _a)`$ the functions $`\Phi (d)`$ are negative and depend linearly on the thickness (with slope equal to the value of the shape anisotropy). In thiner layers ($`d\lambda _a`$) the functions $`\Phi (d)`$ become nonlinear. For $`q<1`$ they are always negative (parallel anisotropy) and reach zero as the thickness decreases to zero. For $`q>1`$, the function $`\Phi (d)`$ changes sign at the critical thickness $`d_0`$ determined by the equation

$$K_0\mathrm{tanh}\left(\frac{d_0}{2\lambda _a}\right)=\frac{2\pi M_0^2d_0}{2\lambda _a}$$
(8)
and remains positive for $`0<d<d_0`$ manifesting the existence of the perpendicular anisotropy in this region. The function $`\Phi (d)`$ reaches the maximum value for the thickness

$$d_{\mathrm{max}}=2\lambda _a\text{arc}\mathrm{cosh}\sqrt{\frac{K_0}{2\pi M_0^2},}$$
(9)
and then monotonically decreases to zero as $`d0`$ (Fig. 4). Experimentally obtained functions $`\Phi (d)`$ are in close accordance with our theoretical results. E.g., nonmonotonic dependencies of $`\Phi (d)`$ have been observed for Co/Au multilayers, Co/Ir multilayers, or for Ni/Cu multilayers and Cu/Ni/Cu sandwiches . Taking experimental values for $`d_0`$ and $`d_{\mathrm{max}}`$ it is possible to calculate the parameters $`q`$ and $`\lambda _a`$ from eqs. (8) and (9). E.g., for Co/Au  films this yields $`q=4.6`$ and $`\lambda _a=1.9`$Å. For Ni/Cu multilayer systems the function$`\Phi (d)`$ in the regions of perpendicular anisotropy has been obtained by three different experimental methods . The results from these experiments show similar functional dependencies for $`\Phi (d)`$. However, the measured values display considerable quantitative differences. Taking the experimental data from  and  we derive the following values for the characteristic parameters $`q=2.1`$ and $`2.2`$, $`\lambda _a=26.4`$Å and $`31.9`$Å. In contrast to the full curves from our theory in Fig. 4 experimental dependencies $`\Phi (d)`$ for these Ni-films become negative again for finite values of the thickness manifesting the reentrance of planar anisotropy in ultrathin films . It should be noted in this respect that for the calculation of $`\Phi (d)`$ (Fig. 4) the phenomenological parameters $`K_0`$ are kept constant in this simple approach. According to experimental observations they may depend on the layer thickness, and they may be sensitive to details of the film preparation, e.g. via surface reconstructions, misfit strain or misfit dislocations etc. In particular, it was established experimentally that in the Ni-films under discussion surface anisotropy decreases with decrease of the layer thickness due to magnetoelastic interactions . For decreasing values of $`K_0`$ the theoretical dependencies $`\Phi (d)`$ become negative in the limit of very thin films (dashed line in Fig. 4). This thickness dependency of the phenomenological parameters $`K_0`$ as well as other phenomenological constants should be considered by introducing appropriate models in future work based on further experimental results and results from ab initio calculations.

The distribution of the induced anisotropy and the magnetization in a magnetic layer can also be investigated within a discretized model which is useful in the limit of ultrathin layers or for comparison with ab initio calculations. In this case the equilibrium values of the induced anisotropy in the magnetic layer consisting of $`N`$ magnetic planes are determined by minimization of the energy

$$\stackrel{~}{W}_A=\underset{n=1}{\overset{N1}{}}\left[\left(1+\frac{2\lambda _a^2}{\mathrm{\Delta }^2}\right)K_n^2\frac{2\lambda _a^2}{\mathrm{\Delta }^2}K_nK_{n+1}+\frac{2b}{a}K_n\right]\mathrm{\Delta }$$
(10)
where $`\mathrm{\Delta }`$ is the distance between two adjacent planes. Equilibrium profiles of $`K_n`$ for symmetric boundary conditions $`K_1=K_N=K_0`$ are shown in Fig. 5. Using these profiles, the orientation of the magnetization vectors in the planes $`\theta _n`$ can be calculated by minimization of the discretized analogue of the energy (5). The characteristic functions for the discretized model $`\Phi _N=2\pi M_0^2N\mathrm{\Delta }+_{n=1}^NK_n`$ shown in Fig. 6 display the same dependence on thickness $`d=N\mathrm{\Delta }`$ as the continuous model (Fig. 4).

Finally, we discuss the relation between our theory and the traditional approach. For large thickness ($`d\lambda _a)`$, the function (4) may be asymptotically written as

$$K(z)=K_0\mathrm{exp}(\frac{d\left|z\right|}{\lambda _a}).$$
(11)
The function $`K(z)`$ has finite values only in the vicinity of the surfaces and equals zero in the main part of the layer. Integrating (11) with respect to $`z`$ yields the anisotropy energy (per area)

$$K_s=K_0_0^{\mathrm{}}\mathrm{exp}(\xi /\lambda _a)𝑑\xi =K_0\lambda _a.$$
(12)
The integral parameter $`K_s`$ is equal to the additional energy connected with surface induced anisotropy. Similar parameters characterize surface tension in liquids and domain wall energy in the theory of magnetic domains. Expanding $`K_{\text{eff}}`$ in (7) for large $`d/\lambda _a`$ and substituting for $`K_s`$ (12) one obtains the Néel anisotropy ansatz (1). In the region of large thickness the surface anisotropy may be transformed into a surface term which should be included into the boundary conditions of the corresponding variation problem for the magnetization. For this case the calculation of the magnetization in magnetic layers firstly has been done by Thiaville and Fert . It is clear that the surface anisotropy $`K_s`$ is meaningful only in the case of a strong confinement of the induced anisotropy to the surface, i.e. when $`d\lambda _a`$. This relation determines the region of the validity for the Néel approach and other approaches that reduce surface effects to an effective surface energy contribution. Beyond this region the function $`\stackrel{~}{K}(z)`$ of eq. (6) has finite values within all the layer and the separation of the induced anisotropy into a surface and a volume contribution becomes meaningless.

In conclusion, our theory considers inhomogeneous distributions of induced magnetic anisotropy. It is able to describe many of the experimental features found for thin ferromagnetic films. In particular, the nonlinear deviations from Néel’s approach for surface anisotropy found in various experiments are consistently explained.

We thank H. Eschrig, R. Hayn, K.-H. Müller, and P.M. Oppeneer for discussion. A.N.B. thanks P. Fulde for hospitality and support. U.K.R. is supported by DFG.

Fig. 1
A.N. Bogdanov
Phenomenological theory of magnetic anisotropy $`\mathrm{}`$

Fig. 2
A.N. Bogdanov
Phenomenological theory of magnetic anisotropy $`\mathrm{}`$

Fig. 3
A.N. Bogdanov
Phenomenological theory of magnetic anisotropy $`\mathrm{}`$

Fig. 4
A.N. Bogdanov
Phenomenological theory of magnetic anisotropy $`\mathrm{}`$

Fig. 5
A.N. Bogdanov
Phenomenological theory of magnetic anisotropy $`\mathrm{}`$

Fig. 6
A.N. Bogdanov
Phenomenological theory of magnetic anisotropy $`\mathrm{}`$


text:
# Search for Critical Phenomena in Pb+Pb Collisions

## Abstract

NA44 uses a 512 channel Si pad array covering $`1.5<\eta <3.3`$ to study charged hadron production in Pb+Pb collisions at the CERN SPS. We apply a multiresolution analysis, based on a Discrete Wavelet Transformation, to probe the texture of particle distributions event-by-event, by simultaneous localization of features in space and scale. Scanning a broad range of multiplicities, we look for a possible critical behaviour in the power spectra of local density fluctuations. The data are compared with detailed simulations of detector response, using heavy ion event generators, and with a reference sample created via event mixing.

An ultrarelativistic collision of heavy ions presents a phenomenon whose most interesting features are conditioned by the large multitude of degrees of freedom involved, and yet offer an opportunity for the fundamental physics of the strong interaction to manifest itself. The very notion of a phase transition in such collisions is inherently of a multiparticle nature. Truly multiparticle observables, defined on event-by-event basis – a few have been constructed so far – therefore attract attention. Recently published event-by-event analyses of the 158 GeV/A $`Pb+Pb`$ data either deal with a small number of events  or analyse properties of a large ensemble of events using a single observable ($`p_T`$) to compare different ensemble averages. In the first case, accumulation of feature information from large data sets remains open. In the second case , one can not establish a scale independency in event textures by observing a logical consequence thereof . Furthermore, an ensemble average on a set of *post-freeze-out* events is not representative of the *pre-freeze-out* history of those events, due to the dramatic non-stationarity of the open system, with a consequent lack of ergodicity.

Here we concentrate on *texture*, or *local fluctuation* observables, where a single event determines its own correlation/fluctuation content, and the scale composition thereof manifests itself in the observables in a positive way. The idea to look at particle distributions in rapidity $`y`$ to search for critical behaviour was proposed  based upon a Ginzburg-Landau type of multihadron production theory , where a random hadronic field $`\varphi (y)`$ plays the role of an order parameter in a hadronization transition. Enhanced large scale correlations of hadrons in $`y`$ at the phase transition would signal critical fluctuations in the order parameter. Stephanov and coworkers  indicated a second order QCD phase transition point which should exist under certain initial conditions, within the reach of today’s experiments.

In our work, a power spectrum analysis of event texture in pseudorapidity $`\eta `$ and azimuthal angle $`\zeta `$, based on a Discrete Wavelet Transformation (DWT), is performed on a number of large event ensembles sampled according to their multiplcity, thereby studying the impact parameter dependence of the observables. DWT quantifies contributions of different $`\zeta `$ and $`\eta `$ scales into the overall event’s texture, thus testing the possible large scale enhancement.

The SPS beam was collimated to a $`1\times 2`$ mm profile. The NA44 Si pad array, installed 10 cm downstream from the target, in the magnetic field of the first dipole, measured ionization energy loss of charged particles in its 512 300 $`\mu `$m thick Si pads. The silicon detector had inner radius 7.2 mm and outer radius 43 mm. The detector was split *radially* into 16 rings of equal $`\eta `$ coverage. Each ring was divided *azimuthally* into 32 sectors of equal angular coverage to form pads. The pads were read out by AMPLEX  chips, one chip per sector. $`\delta `$-electrons, produced by the $`Pb`$ beam traversing the target, were swept away to one side by the dipole magnetic field ($`1.6`$ Tl). Only the $`\delta `$-electron-free side was used in this analysis. Channel pedestals had, on the average, $`FWHM=0.48<dE>`$ of 1 MIP. In the texture analysis, every event was represented by a 2D array of the calibrated digitized amplitudes of the channels ( an *amplitude array*). Empty target runs were used to measure the background. Cross-talk in the detector was evaluated off-line.

DWT formalizes the images of the $`PbPb`$ collision events in pseudorapidity $`\eta `$ and azimuthal angle $`\zeta `$ by expanding them into a set of functions orthogonal with respect to scale and position, and allows one to accumulate the texture information by averaging the power spectra of many events. The simplest DWT basis is the Haar one, built upon the scaling function $`\varphi (x)=1`$ for $`0x<1`$ and 0 otherwise. If the interaction vertex lies on the detector’s symmetry axis, every pad’s acceptance is a rectangle in the $`(\zeta ,\eta )`$ space. Then, the Haar basis is the natural one, as its scaling function in two dimenstions (2D) $`\mathrm{\Phi }(\zeta ,\eta )=\varphi (\zeta )\varphi (\eta )`$ is just a pad’s acceptance (modulo units). We set up a 2D wavelet basis:

$$\mathrm{\Psi }_{m,i,j}^\lambda (\zeta ,\eta )=2^m\mathrm{\Psi }^\lambda (2^m\zeta i,2^m\eta j)$$
(1)
$`\mathrm{\Phi }_{m,i,j}(\zeta ,\eta )`$ is constructed from $`\mathrm{\Phi }(\zeta ,\eta )`$ similarly. Here, $`m`$ is an integer scale fineness index; $`i`$ and $`j`$ index the discrete positions of pad centers in $`\zeta `$ and $`\eta `$ ( $`1m4`$ and $`1i,j16`$ because we use $`16=2^4`$ rings and 16 sectors ). Different values of $`\lambda `$ (denoted as $`\zeta `$, $`\eta `$, and $`\zeta \eta `$) distinguish, respectively, functions with azimuthal, pseudorapidity, and diagonal texture sensitivity:

$$\mathrm{\Psi }^\zeta =\psi (\zeta )\varphi (\eta ),\mathrm{\Psi }^\eta =\varphi (\zeta )\psi (\eta ),\mathrm{\Psi }^{\zeta \eta }=\psi (\zeta )\psi (\eta )$$
(2)
In the Haar basis, for any variable $`x`$

$$\psi (x)=\{+1\text{ for }0x<\frac{1}{2};1\text{ for }\frac{1}{2}x<1;0\text{ otherwise}\}$$
(3)
is the wavelet function. Then, $`\mathrm{\Psi }_{m,i,j}^\lambda `$ with integer $`m`$, $`i`$, and $`j`$ are known  to form an othonormal basis in $`L^2(^2)`$. We adopt the existing  1D DWT power spectrum analysis technique and expand it to 2D. Track density in an individual event is $`\rho (\zeta ,\eta )`$ and its *local* fluctuation in a given event is $`\sigma ^2\rho \overline{\rho },\rho \overline{\rho },`$ where $`\overline{\rho }`$ is the average $`\rho `$ in the given event. Using completeness of the basis, we expand

$$\rho \overline{\rho }=\rho ,\mathrm{\Psi }_{m,i,j}^\lambda \mathrm{\Psi }_{m,i,j}^\lambda \overline{\rho },\mathrm{\Psi }_{m,i,j}^\lambda \mathrm{\Psi }_{m,i,j}^\lambda $$
(4)
Notice that $`\overline{\rho }`$, being constant within detector’s rectangular acceptance, is orthogonal to any $`\mathrm{\Psi }_{m,i,j}^\lambda `$ with $`m1`$. Due to the orthonormality condition $`\mathrm{\Psi }_{m,i,j}^\lambda ,\mathrm{\Psi }_{m^{},i^{},j^{}}^\lambda ^{}=\delta _{\lambda ,\lambda ^{}}\delta _{m,m^{}}\delta _{i,i^{}}\delta _{j,j^{}}`$, the $`\rho \overline{\rho }`$ components for different scales do not form cross-terms in the $`\sigma ^2`$ sum, and the sum contains no cross-terms between $`\rho `$ and $`\overline{\rho }`$ for the four observable scales. Instead of a $`\rho ,\mathrm{\Phi }_{m=5,i,j}`$ set, the amplitude array – its closest experimentally achievable approximation – is used as the DWT input.

The Fourier images of 1D wavelet functions occupy a set of wave numbers whose characteristic broadness grows with scale fineness $`m`$ as $`2^m`$; $`2^{2m}`$ should be used in the 2D case. In 2D, we find it most informative to present the three modes of a power spectrum with different directions of sensitivity $`P^{\zeta \eta }(m)`$, $`P^\zeta (m)`$, $`P^\eta (m)`$ separately. We define *power spectrum* as

$$P^\lambda (m)=\frac{1}{2^{2m}}\underset{i,j}{}\rho ,\mathrm{\Psi }_{m,i,j}^\lambda ^2,$$
(5)
where the denominator gives the meaning of spectral *density* to the observable. So defined, the $`P(m)`$ of a random white noise field (for example) is independent of $`m`$ . In the first approximation, the white noise example provides a base-line case for comparisons in search for non-trivial effects.

We used WAILI  software library to obtain the wavelet expansions.

Figure 1 shows the power spectra for one multiplicity range. The first striking feature is that the power spectra of physical events are indeed enhanced on the coarse scale. The task of the analysis is to quantify and, as much as possible, eliminate “trivial” and experiment-specific reasons for this enhancement.

The average event, formed by summing amplitude arrays of the measured events within a multiplicity range, and dividing by the number of events, has a much reduced texture as fluctuations cancel. However it retains the texture associated with the $`d^2N/d\eta d\zeta `$, with the dead channels and the finiteness of the beam’s geometrical cross-section. A better way to get rid of the “trivial” texture is to use mixed events. The event mixing is done by taking different channels from different events. Therefore, the mixed events preserve the texture associated with the detector position offset, the inherent $`dN/d\eta `$ shape and the dead channels. In order to reproduce the electronics cross-talk effects in the mixed event sample, mixing is done sector-wise, i.e. the sectors constitute the subevents subjected to the event number scrambling. In other words, the mixed events preserve the texture coupled with the channels of the detector. This is *static* texture as it reproduces its pattern event after event; we are interested in *dynamic* texture. We reduce sources of static texture in the power spectra by empty target subtraction and by subtraction of the mixed events power spectra, thus obtaining the *texture correlation* $`P^\lambda (m)_{true}P^\lambda (m)_{mixed}`$. Its multiplicity dependence is plotted on Figure 2. For comparison with models, a MC simulation (done with RQMD ) includes the known static texture effects and undergoes the same elimination procedure. This allows the effects irreducible by the subtraction methods to be taken into account in the comparison. One such example is the finite beam size, which has been shown by the MC studies to cause the RQMD points to rise with $`dN/d\eta `$.

The systematic errors were evaluated by removing the $`Pb`$ target and switching magnetic field polarity to expose the analyzed side of the detector to $`\delta `$-electrons, while minimizing nuclear interactions. All correlations (i.e. deviations of $`P^\lambda (m)_{true}`$ from $`P^\lambda (m)_{mixed}`$) in such events are considered to be systematic errors. Thus this component of the systematic error gets a sign, and the systematic errors are asymmetric. The other component (significant only on the coarsest scale) is the uncertainty of our knowledge of the beam’s geometrical cross-section.

This novel method of event-by-event analysis, applied to the SPS $`PbPb`$ data, does not reveal any evidence of critical phenomena.

The authors thank N.Antoniou, I.Dremin, E.Shuryak, M.Stephanov, and T.Trainor for illuminating discussions. The NA44 Collaboration thanks the staff of the CERN PS-SPS accelerator complex for their excellent work, and the technical staff in the collaborating institutes for their valuable contributions. This work was supported by the Austrian Fonds zur Förderung der Wissenschaftlichen Forschung; the Science Research Council of Denmark; the Japanese Society for the Promotion of Science; the Ministry of Education, Science and Culture, Japan; the Science Research Council of Sweden; the US Department of Energy and the National Science Foundation.