format: text

text:
# TCTL Inevitability Analysis of Dense-time SystemsThe work is partially supported by NSC, Taiwan, ROC under grants NSC 90-2213-E-002-131, NSC 90-2213-E-002-132, and by the Internet protocol verification project of Institute of Applied Science & Engineering Research, Academia Sinica, 2001.

## 1 Introduction

In the research of verification, very often two types of specification properties attract most interest from academia and industry. The first type specifies that ”bad things will never happen” while the second type specifies that ”good things will happen”. In linear temporal logics, the former is captured by modal operator $`\mathrm{}`$ while the latter by $`\mathrm{}`$. Tremendous research effort has been devoted to the efficient analysis of these two types of properties in the framework of linear temporal logics. In the branching temporal logics of (timed) CTL, these two types can be mapped to modal operators $`\mathrm{}`$ and $`\mathrm{}`$ respectively. $`\mathrm{}`$ properties are called safety properties while $`\mathrm{}`$’s are usually called inevitability properties. In the domain of dense-time system verification, people have focused on the efficient analysis of safety properties. Inevitability properties in Timed CTL (TCTL) are comparatively more complex to analyze due to the following reason. In the framework of model-checking, to analyze an inevitability property, say $`\mathrm{}\varphi `$, we actually compute the set of states that satisfy the negation of inevitability, in symbols $`[[\mathrm{}\neg \varphi ]]`$, and then use the intersection emptiness between $`[[\mathrm{}\neg \varphi ]]`$ and the initial states for the answer to the inevitability anlaysis. However, property $`\mathrm{}\neg \varphi `$ in TCTL semantics is only satisfied with non-Zeno computations. (Zeno computations are those counter-intuitive infinite computations whose execution times converge to a finite value.) For example, a specification like

”along all computations, eventually a bus collision will happen in three time units”

can be violated by a Zeno computation whose execution time converges to a finite timepoint, e.g. 2.9. Such requirement on non-Zeno computations may add complexity to the evaluation of inevitability properties. In this work, we present our symbolic TCTL model-checking algorithm which can handle the non-Zeno requirement in the evaluation of greatest fixpoints. The evaluation of inevitability properties in TCTL involves nested reachability analysis and demands much higher complexity than simple safety analysis.

To contain the complexity of TCTL inevitability, it is important to integrate new and old techniques for a performance solution. In this paper, we investigate three approaches. In the first approach, we investigate how to adjust a parameter value in our greatest fixpoint evaluation algorithms for better performance. We have carried out experiments to get insight on this issue.

In the second approach, we present a technique called Early Decision on the Greatest Fixpoint (EDGF). The idea is that, in the evaluation of the greatest fixpoints, we start with a state-space and iteratively pare states from it until we reach a fixpoint. Throughout iterations of the greatest fixpoint evaluations, the state-space is non-increasing. Thus, if in a middle greatest fixpoint evaluation iteration, we find that target states have already been pared from the greatest fixpoint, we can conclude that it is not possible to include these target states in the fixpoint. Through this technique, we can reduce time-complexity irrelevant to the answer of the model-checking. As reported in section 9, significant performance improvement has been observed in several benchmarks.

Our third approach is to use abstraction techniques. We shall focus on a special subclass, TCTL, of TCTL in which every formula can be analyzed with safe abstraction if over-approximation is used in the evaluation of its negation. For example, we may write the following formula in the subclass.

$`\mathrm{}\left(\text{request}\mathrm{}\left(\text{service}\mathrm{}\mathrm{𝚛𝚎𝚚𝚞𝚎𝚜𝚝}\right)\right)`$

This formula says that if a request is responded by a service, then a request will follow the service. This subclass allows for nested modal formulas and we feel that it captures many TCTL inevitability properties.

One challenge in designing safe abstraction techniques in model-checking is making them accurate enough to discern many true properties, while still allowing us to enhance verification performance. In previous research, people have designed many abstraction techniques for reachability analysis, as have we. However, for model-checking formulas in TCTL, abstraction accuracy can be a bigger issue since the inaccuracy in abstraction can be potentially magnified when we use inaccurate evaluation results of nested modal subformulas to evaluate nesting modal subformulas with abstraction techniques. Thus it is important to discern accuracy of previous abstraction techniques in discerning true TCTLformulas.

In this paper, we also discuss another possibility for abstract evaluation of greatest fixpoints, which is to omit the requirement for non-Zeno computations in TCTL semantics. As reported in section 9, many benchmarks are true even without exclusion of Zeno computations.

Finally, we have implemented these ideas in our model-checker/simulator red 4.1. We report here our experiments to observe the effects of parameter values, EDGF, various abstraction techniques, and non-Zeno requirements on our inevitability analysis. We also compare our analysis with Kronos 5.1, which is another model-checker for full TCTL.

Our presentation is ordered as follows. Section 2 discusses several related works. Sections 3 and 4 give brief presentations of our model, timed automata (TA), and TCTL. Section 5 presents our TCTL model-checking algorithm with requirements for non-Zeno computations. Section 6 improves our model-checking algorithm using an EDGF technique. Section 7 gives another version of a greatest fixpoint evaluation algorithm by omitting the requirement of non-Zeno computations. Section 8 identifies the subclass TCTLof TCTL which supersedes many inevitability properties, while allowing for safe abstract model-checking by using over-approximation techniques. Section 9 illustrates our experiment results and helps clarify how various techniques can be used to improve analysis of inevitability properties. Section 10 is the conclusion.

## 2 Related work

The TA model with dense-time clocks was first presented in . Notably, the data-structure of DBM is proposed in  for the representation of convex state-spaces of TA. The theory and algorithm of TCTL model-checking were first given in . The algorithm is based on region graphs and helps manifest the PSPACE-complexity of the TCTL model-checking problem.

In , Henzinger et al proposed an efficient symbolic model-checking algorithm for TCTL. However, the algorithm does not distinguish between Zeno and non-Zeno computations. Instead, the authors proposed to modify TAs with Zeno computations to ones without. In comparison, our greatest fixpoint evaluation algorithm is innately able to quantify over non-Zeno computations.

Several verification tools for TA have been devised and implemented so far. UPPAAL is one of the popular tool with DBM technology. It supports safety (reachability) analysis in forward reasoning techniques. Various state-space abstraction techniques and compact representation techniques have been developed. Recently, Moller has used UPPAAL with abstraction techniques to analyze restricted inevitability properties with no modal-formula nesting. The idea is to make model augmentations to speed up the verification performance. Moller also shows how to extend the idea to analyze TCTL with only universal quantifications. However, no experiment has been reported on the verification of nested modal-formulas.

Kronos is a full TCTL model-checker with DBM technology and both forward and backward reasoning capability. Experiments to demonstrate how to use Kronos to verify several TCTL bounded inevitability properties is demonstrated in . (Bouonded inevitabilities are those inevitabilities specified with a deadline.) But no report has been made on how to enhance the performance of general inevitability analysis. In comparison, we have discussed techniques like EDGF and abstractions which handle both bounded and unbounded inevitabilities.

DDD is a reachability analyzer based on BDD-like data-structures for TA.

SGM is a compositional safety (reachability) analyzer for TA, also based on DBM technology. A newer version also supports partial TCTL model-checking.

CMC is another compositional model-checker. Its specification language is a restricted subclass of $`T_\mu `$ and is capable of specifying bounded inevitabilities.

Our tool red (version 4.1) is a full TCTL model-checker/simulator with a BDD-like data-structure, called CRD (clock-restriction diagram). Previous research with red has focused on enhancing the performance of safety analysis.

Abstraction techniques for analysis have been studied in great depth since the pioneering work of Cousot et al. For TA, convex-hull over-approximation has been a popular choice for DBM technology due to its intuitiveness and effective performance. It is difficult to implement this over-approximation in red since variable-accessing has to observe variable-orderings of BDD-like data-structures. Nevertheless, many over-approximation techniques for TA have been reported in  for BDD-like data-structures and in  specifically for CRD.

Relations between abstraction techniques and subclasses of CTL with only universal (or existential respectively) path quantifiers has been studied in . As mentioned above, the corresponding framework in TCTL is noted in .

## 3 Timed automata (TA)

We use the widely accepted model of timed automata, which is a finite-state automata equipped with a finite set of clocks which can hold nonnegative real-values. At any moment, a timed automata can stay in only one mode (or control location). In its operation, one of the transitions can be triggered when the corresponding triggering condition is satisfied. Upon being triggered, the automata instantaneously transits from one mode to another and resets some clocks to zero. Between transitions, all clocks increase readings at a uniform rate.

For convenience, given a set $`Q`$ of modes and a set $`X`$ of clocks, we use $`B(Q,X)`$ as the set of all Boolean combinations of atoms of the forms $`q`$ and $`xx^{}c`$, where $`qQ`$, $`x,x^{}X\{0\}`$, “$``$” is one of $`,<,=,>,`$, and $`c`$ is an integer constant.

###### Definition 1

timed automata (TA)

A TA $`A`$ is given as a tuple $`X,Q,I,\mu ,T,\tau ,\pi `$ with the following restrictions. $`X`$ is the set of clocks. $`Q`$ is the set of modes. $`IB(Q,X)`$ is the initial condition. $`\mu :QB(\mathrm{},X)`$ defines the invariance condition of each mode. $`TQ\times Q`$ is the set of transitions. $`\tau :TB(\mathrm{},X)`$ and $`\pi :T2^X`$ respectively define the triggering condition and the clock set to reset of each transition.  $`\mathrm{}`$

A valuation of a set is a mapping from the set to another set. Given an $`\eta B(Q,X)`$ and a valuation $`\nu `$ of $`X`$, we say $`\nu `$ satisfies $`\eta `$, in symbols $`\nu \eta `$, iff it is the case that when the variables in $`\eta `$ are interpreted according to $`\nu `$, $`\eta `$ will be evaluated as true.

###### Definition 2

states

A state $`\nu `$ of $`A=X,Q,I,\mu ,T,\tau ,\pi `$ is a valuation of $`XQ`$ such that

* there is a unique $`qQ`$ such that $`\nu (q)=\text{true}`$ and for all $`q^{}q,\nu (q^{})=\text{false}`$;
* for each $`xX`$, $`\nu (x)^+`$ (the set of nonnegative reals) and $`qQ,\nu (q)\nu \mu (q)`$.

Given state $`\nu `$ and $`qQ`$ such that $`\nu (q)=\text{true}`$, we call $`q`$ the mode of $`\nu `$, in symbols $`\nu ^Q`$.  $`\mathrm{}`$

For any $`t^+`$, $`\nu +t`$ is a state identical to $`\nu `$ except that for every clock $`xX`$, $`\nu (x)+t=(\nu +t)(x)`$. Given $`\overline{X}X`$, $`\nu \overline{X}`$ is a new state identical to $`\nu `$ except that for every $`x\overline{X}`$, $`\nu \overline{X}(x)=0`$.

###### Definition 3

runs (computations)

Given a TA $`A=X,Q,I,\mu ,T,\tau ,\pi `$, a run is an (infinite) sequence of state-time pairs, $`(\nu _0,t_0)(\nu _1,t_1)\mathrm{}(\nu _k,t_k)\mathrm{}\mathrm{}`$ , such that $`\nu _0I`$ and $`t_0t_1\mathrm{}t_k\mathrm{}\mathrm{}`$ is a monotonically increasing real-number (time) divergent sequence, and for all $`k0`$,

* invariance conditions are preserved in each interval: that is,
for all $`t[0,t_{k+1}t_k]`$, $`\nu _k+t\mu (\nu _k^Q)`$; and
* either no transition happens at time $`t_k`$, that is, $`\nu _k^Q=\nu _{k+1}^Q`$ and $`\nu _k+(t_{k+1}t_k)=\nu _{k+1}`$; or a transition happens at $`t_k`$, that is,

+ there is such a transition, that is $`(\nu _k^Q,\nu _{k+1}^Q)T`$; and
    + the corresponding triggering condition is satisfied, that is, $`\nu _k+(t_{k+1}t_k)\tau (\nu _k^Q,\nu _{k+1}^Q)`$; and
    + the clocks are reset to zero accordingly, that is, $`(\nu _k+(t_{k+1}t_k))\pi (\nu _k^Q,\nu _{k+1}^Q)=\nu _{k+1}`$.  $`\mathrm{}`$

## 4 TCTL (Timed CTL)

TCTL is a branching temporal logic for the specification of dense-time systems.

###### Definition 4

(Syntax of TCTL formulas):

A TCTL formula $`\varphi `$ has the following syntax rules.

$`\begin{array}{ccc}\hfill \varphi & \hfill ::=& \eta |\neg \varphi _1|\varphi _1\varphi _2|x.\varphi _1|\varphi _1𝒰\varphi _2|\mathrm{}\varphi _1\hfill \end{array}`$

Here $`\eta B(Q,X)`$ and $`\varphi _1`$, $`\varphi _2`$ are TCTL formulas.  $`\mathrm{}`$

The modal operators are intuitively explained in the following.

* $`x.\varphi `$ means that ”if there is a clock $`x`$ with reading zero now, then $`\varphi `$ is satisfied.”
* $``$ means “there exists a run”
* $`\varphi _1𝒰\varphi _2`$ means that along a computation, $`\varphi _1`$ is true until $`\varphi _2`$ becomes true.
* $`\mathrm{}\varphi _1`$ means that along a computation, $`\varphi _1`$ is always true.

Besides the standard shorthands of temporal logics, we adopt the following for TCTL:

$`\begin{array}{cccc}\hfill & \mathrm{}\varphi _1\text{ for }\text{true}𝒰\varphi _1\hfill & \hfill & \mathrm{}\varphi _1\text{ for }\neg \mathrm{}\neg \varphi _1\hfill \\ \hfill & \varphi _1𝒰\varphi _2\text{ for }\neg (((\neg \varphi _2)𝒰\neg (\varphi _1\varphi _2))(\mathrm{}\neg \varphi _2))\hfill & \hfill & \mathrm{}\varphi _1\text{ for }\text{true}𝒰\varphi _1\hfill \end{array}`$

###### Definition 5

(Satisfaction of TCTL formulas):

We write in notations $`A,\nu \varphi `$ to mean that $`\varphi `$ is satisfied at state $`\nu `$ in TA $`A`$. The satisfaction relation is defined inductively as follows.

* When $`\varphi _1B(Q,X)`$, $`A,\nu \varphi _1`$ according to the definition in the beginning of subsection 3.
* $`A,\nu \varphi _1\varphi _2`$ iff either $`A,\nu \varphi _1`$ or $`A,\nu \varphi _2`$
* $`A,\nu \neg \varphi _1`$ iff $`A,\nu \vDash ̸\varphi _1`$
* $`A,\nu x.\varphi _1`$ iff $`A,\nu \{x\}\varphi _1`$.
* $`A,\nu \varphi _1𝒰\varphi _2`$ iff there exists a run $`(\nu _1,t_1)(\nu _2,t_2)\mathrm{}`$ such that $`\nu _1=\nu `$ in $`A`$, and there exist an $`i1`$ and a $`\delta [0,t_{i+1}t_i]`$, s.t.

+ $`A,\nu _i+\delta \varphi _2`$,
    + for all $`j,\delta ^{}`$, if either $`(1j<i)(\delta ^{}[0,t_{j+1}t_j])`$ or $`(j=i)(\delta ^{}[0,\delta ))`$, then $`A,\nu _j+\delta ^{}\varphi _1`$.

In words, $`\nu `$ satisfies $`\varphi _1𝒰\varphi _2`$ iff there exists a run from $`\nu `$ such that along the run, $`\varphi _1`$ is true until $`\varphi _2`$ is true.
* $`A,\nu \mathrm{}\varphi _1`$ iff there exists a run $`(\nu _1,t_1)(\nu _2,t_2)\mathrm{}`$ such that $`\nu _1=\nu `$ in $`A`$, and for every $`i1`$ and $`\delta [0,t_{i+1}t_i]`$, $`A,\nu _i+\delta \varphi _1`$. In other words, $`\nu `$ satisfies $`\mathrm{}\varphi _1`$ iff there exists a run from $`\nu `$ such that $`\varphi _1`$ is always true.

A TA $`A`$ satisfies a TCTL formula $`\varphi `$, in symbols $`A\varphi `$, iff for every state $`\nu _0I`$, $`A,\nu _0\varphi `$.  $`\mathrm{}`$

## 5 TCTL Model-checking algorithm with non-Zeno requirements

Our model-checking algorithm is backward reasoning. We need two basic procedures, one for the computation of the weakest precondition of transitions, and the other for backward time-progression. These two procedures are important in the symbolic construction of backward reachable state-space representations. Various presentations of the two procedures can be found in . Given a state-space representation $`\eta `$ and a transition $`e`$, the first procedure, $`\text{xtion\_bck}(\eta ,e)`$, computes the weakest precondition

* in which, every state satisfies the invariance condition imposed by $`\mu ()`$; and
* from which we can transit to states in $`\eta `$ through $`e`$.

The second procedure, $`\text{time\_bck}(\eta )`$, computes the space representation of states

* from which we can go to states in $`\eta `$ simply by time-passage; and
* every state in the time-passage also satisfies the invariance condition imposed by $`\mu ()`$.

With the two basic procedures, we can construct a symbolic backward reachability procedure as in . We call this procedure $`\text{reachable-bck}(\eta _1,\eta _2)`$ for convenience. Intuitively, $`\text{reachable-bck}(\eta _1,\eta _2)`$ characterizes the backwardly reachable state-space from states in $`\eta _2`$ through runs along which all states satisfy $`\eta _1`$. Computationally, $`\text{reachable-bck}(\eta _1,\eta _2)`$ can be defined as the least fixpoint of the equation $`Y=\eta _2\left(\eta _1\text{time\_bck}(\eta _1_{eT}\text{xtion\_bck}(Y,e))\right)`$, i.e.,

$`\text{reachable-bck}(\eta _1,\eta _2)\text{lfp }Y.\left(\eta _2\left(\eta _1\text{time\_bck}(\eta _1_{eT}\text{xtion\_bck}(Y,e))\right)\right)`$.

Our model-checking algorithm is modified from the classic model-checking algorithm for TCTL. The design of the greatest fixpoint evaluation algorithm with consideration of non-Zeno requirement is based on the following lemma.

###### Lemma 6

Given $`d1`$, $`A,\nu \mathrm{}\eta `$ iff there is a finite run from $`\nu `$ of duration $`d`$ such that along the run every state satisfies $`\eta `$ and the finite run ends at a state satisfying $`\mathrm{}\eta `$.

Proof : Details are omitted due to page-limit. But note that we can construct an infinite and divergent run by concatenating an infinite sequence of finite runs with durations $`d1`$. The existence of infinitely many such concatenable finite runs is assured by the recursive construction of $`\mathrm{}\eta `$.  $`\mathrm{}`$

Then $`\mathrm{}\eta `$ can be defined with the following greatest fixpoint.

$`\mathrm{}\eta \text{gfp }Y.(\text{ZC}.\text{reachable-bck}(\eta ,Y\text{ZC}d))`$

Here clock ZC is used specifically to measure the non-Zeno requirement. The following procedure can construct the greatest fixpoint satisfying $`\mathrm{}\eta `$ with a non-Zeno requirement.   $`\text{gfp}(\eta )`$ /\* d is a static parameter for measuring time-progress \*/ {
$`Y:=\eta `$; $`Y^{}:=\text{true}`$;
repeat until $`Y=Y^{}`$, { (1)
$`Y^{}:=Y`$; $`Y:=Y\text{clock\_eliminate}(\text{ZC}=0\text{reachable-bck}(\eta ,Y\text{ZC}d),\text{ZC})`$; (2)
}
return $`Y`$;
}   Here $`\text{clock\_eliminate}()`$ removes a clock from a state-predicate without losing information on relations among other clocks. Details can be found in appendix A.

Note here that $`d`$ works as a parameter. We can choose the value of $`d1`$ for better performance in the computation of the greatest fixpoint.

Procedure gfp() can be used in the labeling algorithm in  to replace the evaluation of $`\mathrm{}`$-formulas. For completeness of the presentation, please check appendix B to see our complete model-checking algorithm with non-Zeno requirement. The correctness follows from Lemma 6.

## 6 Model-checking with Early Decision on the Greatest Fixpoint (EDGF)

In the evaluation of the greatest fixpoint for formulas like $`\mathrm{}\varphi _1`$, we start from the description, say $`Y`$, for a subspace of $`\varphi _1`$ and iteratively eliminate those subspaces which cannot go to a state in $`Y`$ through finite runs of $`d1`$ time units. Thus, the state-space represented by $`Y`$ shrinks iteratively until it settles at a fixpoint. In practice, this greatest fixpoint usually happens in conjunction with other formulas. For example, we may want to specify $`\text{collision}y.\mathrm{}(y<26\text{idle})`$ meaning that a bus at the collision state, will enter the idle state in 26 time-units. After negation for model-checking, we get $`\text{collision}y.\mathrm{}(y26\neg \text{idle})`$. In evaluating this negated formula, we want to see if the greatest fixpoint for the $`\mathrm{}`$-formula intersects with the state-space for collision. We do not actually have to compute the greatest fixpoint to know if the intersection is empty. Since the value of $`Y`$ iteratively shrinks, we can check if the intersection between $`Y`$ and the state-space for collision becomes empty at each iteration of the greatest fixpoint construction (i.e., the repeat-loop at statement (1) in procedure gfp). If at an iteration, we find the intersection with $`Y`$ is already empty, then there is no need to continue calculating the greatest fixpoint and we can immediately return the current value of $`Y`$ (or false) without affecting the result of the model-checking.

Based on this idea, we rewrite our model-checking algorithm with our Early Deision on the Greatest Fixpoint (EDGF). We introduce a new parameter $`\beta `$ to pass the information of the target states inherited from the scope.   $`\text{Eval-EDGF}(A,\chi ,\beta ,\overline{\varphi })`$
/\* $`\chi `$ is the set of clocks declared in the scope of $`\overline{\varphi }`$ \*/
/\* $`\beta `$ is constraints inhereted in the scope of $`\overline{\varphi }`$ for early decision of gfp \*/ {
switch ($`\overline{\varphi }`$) {
case (false): return false;
case ($`p`$): return $`p_{x\chi }x0`$;
case ($`xyc`$): return $`xyc_{x\chi }x0`$;
case ($`\varphi _1\varphi _2`$): return $`\text{Eval-EDGF}(A,\chi ,\beta ,\varphi _1)\text{Eval-EDGF}(A,\chi ,\beta ,\varphi _2)`$;
case ($`\varphi _1\varphi _2`$):
if $`\varphi _2`$ does not contain modal operator, {
$`\eta _2:=\text{Eval-EDGF}(A,\chi ,\beta ,\varphi _2)`$;
return $`\eta _2\text{Eval-EDGF}(A,\chi ,\beta \eta _2,\varphi _1)`$; (3)
}
else {
$`\eta _1:=\text{Eval-EDGF}(A,\chi ,\beta ,\varphi _1)`$;
return $`\eta _1\text{Eval-EDGF}(A,\chi ,\beta \eta _1,\varphi _2)`$; (4)
}
case ($`\neg \varphi _1`$): return $`\neg \text{Eval-EDGF}(A,\chi ,\text{true},\varphi _1)`$;
case ($`x.\varphi _1`$): return $`\text{clock\_eliminate}(x=0\text{Eval-EDGF}(A,\chi \{x\},\beta ,\varphi _1x0),x)`$;
case ($`\varphi _1𝒰\varphi _2`$):
$`\eta _1:=\text{Eval-EDGF}(A,\chi ,\text{true},\varphi _1)`$; $`\eta _2:=\text{Eval-EDGF}(A,\chi ,\text{true},\varphi _2)`$;
return $`\text{reachable-bck}(\eta _1,\eta _2)`$;
case ($`\mathrm{}\varphi _1`$): return $`\text{gfp\_EDGF}(\text{Eval-EDGF}(\varphi _1),\beta )`$;
}
}
$`\text{gfp\_EDGF}(\eta ,\beta )`$ /\* d is a static parameter for measuring time-progress \*/ {
$`Y:=\eta `$; $`Y^{}:=\text{true}`$;
repeat until $`Y=Y^{}`$ or $`(Y\beta )=\text{false}`$, { (5)
$`Y^{}:=Y`$; $`Y:=Y\text{clock\_eliminate}(\text{ZC}=0\text{reachable-bck}(\eta ,Y\text{ZC}d),\text{ZC})`$; (6)
}
return $`Y`$;
}   To model-check TA $`A`$ against TCTL formula $`\varphi `$, we reply true iff $`\text{Eval-EDGF}(A,\mathrm{},\text{true},\neg \varphi )`$ is false. As can be seen from statement (3) and (4) in the case of conjunction formulas, we strengthen the target-state information. In the evaluation of the greatest fixpoint, we use condition-testing $`(Y\beta )=\text{false}`$ in statement (5) respectively to check for early decision.

## 7 Greatest fixpoint computation by tolerating Zenoness

In practice, the greatest fixpoint computation procedures presented in the last two sections can be costly in computing resources since their characterizations have a least fixpoint nested in a greatest fixpoint. This is necessary to guarantee that only nonZeno computations are considered. In reality, it may happen that, due to well-designed behaviors, systems may still satisfy certain inevitability properties for both Zeno and non-Zeno computation. In this case, we can benefit from a less expensive procedure to compute the greatest fixpoint. For example, we have designed the following procedure which does not rule out Zeno computations in the evaluation of $`\mathrm{}`$-formulas.   $`\text{gfp\_Zeno\_EDGF}(\eta ,\beta )`$ {
$`Y:=\eta ;Y^{}:=\text{true}`$;
repeat until $`Y=Y^{}`$ or $`(Y\beta )=\text{false}`$, { (7)
$`Y^{}:=Y`$; $`Y:=Y\eta \text{time\_bck}\left(\eta _{eT}\text{xtion\_bck}(Y,e)\right)`$;
}
return $`Y`$;
}   Even if the procedure can be imprecise in over-estimation of the greatest fixpoint, it can be much less expensive in the verification of well-designed real-world projects.

## 8 Abstract model-checking with TCTL

We have also experimented with abstraction techniques in the evaluation of greatest fixpoints. Due to page-limit, we shall leave the explanation in appendix C. The corresponding xperiment report is in subsection 9.4.

## 9 Implementation and experiments

We have implemented the ideas in our model-checker/simulator, red version 4.1, for TA. red uses the new BDD-like data-structure, CRD (Clock-Restriction Diagram), and supports both forward and backward analysis, full TCTL model-checking with non-Zeno computations, deadlock detection, and counter-example generation. Users can also declare global and local (to each process) variables of type clock, integer, and pointer (to identifier of processes). Boolean conditions on variables can be tested and variable values can be assigned. The TCTL formulas in red also allow quantification on process identifiers for succinct specification. Interested readers can download red for free from

```
                             http://cc.ee.ntu.edu.tw/~val/

```

We design our experiment in two ways. First, we run red 4.1 with various options and benchmarks to test if our ideas can indeed improve the verification performance of inevitability properties in TCTL. Second, we compare red 4.1 with Kronos 5.2 to check if our implementation remains competitive in regard to other tools. However, we remind the readers that comparison report with other tools should be read carefully since red uses different data-structures from Kronos. Moreover, it is difficult to know what fine-tuning techniques each tool has used. Thus it is difficult to conclude if the techniques presented in this work really contribute to the performance difference between red and Kronos. Nevertheless, we believe it is still an objective measure to roughly estimate how our ideas perform.

In the following section, we shall first discuss the design of our benchmarks, then report our experiments. Data is collected on a Pentium 4 1.7GHz with 256MB memory running LINUX. Execution times are collected for Kronos while times and memory (for data-structure) are collected for red. ”s” means seconds of CPU time, ”k” means kilobytes for memory space for data-structures, ”O/M” means ”out-of-memory.”

### 9.1 Benchmarks

We do not claim that the benchmarks selected here represent the complete spectrum of model-checking tasks. The evaluation of TCTL formulas may incur various complex computations depending on the structures of the timed automata and the specification formulas. But we do carefully choose our benchmarks according to the broad spectrum of combination of models and specifications so that we can gain some insights about performance enhancement of TCTL inevitability analysis. Benchmarks include three different timed automatas and specifications for unbounded inevitability, bounded inevitability, and modal operators with nesting depth zero, one, and two respectively. We identify one important benchmark which can only be verified with non-Zeno computations. The other benchmarks can be (safely) verified without requirement of non-Zeno computations.

Due to page-limit, we leave the description of the benchmarks in appendix D.

### 9.2 Performance w.r.t. parameter for measuring time-progress

In statement (2) of procedure gfp() and statement (6) of procedure gfp\_EDGF(), we use inequality $`\text{ZC}d`$ to check time-progress in non-Zeno computations, where $`d`$ is a parameter $`1`$. We can choose various values for the parameter in our implementations. In our experiment reported in this subsection, we have found that the value of parameter $`d`$ can greatly affect the verification performance.

In this experiment, we shall use various values of parameter $`d`$ ranging from $`1`$ to beyond the biggest timing constants used in the models. For the leader-election benchmark, the biggest timing constant used is $`2`$. For the Pathos benchmark, the biggest timing constant used is equal to the number of processes. For the CSMA/CD benchmarks (A), (B), and (C), the biggest timing constant used is equal to 808.

In fact, we can also use inequality $`\text{ZC}>d`$, with $`d1`$, in statements (2) and (6) of procedures gfp() and gfp\_EDGF() respectively. Due to page-limit, we shall leave the performance data table to appendix E. We have drawn charts to show time-complexity for the benchmarks w.r.t. $`d`$-values in figure 1.

More charts for the space-complexity can be found in appendix F.

As can be seen from the charts, our algorithms may respond with different complexity curves to various model structures and specifications. For benchmarks leader-election and PATHOS, it seems that the bigger the $`d`$-value, the better the performance. For the three CSMA/CD benchmarks, it seems that the best performance happens when $`d`$ is around 80. But one thing common in these charts is that $`d=1`$ always gives the worst performance.

We have to admit that we do not have a theory to analyze or predict the complexity curves w.r.t. various model structures and specifications. More experiments on more benchmarks may be needed in order to get more understanding of the curves. In general, we feel it can be difficult to analyze such complexity curves. After all, our models of TA are still ”programs” in some sense.

Nevertheless, we have still tried hard to look into the execution of our algorithms for explanation of the complexity cuves. Procedures gfp() and gfp\_EDGF() both are constructed with an inner loop (for the least fixpoint evaluation of reachable-bck()) and an outer loop (for the greatest fixpoint evaluation). With bigger $`d`$-values, it seems that the outer loop converges faster while the inner loop converges slower. That is to say, with bigger $`d`$-values, we may need less iterations of the outer-loop and, in the same time, more iterations of the inner loop to compute the greatest fixpoints. The complexity patterns in the charts are thus superpositions between the complexities of the outer loop and the inner loop.

We have used the $`d`$-values with the best performance for the experiments reported in the next few subsections. For benchmarks PATHOS and leader-election, $`d`$ is set to $`>C_{A:\varphi }`$. ($`C_{A:\varphi }`$ is the biggest timing constant used in model $`A`$ and TCTL specification $`\varphi `$.) For the three CSMA/CD benchmarks, $`d`$ is set to $`80`$.

### 9.3 Performance w.r.t. non-Zeno requirement and EDGF

In our first experiment, we observe the performance of our inevitability analysis algorithm w.r.t. the non-Zeno requirement and the EDGF policy. The performance data is in table 1.

In general, we find that with or without EDGF technique, a non-Zeno requirement does add more complexity to the evaluation of the inevitability properties. Especially for the three specifications of CSMA/CD model, exponential blow-ups have been observed.

For PATHOS benchmark, it is a totally different story. The non-Zeno requirement seems to incur much less complexity than without it. After we have carefully traced the execution of our mode-checker, we found that this benchmark incurs very few iterations of loop (5) in gfp\_EDGF() although each iteration can be costly to run. On the other hand, it incurs a significant number of iterations of loop (6) in gfp\_Zeno\_EDGF() although each iteration is not so costly. The accumulative effect of the loop iterations result in performance that contradicts our expectation. This benchmark shows that the evaluation performance of inevitability properties is very involved and depends on many factors.

Futhermore, benchmark CSMA/CD (B) shows that some inevitability properties can only be verified with non-Zeno computations.

As for the performance of the EDGF technique, we find that when the technique fails, it only incurs a small overhead. When it succeeds, it significantly improves performance two to three-fold.

### 9.4 Performance w.r.t. abstraction techniques

In table 2, we report the performance data of our red 4.1 with respect to our three abstraction techniques.

In general, the abstraction techniques give us much better performance. Notably, game-discrete and game-magnitude abstractions seem to have enough accuracy to discern true properties.

It is somewhat surprising that the game-magnitude abstraction incurs excessive complexity for PATHOS benchmark. After carefully examining the traces generated by red, we found that because non-magnitude constraints were eliminated, some of the inconsistent convex state-spaces in the representation became consistent. These spurious convex state-spaces make many more paths in our CRD and greately burden our greatest fixpoint calculation. For instance, the outer-loop (5) of procedure gfp\_EDGF() takes two iterations to reach the fixpoint with the game-magnitude abstraction. It only takes one iteration to do so without the abstraction. In our previous experience, this abstraction has worked efficiently with reachability analysis. It seems that the performance of abstraction techniques for greatest fixpoint evaluation can be subtle.

### 9.5 Performance w.r.t. Kronos

In table 3, we report the performance of Kronos 5.2 w.r.t. the five benchmarks.

For PATHOS and leader election, Kronos did not succeed in constructing the quotient automata. But our red seems to have no problem in this regard with its on-the-fly exploration of the state-space. Of course, the lack of high-level data-variables in Kronos’ modeling language may also acerbate the problem.

As for benchmark CSMA/CD (A), Kronos performs very well. We believe this is because this benchmark uses a bounded inevitability specification. Such properties have already been studied in the literature of Kronos.

On the other hand, benchmarks CSMA/CD (B) and (C) use unbounded inevitability specifications with modal-subformula nesting depths 1 and 2 respectively. Kronos does not scale up to the complexity of concurrency for these two benchmarks. Our red prevails in these two benchmarks.

## 10 Conclusion

How to enhance the performance of TCTL model-checking is a research issue to which people have not paid much attention. The reason may be that the reachability analysis problem for TA is already difficult enough and has absorbed much of our energy. Nevertheless, the issue is still important both theoretically and practically. Hopefully, this work can give us some insight to the complexity of the issue and attract more research interest in this regard. Specifically, the charts reported in subsection 9.2 may imply that further research is needed to investigate how to predict good $`d`$-values in real-world verification tasks. Our implementation shows that the ideas in this paper can be of potential use.

APPENDICES

## Appendix A Procedure clock\_elimiante()

$`\text{clock\_eliminate}(\eta ,x)`$ {
for each $`x_1xc`$ and $`xx_2^{}c^{}`$, if $`\eta x_1xcxx_2^{}c^{}`$ is not empty, {
$`\eta _1:=\eta x_1xcxx_2^{}c`$;
$`\eta :=\eta \neg \eta _1`$; $`\eta :=\eta (\eta _1x_1x_2\text{compose\_upperbound}(,c,^{},c^{}))`$;
}
return $`\eta `$;
}
$`\text{compse\_upperbound}(,c,^{},c^{})`$ {
if $`c=\mathrm{}c^{}=\mathrm{}`$, return ”$`<\mathrm{}`$
else if $`c=\mathrm{}`$,
if $`c^{}0`$, return ”$`<\mathrm{}`$”; else return ”$`<C_{A:\varphi }+c^{}`$;”
else if $`c^{}=\mathrm{}`$,
if $`c0`$, return ”$`<\mathrm{}`$”; else return ”$`<C_{A:\varphi }+c`$;”
$`c_r:=c+c^{}`$;
if $`c_r>C_{A:\varphi }(_r=\text{}<\text{}c_r=C_{A:\varphi })`$, return ”$`<\mathrm{}`$;”
if either $``$ or $`^{}`$ is ”$`<`$”, $`_r`$ is assigned ”$`<`$”, else $`_r`$ is assigned ”$``$
else if $`c_r<C_{A:\varphi }(_r=\text{}<\text{}c_r=C_{A:\varphi })`$, return ”$`<\mathrm{}`$;”
else return ”$`_rc_r`$”;
}

Procedure compose\_upperbound() computes the new upperbound as the result of adding two, up to the absolute bound of $`C_{A:\varphi }`$. For example, with $`C_{A:\varphi }=5`$, $`\text{compose\_upperbound}(<,2,,3)=\text{}<\mathrm{}\text{}`$ and $`\text{compose\_upperbound}(,2,,1)=\text{}3\text{}`$.

## Appendix B Complete model-checking procedure with non-Zeno requirement

$`\text{model-check}(A,\varphi )`$ {
if $`\text{Eval}(A,\mathrm{},\neg \varphi )`$ is false, return true; else return false.
}
$`\text{Eval}(A,\chi ,\overline{\varphi })`$
/\* $`\chi `$ is the set of clocks declared in the scope of $`\overline{\varphi }`$ \*/ {
switch ($`\overline{\varphi }`$) {
case (false): return false;
case ($`p`$): return $`p_{x\chi }x0`$;
case ($`xyc`$): return $`xyc_{x\chi }x0`$;
case ($`\varphi _1\varphi _2`$): return $`\text{Eval}(A,\chi ,\varphi _1)\text{Eval}(A,\chi ,\varphi _2)`$;
case ($`\varphi _1\varphi _2`$):
return $`\text{Eval}(A,\chi ,\varphi _1)\text{Eval}(A,\chi ,\varphi _2)`$;
case ($`\neg \varphi _1`$): return $`\neg \text{Eval}(A,\chi ,\varphi _1)`$;
case ($`x.\varphi _1`$): return $`\text{clock\_eliminate}(x=0\text{Eval}(A,\chi \{x\},\varphi _1x0),x)`$;
case ($`\varphi _1𝒰\varphi _2`$):
$`\eta _1:=\text{Eval}(A,\chi ,\varphi _1)`$; $`\eta _2:=\text{Eval}(A,\chi ,\varphi _2)`$;
return $`\text{reachable-bck}(\eta _1,\eta _2)`$;
case ($`\mathrm{}\varphi _1`$): return $`\text{gfp}(\text{Eval}(A,\varphi _1))`$;
}
}

## Appendix C Abstract model-checking with TCTL

In the application of abstraction techniques, it is important to make them safe. That is to say, when the safe abstraction analyzer says a property is true, the property is indeed true. (But when it says false, we do not know whether the property is true.) There are two types of abstractions: over-approximation and under-approximation. The former means that the abstract state-space is a superset of the concrete state-space, while the latter means that the abstract state-space is a subset of the concrete state-space. To make an abstraction safe means that we should stick to over-approximation while evaluating $`\mathrm{}\neg \varphi `$ (the negation of the inevitability). But this can be difficult to enforce in general since negations deeply nested in formulas can turn over-approximation into under-approximation and thus make abstractions unsafe.

Usually inevitability properties do not occur on their own either. Instead, they are usually nested in other modal-formulas. For example, normally we may specify that

When there is a request, then eventually there is a service.

In TCTL, this can be written as

$`\mathrm{}\left(\text{request}\mathrm{}\text{service}\right)`$ $`(f1)`$

In general, it can be difficult to restrict the negation of over-approximation from happening. But fortunately, formula $`(f1)`$ does not have such a problem. Consider the negation of formula $`(f1)`$, which is the following.

$`\mathrm{}\left(\text{request}\mathrm{}\neg \text{service}\right)`$ $`(f2)`$

Since there is no negation sign before any modal-formula, there is no problem with negation of over-approximation. Thus over-approximation can be applied here by doing over-approximations of the state sets thus satisfying the following subformulas in sequence (from left to right).

$`\neg \text{service},\mathrm{}\neg \text{service},\text{request},\text{request}\mathrm{}\neg \text{service},(f2),I(f2)`$

If the state-set for $`I(f2)`$ is empty, then formula $`(f1)`$ is satisfied; else we do not have a conclusion. Note that the evaluation of state sets for $`\neg \text{service}`$ and request respectively does not need any abstraction.

The reasoning in the last paragraph can be extended to a general subclass of TCTL, TCTL. A formula is in TCTLiff the negation signs only appear before its atoms, and only universal path quantifications are used. For example, we may want to verify the following specification:

$`\mathrm{}\left(\text{request}\mathrm{}\left(\text{service}\mathrm{}\text{request}\right)\right)`$

The formula says that if there is a request followed by a service, then from that service on, there will be a request. The negation of the specification is $`\mathrm{}\left(\text{request}\mathrm{}\left(\text{service}\mathrm{}\neg \text{request}\right)\right)`$, which is in the subclass of TCTL, the subclass of TCTL with negations right before atoms and with only existential path quantifications. Note that the negations of formulas in TCTL fall correctly in TCTL. The following lemma shows that over-approximation techniques with TCTLformulas always yield over-approximation.

###### Lemma 7

: Given a TCTLformula $`\varphi `$, if we evaluate each modal-subformula in $`\varphi `$ with over-approximation, then we still get an over-approximation of the state set satisfying $`\varphi `$.

Proof : This can be done by an inductive analysis on the structure of $`\varphi `$. If $`\varphi `$ is a literal expression of the forms $`p`$ or $`\neg p`$, then the evaluation does not involve any approximation. If $`\varphi `$ is like $`\varphi _1\varphi _2`$ or $`\varphi _1\varphi _2`$, then the evaluation of $`\varphi `$ still yields over-approximation with the inductive hypothesis that the evaluations of $`\varphi _1`$ and $`\varphi _2`$ are both over-approximations. If $`\varphi `$ is like $`\varphi _1𝒰\varphi _2`$, then since the modal-formula is to be evaluated with over-approximation, with the inductive hypothesis, we know that $`\varphi `$ is evaluated with over-approximation. The case for $`\mathrm{}\varphi _1`$ is similar. Thus the lemma is proven.  $`\mathrm{}`$

Our over-approximation technique is directly applied in procedure $`\text{reachable-bck}()`$. In other words, we can extend $`\text{reachable-bck}()`$ with over-approximation techniques as following.

$`\text{reachable-bck}^O(\eta _1,\eta _2)\text{lfp}Y.\text{abs}\left(\eta _2\left(\eta _1\text{time\_bck}(\eta _1_{eT}\text{xtion\_bck}(Y,e))\right)\right)`$.

Here $`\text{abs}()`$ means a generic over-approximation procedure. Thus procedure $`\text{reachable-bck}^O()`$ can be used in place of $`\text{reachable-bck}()`$ in procedures gfp(), Eval-EDGF(), and gfp\_EDGF().

In our tool red 4.1, we have implemented a series of game-based abstraction procedures suitable for BDD-like data-structures and concurrent systems. We use the term ”game” here because we envision the concurrent system operation as a game. Those processes, which we want to verify, are treated as players while the other processes are treated as opponents. In the game, the players try to win (maintain the specification property) under the worst (i.e., minimal) assumption on their opponents. A process is a player iff its local variables appear in the inevitability properties. The other processes are called opponents. According to the well-observed discipline of modular programming, the behavioral correctness of a functional module should be based on minimal assumption on the environment. These game-based abstraction procedures omit opponents’ state-information to make abstractions.

* Game-abstraction: The game abstraction procedure will eliminate the state information of the opponents from its argument state-predicate.
* Game-discrete-abstraction: This abstraction procedure will eliminate all clock constraints for the opponents in the argument state-predicate.
* Game-magnitude-abstraction: A clock constraint like $`xx^{}c`$ is called a magnitude constraint iff either $`x`$ or $`x^{}`$ is zero itself (i.e. the constraint is either $`xc`$ or $`x^{}c`$). This abstraction procedure will erase all non-magnitude constraints of the opponents in the argument state-predicate.

Details can be found in .

## Appendix D Benchmarks

We use the following benchmarks to test our ideas and implementations. The specifications for the benchmarks fall in TCTL. Thus we can also carry out experiments with our abstraction techniques.

* PATHOS real-time operating system scheduling specification:
In the system, each process runs with a distinct priority in a period equal to the number of processes. The biggest timing constant used is equal to the number of processes. The unbounded inevitability property we want to evaluate is that ”if the process with lowest priority is in the pending state, then inevitably it will enter the running state thereafter.” For a system with three processes, this property is as follows.

$`\mathrm{}\left(\text{pending}_3\mathrm{}\text{running}_3\right)`$

The nesting depth of the modal-operators is one.
* Leader election specification:
Each process has a local pointer parent and a local clock. All processes initially come with its parent = NULL. Then a process with its parent = NULL may broadcast its request to be adopted by a parent. Another process with its parent = NULL may respond. The process with the smaller identifier will become the parent of the other process in the requester-responder pair. The biggest timing constant used is $`2`$. The unbounded inevitability we want to verify is that eventually, the algorithm will finish with a unique leader elected. That is

$`\mathrm{}(\text{parent}_1=\text{NULL}i:i1,(\text{parent}_i\text{NULL}\text{parent}_i<i))`$

There is no nested modal-operators. To guarantee the inevitability, we assume that a process with parent = NULL will finish an iteration of the algorithm in 2 time units.
* CSMA/CD protocol :
Basically, this is the ethernet bus arbitration protocol with the idea of collision-and-retry. The timing constants used are 26, 52, and 808. We have used three TCTL specifications for these benchmarks.

+ The first one requires that, when two processes are simultaneously in the transmission mode, then in 26 time units, the bus will inevitably go back to the idle state. The property is a bounded inevitability and can be written as follows:

    $`\mathrm{}((\text{transm}_1\text{transm}_2)x.\mathrm{}(x<26\text{bus\_idle}))`$ (A)

    Note that the inevitability is timed to happen in 26 time units. This experiment allows us to observe how our techniques perform with bounded inevitability.
    + The second specification requires that if sender 1 is in its transmission mode for no less than 52 time units, then it will inevitably enter the wait mode.

    $`\mathrm{}((\text{transm}_1x_152)\mathrm{}\text{wait}_1)`$ (B)

    Specially, this specification can only be verified by quantifying only on non-Zeno computations.
    + The third specification requires that if the bus is in the idle mode and later enters the collision mode, then it will inevitably go back to the idle mode. This unbounded inevitability property is as follows.

    $`\mathrm{}\left(\text{bus\_idle}\mathrm{}\left(\text{bus\_collision}\mathrm{}\text{bus\_idle}\right)\right)`$ (C)

    This property is special in that the nesting depth of modal-operator is two and can give us some insight on how our abstraction techniques scale to the inductive structure of specifications.

## Appendix E Performance data w.r.t. $`d`$-values

The performance data with various $`d`$-value setting in the evaluation of greatest fixpoint with non-Zeno computation requirement is in table 4.

## Appendix F Space complexity charts w.r.t. $`d`$-values

The charts for memory complexity with various $`d`$-parameter values is in figure 2.


text:
# The easy axis of the magnetic response of the conduction electrons of yttrium and rare earths

## Abstract

We describe a scheme for first-principles calculations of the static, paramagnetic, spin susceptibility of metals with relativistic effects such as spin-orbit coupling included. This gives the direction for an applied modulated magnetic field to maximise its response. This easy axis depends on the modulation’s wave-vector $`𝐪`$. For h.c.p. yttrium we find the peak response at a $`𝐪=(0,0,0.57)\pi /c`$, coincident with a Fermi surface nesting vector, to have an easy axis perpendicular to $`𝐪`$. This is consistent with the helical anti-ferromagnetic order found in many dilute rare earth-$`Y`$ alloys. Conversely, the easy axis for the response to a uniform magnetic field lies along the c-axis. The conduction electrons’ role in the canting of magnetic moments in $`GdY`$ alloys and other rare earth materials is mooted.

Conduction electrons play an important role in the magnetism of rare earths. Specifically they determine the Ruderman-Kittel-Kasuya-Yosida (RKKY) indirect exchange interactions which connect the localised 4f magnetic moments of the materials. RKKY  The magnetic response of the conduction electrons, the wave-vector dependent susceptibility $`\chi (𝐪)`$, governs the ordering of the moments. If the greatest value of $`\chi (𝐪)`$ is for $`𝐪=(0,0,0)`$, the moments order ferromagnetically. A peak for $`\chi (𝐪)`$ at some finite $`𝐪`$, on the other hand, signifies a more complicated magnetic structure with a modulation wave-vector equal to $`𝐪`$. In many of the heavier rare earths, e.g. Tb, Dy, Ho and Er the incommensurate anti-ferromagnetic (AF) ordered states, modulated along the c-axes of these h.c.p. crystalline elements, depend strongly on the Fermi surface topology of the conduction electrons. The salient part of the Fermi surface (FS) that causes this anisotropy is the webbing feature that contains flat parallel sheets perpendicular to the $`c`$-axis forming a strong nesting effect Keeton\_Loucks ; Loucks .

The link between FS nesting and the $`𝐪`$-dependence of $`\chi (𝐪)`$ is well-known and has been explored by electronic structure calculations. Liu+Yttium\_chi  There is however another anisotropy in the interaction between the rare earth (RE) magnetic moments that we suggest is also traceable to $`\chi (𝐪)`$. This comes from spin-orbit coupling effects on the conduction electrons and affects whether the localised moments lie perpendicular to the c-axis and thus form helical AF magnetic states or instead have components along the c-axis and order into b-axis modulated AF phases. Of course the anisotropic crystal field acting on the 4f-shell plays an important role in this aspect but the similarity between Gd alloys in which the crystal field effects are small and alloys of Dy,Tb, etc suggests a generic feature from the conduction electrons. This is the issue we address in this letter. We describe a ‘first-principles’ theoretical formalism for the magnetic response of paramagnetic metals in which all relativistic effects such as spin-orbit coupling are included. For the first time the easy axis and its dependence upon wave-vector $`𝐪`$ for the magnetic response can be calculated. In effect we elucidate a relativistic RKKY interaction. jbs-rkky

We present specific calculations of the magnetic response of the transition metal yttrium. With the same h.c.p. crystal structure and electronic configuration as the rare earths, apart from the f-electrons, it provides an excellent model for their conduction electrons. Much is being learnt about the magnetism of rare earths from the varied and complex magnetic structures that they form in multilayers with yttrium. Majkrzak  These range from canted, antiphase domain structure in $`Gd/Y`$ multilayers to coherent, incommensurate magnetic helices found in $`Dy/Y`$ and $`Ho/Y`$ structures. The mediating role of the conduction electrons of $`Y`$ for these rare earth magnetic interactions is crucial in these tailored systems. Yafet  High quality samples of yttrium are more easily come by than the heavy rare earths and recent experiments Dugdale  have measured directly the ‘webbing’ feature of its Fermi surface that it has in common with several rare earths. Moreover the nature of the magnetic ordering of low concentrations of rare earth dopants in $`Y`$ can be directly related to this study. YGd ; Child ; Bates

We begin by considering a paramagnetic metal subjected to a small, external, inhomogeneous magnetic field, $`\delta 𝐛^{ext}(𝐫)`$, and obtain an expression for the induced magnetisation $`\delta 𝐦(𝐫)`$. We use relativistic density functional theory (RDFT) Rajg+Call  to treat the interacting electrons of the system and derive an expression via a variational linear response approach Vosko+Perdew ; Applied\_to\_TM ; staun1 . Although there are a number of non-relativistic studies of this type SW+Sav  here we include relativistic effects and pay particular attention to the magnetic anisotropy of the response. From our RDFT starting point we make a Gordon decomposition of the current density Rajg+Call  and retain the spin-only part of the current, namely the spin magnetisation $`𝐦(𝐫)`$. This results in a spin-only version of RDFT. Rajg+Call  in which the self-consistent solution of Kohn-Sham-Dirac equations is sought, i.e.

$`[c\stackrel{\mathbf{~}}{𝜶}.\widehat{𝐩}+\stackrel{~}{\beta }mc^2+\stackrel{~}{\mathrm{𝟏}}V^{eff}[\rho ,𝐦]\stackrel{~}{\beta }\stackrel{\mathbf{~}}{𝝈}.𝐛^{eff}[\rho ,𝐦]\epsilon ]`$
$`\times G(𝐫,𝐫^{};\epsilon )=\stackrel{~}{\mathrm{𝟏}}\delta (𝐫𝐫^{})`$ (1)
which describes the motion of a single electron through effective fields and $`\stackrel{\mathbf{~}}{𝜶}`$ and $`\stackrel{~}{\beta }`$ are Dirac $`4\times 4`$ matrices. $`G(𝐫,𝐫^{};\epsilon )`$ is the one electron Greens function and the charge $`\rho (𝐫)`$ and magnetisation densities $`𝐦(𝐫)`$ can be written in terms of it i.e.

$$\rho (𝐫)=Tr.𝑑\epsilon f(\epsilon ,\mu ,T)\frac{Im}{\pi }G(𝐫,𝐫;\epsilon )$$

$$𝐦(𝐫)=Tr.\stackrel{~}{\beta }\stackrel{\mathbf{~}}{𝝈}𝑑\epsilon f(\epsilon ,\mu ,T)\frac{Im}{\pi }G(𝐫,𝐫;\epsilon )$$
(2)
where $`\mu `$ is the chemical potential, $`T`$ the temperature and $`f(\epsilon ,\mu ,T)`$ the Fermi-Dirac function. These expressions can be converted into sums over fermionic Matsubara frequencies $`\omega _n=i(2n+1)\pi k_BT`$ FW . The effective potential $`V^{eff}[\rho ,𝐦]`$ consists of the usual combination of external potential (from the lattice of nuclei), the Hartree potential and functional derivative of exchange-correlation energy $`E_{xc}[\rho ,𝐦]`$ with respect to $`\rho `$ whilst the effective magnetic field $`𝐛^{eff}[\rho ,𝐦]`$ is the sum of any external magnetic field (from magnetic impurities, for example) and the functional derivative of $`E_{xc}`$ with respect to magnetisation. We use the local density approximation (LDA) vBH  for $`E_{xc}`$. The leading relativistic effects contained in the Kohn-Sham-Dirac hamiltonian of eq.(1) are the well-known mass-velocity, Darwin and spin-orbit coupling effects.

If a small external field $`\delta 𝐛^{ext}`$ is applied along a direction $`\widehat{𝐧}`$ with respect to the crystal axes of a paramagnetic system, a small magnetisation, $`\delta 𝐦(𝐫)`$, and effective magnetic field, $`\delta 𝐛^{eff}`$ are set up. The effective magnetic field is given by $`\delta 𝐛^{eff}[\rho (𝐫),𝐦(𝐫)]=\delta 𝐛^{ext}(𝐫)+I_{xc}(𝐫)\delta 𝐦(𝐫)`$ where $`I_{xc}(𝐫)`$ is the functional derivative of the effective exchange and correlation magnetic field (within the LDA) with respect to the induced magnetisation density. The Green’s function satisfying equation (1) can be expanded in a Dyson equation in terms of the unperturbed Green’s function, $`G_o(𝐫,𝐫^{};\epsilon )`$ of the paramagnetic system ($`\delta 𝐛^{eff}`$ = 0) and perturbation $`\stackrel{~}{\beta }\stackrel{\mathbf{~}}{𝝈}.\delta 𝐛^{eff}`$ and the first order terms enable the magnetic response function to be obtained.

For a general crystal lattice with $`N_s`$ atoms located at positions $`𝐚_l`$ ($`l`$ = 1,..,$`N_s`$) in each unit cell, a lattice Fourier transform can be carried out over lattice vectors $`\{𝐑_i\}`$. This can be written

$`\chi ^{\widehat{𝐧}}(𝐱_l,𝐱_l^{}^{},𝐪)=\chi _o^{\widehat{𝐧}}(𝐱_l,𝐱_l^{}^{},𝐪)`$
$`+{\displaystyle \underset{l^{\prime \prime }}{\overset{N_s}{}}}{\displaystyle \chi _o^{\widehat{𝐧}}(𝐱_l,𝐱_{l^{\prime \prime }}^{\prime \prime },𝐪)I_{xc}(𝐱_{l^{\prime \prime }}^{\prime \prime })}`$
$`\times \chi ^{\widehat{𝐧}}(𝐱_{l^{\prime \prime }}^{\prime \prime },𝐱_l^{}^{},𝐪)d𝐱_{l^{\prime \prime }}^{\prime \prime }`$ (3)
where the $`𝐱_l`$ are measured relative to the positions of atoms centred on $`𝐚_l`$. The non-interacting susceptibility of the static unperturbed system is given by

$`\chi _o^{\widehat{𝐧}}(𝐱_l,𝐱_l^{}^{},𝐪)=`$
$`(k_BT)Tr\stackrel{~}{\beta }\stackrel{\mathbf{~}}{𝝈}\widehat{𝐧}{\displaystyle \underset{n}{}}{\displaystyle \frac{d𝐤}{\nu _{BZ}}G_o(𝐱_l,𝐱_l^{}^{},𝐤,\mu +i\omega _n)}`$
$`\times \stackrel{~}{\beta }\stackrel{\mathbf{~}}{𝝈}\widehat{𝐧}G_o(𝐱_l^{}^{},𝐱_l,𝐤+𝐪,\mu +i\omega _n)`$ (4)
The integral is over the Brillouin zone with wave vectors $`𝐤`$, $`𝐪`$ and $`𝐤+𝐪`$ within the Brillouin zone of volume $`\nu _{BZ}`$. The sum is over the fermionic Matsubara frequencies. The Green’s function for the unperturbed, paramagnetic system containing the band structure effects is obtained via relativistic multiple scattering (Korringa-Kohn-Rostoker, KKR) theory Faulkner+Stocks . We solve equation (3) using a direct method of matrix inversion. The full Fourier transform is then generated

$`\chi ^{\widehat{𝐧}}(𝐪)=(1/V){\displaystyle \underset{l}{}}{\displaystyle \underset{l^{}}{}}e^{i𝐪.(𝐚_l𝐚_l^{}^{})}`$
$`\times {\displaystyle }d𝐱_l{\displaystyle }d𝐱_l^{}^{}e^{i𝐪.(𝐱_l𝐱_l^{}^{})}\chi ^{\widehat{𝐧}}(𝐱_l,𝐱_l^{}^{},𝐪,\theta ,\phi )`$ (5)
where $`V`$ is the volume of the unit cell. Some aspects of the numerical methods used to evaluate equations (3-5) of this type can be found in staun1 . Note that this expression for the non-interacting susceptibility can be shown to be formally equivalent to one of the familiar type

$`\chi _o^{\widehat{𝐧}}(𝐪){\displaystyle 𝑑𝐤\underset{j,j^{}}{}\frac{|M(𝐤,𝐤+𝐪,\widehat{𝐧})|^2f_{𝐤,j}(1f_{𝐤+𝐪,j^{}})}{\epsilon _j^{}(𝐤+𝐪)\epsilon _j(𝐤)}}`$ (6)
where $`j`$ is an electronic band index, $`\epsilon _j(𝐤)`$ a single electron energy, $`M`$ is a matrix element and $`f_{𝐤,j}=f(\epsilon _j(𝐤),\mu ,T)`$, the Fermi Dirac function.

The important feature of the response function (equations 3-5), is its dependence on the direction of the magnetic field, $`\widehat{𝐧}`$, which vanishes when relativistic, spin-orbit coupling effects are omitted. The direction, $`\widehat{𝐧}`$=(sin$`\theta `$cos$`\phi `$,sin$`\theta `$sin$`\phi `$,cos$`\theta `$) is defined by polar and azimuth angles ($`\theta ,\phi `$). In an h.c.p system such as $`Y`$ or $`Sc`$, if an external magnetic field is applied along the $`\widehat{𝐧}=(0,0,1)`$ direction, i.e. the c-axis, where ($`\theta ,\phi `$$``$ (0,0), $`\chi _o^z`$ is produced. On the other hand $`\chi _o^x`$ is the response of the system when the field is applied in the $`ab`$-plane, $`\widehat{𝐧}=(1,0,0)`$, and ($`\theta ,\phi `$$``$ ($`\pi /2`$,0). We obtain an anisotropy as the difference in the non-interacting susceptibility when an external magnetic field is applied in two directions with respect to the crystal axes, i.e. ($`\chi _o^{\widehat{𝐧}}\chi _o^z`$). This is enhanced by exchange and correlation effects (eq.3). The approach presented here is applicable to ordered compounds and elemental metals and can be modified to study disordered alloys staun1  owing to its KKR multiple-scattering framework. In order to gauge the importance of these relativistic effects with atomic number, we compare our calculations of the magnetic response of $`Y`$ with its lighter 3d counterpart $`Sc`$.

We use atomic sphere approximation (ASA), effective one-electron potentials and charge densities in the calculations for $`Y`$ and $`Sc`$ with experimental lattice constants $`a=6.89`$, $`c=10.83`$ and $`a=6.24`$, $`c=9.91`$ respectively, in atomic units, $`a_0`$Pearson  The details of the electronic structures using a fully relativistic KKR method compare well with those from full potential calculations. Wien97  The Fermi surface for $`Y`$ in the H-L-M-K plane is shown in Fig. 1(a). It shows two relatively flat parallel sheets. The nesting vector $`𝐪_{inc}=\pi /c(0,0,0.57)`$ is indicated by the arrow. Fig. 1(a) is in very good agreement with previous calculations Loucks  and experiments Dugdale ; Vinokurova . Dugdale et alDugdale  have recently carried out positron annihilation fermiology experiments and measured $`𝐪_{inc}=\pi /c(0,0,0.55\pm 0.02)`$ for $`Y`$. Also, Vinokurova et alVinokurova , measured $`𝐪_{inc}\pi /c(0,0,0.58)`$. We find a similar Fermi surface for Sc with the same nesting vector $`𝐪_{inc}=\pi /c(0,0,0.57)`$, also in good agreement with earlier calculations FlemLoucks .

Our calculations of the enhanced static susceptibility, defined by equation (3), show a peak at this same wave-vector $`𝐪_{inc}=\pi /c(0,0,0.57)`$ for both $`Y`$ and $`Sc`$. This is shown in Fig. 1(b), where we probe wave-vectors along the $`c`$-axis from $`\mathrm{\Gamma }`$ to A. (The special point A, for h.c.p. crystal structures is $`(0,0,\pi /c)`$.) The temperature is set at 100K for these calculations.

Fig. 2(a) shows the anisotropy $`\chi _o^{\widehat{𝐧}}(𝐪,\theta ,\phi )\chi _o^z(𝐪)`$ as a function of $`\theta `$ for $`Y`$ and $`Sc`$ at $`𝐪_{inc}`$ and $`𝐪_o(0,0,0)`$. As expected on symmetry grounds, we find the anisotropy to be invariant in $`\phi `$. This means that the magnetic response is insensitive to the direction along which an external magnetic field is applied in the ab-plane. As $`\theta `$ is increased, an anisotropy is observed which reaches a maximum at $`\theta =\pi /2`$ as shown in Fig. 2(a). In fact, we observe that the anisotropy takes the rather simple form

$$\chi _o^{\widehat{𝐧}}(𝐪,\theta ,\phi )\chi _o^z(𝐪)=(\chi _o^x(𝐪)\chi _o^z(𝐪))|\widehat{𝐪}\times \widehat{𝐧}|^2$$
(7)
for h.c.p. crystal structures Remark  ($`|\widehat{𝐪}\times \widehat{𝐧}|=\mathrm{sin}\theta `$). It is therefore sufficient to show $`\chi _o^x(𝐪)\chi _o^z(𝐪)`$, to determine the easy axes of the magnetic response. Fig. 2(b) shows the weak temperature dependence of the anisotropy of $`Y`$ and $`Sc`$ at both the nesting vector $`𝐪_{inc}`$ and also at $`𝐪_o(0,0,0)`$. As expected, the anisotropy is an order of magnitude larger for $`Y`$ than for $`Sc`$. This difference is a result of spin-orbit coupling being more pronounced in the heavier 4d metal $`Y`$ than in the 3d $`Sc`$. We infer that a still greater but similar anisotropy should be evident in the magnetic response of the conduction electrons in the heavier still RE materials.

It is apparent from Figs.1 and 2 that $`Y`$ shows its strongest response at the wave-vector $`𝐪_{inc}`$ which corresponds to the FS nesting and that the response here is strongest to magnetic fields directed in the basal $`ab`$-plane. Here, then, is an explanation as why $`Gd`$ impurities, denoted typically as S ions (zero angular momentum) and thus free from crystal field effects, order into helical incommensurate AF states. YGd  It is also a strong factor for the similar magnetic order found in other dilute rare-earth alloys such as $`\underset{¯}{Y}Tb`$, $`\underset{¯}{Y}Dy`$ and $`\underset{¯}{Y}Ho`$Child

Fig.1(b) shows that yttrium’s responses to ferromagnetic ($`\chi (𝐪)=(0,0,0))`$) and $`𝐪_{inc}`$-AF modulations are close in strength. Figs.2(a) and 2(b) demonstrate, however, that the easy axis for a FM modulation lies along the c-axis and therefore contrary to that of the AF one. These features are pertinent to the magnetic structure of $`GdY`$ alloys. Pure gadolinium has a ferromagnetic phase with the easy axis along the c-axis so that in the paramagnetic phase at higher temperatures $`\chi (𝐪)`$ peaks at $`𝐪=(0,0,0)`$. Moreover its FS here has been observed not to possess a webbing feature Fretwell . Some recent experiments Fretwell  have shown that adding more than around 30$`\%`$ $`Y`$ to $`Gd`$ changes the topology of the Fermi surface of the paramagnetic phase and the webbing is observed again. This is coincident with the alloys’ forming helical AF states at lower temperatures. Moreover, in this concentration range, neutron-diffraction Bates ; Fretwell  has shown that application of a modest uniform magnetic field along the c-axis leads to a ferromagnetic alignment of the $`Gd`$-moments along the c-axis. Once the magnetic field is switched off the system reverts to its helical AF state. This suggests that the alloys’ $`\chi (𝐪)`$’s and anisotropies are similar to those shown in Fig.1(b) and Fig.2(a) but with the relative peak heights of $`\chi (𝐪)`$ at $`𝐪=(0,0,0)`$ and $`𝐪_{inc}`$ to be finely balanced.

In conclusion, we have described an ab initio theoretical formalism to calculate the relativistic static paramagnetic spin susceptibility for metals at finite temperatures. Since relativistic effects such as spin-orbit coupling are included we can identify the anisotropy or easy axes of the magnetic response. We applied this formalism to the 4d metal $`Y`$. Its enhanced susceptibility displays a peak at the incommensurate wave-vector $`𝐪_{inc}=(0,0,0.57\pi /c)`$ traceable to a FS nesting feature found in the calculated electronic structure and fermiology experiments Dugdale ; Vinokurova . The explanation of magnetic order in rare-earth-$`Y`$ alloys and superstructures in terms of the FS topology of $`Y`$ is well established both on theoretical Liu+Yttium\_chi  and, more recently, experimental grounds. Dugdale ; Fretwell  Results in this letter are fully consistent with these findings but we have also shown how the electronic structure of $`Y`$ influences the canting of the rare earth moments with respect to the crystal axes. The finding that the incommensurate AF order of these moments has a helical structure whilst any potential ferromagnetism has an easy-axis parallel to the c-axis is suggested as a more general property of the conduction electrons of h.c.p. rare-earth materials and has a role in the magnetism of RE/$`Y`$ multilayers.

We acknowledge support from the EPSRC (UK).