Stephen Wright

Kerstin Eder

Using Event-B to Construct Instruction Set Architectures

The Instruction Set Architecture (ISA) of a computing machine is the definition of the binary instructions, registers, and memory space visible to an executable binary image. ISAs are typically implemented in hardware as microprocessors, but also in software running on a host processor, i.e. Virtual Machines (VMs). Despite there being many ISAs in existence, all share a set of core properties which have been tailored to their particular applications. An abstract model may capture these generic properties and be subsequently refined to a particular machine, providing a reusable template for development of robust ISAs by the formal construction of all normal and exception conditions for each instruction. This is a task to which the Event-B [MAV05, Sch01] formal notation is well suited. This paper describes a project to use the Rodin tool-set [ABH06] to perform such a process, ultimately producing two variants of the MIDAS (Microprocessor Instruction and Data Abstraction System) ISA [Wri08,Wri09/1] as VMs. The abstract model is incrementally refined to variant models capable of automatic translation to C source code, which this is compiled to create useable VMs. These are capable of running binary executables compiled from high-level languages such as C [KR88], and compilers targeted to each variant allow demonstration programs to be executed on them.

Introduction

One of the most fascinating aspects of a computing machine is its instructions: a relatively small set of obscure functions that, when combined in large sequences, forms the emergent properties of a computer program. Thus, the same machine may be used to calculate a tax return, play a video game, or keep an otherwise uncontrollable aircraft in the air, only by changing the seemingly unfathomable binary sequence loaded into its memory. However, the lack of intuitive connection between the detail of a computing machine's instructions and the large-scale functionality of a program presents a problem in attempting to define those instructions: how to ensure the machine reliably handles all combinations of instructions, whether functionally correct or erroneous. This problem has only been exacerbated with the appearance of Reduced Instruction Set Computer machines, in which the instruction set is designed purely for efficiency and simplicity, with no regard to easing the burden of a programmer attempting to write programs or understand the rationale of an instruction sequence. This has motivated a method for formally managing this comprehension gap through the use of abstraction (the simplification of a system description to an understandable level), and subsequent refinement (the incremental additional of detail in understandable steps).

Related Work

Since the 1980s there has been much activity in the field of automatic verification of processor microarchitecture (i.e. the hardware used to implement the ISA). Two approaches have been taken: definition of simplified academic examples capable of being studied in their entirety [START_REF] Graham | Formalising the design of an SECD chip[END_REF][START_REF] Hunt | A Verified Microprocessor[END_REF], and verification of fragments of pre-existing commercial processors [START_REF] Brock | Report on the Formal Specification and Partial Verification of the VIPER Microprocessor[END_REF][START_REF] Fox | Formal Specification and Verification of ARM6[END_REF][START_REF] Srivas | Formal Verification of an Avionics Microprocessor[END_REF]. A common pattern emerges throughout these projects. Instruction sets are defined as a series of individual specifications defining the behavior of some or all of the entire instruction set. These target-specific specifications are then used to manually or semi-automatically prove theorems stating the assumed adherence to these behavioural descriptions by the microarchitecture, usually down to the level of gate logic.

In the domain of VMs, the ubiquitous Java Virtual Machine (JVM) has attracted much interest in its formal verification. JVM implementations include an instruction interpreter (which performs the functions of a microprocessor's ISA) as well as other components. Architectural models have been constructed [START_REF] Stark | Java and the Java Virtual Machine[END_REF], and other components formally specified and verified, such as the byte-code verifier [START_REF] Caset | Formal Development of an Embedded Verifier for Java Card Byte Code[END_REF][START_REF] Klein | Verified bytecode verifiers[END_REF] and parts of the instruction interpreter [START_REF] Qian | A Formal Specification of Java Virtual Machine Instructions for Objects, Methods and Subroutines[END_REF]. At implementation level, Event-B has been used to construct a specification of some of the JVM instruction set in order verify an existing microcoded interpreter implementation [START_REF] Evans | Towards the Formal Verification of a Java Processor in Event-B[END_REF].

Such work highlights the need for systematic construction of the specification itself: instruction sets have been defined as a series of individual specifications for each instruction, making any verification effort vulnerable to errors at this point. Such errors are generally only discovered by manual review or mismatch with the very implementations they are intended to check. In [START_REF] Brock | Report on the Formal Specification and Partial Verification of the VIPER Microprocessor[END_REF] "a few" errors in the specification are reported and two discussed in detail. One is a simple type mismatch, but another is more serious in which handling of a possible invalid instruction during instruction fetching is omitted. In [START_REF] Srivas | Formal Verification of an Avionics Microprocessor[END_REF] twenty-eight bugs in a formal specification were found by manual review, in comparison with the two bugs finally found in the implementation

Event-B and Rodin

Event-B [START_REF] Metayer | Rodin Deliverable 3.2 Event-B Language[END_REF] is a formal method that combines mathematical techniques, set theory, and first-order logic. It is an evolution of B-Method [START_REF] Abrial | The B-Book: Assigning Programs to Meanings[END_REF] (now often referred to as Classic-B), which was itself derived from Z notation [START_REF] Spivey | The Z Notation: A Reference Manual[END_REF]. Event-B is a simplification of Classic-B, decomposing machines into small, discrete events explicitly linked to their abstractions. This allows more convenient practical development by encouraging more incremental refinement and allowing easier verification of logical proofs. Event-B's other advantage is its flexibility, both in the notation itself and its supporting tools. For example, Event-B does not strictly define a notation, but a methodology for the construction and analysis of linked logical objects, and the notation seen when using support tools is actually an arbitrary front-end for users familiar with Classic-B and Z notation. Event-B was therefore developed closely alongside a supporting tool-set, Rodin [START_REF] Abrial | An Open Extensible Tool Environment for Event-B[END_REF], which allows flexibility in the presentation and manipulation of an underlying linked database representing a formal model. Rodin has a modular architecture based on the Eclipse framework [Ecl09] with clearly defined Application Programming Interfaces exposing the database and supporting expansion via the addition of Eclipse plug-ins.

Common Properties of Instruction Set Architectures

In spite of a vast number of general-purpose microprocessor ISAs in existence, all share a set of core properties [START_REF] Hennessy | Computer Architecture, A Quantitive Approach[END_REF]. Programs are stored within the machine as a contiguous array of binary values, each locatable via an integer address stored in a special register: the program counter (PC). The selection of the next instruction to be executed ("control flow") is achieved by simple incrementing of the PC, or its overwriting with a calculated value to perform a jump operation. Taking one or two elements of a bank of fixed-size registers (the "register file"), and outputting a result back to it performs calculations. In the case of a destructive operation, the output register may be one of the input registers. Certain ISAs implement specialist instructions with more than two inputs [START_REF]AMD Inc "128-Bit SSE5 Instruction Set[END_REF], but are not targeted at general purpose applications. Calculations consist of bit-manipulations, typically implementing the basic arithmetic operations (i.e. addition, subtraction, multiplication and division) plus a number of applicationoriented operations. Data may be loaded into the register file via "immediate mode" addressing, in which a constant value is incorporated in the executed instruction, moved between registers via "register direct mode" addressing, or transferred from a specified location in a read-writable data memory via "register indirect mode". Register indirect mode may also be used to copy data back to this memory from the register file. For all addressing modes source and destination registers are specified by index values incorporated within the instruction. Variations on this general paradigm exist, such as PC-relative mode [Hit98], which may be considered a form of indirect addressing. The data memory consists of a contiguous array of binary values: large register values may be stored across two or more contiguous memory locations.

In order to execute binary images compiled from the C high level language using existing tools, particularly the GNU Compiler Collection (GCC), the ISA must be further specified to implement multiple thirty-two bit general registers, eight bit wide data memory and all instruction and data memory addressing from a single register [START_REF] Stallman | Using and Porting the GNU Compiler Collection[END_REF].

The MIDAS VM

In order to demonstrate the completeness and usefulness of the modeling technique, a working ISA was developed. Implementation of an existing instruction set and machine architecture was considered: such an approach would allow use of existing support tools and benefit from existing development. Two example machines were considered: an existing Virtual Machine standard, the Java Virtual Machine (JVM), and a typical hardware microprocessor targeted at embedded systems, the Hitachi SH4. Existing architectures and instruction sets contain complexities related to achieving particular requirements not relevant to this application. For example, the JVM contains various features to allow efficient support for the Java high level language [START_REF] Lindholm | The Java Virtual Machine Specification, Second Edition[END_REF], and the SH4 instruction set's use of delay-slotted branch instructions [Hit98], enhances performance in a hardware implementation but increases processor and support tool complexity. Full control over the design allows the instruction set to be reduced to a minimum required for program execution and instruction code formats to be selected that are logically based on the formal model.

The MIDAS (Microprocessor Instruction & Data Abstraction System) specification [Wri09/1] describes two variants of a Modified Harvard Architecture, thirty two bit register ISA with a total of thirty-five instructions in eight orthogonal groups [START_REF] Hennessy | Computer Architecture, A Quantitive Approach[END_REF], and a little-endian memory system. The two variants are a stack-based machine and randomly accessible register array machine [START_REF] Hennessy | Computer Architecture, A Quantitive Approach[END_REF], employing the same numerical values to implement similar instruction functionality, the differences being limited to Register File behavior. Randomly accessible register array machines are now the most common form of hardware-implemented microprocessors [START_REF] Hennessy | Computer Architecture, A Quantitive Approach[END_REF], but stack machines are still used, particularly in the field of true VMs not intended for hardware implementation [START_REF] Lindholm | The Java Virtual Machine Specification, Second Edition[END_REF]. Therefore, demonstration of both architectures is desirable.

The MIDAS instruction groups implement no-operation (one instruction), register-fetch (seven instructions), register-store (four instructions), single-operand calculations (two instructions), dualoperand calculations (thirteen instructions), operand compare (five instructions), control-flow jump (two instructions), and machine halt (one instruction). A Harvard Architecture was selected for its guaranteed prevention of executable corruption by bad data accesses, increasing integrity in safety-critical applications. In order to reduce ISA size and complexity the MIDAS ISA is not optimized for performance: for example the single operand compare-to-zero instruction provided by many typical ISAs is not implemented [Hit98]. The basic instruction is an eight bit field, treated as two four bit sub-fields (nibbles). The instruction group is specified by the most significant nibble (MSN) and the precise operation given by a modifier in the least significant nibble (LSN). Instruction groups and modifiers are derived from the groupings constructed in the formal model, allowing for an efficient decoding scheme to be implemented, as events applying to whole instruction groups need only decode the MSN.

Model Description

High-Level Structure

The Event-B model addresses only those aspects of a microprocessor's functionality needed to specify an ISA. Therefore, only a functional description of the instruction set and its associated aspects are included, such as those parts of the memory system visible to the instruction stream. No consideration of the mechanisms within a typical microprocessor needed to efficiently implement an ISA is necessary. For example the model does not address instruction pipelining, memory caching or speculative execution mechanisms [START_REF] Hennessy | Computer Architecture, A Quantitive Approach[END_REF]. The model consists of a common refinement chain, which is then split into separate refinements for two variants of an example ISA constructed for demonstration of the technique. The refinement structure is summarized in Figure 1. The model consists of twenty-eight stages of common refinement, which may be summarized into seven distinct layers: StateMch, ControlFlow, RegArrayMch, MemMch, ConFlagMch, CalcMch and MidasMch. Each layer confers particular properties on the model, described individually. The MidasStk refinement branch introduces a further five layers of refinement: the MidasReg branch another two. Consideration is given to the order of model refinement in order to allow the model to be re-used with minimal modifications. Some layers of the model provide essential definitions necessary to facilitate later refinements, and therefore these must be positioned earlier in the refinement process. For example, the

StateMch

ControlFlow

RegArrayMch

MemMch

ConFlagMch

CalcMch

MidasMch

MidasReg MidasStk

StateMch layer is required to define the exception condition concept and its triggering of machine deadlock, before any specific exception conditions, such as program counter out-of-range detection, may be defined in the ControlFlow. Other refinements may be performed at any point in the modeling process, and are therefore postponed to as late as possible in order to maximize the flexibility of the model. For example, the refinements specifying the exact calculations supported by the MIDAS VM in CalcMch are positioned prior to only the final MidasMch layer. The decision to place the CalcMch layer after the ConFlagMch layer is made on the assumption that supported calculation operations will be modified more regularly than the mechanism used to implement conditionality. Refinements describing the higher level architectural features of the machine, such as the register and memory address spaces, may be easily identified and placed before implementation-specific features, such as exact instruction code values. Decisions on the relative importance of separate architectural issues are harder to make, and judgment based on projected applications is required. However, such an approach does not consider the consequences of a large number of events being constructed early in the refinement process, and therefore propagating the associated management and proving burden to all subsequent refinements.

The separate modeling of executable and data memories allows its applicability to Harvard or Von Neumann memory architectures (i.e. whether the loaded binary program is visible or write-able in the machine's memory system) [START_REF] Lapsley | DSP Processor Fundamentals[END_REF][START_REF] Lee | Programmable DSP Processors part I and II[END_REF]. In the case of the MIDAS demonstrator, a Modified Harvard Architecture, in which the read-only region is visible but the binary program is not, is selected at the MidasMch layer.

The StateMch Layer

StateMch initially abstracts the ISA to the most trivial possible form by defining a single unguarded event acting on only two state variables: an instruction and a machine status, shown in the simplified event in Equation [START_REF] Wright | Using EventB to Create a Virtual Machine Instruction Set Architecture[END_REF].

Iterate ≙ BEGIN act1: inst :∈ INST act2: status :∈ STATE END (1)
The instruction is a member of INST, defined as an abstract SET in the accompanying context, which represents the complete instruction space of the machine (i.e. both valid and invalid instructions). Thus all aspects of the actual instruction, including its symbolic representation and method of extraction, are abstracted. The StateMch layer then uses four more refinement steps to construct the state machine illustrated in Figure 2, defining the modes of operation of the machine and their possible transitions.

Figure 2: StateMch State Transitions

These modes are LOADING (of the program to be executed), RUNNING (of the loaded program), HALTED (i.e. commanded shutdown of the machine) and FAILED (i.e. detection of an error in the program). The model is initialized to LOADING and may continue in that state until loading is complete, when RUNNING is entered. From RUNNING, RUNNING is re-entered on successful execution of a valid instruction, HALTED is entered by the execution of a valid halt instruction, or FAILED entered by execution of an invalid instruction or unsuccessful execution of a valid instruction. Execution results in the modification of inst by unspecified means. Once HALTED or FAILED is entered, the machine enters explicit deadlock by the definition of guarded events with no actions defined. Thus a mechanism for the refining of execution exceptions and the guaranteed stopping of the machine in this event is defined.

Successive division of the INST SET defines instruction groups affecting the state machine. Statements are introduced to explicitly state that the derived sub-sets are inclusive of the entire parent set (in order to ensure that the entire instruction space is decoded) and that the constructed sub-sets are mutually exclusive of each other (to ensure that instruction sub-sets take only one set of properties). This technique is used throughout the model to hierarchically derive all instruction sub-sets. For example, the Event-B fragment shown in Equation (2) describes the initial decomposition of INST into the valid and invalid instructions within the instruction space. In the example the ValidInst sub-set is explicitly stated as non-empty, whilst the InvalidInst sub-set is not. This represents the possibility of a machine in which all possible instruction values are populated with decoded instructions, but a machine with no instructions cannot occur. The actual numerical values of the instruction sub-sets is a feature of a particular ISA, and is therefore postponed to the last refinement stages of the model. (3)

The ControlFlow Layer

In StateMch, the selection of the instruction is abstracted by the apparent modification of a single instruction variable inst. ControlFlow refines this abstraction to a mechanism common to most computing machines: selection of the instruction from a fixed array instArray via a modifiable indexing variable instPtr, representing the PC. The invariants establishing this refinement are given in Equation (4).

ControlFlow uses four more refinement steps to construct the mechanisms by which the PC may be modified (shown in Figure 3), or its attempted modification to illegal values detected. Two instruction groups are defined by their action on the PC: IncrInst which increments the PC by a fixed size, and FlowInst which allows it to be overwritten by a computed value (i.e. a jump or branch). Events guarded by FlowInst are further refined to introduce conditionality and separate PC-relative or absolute calculation of the overwriting vector. (5) Definition of the registers as an array allows the greatest flexibility for the construction of different register architectures in later refinements. Conventional randomly accessible arrays may be constructed by defining indexes free to select any value in the domain RegArrayDom. Fixed, specialist registers may also be constructed by the selection of an element via a fixed index and subsequently refined to a separate data element if desired. Stack architectures may be constructed by the refinement of the selecting index to a global stack-pointer variable.

An instruction sub-set allowing the modification of this newly introduced state is initially constructed by dividing the normal, PC-incrementing instruction using the previously described method, creating NullInst and RegWriteInst sub-sets. This method of division of NullInst is used throughout the model to construct instructions capable of modifying newly introduced state variables. Further refinement steps within RegArrayMch sub-divide RegWriteInst into different operation types, summarized in Figure 4. FetchInst imports data from outside the register file. MoveInst transfers data between locations in the register file unchanged. SingleOpInst populates a member of the register file with the result of a function that takes a single member of the file as input. DualOpInst populates a member of the register file with the result of a function that takes two members of the file as input.

Figure 4: RegArrayMch Data Movements

FetchInst is further refined to specify three possible external sources of data: the numerical value of the PC represented as data, an item of immediate data (i.e. data embedded within the instruction itself), and data from the as yet undefined memory system.

The MemMch Layer

MemMch introduces another new state variable memArray, representing the machine's memory system. The memory is again defined as an indexed array of indeterminate type DATA by the invariant shown in Equation (6).

memArray ∈ MemArrayDom → DATA (6)
Division of NullInst is again used to construct a modifying instruction MemWriteInst. Unlike the multiple refinements anticipated for the regArray state, refinement of memArray to only a single paradigm is expected: a randomly accessible data array, indexed from a selected element of the register file (see Figure 5). This reflects the universality of this pattern in most computing machines [START_REF] Hennessy | Computer Architecture, A Quantitive Approach[END_REF] and is incrementally constructed by the first five refinement steps of MemMch.

The final four refinement steps construct a function performed by most computing machines: the storage of thirty-two and sixteen bit data structures across multiple locations of an eight bit memory system [START_REF] Hennessy | Computer Architecture, A Quantitive Approach[END_REF]. The details of this refinement are given in Section 7.4.

The ConFlagMch Layer

ConFlagMch introduces a new boolean state variable conFlag to refine the abstract conditional defined in the jump instruction constructed in the ControlFlow layer. Use of a separately modified boolean flag is representative of some ISAs [Hit98] but not all [START_REF] Patterson | Computer Organization and Design: The Hardware/Software Interface[END_REF]. Therefore, this may be regarded as the point at which the refinement process begins to specify a particular ISA from the much more general model constructed in the previous layers.

Division of NullInst is again used to construct an abstract instruction capable of modifying the flag, which is then split into instructions performing comparisons of two register file elements: equal, integergreater-than, integer-less-than, float-greater-than and float -less-than.

The CalcMch Layer

CalcMch continues the specification of an exact ISA by the splitting of SingleOpInst and DualOpInst defined in RegArrayMch into the actual functions supported by the ISA. SingleOpInst is refined to integer-to-floating-point cast and floating-point-to-integer cast. DualOpInst is refined to addition, subtraction, multiplication and division for both integer and floating point, bit-wise OR, AND, XOR, bitmap shift-left and shift-right.

The MidasMch Layer

MidasMch specifies the remaining details particular to the MIDAS ISA, excepting the register file differences between the stack and register variants. Jump vector sources are specified as register elements. Instruction sizes are fixed at constants (although still undefined). Instructions not implemented by MIDAS (i.e. sixteen bit fetch-immediate) are eliminated. Precise numerical values are assigned to each instruction via a mapping function between the INST SET and a range of the natural numbers, shown in Equation (7).

The separate arrays instArray and memArray are merged into a single data array memByte using the gluing invariants shown in Equation (8).

INVARIANTS memByte ∈ MemDom → INST ∀x.x∈InstArrayDom ⇒ instArray(x)=(DataByte2Int;Int2Inst)(memByte(x)) ∀x.x∈MemArrayDom ⇒ memArray(x)=memByte(x) (8)
Thus concrete numerical values may be allocated to the domains of instArray and memArray within the single domain of memByte, defining the memory map of the machine.

The event associated with the LOADING mode is refined to a state machine capable of sequentially loading the instruction and read-only regions of the new contiguous memory space. Memory-write events are refined to define a small region of writable IO within the memory region.

The MidasStkMch Layer

MidasStkMch introduces all refinements specific to the stack variant of the MIDAS register file. The register specifier index for direct mode addressing is initially refined to an offset relative to a datum index, which is subsequently refined to the top value of the stack. A new variable stkSize is introduced and register access events are refined to access regArray relative to it.

The stkSize variable is type defined by the INVARIANT shown in Equation (9), thus enforcing the introduction of guards defending against stack overflow and underflow, and the corresponding error events if these conditions are entered. stkSize ∈ 0..MaxRegArraySize (9)

The MidasRegMch Layer

MidasRegMch introduces all refinements specific to the random access variant of the MIDAS register file. All register indexes are trivially replaced with references to the specifier byte fields within the MIDAS instruction. As the value-range expressible by a byte field exceeds the register domain, events detecting out-of-domain accesses are used. This is in contrast to many microprocessor ISAs in which the instruction field may only express the register domain [Hit98].

A refined event constructing a successful NOP instruction from this layer is subsequently shown in Figure 7 of Section 8.

Model Summary

A list of the refinement stages within the general layers described previously is given in Table 1. Proof Obligations (PO) for all refinement stages are discharged using all of the Rodin proving tools. The increasing necessity for manual interaction during later refinement stages reflects the large number of hypotheses visible to the Rodin automatic proving tools at this point, thus requiring selection by the developer in order to achieve proof. A summary of the proof obligations and the method of discharge is given in Table 2.

Modeling of Data Elements

Modeling of an ISA within a formal environment highlights two common aspects of microprocessor design that greatly increase specification complexity: multiple interpretations of stored data and the fragmentation of large data elements across multiple smaller elements. For example the same thirty two bit long-word may be treated as an integer during an addition operation or a simple bit field during shift or bit-wise AND operations. This long-word will also be stored in memory as four contiguous 8-bit bytes, as illustrated by Figure 6.

Modeling of Data Meaning

The possibility of multiple meanings of a stored data element make the basic Event-B types insufficient, and all data stored within the modeled machine is initially abstracted to a generic SET DATA. This allows simple data transfers within the machine to be modeled at an abstract level. For example, the event shown in Equation ( 10) specifies the transfer of data from memory to the register file during an abstract memoryfetch operation.

However, if a stored data element is used to affect the behavior of the machine within the model, type translation functions must be introduced to explicitly state the relationship between DATA and its local meaning. In hardware design, and in C applications, this translation is implicit, potentially leading to ambiguity of behavior. A typical example is the translation between DATA and the integer set during storage and extraction of the PC, requiring the introduction of the mapping functions shown in Equation (11).

DataInt ⊆ ℤ Data2Int ∈ DATA → DataInt Int2Data ∈ DataInt → DATA (11)
An example showing storage of the PC to the register file using this method is shown in Equation (12). act1: regArray(regIndex) := Int2Data(instPtr) (12)

Modeling of Data Size

Beyond the definition of the abstract DATA type, refinements must be introduced to capture the different data sizes used within a machine: 8-bits (byte), 16-bits (short) and thirty-two bits (long). Therefore DATA is recursively refined into three sub-sets, representing the three sizes:

Modeling of Data Fragmentation

The construction of a thirty two bit DataLong from four eight bit DataByte elements is initially specified by the following abstract mapping function:

DataBytes2DataLong ∈ DataByte × DataByte × DataByte × DataByte → DataLong (15)
Conversely, abstract record accessor mapping functions [START_REF] Evans | A Proposal for Records in Event-B[END_REF] are defined to specify the re-extraction of these DataByte elements:

DataLong2DataByte0 ∈ DataLong → DataByte DataLong2DataByte1 ∈ DataLong → DataByte DataLong2DataByte2 ∈ DataLong → DataByte DataLong2DataByte3 ∈ DataLong → DataByte (16)
The location of each DataByte within a DataLong must then be defined. In the AXIOM given in Equation ( 17), the DataLong implicitly constructed by DataBytes2DataLong is specified as four contiguous DataBytes within a generic array:

∀A,m,n,i,lA∈m‥n+3→DataByte ∧ i∈m‥n ⇒ l=DataBytes2DataLong(A(i)֏A(i+1)֏A(i+2)֏A(i+3)) (17)
This AXIOM is complemented by that given in Equation ( 18), defining the DataByte accessed by each record accessor function:

∀a,b,c,d,l a∈DataByte ∧ b∈DataByte ∧ c∈DataByte ∧ d∈DataByte ∧ l∈DataLong ∧ l=DataBytes2DataLong(a֏b֏c֏d) ⇒ a=DataLong2DataByte0(l) ∧ b=DataLong2DataByte1(l) ∧ c=DataLong2DataByte2(l) ∧ d=DataLong2DataByte3(l) (18)
Collectively these statements define the machine's byte-order, or endianess [START_REF] Hennessy | Computer Architecture, A Quantitive Approach[END_REF], as little-endian.

The arrangement of DataBytes within a DataLong allows similar statements to be derived as THEOREMS for the DataByte ordering within a DataShort:

∀a,b,s a∈DataByte ∧ b∈DataByte ∧ s∈DataShort ∧ s=DataBytes2DataShort(a֏b) ⇒ a=DataLong2DataByte0(s) ∧ b=DataLong2DataByte1(s) (19)

Modeling of Memory-Mapped Data

The machine's memory system is initially modeled as a simple array mapping between an address range MemDom and the generic DATA set:

memArray ∈ MemDom → DATA (20)
Refinement to a more precise specification of the memory system is performed in two stages. The memArray is initially refined to a new array of DataLongs: Proof Obligations generated to prove the maintenance of this gluing invariant during memory-write operations are discharged by instantiation of (18) or ( 19) for all memory locations modified by that write.

memDataLongArray ∈ MemDom → DataLong

Implementation Generation

Event-B and the Rodin tool are intended to support automatic generation of executable source code from sufficiently refined models [START_REF] Butler | Rodin Deliverable D16 Prototype Plug-in Tools[END_REF]. This functionality is not yet part of the current Rodin functionality, and therefore a plug-in extension [START_REF] Shavor | The Java Developer's Guide to Eclipse[END_REF] was developed to support a sufficient subset of Event-B to support the VM project, translating to the C language [Wri09/2]. An example showing the final refinement of the MIDAS NOP instruction, and its translated C implementation is given in Figure 7. Each event is translated to a separate C function returning a boolean signifying whether the event has been triggered. A function is generated to call all event functions in turn until an event is triggered, or signal if no event has been triggered at the end of the machine iteration (i.e. deadlock has occurred).

NopOk

The calling function implicitly introduces determinism into models containing non-deterministic event triggering, as events are run in the same order, defined by their position in the Event-B model, Therefore in the case of multiple events being enabled precedence is always given to earlier events. The translator requires that precise values are assigned to all ranges and codes, and set membership is reduced to direct comparison operations. Range checking is performed to ensure that the implementing C type may contain numerical values and ranges. State variables are disallowed from the right side of actions, in order to prevent use after modification by preceding action-derived statements.

Guard statements are automatically evaluated for one of three possible interpretations: type definition of local parameters, assignments to local parameters, or conditional evaluations. Conditionals are implemented as negations of the basic comparisons enabling early returns from a function, and local parameter assignments are only calculated immediately prior to use. Thus execution is optimized and assignments are only evaluated in a valid context. Comments and instrumentation is inserted to provide traceability between the model and implementation.

MIDAS Demonstration

In order to test the constructed MIDAS VMs, a support environment is provided using conventional coding techniques. The environment provides a mechanism for download of binary images into MIDAS executable memory, and text output via a virtual console. These facilities are integrated with the automatically generated MIDAS implementations using conventional C development tools. Each MIDAS variant was built into 2 different calling environments: a standard environment used to run an error-free binary compiled from C, and a test environment used to execute a series of binary images coded in C, assembler or machine-code to exercise all the events-derived functions of the VM.

C compilers targeted at each MIDAS variant were provided via an appropriate GCC assembler and compiler back-end [START_REF] Stallman | Using and Porting the GNU Compiler Collection[END_REF]. Thus a hand-coded C test-suite could be compiled, loaded and executed, demonstrating the suitability of the MIDAS VM for supporting C programs. The test-suite, which was not developed using Formal Methods, is not a complete test of the C language, but includes the most common constructs and invokes the major executable fragments implemented by the compiler. An assemblercoded bootstrap performs segment initialization and machine shutdown.

Lessons Learned

Experience gained during the development of the generic VM model and MIDAS demonstrators allow future improvements to both Event-B and the Rodin tool to be suggested. Model development using multiple refinement steps requires considerable repetition of event guard and action sections. In larger models refinement can lead to small changes being difficult to identify amongst previously constructed functionality, and as the principle of guard strengthening allows additional guards to be legitimately introduced, erroneous guards may be introduced, potentially leading to unintended deadlock conditions. Such issues may be prevented and detected by improved tool support: more concise statement of model refinements (as now provided by the Event-B EXTENDS statement) and enabledness checking tools (as currently proposed).

Within the Rodin tool, the concept of a stored model database, indirectly editable via specialized tools, allows for efficient storage and traversing of model logic by tool extensions. Improvements in the ergonomics of the editing interface are required to achieve productivity similar to that of conventional text editing interfaces. Search, replace and pasting capabilities within a single view would allow faster, less error-prone development. Early versions of these tools are now available and should greatly enhance the development process.

As shown in Section 6.11, current Rodin automatic proving tools are vulnerable to scaling issues. A particular problem is the inability to select pertinent hypotheses from the large lists generated for large models, requiring developer intervention to select manageable sub-sets for PO discharge.

In common with other logical proving tools, the Rodin provers are susceptible to incorrect discharge of proof obligations due to accidental introduction of contradictory predicates. It is suggested that vacuitychecking techniques be introduced to defend against such conditions, as incorporated in other commercially available tools [START_REF] Beer | RuleBase: Model checking at IBM[END_REF].

Future Work 11.1 Model Checking

The VM model has been developed using incremental refinement techniques, with all generated POs discharged (as summarized in Table 2) in order to prove consistency between refinement steps. However, this analysis does not guarantee correctness of the model itself. In the case of the VM, checking against implicit deadlock states is of particular importance. Tools exist for the checking of model correctness [START_REF] Leuschel | ProB: A Model Checker for B[END_REF], and use of these for checking of the VM model is desirable.

Architecture Management Tools

The construction of larger models would be enhanced by the development of high-level tools to assist in the understanding and navigation of deeply refined models containing large numbers of events, such as MIDAS. For example graphical display of a model's refinement hierarchy allowing the browsing and selection of particular events would greatly enhance productivity.

Prover Scaling Support

Relatively superficial enhancements to existing proving tools could be implemented to overcome the scaling issues already discussed and enhance automatic proving rates. For example, the ability for a developer to add programmable tactics, allowing pattern searching of available hypotheses and automatic addition of new hypotheses would be a powerful enhancement. Such techniques are particularly relevant to applications such as MIDAS, which contain many similar POs requiring repetitive application of relatively simple (currently manual) tactics.

Proof Inspection Tools

Current Rodin proving tools enable the automatic and semi-automatic discharge of the hundreds of POs generated by a model such as MIDAS: an essential capability if such Formal Methods are to be practical in this scale. However, only limited ability to review a discharged proof is currently available, although provided by other tools [START_REF]The B-Toolkit User Manual[END_REF]. This capability would allow the manual checking of proofs if required, and allow a developer to plan proof tactics during manual proving, by inspection of similar existing proofs. Tools listing the hypotheses and axioms used by a proving tool to discharge a PO could provide a useful first step towards a more complete proof description.

Automatic Test Generation

Formally derived models have been recognized as a possible basis for generation of testing criteria [START_REF] Utting | B Practical Model-Based Testing -A Tools Approach[END_REF]. The VM application allows the opportunity for testing to be performed on a deployed machine, test inputs being provided by the loading of appropriate binary executables.

Full C Testing

The MIDAS ISA and compiler have been demonstrated to support a hand-coded example application: testing of the VM against a complete C test-suite is required to fully demonstrate the VM [START_REF] Sherridan | Practical Testing of a C99 Compiler Using Output Comparison[END_REF]. The MIDAS ISA has been specified to demonstrate the Event-B modeling technique without regard to other metrics such as performance: expansion of the ISA to include additional performance-enhancing instructions and features, within the Event-B modeling paradigm discussed here, is possible.

Deterministic Model

The VM model currently allows for non-determinism under certain conditions. For example, two separate events are constructed to raise separate exceptions if either the source or destination are unavailable on an attempted data transfer, but the event triggered in the case of both being unavailable is not defined. Expansion of the model to deterministically enumerate all such conditions would yield a more precise specification of VM behavior under such error conditions, albeit at the cost of greatly increased model size.

Conclusions

Event-B allows the generic properties of binary Instruction Set Architectures to be captured in an abstract model, thus providing a re-usable template for the development of Formally Proved computing machines. The Event-B refinement process allows an incremental structure in this abstract model, maximizing its reusability, and its concretization to a level sufficient for automatic conversion to a usable implementation. Constructed relationships derived within the model may also be used to guide specification of new ISAs.

The Rodin tool enables the management of the multiple refinement stages and Formal Proof analysis required for such a technique, and provides the capability for necessary implementation generation tools to be developed.

The technique has been demonstrated by the construction of such a model, and its refinement to multiple implementations in the form of Virtual Machines capable of running compiled binary images.

Figure 1 :

 1 Figure 1: Model Refinement Structure

  CONSTANTSValidInst // Valid instructions InvalidInst // Invalid instructions AXIOMS ValidInst ⊆ INST // Subset of all instructions ValidInst ≠ ∅ // Some valid instructions must exist InvalidInst ⊆ INST // Subset of all instructions ValidInst ⋂ InvalidInst = ∅ // Mutually exclusive ValidInst ⋃ InvalidInst = INST // Complete coverage(2)An example refinement is shown in Equation (3), in which the detection of an invalid instruction leading the FAILED state being entered is constructed.

  INVARIANTS typedef1: instArray ∈ InstArrayDom → INST typedef2: instPtr ∈ InstArrayDom gluing: inst = instArray(instPtr)

Figure 3 :

 3 Figure 3: ControlFlow Actions

Figure 5 :

 5 Figure 5: MemMch location selection

Figure 6 :

 6 Figure 6: MemMch DATA/byte mappingFigure6illustrates a significant side effect of this storage method for formal specification: modification of such a data element implies modification of the data elements mapped at subsequent locations. Such dualities are exposed to the programmer by the weak typing of the C language [KR88], and must therefore be captured by a model.

  (regIndex) := memArray(memIndex)

  can be contained by DATA DataShort ⊆ DataLong // Short can be contained by Long DataByte ⊆ DataShort // Byte can be contained by Short (13) Construction of these refinements is accompanied by refinement of the type-translation mappings discussed in 7.1:

  establishes equivalence across the entire domain of the abstract array:∀x • x∈MemDom ⇒ memArray(x) = memDataLongArray(x)(22)A second refinement step refines the mapping of each address to a DataLong by refining the memory itself to an array of DataByte in Equation (23), and a gluing INVARIANT in Equation (24) establishing the DataLong/DataBytes mapping described in 7.3.

Figure 7 :

 7 Figure 7: Event-B event and derived C

Table 1 .

 1 Refinement summary.

	Layer	Events POs	Description
	StateMch	2	4	Define top-level single event
	StateMchR1	5	2	Split event by status value
	StateMchR2	6	1	Split RUNNING events into behavior for valid and
				invalid insts
	StateMchR3	7	5	Split good-inst into exec and halt instructions
	StateMchR4	8	2	Split exec-inst into success and failure outcomes
	ControlFlow	9	11 Refine inst to instArray/instPtr. Split exec-inst into
				OK/fail/bad-PC
	ControlFlowR1	12	15 Split exec-inst into increment and control-flow
	ControlFlowR2	13	3	Split flow-inst into conditional/vector available/not
				-available
	ControlFlowR3	15	5	Split flow-inst into conditional true/false
	ControlFlowR4	17	31 Split flow-inst-true into relative and absolute
	RegArrayMch	20	17 Introduce data array. Split incr-inst into null-inst
				and those acting on data array
	RegArrayMchR1	21	4	Refine operation-ok in reg-write to src-OK and
				destination-OK
	RegArrayMchR2	35	70 Refine reg-write into fetch, move, single-op and
				dual-op instructions
	RegArrayMchR3	35	18 Define reg indexes to be immediate data in inst
	RegArrayMchR4	42	47 Refine fetch data into inst-ptr, immediate or
				memory
	MemMch	45	24 Introduce memory. Split null-inst into null and insts
				acting on memory
	MemMchR01	46	5	Fetch-mem insts refined to select element from
				mem via mem-index, guarded by readable flag
	MemMchR02	48	6	Mem-write insts refined to select element from
				mem via mem-index, guarded by readable flag
	MemMchR03	48	17 Refine mem-index to address stored in reg-array.
	MemMchR1	48	20 Refine reg-indexes for mem-indexes as immediate
				data in instruction
	MemMchR2	48	113 Refine all other reg-indexes as immediate data in
				instruction
	MemMchR3	52	129 Refine Data-type elements to DataLong-type. Split
				fetch-immediate to byte/short/long
	MemMchR4	56	57 Split mem-access insts into byte/short/long
	MemMchR5	56	89 Refine mem to DataBytes. Define multi-byte data
				mappings for short/long accesses
	ConFlagMch	59	32 Introduce con-flag. Split null to to null and con-
				flag-write. Refine flow conditionals to con-flag
	ConFlagMchR1	59	16 Refine con-flag write to 2-input compare
	ConFlagMchR2	59	4	Define compare src indexes as immediate data in
				instruction
	ConFlagMchR3	63	5	Split compares into equals/greater-than/less-than
	CalcMch	82	207 Split dual-op insts into add/subtract/multiply etc.
	CalcMchR1	83	24 Split single-ops insts into casts
	MidasMch	83	29 Define flow vector as reg data.

Table 2 .

 2 Proof Obligation summary.

	Layer	Automatic	Manual	Total
		Discharge	Discharge
	StateMch	4	0	4
	StateMchR1	2	0	2
	StateMchR2	1	0	1
	StateMchR3	5	0	5
	StateMchR4	2	0	2
	ControlFlow	25	6	31
	ControlFlowR1	15	0	15
	ControlFlowR2	3	0	3
	ControlFlowR3	5	0	5
	ControlFlowR4	11	0	11
	RegArrayMch	17	4	17
	RegArrayMchR1	4	0	4
	RegArrayMchR2	69	1	70
	RegArrayMchR3	18	0	18
	RegArrayMchR4	45	2	47
	MemMch	21	3	24
	MemMchR01	5	0	5
	MemMchR02	6	0	6
	MemMchR03	17	0	17
	MemMchR1	20	0	20
	MemMchR2	111	2	113
	MemMchR3	33	96	129
	MemMchR4	29	28	57

  String initialization demonstrates the correct use of the read-only data region, length calculations test integer-based looping, and output to the virtual console tests character manipulation. Integer-to-string conversions test integer arithmetic and integer digit display via C switch statements test use of dispatch tables. Passing of function arguments and results test variable passing via stack pushes. Explicit tests exercise floating point arithmetic and casting. Integer field extraction functions test bit-wise shift and masking functions. Nested for and if statements test in-function nesting. Nested calls test deeply stacked function calls and returns.


V A Kudryavtsev
email: v.kudryavtsev@sheffield.ac.uk

N J C S P O O N Er

The expected background spectrum in NaI dark matter detectors and the D A M A result

Keywords: D ark m atter, W IM P s, B ack g ro un d radiatio n, R adio ac tivity, D AM

Detailed Monte Carlo simulations of the expected radioactiv e b ack g round rates and spectra in N aI cry stals are presented. T he ob tained spectra are then compared to those measured in the DAMA/ N aI and DAMA/ L IB R A experiments. T he simulations can b e made consistent w ith the measured DAMA spectrum only b y assuming hig her than reported concentrations of some isotopes and ev en so leav e v ery little room for the dark matter sig nal. W e conclude that any interpretation of the annual modulation of the ev ent rate ob serv ed b y DAMA as a dark matter sig nal, should include full consideration of the b ack g round spectrum. T his w ould sig nifi cantly restrict the rang e of dark matter models capab le of explaining the modulation eff ect.

Introduction

The DAMA group has reported a possible signal from dark matter particles. The fi rst evidence came from the DAMA/ NaI experiment (100 kg of NaI), published in Refs. [1,2]. The most recent papers [3,[START_REF] Bernabei | Instru m. & M eth[END_REF] based on DAMA/ LIBRA measurements using 250 kg of NaI combined w ith the DAMA/ NaI results presented improved statistical signifi cance of the signal and evaluation of possible background sources and systematic uncertainties. The positive identifi cation of the eff ect from dark matter particles is based on observation of the annual modulation of the event rate at low energies w ith period (one year) and phase consistent w ith those expected from dark matter halo models [5,6].

U nlike many other experiments, DAMA/ NaI and DAMA/ LIBRA did not use pulse shape discrimination betw een electron and nuclear recoils to eliminate the background from gamma-rays for their high statistics exposure, relying only on the annual modulation of the total event rate and energy spectrum. In the presence of a non-z ero signal, the event rate in any particular energy bin can be described by the eq uation:

R(E, t) = b(E) + S 0 (E) + S m (E) cos (ω(t -t 0 )) , (1)
w here R(E, t) is the time-dependent event rate in a particular energy bin, b(E) is the time-independent rate of background events, S 0 (E) is the time-independent nonmodulated part of the dark matter signal (average event rate), S m (E) is the amplitude of the modulated part of the dark matter signal, ω is the freq uency of modulation (ω = 2π / T , w here T = 1 year) and t 0 is the time in years for the maximum of the signal. A similar eq uation has been used in the analysis of the DAMA data [2,3] but the time-independent terms b and S 0 of the eq uation w ere combined together. Such an approach, how ever, neglects the correlation betw een the modulated and non-modulated parts of the signal. Interpretation of the DAMA signal in terms of dark matter models carried out by the DAMA group [2,3] and in many other papers (see, for instance [7 , 8, 9] for recent analyses) w ere limited to fi ts to the modulated signal and, in some cases, to reconstruction of the average energy spectrum of the signal S 0 (E).

According to E q . (1), the total average rate of events in each energy bin is the sum of the background rate b(E) and the average rate of events from dark matter particles S 0 (E), the third term being cancelled out w hen averaged over a period of a year or several years. This means that if the background spectrum w ere know n w e w ould be able to check w hether the sum of the reconstructed spectrum of the signal (S 0 (E)) and background (b(E)) matches the measured distribution. It turns out, how ever, that most eff orts so far w ere directed to calculation of the signal (mainly S m (E)) in various models helping w ith interpretation of the measured annual modulation eff ect, but totally neglecting the background contribution to the event spectrum and the information about the dark matter models that can be extracted from the measured spectrum. The background event rate and spectrum can be simulated w ith a reasonable degree of accuracy if the abundances of radioactive isotopes in detector materials are know n. The DAMA Collaboration should have included the evaluation of this spectrum in their results to support their claim of dark matter discovery.

Limiting the data analysis to the modulated part of the spectrum only, and ignoring information from the non-modulated component of the signal and background, severely weakens the conclusions that may be derived from the experiment. Self consistent analysis of the experimental data should include all components of the measured spectrum: background, modulated and non-modulated parts of the signal. In this way no information is missed and a proper fit of the measured spectrum to the sum of predicted distributions is possible. In practice, accurate calculation of the background spectrum in terms of the absolute event rate (events/kg/day/keV) can only be done with high accuracy if concentrations of radioactive isotopes in all materials in and around the detector are precisely known. H owever, any particular source of background radiation has a characteristic spectrum and a combination of spectra from different sources with certain weights should match observations with or without signal. H ence the difference between the measured and background spectra should be equal to the spectrum of dark matter events. Analysing the measured spectrum and simulated background would allow us to constrain the dark matter signal.

We have used our expertise in running NaI(Tl) dark matter experiment NAIAD [START_REF] Alner | The UK Dark Matter[END_REF] and in modelling radioactive background in different types of detectors [11,[START_REF] Carson | Instru m. & M eth[END_REF]13] to simulate gamma-ray production, transport and detection in the DAMA/LIBRA experiment. In this paper we present detailed Monte Carlo simulations of the background spectrum in a modelled setup of the DAMA/LIBRA NaI(Tl) detector. We then discuss the implications that these results may have on interpretation of the signal claimed by the DAMA group. We demonstrate that the energy spectrum measured in the DAMA/LIBRA experiment can hardly include any dark matter signal by (i) subtracting some possible WIMP signals from the measured spectrum and comparing this to the simulated background and (ii) subtracting the simulated background spectrum from the measured one and showing that the difference spectrum has a minimum at low energies which would not be a usual feature for dark matter signals. We also show that some inconsistencies exist between low-energy and high-energy parts of the measured DAMA/LIBRA spectrum.

A possible problem with the background spectrum in the DAMA/NaI experiment in the presence of dark matter signal was pointed out also in Ref. [START_REF] Drees | Dark Matter[END_REF]. Several recent analyses (see, for instance, [7,9]) used the measured spectrum to constrain theoretical models but did not take into account the spectrum of background events that should contribute to the observed rate.

Modelling of the background radiation

To model the DAMA/LIBRA background we have considered three locations of radioactive sources: (i) external sources producing gamma-rays and other particles outside the NaI crystals; ii) internal sources of radiation evenly distributed within the volume of the crystal; and iii) surface sources with radioactive isotopes concentrated within the thin surface layer of the crystal (assumed here to be 50 µm thick).

External background sources

The setup of the DAMA/LIBRA experiment comprising 25 NaI(Tl) crystals, ∼10 kg each, was modelled using the published geometry [3,[START_REF] Bernabei | Instru m. & M eth[END_REF]. F or case (i) the source of background radiation was assumed to be in the 100 g envelopes of the photomultiplier tubes (PMTs) attached to the light guides connected to the two fl at surfaces of the crystals. The exact position of the radiation source (PMTs, light guides or similar) is not important since only the high-energy gamma-rays can reach the crystal and the spectrum of events at low energies (below 30 keV) is fully determined by Compton electrons. The PMT envelopes were populated with radioactive isotopes of 238 U, 232 Th and 4 0 K. The isotopes were allowed to decay in our Monte Carlo code, based on the G EANT4 toolkit (version 9.2) [START_REF] Agostinelli | N u cl. Instru m. & M eth[END_REF] . All particles were generated according to G EANT4 library. All isotopes in the decay chains of 238 U and 232 Th decayed in our simulations in secular equilibrium (the spectra of gamma-rays from uranium and thorium decay chains in secular equilibrium have also been obtained in Ref. [16]). O bviously, only high-energy gamma-rays can reach NaI from large distances, with X -rays, electrons and alphas being absorbed in the surrounding materials.

We have also considered 6 0 Co in copper around the detectors and in a shield, as a possible source of background events. Simulations of 6 0 Co induced events have been carried out in the same way as for other isotopes.

F igure 1 shows the spectra of energy depositions from electron recoils in the crystal with PMTs contaminated with 238 U, 232 Th and 4 0 K. F or normalisation we have used typical concentrations of radioactive isotopes in ultra-low background PMTs produced by ETL and used by DAMA: 30 ppb of uranium and thorium, 60 ppm of potassium. Events produced by 6 0 Co decays in copper are also shown in F igure 1 assuming a decay rate of 6 0 Co of 10 mBq/kg. This rate is significantly higher than typically measured in low-background copper. O nly events in which an energy deposition above 0.5 keV is detected in a single crystal (close to the energy threshold of the DAMA experiment), were included. Multiple scattering events in two or more crystals with energy depositions exceeding 0.5 keV were excluded. If a photon scattered two or more times in one crystal, all energy depositions were summed together. The energy deposition spectra were convolved with G aussian distributions describing the energy resolution of the DAMA/LIBRA detectors [START_REF] Bernabei | Instru m. & M eth[END_REF].

As can be seen from F igure 1 the energy deposition spectra from all decay chains at low energies from external sources are essentially fl at, due to Compton scattering of highenergy gamma-rays in the crystal. Spectra from any other decay outside the crystal have very similar shape at low energies dominated by Compton scattered electrons.

The assumption of secular equilibrium does not affect the shape of the spectrum at low energies: the spectrum from any particular decay has the same fl at shape below 30 keV. The broad peaks at about 150 keV are the well-known back-scatter peaks due to the scattering of photons in materials around the crystal prior to entering the sensitive volume.

Internal background sources

For case (ii) the sources of background radiation were distributed evenly in all the crystals. Figure 2 shows the energy deposition spectra from 238 U, 232 Th and 40 K decay chains, as well as from the decays of some other isotopes ( 129 I and 22 Na) identified as possible sources of background in the DAMA/LIBRA experiment [3,[START_REF] Bernabei | Instru m. & M eth[END_REF]. The spectra were convolved with the energy resolution function. For energy depositions from αparticles we have assumed a quenching factor as measured by DAMA [START_REF] Bernabei | Instru m. & M eth[END_REF].

The calculated spectra were normalised using the measured concentrations of radioactive isotopes or their decay rates as reported by DAMA [START_REF] Bernabei | Instru m. & M eth[END_REF]. The main problem here is that only a range of concentrations in the different crystals can be found in Ref. [START_REF] Bernabei | Instru m. & M eth[END_REF], together with a value determined for one crystal as an example. The DAMA Collaboration reported [START_REF] Bernabei | Instru m. & M eth[END_REF] the ranges for decay rates of the isotopes 40 K, 129 I and 22 Na, as well as for some isotopes from the 238 U and 232 Th decay chains. We used the reported typical concentrations of these isotopes in the crystals to normalise the spectra shown in Figure 2: 5 ppt of 238 U (0.7 to 10 ppt was reported in Ref. [START_REF] Bernabei | Instru m. & M eth[END_REF] assuming equilibrium), 5 ppt of 232 Th (0.5 to 7.5 ppt was reported), 10 ppb of natural potassium (≤20 ppb was reported), 0.2 ppt of 129 I (a similar value was reported but for two crystals only) and 6.46 × 10 -14 ppb of 22 Na. The 232 Th decay chain was reported to be in equilibrium [START_REF] Bernabei | Instru m. & M eth[END_REF] unlike 238 U for which the equilibrium was found to be broken: a factor of 5.5 higher rate was found at the end of the chain indicating higher concentration of 210 Pb than expected from the equilibrium. To plot the uranium spectrum in Figure 2 we assumed equilibrium in the uranium decay chain. Since the spectrum from 238 U chain does not have any features below 30 keV, the assumption of equilibrium should not affect the shape of the spectrum at low energies. Assuming higher concentration of 210 Pb compared to the equilibrium would result in a higher peak at 40-60 keV, inconsistent with the measurements.

Present simulations were carried out using the most recent version 9.2 of GEANT4.

The spectra are slightly different from those obtained with the previous version 9.1p3. We have checked the production of the low-energy X-rays, conversion electrons and Auger electrons in GEANT4 and have found that, despite significant improvements, some low-energy photons and electrons are missing resulting in an energy imbalance and spectrum distortions if the radioactive decays happened inside the crystal. We have corrected this by requiring the conservation of energy, with a missing energy being deposited at the point of the decay. The correction becomes significant when a low-energy beta-decay is accompanied by the emission of low-energy X-rays or electrons that are missing from the GEANT4 library, in particular for 129 I decay.

For case (iii) (radioactive sources distributed in the surface layer), the spectra (not shown here) are found to be similar to case (ii) with some enhancement in the number of events at low energies due to particles escaping from the crystal. This effect, however, does not significantly change the shape of spectra at low energies. The assumption that radioactive sources are located in a thin layer at the crystal surface requires their concentrations to be very high, the enhancement over the bulk concentration being equal to the ratio of bulk to surface volumes. Such high concentration looks unrealistic in a well-controlled low-background dark matter experiment.

3 Implications for DAMA results

DAMA spectrum

In this section we will discuss how the measurements and simulations of the background radiation may help with interpretation of the effect reported by DAMA. Figure 3 shows the energy spectra reported by the DAMA/NaI (open circles) and DAMA/LIBRA (filled circles) experiments. The spectra, as reported [2,3], were corrected for all detection and reconstruction effi ciencies. The energy spectra from the DAMA/NaI and DAMA/LIBRA experiments are slightly different. This may be due to different cuts to remove PMT noise and different effi ciencies in the two data sets at small energies where the rate is dominated by the PMT noise. The difference looks to be significant in the position of the peak at about 3 keV that is due to 40 K decay [START_REF] Bernabei | Instru m. & M eth[END_REF]. The energy spectra measured in other NaI detectors at low energies are different from those observed by the DAMA/NaI and DAMA/LIBRA experiments. The filled triangles show the spectrum of energy depositions as measured in NaI by the NAIAD experiment [START_REF] Alner | The UK Dark Matter[END_REF]. The DAMA spectra show the minimum at about 2 keV that was not seen in NAIAD. However, the absolute rate of events was higher in NAIAD than in DAMA/LIBRA. Event rate of about 2 events/kg/day/keV above 5 keV (similar to DAMA) was observed in Ref. [17] but the rate in the 3 keV peak was higher than in DAMA.

The interaction of WIMPs with the target (Na and I) material is one interpretation of the measured annual modulation effect seen by DAMA. The curve in Figure 3 shows the energy spectrum of events from spin-independent WIMP interactions with parameters M W IM P = 60 GeV and σ S I = 7 × 10 -6 pb in the isothermal halo model. These parameters are close to the best fit to the measured signal as reported in Ref. [2]. Filled squares show the difference between the measured spectrum and the spectrum of the WIMP signal described above. So the pure background spectrum should look like this curve. This 'background spectrum' shows a deep minimum at about 1.5-2.0 keV. The depth and exact position of the minimum depends on the dark matter particle and halo model but the minimum should be present in all models, most of them requiring the ratio of non-modulated to modulated parts of the signal to be at least 20-30 at 2-4 keV. The rise of the measured event rate below 1.5 keV is most likely due to the PMT noise that was not included in our simulations.

The rate at 4-10 keV in the DAMA/LIBRA experiment was reported to be about 1 event/kg/day/keV [3,[START_REF] Bernabei | Instru m. & M eth[END_REF]. The peak at about 3 keV is due to the intrinsic contamination of the crystals by 40 K [START_REF] Bernabei | Instru m. & M eth[END_REF]. The reported rate at 20-30 keV is about 0.5 events/kg/day/keV. It is not clear, however, at what energy and why the event rate drops from about 1 events/kg/day/keV as observed at 4-10 keV to about 0.5 events/kg/day/keV as observed at 20-30 keV. It is diffi cult to explain such a decrease in rate with our simulations. Note, however, that the DAMA Collaboration reported the spectra at 20-80 keV from two individual crystals only, whereas the spectrum below 10 keV is averaged over all 25 crystals. If the majority of the crystals show much higher rate at 20-30 keV, then the rates below 10 keV and above 20 keV might be consistent. Assuming typical concentrations of uranium, thorium and potassium in PMTs, this external source of background does not contribute substantially to the total spectrum. The simulated absolute rate below 10 keV (solid curve in Figure 4) is significantly smaller than observed by DAMA (∼1 event/kg/day/keV, see Figure 3 above and Figure 27 in Ref. [START_REF] Bernabei | Instru m. & M eth[END_REF]). If the copper of the crystal case or shielding, or other materials are contaminated by 60 Co with a decay rate about 500 mBq/kg, then this source of radiation can explain the rate of events at 4-10 keV as measured by DAMA. The broad peak at 40-100 keV is due to the decay of 129 I. The event rate in the peak is consistent with that reported by the DAMA Collaboration [START_REF] Bernabei | Instru m. & M eth[END_REF].

Combined background spectra from simulations and comparison w ith measurements

In an attempt to make the simulated spectrum compatible with the measurements we examined an effect of assuming higher concentration of some radioactive isotopes compared to the reported values. The flat part of the spectrum at low energies could be due to either intrinsic 238 U (with or without equilibrium) or any external source (decay chains of 238 U or 232 Th, 40 K, 60 Co). Each spectrum would result in characteristic features at higher energies.

We have first assumed that the internal contamination of the crystals with radioisotopes is higher than reported by DAMA. The dashed curve in Figure 4 shows the spectrum calculated assuming higher concentrations of certain isotopes inside the crystals compared to the values used in Figure 2: twice as much of natural potassium (20 ppb) to bring the 3 keV peak to the measured value, 8 times more 238 U (40 ppt) to make the spectrum flat at 4-10 keV with a measured rate of about 1 event/kg/day/keV, and 4 times more 232 Th (20 ppt) to achieve a better match at 4-10 keV (see also Figure 5 below). An important feature of this spectrum is the presence of a peak from 210 Pb decay at ∼40-50 keV (β-decay accompanied by gamma-rays and X-rays) which dominates over 129 I decay. If the equilibrium is broken in the uranium decay chain and the decay rate of 210 Pb is higher than that of the parent isotopes as reported by DAMA [START_REF] Bernabei | Instru m. & M eth[END_REF], then the peak from 210 Pb will be even more pronounced making the predicted spectrum inconsistent with the measurements (see Figure 11 of Ref. [START_REF] Bernabei | Instru m. & M eth[END_REF]).

Then we have looked at the possibility of having higher background rate due to external sources. Any external source of background radiation would produce a broad back-scatter peak centered at about 150 keV. The dotted curve in Figure 4 shows the spectrum calculated using higher contaminations of potassium in the PMTs (12000 ppm or 1.2% ) and in the crystals (20 ppb) to account for the flat part of the measured spectrum and absolute rate at low energies (including the peak at about 3 keV). Such a high concentration of potassium (or similarly very high concentrations of uranium or thorium) in low background PMTs looks unrealistic but is used to illustrate the effect of possible external source of background. In fact, such a source could be much closer to the crystal and not shielded by any light guides or copper. In this case the concentration needed to obtain the flat spectrum and the rate consistent with measurements at low energies might be smaller. For instance, as mentioned above, a high decay rate (about 0.5 Bq/kg) of 60 Co in copper around the detectors (case and shielding) could explain the observed rate at 4-10 keV.

An important feature of any spectrum dominated by external source, is the presence of a broad back-scatter peak at about 150 keV. Figure 11 in Ref. [START_REF] Bernabei | Instru m. & M eth[END_REF] shows, however, that the measured spectrum at 40-80 keV (for two individual crystals) can fully be explained by the sum of the 129 I and 210 Pb decays leaving too little room for the back-scattered gamma-rays. A more detailed investigation of the background could be performed if the measured spectrum averaged over all crystals for a broad energy range (from 1 keV to a few MeV) were to be published by the DAMA Collaboration.

The minimum at 1.5-2 keV that should be present in the background spectrum (Figure 3) due to the contribution of the WIMP signal to the observed event rate, is not seen in the simulated spectra (Figure 4) independently of the assumptions made about radioactive contaminations in different materials in and around the detectors.

Figure 5 shows simulated spectra below 10 keV with 'optimised' concentrations to match the measured spectrum at low energies: external background sources were assumed to be in the PMTs with concentrations specified above for Figure 1; internal sources are as for the dashed curve of Figure 2 -potassium at 20 ppb, 129 I at 0.2 ppt, 238 U at 40 ppt and 232 Th at 20 ppt. The simulated spectrum is compared with the measurements without signal subtraction (open squares) and with signal subtraction.

For the latter case we assumed 60 GeV mass WIMPs with σ SI = 7 × 10 -6 pb (filled squares) and σ SI = 2 × 10 -6 pb (open circles). A similar simulated spectrum but over a wide energy range is shown in Figure 4 by the dashed curve. Even for a model with small cross-section, the measured spectrum minus signal (open circles in Figure 5) does not match the predictions at 1.5-2.5 keV. A better match between the simulated background spectrum and the measurements minus signal would require a very small ratio of non-modulated to modulated parts of the signal spectrum (certainly less than a factor of 10 at 2-4 keV).

Our 'optimised' background spectrum, although being in reasonable agreement with the measured one, does not satisfy a crucial requirement in the presence of a signal: the background event rate should not exceed the measured one in any energy bin since a contribution from a signal is expected too. In some energy bins in Figure 5, in particular at 2-3 keV, the calculated background rate exceeds the measured one. This can partly be corrected if we assume an uncertainty in the energy calibration of DAMA/LIBRA. For instance, shifting the experimental data by about 0.25 keV towards low energies can improve the agreement with simulations and opens the way to seeing what the signal spectrum could look like after background being subtracted from the measured rate.

Figure 6 shows the measured DAMA/LIBRA spectrum shifted by 0.25 keV towards low energies. (Note that this is only an assumption and other ways of improving the agreement between the measured spectrum and the background taking into account the uncertainties in energy calibration, resolution, etc are possible.) Let us now require that the expected background spectrum should be less than the measured rate in any energy bin. The red solid curve shows the background spectrum with internal concentrations of 34 ppt of 238 U, 20 ppt of 232 Th and 20 ppb of natural potassium.

The reduction of uranium concentration from 40 ppt to 34 ppt is the consequence of the above requirement. The background spectrum still looks reasonably close to the measured one and leaves practically no room for the signal, the difference between the measured spectrum and the background being shown by filled squares.

To allow for a signal to be present we should now reduce the concentrations of radioisotopes in our calculations. Reducing the internal concentration of uranium to 20 ppt (dashed blue curve in Figure 6 corresponds to the sum spectrum) results in a 'signal' spectrum shown by open circles which is obviously inconsistent with any dark matter signal. Internal concentration of thorium does not affect the spectrum below 6 keV at all (see Figure 5). Lowering the internal concentration of potassium to 10 ppb reduces the background event rate at the peak only and does not affect energies above 4.5 keV.

The difference between the measured and expected background spectra that can be associated with the signal, looks like a peak at about 3 keV (filled triangles) -copying the peak from 40 K. If realistic dark matter models exist with such a feature and with a low ratio (about 10) of non-modulated to modulated parts of the signal spectrum, they should be carefully analysed. There may be, however, a simpler explanation of the annual modulation around 3 keV peak related to the selection of single hit events at these energies which are dominated by 40 K (see also Ref. [START_REF] Gerbier | Talk at the 16th International Conference on Supersymmetry and the Unification of Fundamental Interactions -SUSY 08[END_REF] for discussion). This can be checked by looking at the time behaviour of the rate around 1460 keV and of the coincident events with a 3 keV energy deposition in one crystal and 1460 keV deposition in another one. The problem with a smaller measured rate at 40-60 keV compared to the required one from a higher concentration of 238 U is still present in all these scenarios.

If we do not assume high concentrations of some isotopes as suggested above, the simulated event rate at low energies is much smaller than the DAMA measurements and no dark matter signal can explain the difference above 6 keV. At present it is difficult to explain the apparent inconsistency between the low-energy part (requiring high concentrations of uranium and possibly thorium in the crystals or high rate from an external source) and the high-energy part (not showing associated features) of the DAMA spectrum.

Several fits of the measured modulated part of the signal to the predictions of different dark matter models were shown in Ref. [3]. Unfortunately, no comparison between the measured spectrum and the predictions including the background and the time-averaged signal was presented, although such a comparison would provide a very important test of the interpretation of the observed effect. For instance, most WIMP models that fit the DAMA/NaI data [2,7,8] (some of them are already excluded by other experiments) are not compatible with the combination of background and signal spectra as shown in the figures. Most models require the non-modulated part of the signal to be at least 20-30 times bigger than the modulated part at 2-4 keV resulting in a deep minimum after the signal spectrum is subtracted from the measured one (see Figure 5). Independently of the relative contributions of external sources of different isotopes (no accurate measurements for the DAMA experiment were found), and even varying intrinsic contaminations of the crystals, it is hard to explain such a minimum in the overall background spectrum at 1.5-2 keV, the measured rate in the peak at 3 keV and almost flat rate of 1 event/kg/day/keV at 4-10 keV as seen in Figure 3, with our simulations.

Accurate comparison of the simulations with the measured background spectrum in the DAMA crystals is difficult due to the absence of published data on the spectrum at all energies (from 1 keV to a few MeV). We stress, however, that any analysis of the DAMA signal should include predictions of the background spectrum and its comparison with the measured spectrum minus the averaged dark matter signal extracted from the annual modulation effect.

Our results are practically independent of the assumption about the channeling effect in NaI crystals (see [3] and references therein). Our conclusions assume that dark matter models predict a high ratio of the non-modulated to modulated parts of the spectrum (only a small fraction of spectrum is modulated) which does not depend much on the channeling effect.

Conclusions

We have presented Monte Carlo simulations of the background rates and spectra due to different sources expected in the DAMA/LIBRA NaI crystals. Any specific source of background radiation has a characteristic spectrum and, ideally, the measured spectrum should be consistent with a combination of predicted spectra from the different sources. We have found that it is hard to match the measured spectrum of DAMA/LIBRA at low energies (even without signal) with the simulated background assuming measured concentrations of radioactive isotopes in NaI crystals. Some features (drop in observed intensity between 10 and 20 keV, low total rate from simulations) are difficult to explain and publication of the measured energy spectrum up to MeV energies by DAMA is necessary to clarify this point. If a signal is present, then the background (measured spectrum minus signal) should have a minimum at about 1.5-2 keV. Its position and depth depend on the model but it is very hard to obtain such a minimum with simulations. Significantly enhanced concentrations of certain isotopes in the crystals help to bring the predicted rate and shape of the spectrum to agreement with the measured spectrum at low energies but leave too little room for the dark matter signal. Any dark mater model that predicts a large ratio (above 10 at 2-4 keV) of non-modulated (average) to modulated parts of the spectrum can certainly be excluded. The assumption of higher concentrations of certain isotopes, although improving agreement with the low-energy part of the spectrum, results in some features that are not seen by DAMA at high energies. We do emphasise that any interpretation of the annual modulation of the event rate observed by DAMA, as dark matter signal, is incomplete without clear understanding of the background spectrum. Such integrated analysis would significantly restrict the dark matter models capable of explaining the modulation effect.

Acknowledgments

This work has been supported by the ILIAS integrating activity (Contract No. RII3-CT-2004-506222) as part of the EU FP6 programme in Astroparticle Physics. We acknowledge the financial support from the Science and Technology Facility Council Energy, keV Event rate, events/kg/day/keV

Figure 4 (

 4 Figure4(solid curve) shows the sum spectrum of simulated events in NaI crystals due to internal and external sources of background as plotted in Figures 1 (except 60 Co) and 2. The aforementioned concentrations of the radioactive isotopes were used. Assuming typical concentrations of uranium, thorium and potassium in PMTs, this external source of background does not contribute substantially to the total spectrum. The simulated absolute rate below 10 keV (solid curve in Figure4) is significantly smaller than observed by DAMA (∼1 event/kg/day/keV, see Figure3above and Figure27in Ref.[START_REF] Bernabei | Instru m. & M eth[END_REF]). If the copper of the crystal case or shielding, or other materials are contaminated by60 Co with a decay rate about 500 mBq/kg, then this source of radiation can explain the rate of events at 4-10 keV as measured by DAMA. The broad peak at 40-100 keV is due to the decay of 129 I. The event rate in the peak is consistent with that reported by the DAMA Collaboration[START_REF] Bernabei | Instru m. & M eth[END_REF].

Figure 1 :Figure 2 :Figure 4 :

 124 Figure1: Spectra of energy depositions from electron recoils in the NaI crystals from 238 U, 232 Th and 40 K decay chains in secular equilibrium. The source of radiation was the PMT envelopes (100 g each) attached to the light guides connected to the crystals. Also shown is the spectrum of 60 Co events from Cu. Only events in which a single crystal was hit, are included. See text for details.

Figure 5 :

 5 Figure 5: Simulated spectra of energy depositions in the NaI crystals from internal sources at low energies assuming 'optimised' concentrations of isotopes in the crystals (solid curve): 40 ppt of 238 U (pink dashed-dotted curve), 20 ppt of 232 Th (green solid curve), 20 ppb of natural potassium (blue dashed curve), 0.2 ppt of 129 I (not seen on the graph). The measured spectrum of DAMA/LIBRA without signal subtraction (open squares) and with signal subtraction are also shown. For the latter case we assumed a signal from 60 GeV mass WIMPs with σ SI = 7 × 10 -6 pb (filled squares) and σ SI = 2 × 10 -6 pb (open circles) (see text for details).

Figure 6 :

 6 Figure 6: Simulated spectra of energy depositions in the NaI crystals from internal sources at low energies assuming 'optimised' concentrations of isotopes in the crystals (solid curve): 34 ppt of 238 U, 20 ppt of 232 Th and 20 ppb of natural potassium. The measured spectrum of DAMA/LIBRA (open squares) is shifted by 0.25 keV towards low energies to achieve better agreement with simulations around the 3 keV peak. The difference between the two above spectra is shown by filled squares. The dashed blue line is the calculated spectrum assuming 20 ppt of 238 U, 20 ppt of 232 Th and 20 ppb of natural potassium with the corresponding difference between the measurements and calculations shown by open circles. The difference between the measured spectrum and the simulated one assuming 34 ppt of 238 U, 20 ppt of 232 Th and 10 ppb of natural potassium is denoted by filled triangles (see text for details).

ACCEPTED MANUSCRIPT

-STFC). We thank Prof. G. Gerbier and Dr. R. Lemrani (CEA-Saclay, France) for fruitful discussions and Dr. L. Pandola (LNGS-INFN) for helping us with generating radioactive decays with GEANT4. We are grateful to anonymous referees for useful comments.

 [2](open circles) and DAMA/LIBRA [3] (filled circles) experiments. The spectrum of events expected from 60 GeV WIMP interactions with the spin-independent cross-section of 7×10 -6 pb in the isothermal halo model is shown as example by the solid curve (labeled as '60 GeV WIMPs'). The difference between the measured DAMA/LIBRA spectrum and the WIMP signal is plotted as filled squares (labeled as 'LIBRA-WIMPs'). An example spectrum from one of the NAIAD crystals is shown by filled triangles.