Operating theatre scheduling with patient recovery in both operating rooms and recovery beds

Keywords: operating theatre, surgery scheduling, Lagrangian relaxation, recovery bed

This paper investigates the impact of allowing patient recovery in the operating room when no recovery bed is available. Three types of identical resources are considered : transporters, operating rooms and recovery beds. A fixed number of patients must be planned over a term horizon, usually one or two weeks. The surgery process is modelled as follows : each patient is transported from the ward to the operating theatre. Then the patient visits an operating room for surgery operation and is transferred to the recovery room. If no recovery bed is available, the patient wakes up in the operating room until a bed becomes available. The operating room needs to be cleaned after the patient's departure, before starting another operation. Finally, the patient is transported back to the ward after his recovery. We consider several criteria based on patients' completion times. We propose a Lagrangian relaxation-based method to solve this operating theatre scheduling problem. The efficiency of this method is then validated by numerical experiments. A comprehensive numerical experiment is then performed to quantify the benefit of allowing patient recovery in operating rooms. We show that the benefit is high when the workload of the recovery beds is high.

1 Introduction

Context

Application of operational research methods in health care systems has proven effective for design and operation of the operating theatres. New organisations of health care systems and new health care delivery technologies emerge to better cope with cost pressure while keeping the quality of service and care high. [START_REF] Guinet | Operating theatre planning[END_REF] [START_REF] Guinet | Operating theatre planning[END_REF] described an hospital as a multi-service production system which is constrained by limited material and human resources. Optimal usage of expensive resources is highly important in cost containment. Surgeries are critical processes in hospital operation, not only because of their high costs but also because of their direct impact on patient's health and on quality perception. Therefore the operating theatre scheduling has been widely studied in order to maximize the utilization of operating rooms and enhance the patient's satisfaction. Surgery procedures involve about 9% of total costs of the hospital [START_REF] Gordon | Surgical unit time utilization review : resource utilization and management implications[END_REF] [START_REF] Gordon | Surgical unit time utilization review : resource utilization and management implications[END_REF]. [START_REF] Dexter | Cost implications of various operating room scheduling strategies[END_REF] [START_REF] Dexter | Cost implications of various operating room scheduling strategies[END_REF] showed that one hour overtime cost of any operating room is about 1.75 time of its usual cost.

In an operating theatre, when a surgery operation ends, the patient recovers slowly under appropriate monitoring of competent personal and with necessary equipment. Patient recovery usually takes places in a bed in the recovery room. As a result, an operated patient must be transferred in a recovery bed after the end of the surgery. However recovery process can be longer than expected, and situations where no recovery beds are available occur sometimes. Instead of delaying surgery operations, letting the recovery occurs in the operating room could be a cleaver solution.

Literatury review

The purpose of this paper is to investigate the impact of allowing patient recovery in operating rooms. This strategy has been scarcely addressed in the literature.

The operating theatre scheduling problem consists in scheduling surgery operations over a certain time period (usually a week), taking into account available resources. Operating theatres can be classified into several categories depending on medical equipments, operating rooms features and localisation in the hospital (centralised theatre or dedicated theatre). Depending of the operating theatre type, there are two surgery scheduling strategies in the literature : « block scheduling » and « open scheduling » [START_REF] Jebali | Vers un outil d'aide à la planification et à l'ordonnancement des ressources dans les services de soins[END_REF].

In the block scheduling policy, a preliminary time table called master schedule is established in order to allocate time slots to surgeons, groups of surgeons or medical specialities. Determination of an optimal schedule is a difficult problem [START_REF] Pinedo | Scheduling theory, algorithms and systems[END_REF]. [START_REF] Guinet | Operating theatre planning[END_REF] [START_REF] Guinet | Operating theatre planning[END_REF] proposed a two-step approach where patients are first assigned to operating rooms, and then each operating room is scheduled using an extension of the Hungarian method. [START_REF] Jebali | Operating rooms scheduling[END_REF] [START_REF] Jebali | Operating rooms scheduling[END_REF] studied the same problem by focusing on operating rooms and provided a mixed integer linear programming formulation that they solved with commercial software package. Fei et al. (2006b) [START_REF] Fei | An operating theatre planning and scheduling problem in the case of a "block scheduling" strategy[END_REF] also proposed a two-step approach : a weekly schedule is created using column generation, then the daily scheduling problem is solved using an hybrid genetic algorithm. Finally, [START_REF] Kharraja | Open scheduling : A block scheduling strategybased operating theatre planning[END_REF] [START_REF] Kharraja | Open scheduling : A block scheduling strategybased operating theatre planning[END_REF] studied using simulation a problem similar to ours where patient recovery in operating room was allowed but patient transportation was not taken into account.

In the open scheduling policy, no pre allocation is used, patients are scheduled without any specialty related restriction. In this context, operating theatre scheduling with constraints of patient recovery and transportation is scarcely addressed. Most existing models do not consider patient recovery in operating rooms. The work most related to this paper is Fei et al. (2006a) [START_REF] Fei | Ordonnancement journalier dans un bloc opératoire dans le cadre d'une stratégie open scheduling[END_REF], where a daily surgery schedule is built for an operating theatre. Operating rooms and recovery beds are identical. Patient recovery in operating rooms is allowed. No transporters are considered. A hybrid genetic algorithm was proposed to solve the problem. [START_REF] Marcon | Determining the number of beds in the postanesthesia care unit : A computer simulation flow approach[END_REF] [START_REF] Marcon | Determining the number of beds in the postanesthesia care unit : A computer simulation flow approach[END_REF] developed a simulation model to find an optimal scheduling strategy, taking into account the whole surgery process. A mathematical model is proposed to minimize the overtime for operating rooms. The simulation model is used to model the patient pathway throughout the operating theatre. Bottleneck resources of the system are the transporters. Emergencies have been considered in the operating theatre by [START_REF] Lamiri | A stochastic model for operating room planning with elective and emergency demand for surgery[END_REF] in [START_REF] Lamiri | A stochastic model for operating room planning with elective and emergency demand for surgery[END_REF] : a stochastic scheduling model has been presented in this article. In [START_REF] Lamiri | Column generation for operating theatre planning with elective and emergency patients[END_REF] the same problem is solved using a column generation procedure.

The modified block scheduling policy is an intermediate strategy which combines the open and block schedulin strategies. First the weekly planning is divided into time slots. Two strategies are commonly used : (i) some time slots are allocated to surgeons, surgeons teams or medical specialities, and the remaining time slots are common (unassigned block) ; (ii) after a given time (block release time), unaffected time slots are re-assigned in order to maximize the operating theatre utilisation [START_REF] Marcon | Survol des différentes stratégies de programmation opératoire : Analyse de leurs performances et des outils informatiques de leur mise en oeuvre[END_REF].

Scientific contribution

In this paper we focus on the open scheduling policy. The operating theatre scheduling problem considered here is similar to the one considered in (anonymous ref.

) and both take into account not only availability of operating rooms but also availability of porters and recovery beds. In contrast to the previous work which does not allow recovery in the operating room, patient recovery in the operating room is allowed in this paper when no recovery bed is available at the end of the surgery. This situation occurs when the when workload of the recovery beds is high with respect to the workload of operating rooms, i.e. when the operating theatre is crowded : a large number of emergency patients have to be operated in regular operating rooms could be a reason to allow recovery in operating rooms if the recovery room is full. Furthermore allowing patient recovery in operating rooms could be a cleaver strategy when recovery process times grow high, i.e. when workload of recovery beds becomes greater than that of operating rooms. We consider here critical situations in the operating theatres that often occur according to surgeons' expertise. A Lagrangian relaxation method is proposed. A comprehensive numerical experiment is performed to investigate the benefit of allowing patient recovery in operating rooms. In particular, we show that the benefit is significant when the workload of the recovery beds is high with respect to the workload of operating rooms.

Notice that cost of operating room utilization during one hour is two to five times greater than cost of recovery bed utilization during the same amount of time, depending on surgery type. From this point of view, allowing recovery in operating rooms during patient scheduling may increase surgery costs. However we focus in this problem on makespan criterion (minimize operating theatre closing time), in order to process as much patients as possible during a short amount of time. Taking into account costs comparison may be examined in perspectives of this work.

The rest of the paper is organized as follows. The problem setting is introduced in Section 2. Section 3 presents a mathematical model and Section 4 describes the Lagrangian relaxation approach. Numerical results are given and analyzed in Section 5.

Problem setting

The surgery process followed by a patient can be divided in several steps ; we consider that a patient must go through the surgery process on a fixed day, at a fixed hour.

1. The patient is transported from the ward to the operating theatre by a porter team.

2. If no operating room is available, the patient must wait in a waiting room. We consider the surgeon, the nurses, the medical equipment and the operating room as an unique resource.

3. When an operating room becomes available, the surgery operation starts.

4. At the end of this operation, the patient is transferred to a recovery bed in the recovery room.

If no recovery bed is available, the patient recovery begins in the operating room which becomes blocked. The patient leaves the operating room when the recovery is completed or when a recovery bed becomes available. In the latter case, the patient is immediately transferred to the recovery room.

5. The operating room is cleaned as soon as the patient leaves.

6. When the patient has recovered from his surgery operation, he must wait for a porter team to be transferred to the ward. The entire process is illustrated in figure 1. Taking into consideration the processing time of each activity, the surgery and recovery operations are the most critical activities of the operation process. In this paper, we assume deterministic processing times for all operations although surgery time for example is hard to predict and directly affects the utilization of the most critical resources, the operating rooms. Operating time estimation depends on surgeon expertise and on patient pathology [START_REF] Wright | Statistical modeling to predict elective surgery time[END_REF] [START_REF] Wright | Statistical modeling to predict elective surgery time[END_REF].

The problem studied in this paper can be considered as a four-stage hybrid flowshop, where each resources (transporters, operating rooms and recovery beds) are identical parallel machines. The four stages are : (i) the transport from the ward to the operating theatre, (ii) the surgery and recovery, (iii) the cleaning task and (iv) the transport from the operating theatre to the ward. We also have a blocking constraint between the surgery operation and the recovery operation for each patient. A finite time horizon is given for scheduling all patients. The goal is to minimize a regular criterion of patients' completion times in order to minimize the total cost of the operating theatre activity.

The hybrid flowshop scheduling problem consists in optimizing a criterion (usually minimizing the makespan) in a serial multiple-stage manufacturing system, where each stage consists of parallel identical machines. During the last decade, both hybrid flowshop scheduling problems and its numerous variations alike have been intensely investigated in the scheduling literature. The reader is referred to [START_REF] Linn | Hybrid flow shop scheduling : A survey[END_REF] and [START_REF] Kis | A review of exact solution methods for the non-preemptive multiprocessor flowshop problem[END_REF] for comprehensive surveys.

Mathematical model

In this paper, N surgery operations must be scheduled in the operation theatre during H time units, for example one day. Let p ij be the processing time of task j for patient i. The surgery of each patient i involves five different tasks : transport from the ward to the operating theatre (p i1 ), surgery operation (p i2 ), recovery in a recovery bed (p i3 ), transport from the operating theatre to the ward (p i4 ) and operating room cleaning (p i5 ).

The following assumptions are made. (H 1 ) Surgeries from the emergency department are not taken into account : in reality, emergency surgeries are often treated in a dedicated operating room which is not used for elective interventions. In addition, there is a surgeon in charge of emergency cases. This is why we do not consider emergency surgeries in our scheduling model. (H 2 ) Surgery team resource constraints are not taken into account : we consider all needed resources (operating room, surgeon, nurses, etc.) for surgery as only one resource called "operating room". (H 3 ) The recovery operation starts immediately after the surgery operation, and can occur in a recovery bed, in an operating room, or in both places.

Considering surgery team resources individually could be possible : we would have to take into account a certain number of additional resources, with various availability time constraints. Complexity of the problem would be much greater. Moreover, considering the operation room as several individual resources is not really relevant in this case because surgeon, nurses, anaesthesiologist, etc. are considered as a whole team on the field. Some studies focused on taking into account surgeon's availability, which is still a problem in operating theatre planning. Considering this problem is part of our perspectives.

The following resources are involved in the problem.

-P t : Number of identical porter teams at time t.

-R t : Number of identical operating rooms at time t.

-B t : Number of identical recovery beds at time t. Decision variables are s ij , the starting time for task j of patient i and α i , the recovery time patient i spent in an operating room. We can now state the mathematical model of the problem : Minimize :

J = i f (c i ) (1)
Subject to the following constraints :

s i1 + p i1 ≤ s i2 ∀i ∈ {1 . . . n} (2a)
s i2 + p i2 + α i = s i3 ∀i ∈ {1 . . . n} (2b)
s i2 + p i2 + α i = s i5 ∀i ∈ {1 . . . n} (2c)
s i3 + p i3 -α i ≤ s i4 ∀i ∈ {1 . . . n} (2d)
c ij = s ij + p ij ∀i ∈ {1 . . . n} (2e)
α i ≤ p i3 ∀i ∈ {1 . . . n} (2f) i (δ i1t + δ i4t ) ≤ P t ∀t ∈ {1 . . . H} (2g) i (δ i2t + δ i5t ) ≤ R t ∀t ∈ {1 . . . H} (2h) i δ i3t ≤ B t ∀t ∈ {1 . . . H} (2i)
where

δ ijt = 1 s ij ≤ t ≤ s ij + p ij -1 ∀j ∈ {1, 4, 5} (3a)
δ i2t = 1 s i2 ≤ t ≤ s i2 + p i2 + α i -1 (3b) δ i3t = 1 s i3 ≤ t ≤ s i3 + p i3 -α i -1 (3c)
The criterion J to optimize is the sum of non-decreasing functions of c i , where c i is the completion time for patient i, i.e. c i = c i4 . Therefore, we can have f

(c i ) = c i , f (c i ) = c 2 i or f (c i ) = c 3 i , for example. A criterion like f (c i ) = c 3 i or f (c i ) = c 4
i allows us achieve similar results as if we use the makespan as the assessment measure.

Constraints (2a-2d) are the precedence constraints. Constraint (2e) defines the completion time c ij of the task j for the patient i. Constraints (2g-2i) are the capacity constraints for the porter teams, the operating rooms and the recovery beds, respectively. Constraints (3a-3c) indicate whether task j of patient i is active in discrete time period t.

We can deduce from (2b) and (2c) that s i3 = s i5 . Therefore variable s i5 and constraint (2c) can be removed from the model.

Lagrangian relaxation heuristic 4.1 Lagrangian relaxation of the problem

To apply the method, the capacity constraints (2g-2i) are relaxed using the Lagrangian multipliers, which are interpreted as utilization costs of resources. In this case the multipliers λ t , β t and γ t act as prices for using transporters, operating rooms and recovery beds at time t, respectively. This way we obtain the following relaxed problem :

RP = L(λ t , β t , γ t ) (4)
where :

L(λ t , β t , γ t ) = min i f (c i ) + t λ t i (δ i1t + δ i4t ) -P t + t β t i (δ i2t + δ i5t ) -R t + t γ t i δ i3t -B t (5)
subject to :

(2a -3c)
Problem ( 5) can be further expressed as follows :

L(λ t , β t , γ t ) = - t λ t P t - t β t R t - t γ t B t + i SP i (λ t , β t , γ t ) (6)
where :

SP i (λ t , β t , γ t ) = min f (c i ) + t λ t δ i1t + t λ t δ i4t + t β t (δ i2t + δ i5t ) + t γ t δ i3t (7)
which can be simplified as :

SP i (λ t , β t , γ t ) = min f (c i ) + si1+pi1-1 t=si1 λ t + si4+pi4-1 t=si4 λ t + si2+pi2+αi+pi5-1 t=si2 β t + si2+pi2+pi3-1 t=si2+pi2+αi γ t (8)
subject to : (2a -3c)

Hence we have to solve N sub-problems SP i in order to solve the relaxed problem RP . This can be achieved using dynamic programming.

Solving the sub-problems

The N sub-problems (8) can be solved separately using dynamic programming. We will therefore obtain the optimal starting times s RP ij for each patient i. Let SP P 4i , SP P 2i and SP P 1i be the following sub-problems :

SP P 4i (T ) = min f (c i ) + si4+pi4-1 t=si4 λ t (9)
subject to :

(2a -3c) (10a)
s i4 + p i4 -1 ≤ H (10b) s i4 = T ( 10c
)
where T is the decision variable, the starting time of the transport from the operating theatre to the ward, hence s i4 . Recall that H is the fixed time horizon. Constraint (10b) ensures that the transport operation does not exceed the time horizon H.

SP P 2i (T, A) = min f (c i ) + si4+pi4-1 t=si4 λ t + si2+pi2+A+pi5-1 t=si2 β t + si2+pi2+pi3-1 t=si2+pi2+A γ t (11)
subject to :

(2a -3c) (12a)

s i2 + p i2 + α i + p i3 + p i4 -1 ≤ H (12b) s i2 = T ( 12c
)
α i = A (12d)
where T and A are the decision variables, respectively the starting time of the surgery operation (s i2 ) and the recovery time spent in the operating room (α i ). Constraint (12b) ensures that s i2 does not cause a time horizon violation.

SP P 1i (T ) = min f (c i ) + si4+pi4-1 t=si4 λ t + si2+pi2+αi+pi5-1 t=si2 β t + si2+pi2+pi3-1 t=si2+pi2+αi γ t + si1+pi1-1 t=si1 λ t (13)
subject to :

(2a -3c) (14a)

s i1 + p i1 + p i2 + α i + p i3 + p i4 -1 ≤ H ( 14b
)
s i1 = T (14c) ( 14d
)
where T is the decision variable, the starting time of the transport from the ward to the operating theatre, hence s i1 . Constraint (14b) ensures that s i1 does not cause a time horizon violation. These sub-problems can be rewritten as follow :

SP P 2i (T, A) = T +pi2+A+pi5-1 t=T β t + T +pi2+pi3-1 t=T +pi2+A γ t + min si2+pi2+pi3≤t≤H SP P 4i (t) s.t. (2a -3c), s i2 + p i2 + α i + p i3 + p i4 -1 ≤ H, s i2 = T, α i = A (15)
and :

SP P 1i (T ) = min T +pi1-1 t=T λ t + min si1+pi1≤t≤H SP P 2i (t) s.t. (2a -3c), s i1 + p i1 + p i2 + α i + p i3 + p i4 -1 ≤ H, s i1 = T (16)
where solution of SP P 4i (T ) is obvious. The sub-problem can thus be solved using dynamic programming.

Once the sub-problem is solved, the solution for the relaxed problem can be easily found. The solution of the relaxed-dual problem gives a lower bound of the initial problem J, namely L.

Building a feasible surgery schedule

The solution of the relaxed problem is in general not feasible. The schedule built with the starting times given by the solution of the dual problem might violate the capacity constraints (2g-2i) in some time periods. This subsection presents two algorithms to transform the infeasible solution into a feasible one. These algorithms take as input the infeasible solution given by the resolution of the dual problem and return a feasible solution.

Algorithm 1

Patients are sorted according to the starting time of their surgery operation and scheduled myopically one after another : the schedule is built without tanking into account the future, i.e. the remaining patients to be scheduled. For each patient, all tasks are scheduled one after another at the earliest date : the first transport is scheduled, then the surgery. At the end of the surgery process, the recovery begins in a recovery bed if available, otherwise in the operating room. If a recovery bed becomes available during this task, the patient is transferred and the operating room is released. The last two tasks (operating room cleaning and transport to the ward) are then scheduled. Patients are assigned as soon as possible when resources are available.

The main advantage of this algorithm is its simplicity and therefore its execution speed. Its complexity is O(2N ). It allows the generation of a surgery planning which is not too far from the original one. Critical tasks are the surgery and recovery operation (with a part of recovery eventually in the operating room) : keeping the order of these tasks allow us to generate a planning which have the same structure as the original infeasible solution. Transport tasks are short and can be scheduled fairly easily once the central « surgery-recovery » block is scheduled.

Algorithm 2

This algorithm takes as input the whole schedule given by the infeasible solution. The process is the following : each patient is picked from the planning and rescheduled as soon as possible, without violating the capacity constraints. The objective function value is evaluated after rescheduling each patient and its value is recorded. The patient is then replaced where he/she was, and we test the next patient. At the end of the process (i.e. when all patient has been rescheduled one time), the patient who gives the best objective function value is definitely rescheduled and the schedule is saved. The process restarts for the remaining patients who violate the capacity constraints.

The insertion of the surgery and recovery tasks may be an issue. First, the surgery task is inserted as soon as possible in the planning. Several cases appear when inserting the recovery task :

-A recovery bed is available right after the surgery until the end of the task : no problem here, the recovery task is inserted in the schedule. -A recovery bed is not available right after the surgery, but becomes available some time after until the end of the task : we have to test here if the operating room is still available after the end of the surgery until the recovery bed becomes available. Otherwise, we have to remove the surgery and find another place. -All other cases are infeasible : we have to remove the surgery and find another place. If the recovery task could be inserted respecting the capacity constraints, the cleaning task must be also inserted : if there is no sufficient time when the operating room is released, we have to remove the surgery and recovery tasks and find another place. Finally, algorithm 2 can be summarized as follow :

1. do 2.

for all patient i do 3.

remove patient i from the planning 4.

insert patient i as soon as possible, respecting capacity constraints 5. evaluate objective function value J i 6.

replace patient i where he was 7.

end for 8.

reschedule patient i who maximize J i definitely 9. until planning is feasible and solution cannot be improved anymore This algorithm is more complex than the first one : its complexity is O(N 2 ), and we will see in the results that computation times are higher. Nevertheless, this algorithm takes into account the entire infeasible schedule here. The way we build the feasible solution assure that it will be very close of the original one.

Solving the dual problem

The goal of the Lagrangian relaxation method is to obtain the best lower bound for the optimal solution. The subgradiant method is iteratively applied to solve the dual problem. For iteration n + 1, the multipliers are then updated as follows :

u (n+1) = u (n) + S (n) ∇L(u (n) ) ( 17
)
where u = (λ; β; γ). The following converging step size S (n) is used :

S (n) = ǫ L * -L(u (n) ) ||∇L(u (n) )|| ( 18
)
where L * is the optimal dual, L n is the value of the transformed relaxed problem at the n th iteration and ∇L(u (n) ) is the subgradient of the relaxed problem L with respect to the multiplier. ǫ chosen between 0 and 2 and initialized at 0.5, is a variable factor used to adjust the convergence.

Finally the Lagrangian relaxation heuristic can be expressed as follow :

1. Initialize λ t , β t and γ t with λ

(0) t , β (0) t and γ (0)
t . n ← 0. 2. Solve RP = L(λ t , β t , γ t ). We obtain the s N F ij and L.

3. Deduce a feasible solution. We obtain the s F ij and J. 4. Evaluate the gap :

gap = J -L * L * 5. Update the multipliers. n ← n + 1. 6. Go to step (2) until gap < GAP or n > M AX IT .
Values of GAP and M AX IT are fixed. The duality gap gives a measure of optimality of the feasible schedule calculated in step [START_REF] Fei | An operating theatre planning and scheduling problem in the case of a "block scheduling" strategy[END_REF]. Notice that the optimal value L * is needed to evaluate the duality gap. As we do not know its value, we replace L * by the best L found until now.

Note that, during the Lagrangian relaxation-based heuristic, we also record the schedule that minimizes the makespan, i.e. the schedule that minimizes max(c i ) among all feasible schedules encountered in the Lagrangian relaxation-based heuristic.

Numerical results

Two types of numerical results are presented in this section : (i) problem instances were generated and tested to evaluate the efficiency of the method ; then (ii) a comparison study between two scheduling strategies (allowing recovery in the operating room or not) has been done. We were not able to get data set from a real operating theatre because processing times for surgery operations or recovery are rarely recorded. However, we generate data set as near as possible to the reality, according to our observations.

Efficiency of the Lagrangian relaxation heuristic

In order to test the efficiency the method, eight classes of problem instances were generated and tested (fifteen instances each). Each problem instance class has a different size and is described in table 1. The first column of the table indicates the problem instance number, the second one the number of surgeries cases to perform during the time horizon, the next three columns respectively show the number of porter teams, operating rooms and recovery beds available. The last four columns present the ranges of the processing times of the activities involved in the operation process, which were described at the beginning of Section 3. Data are generated using a uniform distribution. The data sets are very specific for our problem in order to make the recovery in operating room occurs. In reality the number of recovery beds is always greater than the number of operating rooms.

Moreover we observed that the number of porter teams is often very low compared to operating rooms number in small operating theatres of local hospitals. In particular, porter teams are sometimes replaced by nurses in order to save some time. In the state of the art we saw that porter teams are often considered as a bottle-neck resource. This is why we found relevant to test our method with a wide range for the number of porter teams : indeed when we do not have enough porter teams, allowing recovery in the operating rooms could be a cleaver solution to avoid operating theatre overtime or to minimize waiting time of scheduled patients.

Tab. 1 -Problem instances description

In order to test our method the time was discretized into time units of ten minutes. The time horizon H is fixed between 100 and 200 time units depending of the size of the problems. The maximum number of iterations (M AX IT ) is 3000, the minimum gap value (GAP ) is 0.01% and the step size for updating multipliers is 1.9. The following four different criteria of the c i function are tested and compared :

-f 1 (c i ) = c i ; -f 2 (c i ) = c 2 i ; -f 3 (c i ) = c 3 i ; -f 4 (c i ) = (c i -LB) + .
LB is a lower bound for the makespan, which is calculated by making the sum of surgery times and cleaning times over all patients, and the result is divided by the total number of available operating rooms. As a result, we have :

LB = i (p i2 + p i5 ) max t (R t ) (19)
When used with the Lagrangian relaxation, f 4 is the criterion which gives best results for makespan : LB is a rough lower bound which only takes into account usage of the operating rooms where processing times are high. As described in the next sub-section, feasible solutions are farther from the optimal solution with f 3 or f 4 criteria, however makespans are quite good compared to solutions found with f 1 and f 2 criteria.

Results presented in tables 2 and 3 are average values over the 15 cases. All processing times for each data input file are generated randomly. Tables 2 and3 present the duality gap expressed in percentage for the different assessment criteria tested as well as the computation time between parenthesis in seconds to achieve the solution for algorithm 1 and 2 respectively. The heuristic was coded using C++ and tested with a Pentium 4 at 3.2 GHz.

A first overview of gap evaluations for the two algorithms shows that results are quite good in general. For algorithm 1, the smallest and greatest gaps are 0.04% (class 8) and 16.5% (class 4), and for algorithm 2 the smallest is 0.01% (class 8) and the greatest 31.25%. For each algorithm, the best gaps are reached for criterion f 1 : even for a large size problem as class 5 with 30 patients, the gap is only 1.62%, which suggests that larger problems could be successfully solved with this method. Duality gaps are larger for other criteria. Indeed, the two last criteria (f 3 and f 4 ) show poor gap performances for large problems (classes 4 and 5). It is worth noticing that feasible solutions of the last two criteria have smaller makespan (not shown here). Finally, computation times remain very high for almost all instances.

Tab. 2 -Gap and time performance (algorithm 1)

Tab. 3 -Gap and time performance (algorithm 2)

With respect to the performances of the two algorithms, it can be observed that the second one gives better results for small size problems (classes 1, 2, 7 and 8), whereas the first one achieves smaller gaps for large size problems (classes 4 and 5). Although the second algorithm should be chosen for solving small size problems, the benefit in this case is not really significant. For large size problems the first algorithm works far better with gaps less than one third of the ones obtained with the second algorithm. Computation times are also higher for algorithm 2, which can be explained by its complexity, as said in Section 4.4.2. Hence, we recommend the use of first algorithm for real size problems.

Benefit of patient recovery in operating rooms

In a previous paper (anonymous ref.) we studied a similar problem where recovery is not allowed in the operating room : we had to build a feasible schedule with a no-wait constraint between the surgery beds/operating rooms is lower than 3/2. This kind of situation occurs when the operating theatre is crowded : a large number of emergency patients have to be operated in regular operating rooms could be a reason to allow recovery in operating rooms if the recovery room is full. On the other hand, allowing recovery in the operating room could be a good strategy when recovery processing times grow high, i.e. when workload of recovery beds becomes greater than that of operating rooms. This kind of situation occurs when operated patient begin to bleed after a surgery and must be re-operated immediately, increasing the duration of recovery. Notice that transporters represent also a bottleneck resource : patients waiting for transporters in recovery room consume resources uselessly and prevent patients in operating rooms to be transferred into recovery room.

Conclusions and perspectives

In this paper, we presented a new operating theatre scheduling problem where patient recovery is allowed in operating rooms. Recovery in operating rooms is allowed when the recovery room is crowded, which occurs in specific cases (large number of emergency surgeries, recovery time longer than expected, etc.). The problem has been modelled as a 4-stage hybrid flowshop problem with blocking constraint, and successfully solved thanks to the Lagrangian relaxation method. By evaluating the duality gap, we were able to measure the quality of the results obtained. Several data sets were generated in order to test the efficiency of the method on real size problems.

We also compared two scheduling strategies : allowing the recovery in the operating room when to recovery bed is available or not. Using a design of experiments, we were able to determine in which case one strategy should be applied to the operating theatre : a load ratio limit has been found when resources capacity/processing times vary.

Surgeons were enthousiastic about the creation of an automatic surgery planning generation which consider recovery in operating room, as medical teams have to take into account this possibility. Unfortunately, we were not able to apply our method in a real opearting theatre yet.

Further experiments are needed to test the robustness of the model against risks. The randomness of surgery or recovery times is likely to lengthen during the process itself. Allowing recovery in the operating room will certainly lead to better results and further optimize resources utilisation, offering a more flexible schedule. Another future research direction is the introduction of resource costs into our model in order to compute accurate criterion regarding health care staff demands : opening an operating room during an hour is more expensive than using a recovery bed during the same time. Introducing over-time costs is also part of our perspectives. Finally, considering emergent surgeries in the operating theatre would constitute an interesting study.

Tab. This NP hard problem is solved using a Lagrangian relaxation approach. A comprehensive numerical experiment is then performed to quantify the benefit of allowing patient recovery in operating rooms.

Health care system environments are evolving. New policies of health care budget control have been introduced in order to improve the whole system. The operating theatre is the most critical ressource of the hospital, therefore surgery operations scheduling and human ressources planning is operating theatre are crucial issues both for patient's health and resources optimization. This study deals with a concrete problem of surgery scheduling in the operating theatre in order to reduce costs. The Lagragian relaxation method is used . We strongly believe the contribution of this study warrants its publication in the Computers & Industrial Engineering journal.

Finally, this paper is our original unpublished work and it has not been submitted to any other journal for reviews. Some results have been improved since the first paper for the IESM'07 international conference, especially about the benefit of patient recovery in operating room.

Fig. 1 -

 1 Fig. 1 -Surgery Process

Fig. 1 -

 1 Fig. 1 -Surgery Process

Table 1

 1

Table 2

 2 Tab. 1 -Gap and time performance (algorithm 2)

	Class		Criterion
		f 1	f 2	f 3	f 4
	1	0.05 (36.4)	0.12 (40.4)	0.75 (41.2)	1.27 (39.3)
	2	0.06 (45.1)	0.22 (54.8)	0.31 (42)	1.01 (49.2)
	3	0.77 (169.1) 3.30 (142.2) 12.44 (172.4) 5.40 (169.2)
	4	2.87 (376.2) 9.09 (393.5) 17.03 (393.1) 22.5 (396.4)
	5	7.16 (545.2) 16.96 (545.4) 24.49 (549.2) 31.25 (561)
	6	0.85 (136)	3.69 (142.2)	7.21 (129.5)	10.64 (154)
	7	0.02 (93.5)	0.09 (109.8)	0.22 (147)	0.22 (159.6)
	8	0.01 (29.5)	0.18 (37.1)	1.76 (49.9)	0.38 (57.4)

Table 3

 3

				Tab. 1 -Experiments 1
	Class N M Resource load	Improvement of solution (%)
			OR	Beds	f 1	f 2	f 3	f 4	Cmax
	1	4	2 3.88	7.5	10.28 30.45 65.34 44.61	34
	2	4	4 3.88	3.75	0.03	0.19	0.64	1.74	2.16
	3	4	6 3.88	2.5	0	0	0.1	0	0.4
	4	4	8 3.88	1.88	0	0	0	-0.1	0
	5	6	4 2.58	3.75	2.01	6.63 15.06 5.66	12.46
	6	6	6 2.58	2.5	0.25	1.08	2.36	0.58	1.8
	7	6	8 2.58	1.88	0	0	0.04	-0.06	0.43
	8	6 10 2.58	1.5	0	0	0	0	0
	9	8	4 1.94	3.75	9.08 24.25 48.59 19.29	16.4
	10	8	6 1.94	2.5	0.91	2.97	6.67	1.82	2.6
	11	8	8 1.94	1.88	0	0.19	0.87	0	0.1

Table 4

 4

				Tab. 1 -Experiments 2
	Class	p i,2	p i,3	Resource load		Improvement of solution (%)
				OR	Beds	f 1	f 2	f 3	f 4	Cmax
	12	{2, 4}	{2, 4}	1.38	0.5	0.04 0.04	-0.32 -0.07	0.33
	13	{2, 4}	{4, 8}	1.38	1	0	0.06	0.02	-0.14	0.53
	14	{2, 4}	{6, 12} 1.38	1.5	0.53 2.11	4.84	2.61	2.16
	15	{2, 4}	{10, 14} 1.38	2	4.47 11.34 19.82 11.14	6.85
	16	{2, 4}	{12, 18} 1.38	2.5	9.18 23.09 42.59 21.13 17.94
	17	{4, 8}	{4, 8}	2.13	1	0	0	0.02	-0.37	0
	18	{4, 8}	{6, 12} 2.13	1.5	0.29 0.65	0.62	-0.86	0.42
	19	{4, 8}	{10, 14} 2.13	2	0.05 0.09	0.43	0.24	1.17
	20	{4, 8}	{12, 18} 2.13	2.5	1.88 4.91	9.26	5.61	4.73
	21	{10, 14} {10, 14} 3.63	2	0	0.01	0	0	-0.12
	22	{10, 14} {12, 18} 3.63	2.5	0	0.05	0.05	0.08	0
	23	{10, 14} {16, 20} 3.63	3	0.07 0.21	0.44	0.7	0.11
	24	{10, 14} {20, 22} 3.63	3.5	0.74 1.48	2.12	3.1	0

Table 5

 5 December 1st, 2007 Dear Editor of the Computers & Industrial Engineering journal, Enclosed is a paper, entitled "Operating theatre scheduling with patient recovery in both operating rooms and recovery beds". Please accept it as a candidate for publication in the Computers & Industrial Engineering journal. This paper addresses the surgery operation scheduling problem under the open scheduling policy.

operation and the recovery/cleaning operations. In order to decide which method is better depending of the case study, we designed several experiments presented tables 4 and 5. The number of patients is fixed and equals to ten. Comparisons are made using algorithm 1. We define the improvement of the solution in percentage as follows :

The first experiment (table 4) tests the effect of resources capacity variation. Processing times are fixed to the same values as class 1 in the last section. Classes 1 to 4 display results for an operating theatre with 4 operating rooms and 2 to 8 recovery beds. Classes 5 to 8 display results for 6 operating rooms and 4 to 10 recovery beds. At last, classes 9 to 11 display results for 8 operating rooms and 4 to 8 recovery beds. The numbers of transporters is fixed and will not effect the performances of the two strategies.

The second experiment (table 5) describes the influence of resources workload. There are 4 operating rooms and 6 recovery beds, which is the common ratio 3/2 for operating theatres. Two porter teams are also dedicated to the system. Instances of classes 12 to 16 display results for the case with average operating room workload of 1.38 and average recovery bed workload from 0.5 to 2.5. Classes 17 to 20 display results for operating rooms workload of 2.13 and recovery beds workloads from 1 to 2.5. Finally, classes 21 to 24 display results for operating rooms workload of 3.63 and recovery beds workloads from 2 to 3.5.

Tab. 4 -Experiments 1

Tab. 5 -Experiments 2

As expected, the benefit of patient recovery in operating rooms depends on the workload ratio of operating rooms/recovery beds. We observe in tables 4 and 5 that all criteria are improved when the recovery beds load becomes greater than the operating rooms load. In experiment 1 (see table 4), classes 1, 5 and 9 where the recovery beds are heavily loaded with respect to the operating rooms, the makespan is improved by at least 12%. In classes 2, 6 and 11 where the recovery beds and the operating rooms are equally loaded we manage to get makespan improvements from 0.1% to 2.16% for smaller operating theatres. In experiment 2 (see table 5), the makespan is improved by 17.94% for class 16 and 4.73% for class 20. We did not compare the methods with greater loads for recovery beds in order to keep consistency with reality.

If we compare results among the four criteria f 1 , f 2 , f 3 and f 4 , best improvements are reached for criterion f 3 with 65.34% in class 1, 48.59% in class 9 and 42.59% in class 16. These results are obtained for heavily loaded recovery beds, even for a capacity-compliant operating theatre. To the contrary, criterion f 4 is only improved by 20% in classes 9 and 16. Allowing recovery in the operating rooms can even lead to slightly worse performance as we can see in class 12, 17 and 21 where recovery beds load is equal to operating rooms load. This is mainly due to the heuristic nature of the Lagrangian relaxation heuristics. Note however that any feasible schedule not allowing patient recovery in operating rooms is a feasible solution in our case.

The experiments presented in this section allow us to know what is the best strategy depending of the operating theatre configuration and/or the type of surgeries. Results presented in table 4 show that allowing the recovery in the operating room improve performances as soon as the ratio recovery


Y Liu

M Durand

S Chen

A Houard

B Prade

B Forestier

A Mysyrowicz
email: andre.mysyrowicz@ensta.fr.

Energy Exchange between Femtosecond Laser Filaments in Air

Keywords: numbers: 52, 38, Hb, 42, 65, Jx, 42, 65, Sf

We report on the energy exchange between femtosecond laser filaments in air. A traveling plasma grating formed at the intersection of the filaments is proposed to explain the energy transfer. In this moving plasma grating mediated mechanism the laser energy transfers from the lower frequency pulse to its higher frequency partner, while in a traditional Kerr nonlinearity mediated grating, energy transfer occurs in the opposite manner. Based on this mechanism, we demonstrate an energy transfer ratio as high as 50 percent for pulses with energy of several millijoules.

Propagation of powerful femtosecond (10 -15 s) laser pulses in air gives rise to a series of spectacular effects. As it propagates, the light pulse self focuses due to the change of refractive index it induces in the medium. For a pulse power exceeding a critical value P cr (~ 5 GW in air for = 800 nm), diffraction is unable to arrest the beam self-focusing. The intensity in the beam center increases progressively with distance until it starts ionizing air molecules. Once created, the plasma induces a negative change of refractive index, an effect leading to beam defocusing. The ensuing dynamic competition between self-focusing and plasma defocusing is responsible for filamentation, a process where part of the beam maintains a very high intensity over long distances leaving in its trail a column of weakly ionized plasma [1,2].

A large part of the fascination with the process of filamentation is due to its potential for applications. The use of filaments has been proposed and/or demonstrated for various applications ranging from remote sensing of atmospheric pollution [3], recognition of distant targets by Laser Induced Breakdown Spectroscopy [4,5], stand off THz illumination [6,7], to triggering and guiding of electric discharges and even lightning [8]. For all these applications, control and optimization of the filament characteristics are desired. This can be done to a certain extent by modifying the laser beam characteristics before propagation [9,10]. In the last two years, more subtle ways to control filaments have been discovered. One concept relies on molecular lensing: the spontaneous revivals of the refractive index in air induced by a precedent femtosecond pulse can deviate or destroy subsequent filaments [11,12]. In another approach, the use of an intense Airy beam leads to a curved trajectory of the filament [13].

Here, we report on another approach to control filament properties, namely the possibility of exchanging energy between two femtosecond laser filaments in air. This effect should allow their regeneration and increase their length. Energy exchange between powerful femtosecond laser beams in air has been pioneered by Bernstein et al. who showed that two beams crossing at grazing angle could exchange up to 7% of their energy [14]. The exchange reported here, which has a different origin, has an efficiency up to 50% and is also observed at large crossing angle.

The physics of the process is captured in figure 1. Two laser pulses with power P > P cr of slightly different frequencies, form two intersecting plasma filaments in air (see Fig. 1 a). Fig. 1 b is actually a photo of the characteristic luminescence emanating from two such crossing plasma channels [15]. The fields of the optical pulses with wave numbers 1 K  and 2 K  interfere and give rise to a spatially periodic pattern

    2 sin 2    d , see Fig. 1 c. With a frequency difference of 1 2       , this interference pattern moves in the direction of 2 1 K K K g      with phase velocity p g v K   



. The interference field leads in turn to the formation of traveling gratings via the nonlinear response of the medium.

In general, the nonlinear response of a transparent medium to an intense laser field consists of an instantaneous and/or a retarded change of the refractive index. It is only the retarded response that can result in energy exchange because it induces a spatial phase shift between the laser intensity interference and the refractive index modulation [START_REF] Boyd | NONLINEAR OPTICS[END_REF]. In the case where the nonlinear part of the refractive index follows the Debye relaxation, the nonlinear contribution of the refractive index reads

    ' ' , ' exp , 2 dt t r I t t n t r n t                
. From the wave propagation equation, the intensity variation of beam 1 (or 2) can be expressed as:

  2 2 1 2        c I I n dz dI j i i  (i, j = 1, 2 and i ≠ j)
, where  and n 2 are the nonlinear response time and the nonlinear refractive index respectively [START_REF] Boyd | NONLINEAR OPTICS[END_REF]. In this case, two features can be found from above: (a) a positive refractive index response (n 2 ) will induce laser energy flow from higher frequency to lower frequency pulses, while a negative n 2 will induce energy transfer in the opposite manner. (b) To achieve energy exchange, the interacting pulses must have different frequencies. The frequency difference dictates the direction of energy flow. Now, we consider more specifically two retarded nonlinear responses in air, the change of refractive index due to the rotational Raman effect, and that induced by the formation of a plasma. For the Raman response of the molecule, the nonlinearity is related to the laser intensity as

      2 , ' , ' ' t Raman n r t n R t t I r t dt        .
The function R(t) mimics the molecule response with a characteristic time  and frequency  R :

      0 exp sin , R R t R t t    where     2 2 0 R R R t     
. In air for a pulse at 800 nm, 1   = 70 fs and  R = 1.6×10 13 Hz [1,11]. In Fig. 2 a and b, we plot the calculated temporal evolution of the laser intensity and the nonlinear index modulation as a function of position x at propagation distance z = 0 (see Fig. 1 c). In the simulation, the following realistic parameters were chosen: laser intensity in the filament 5×10 13 W/cm 2 , frequency difference Ä = 10 THz (see below) and crossing angle 2º. The interference pattern and index grating both move in the direction of g K  . In Fig. 2 d, the intensity (blue line) and the refractive index modulations (black line) at time = -20 fs are plotted. We note a spatial phase shift between the two modulations (of about 1.2 m), which is essential for possible energy exchange.

For the plasma generation process, the retarded refractive index change is related to the plasma density by

  r c plasma t r n   2 ,    
, where

        t mol dt t r I t r ' ' , , 8 8     
is the plasma density [1,2]. The terms 8

 and  mol are the multiple photon ionization cross section and the density of the neutral oxygen molecules. The calculated evolution of the moving refractive index grating at z = 0 is presented in Fig. 2 c. The grating corresponds here to a free electron density modulation [START_REF] Sheng | [END_REF]18]. At the end of the pulse overlap, the plasma grating stays at rest and remains unchanged because its decay time is much larger than the time scale considered here. Fig. 2 d shows the spatial distribution of the plasma grating at time = -20 fs (red line). A spatial phase shift of nearly 180° with respect to the driving interference pattern is apparent, due to the negative refractive index change. As a result, the energy is expected to flow now from the low frequency pulse to its high frequency partner.

In our experiment, we use a laser delivering pulses of  0 = 45 fs duration centered at 800 nm with energy of up to 15 mJ. To obtain two laser pulses with controllable difference in their instantaneous frequency, a linear chirp is impressed to the laser pulse before splitting it into two replicas. If pulse 1 is delayed with respect to pulse 2 by  d , the difference of the instantaneous frequencies is

2 2 p d C       [19].
Here,  p is the chirped pulse duration with the chirp parameter

  1 2 0      p C
. At a given time delay  d the two pulses possess a constant frequency difference during their time overlap. Figure 3 shows results obtained in air and in Argon at 1 bar pressure, using a crossing angle of  9 . 0   . A linear positive chirp (i.e. an increase of instantaneous frequency from the leading to trailing edge of the pulse) C = +5 was impressed to the initial pulse yielding a pulse duration of 220 p fs   . To reduce the distance required for beam collapse, the two chirped pulse replica, both vertically polarized, were focused by two lens of f = 2000 mm. Above P cr , filaments were formed with a length of the plasma column up to 150 mm at P = 15 GW. Turning the attention to air first, we observe two different regimes. At lower power P = 0.39 P cr (Fig. 3 a), transfer of energy is always directed from the high frequency beam to the low frequency beam. This is similar to the observation by Bernstein et al. and can be assigned to a grating induced by the rotational molecular response of the optical Kerr effect [14]. It would be interesting to explore this regime of energy exchange with very short pulses. One then expects an additional contribution from the vibrational Raman response of the molecules [20]. At higher incident laser power P = 3.5 P cr (Fig. 3 b), the energy transfer changes direction. This new regime is a genuine energy exchange between filaments in the sense that it occurs once plasma channels are formed along the paths of beams 1 and 2. In Fig. 3 c, the maximum energy exchange ratio defined as S = (E 1 -E 2 )/(E 1 +E 2 ) at positive  d is plotted as a function of incident power, where E 1 and E 2 are the energy of the first and second pulses after exchange. The turning point between the two regimes occurs at a power close to P cr . Also, we point out that with negatively chirped pulses in both regimes the energy transfer is reversed between pulse 1 and 2 (for example, in fig. 3(a) for  d < 0 the laser energy now transfers from pulse 2 to pulse 1).

With argon (see Fig. 3 d,e, and f), we notice the absence of energy exchange until the laser power reaches P cr . This is easily understood since there is no rotational Raman effect in an atomic gas. As expected, we find that energy exchange always occurs with a negative ratio S (energy transfer from the lower frequency to the higher frequency pulse). With the increase of the laser power, S increases almost linearly up to 12%.

Because the plasma grating is composed of a weakly ionized gas, a higher refractive index contrast can be expected with a higher intensity in the filament, allowing for more efficient energy exchange. It is possible to achieve this in short filaments by using more converging initial beams [21]. In Fig. 4 a, the exchange ratio S is presented for focal lens of f = 125 mm. It is seen that the ratio S can reach 50% for pulses with initial power P ~ 3 P cr . In that case, the plasma column length is about 5 mm. As an example, the variation of beam energy is shown in Fig. 4 b as a function of pulse delay (frequency difference ) for P = 1.5 P cr .

Finally, the exchange ratio was also studied as a function of the effective interaction length, which was varied by changing the crossing angle while keeping the same focal length f = 300 mm. The results are presented in Fig. 4 c, with the numbers on the top of each point indicating the corresponding crossing angle. As expected, the efficiency increases for longer interaction lengths.

In conclusion, we have shown that it is possible to exchange efficiently laser energy between filaments in air. The source of this novel energy transfer is a traveling plasma grating, which is applicable not only to Raman active molecular gases but also to monatomic noble gases. The magnitude and direction of energy flow from one beam to the other can be controlled by detuning the frequency of the two filamentation pulses. We believe this mechanism can be useful to remotely replenish filaments. refractive index modulation at z = 0 induced by the rotational Raman response. The refractive index modulation also moves downward. A temporal retardation with respect to the laser interference pattern is apparent at the beginning and end of the pulse overlap. (c), Temporal evolution of the plasma induced refractive index modulation at z = 0. This plasma grating travels downward during its production process (-80 fs to 40 fs). After 40 fs, the plasma density is constant because no decay mechanism is considered in the model. (d), Laser intensity (blue line, left ordinate), refractive index modulation due to Raman response (black line, right ordinate) and plasma (red line, right ordinate) at time = -20 fs. The spatial phase shift between the laser intensity and the two index modulations is 1.2 m and 0.5 m, respectively.



  is the critical plasma density for 800 nm pulse and

Fig. 1 .

 1 Fig. 1. Schematic experimental configuration with two interacting filaments. (a) Two focused pulse replica propagate from left to right and form two interacting filaments crossing under angle . The laser polarizations are orthogonal to the plane formed by the two beams. The pulse energies after interaction are measured by detectors d1 and d2. (b) Photo of two interacting filaments, displaying the characteristic near UV luminescence of N 2 in air following ionization. (c) Schematic representation of the magnified area where the two plasma columns overlap, showing the interference pattern of the two laser fields traveling with a velocity p  . The

Fig. 2 .

 2 Fig. 2. Temporal and spatial distribution of laser intensity and refractive index modulation in the interaction region. (a), Temporal evolution of the laser interference intensity at z = 0. Chirped pulse 2 of duration 200 fs is delayed by 100 fs with respect to chirped pulse 1. This corresponds to  = 10 THz, see text. The origin of the time corresponds to the moment when the center of the overlapped pulses reaches z = 0. The laser intensity

Fig. 3 .Fig. 4 .

 34 Fig. 3. Energy exchange of filaments in air (a, b, c) and argon (d, e, f). (a), Energy exchange of filaments in air as a function of delay (frequency difference ) between the chirped pulses with power P = 0.39 P cr and duration  p = 220 fs. Each point corresponds to the average from 100 measurements. (b), Same as in (a) except for P = 3.5 P cr . (c), Maximum Energy exchange ratio S recorded at positive delay  d as a function of incident laser power in air. (d), (e), and (f), similar as (a), (b), (c), except air is replaced by Argon. In all cases, the focal lengths for both pulses are 2000 mm and the crossing angle is 0.9°.

We acknowledge fruitful discussions with Professor V. Tikhonchuk and Dr C. Arnold and B. Zhou. This work has been partially supported by DGA. S. C. appreciates support from the National Science Foundation of China through Grant 10874024.

REFERENCES:


R F Lang

G Angloher

M Bauer

I Bavykina

A Bento

A Brown

C Bucci

C Ciemniak

C Coppi

G Deuter

F Von Feilitzsch

D Hauff

S Henry

P Huff

J Imber

S Ingleby

C Isaila

J Jochum

M Kiefer

M Kimmerle

H Kraus

J.-C Lanfranchi

B Majorovits

M Malek

R Mcgowan

V B Mikhailik

E Pantic

F Petricca

S Pfister

W Potzel

F Pröbst

S Roth

K Rottler

C Sailer

K Schäffner

J Schmaler

S Scholl

W Seidel

L Stodolsky

A J B Tolhurst

I Usherov

W Westphal

Electron and Gamma Background in CRESST Detectors

Keywords: CRESST, CaWO 4, Low Background, Calcium-41, Dark Matter PACS: 29.40.Vj, 95.35.+d

The CRESST experiment monitors 300 g CaWO 4 crystals as targets for particle interactions in an ultra low background environment. In this paper, we analyze the background spectra that are recorded by three detectors over many weeks of data taking. Understanding these spectra is mandatory if one wants to further reduce the background level, and allows us to cross-check the calibration of the detectors. We identify a variety of sources, such as intrinsic contaminations due to primordial radioisotopes and cosmogenic activation of the target material. In particular, we detect a 3.6 keV X-ray line from the decay of 41 Ca with an activity of (26 ± 4) µBq, corresponding to a ratio 41 Ca/ 40 Ca = (2.2 ± 0.3) × 10 -16 .

Introduction

The Cryogenic Rare Event Search with Superconducting Thermometers CRESST [START_REF] Angloher | The CRESST Collaboration)[END_REF]2] aims to detect rare nuclear recoils from elastic scattering of Dark Matter particles [3]. Since less than one such recoil is expected per kilogram of target mass and week of exposure (e.g. [2]), efficient shielding of the target against ambient radioactivity is mandatory.

To this end, the CRESST experiment is located in the Laboratori Nazionali del Gran Sasso under an average of 1400 m of rock overburden, where the cosmic muon flux is reduced by about six orders of magnitude with respect to the surface. The target is surrounded (from the outside to the inside) by 45 cm of polyethylene to moderate the neutron flux, a radon box that is constantly flushed with nitrogen gas to reduce the concentration of radioactive Email address: rafael.lang@mpp.mpg.de (R. F. Lang) radon and its daughters inside the shielding, as well as 20 cm of lead and 14 cm of pure copper against the external gamma background [2]. In addition, the expected WIMP-induced nuclear recoil signal can be discriminated from the remaining electron and gamma background utilizing the scintillation light of the CaWO 4 crystals [START_REF] Angloher | The CRESST Collaboration)[END_REF]. The understanding of the remaining radioactive background present in the experiment prior to this discrimination is subject of this paper. It demonstrates the excellent performance of CRESST detectors and is a necessary step to improve the sensitivity of the experiment to rare processes.

Setup

The target consists of individual 300 g CaWO 4 crystals that are operated as calorimeters [START_REF] Angloher | The CRESST Collaboration)[END_REF]. Up to 33 such detectors can be arranged in a com-pact layout [2]. To enhance the temperature signal following an interaction, the crystals are cooled to ∼ 15 mK where heat capacities are low. The high frequency phonons that are created in interactions in the crystal are collected by a thermometer, which is a thin tungsten film that is evaporated directly onto the crystals. This film becomes superconducting at these temperatures and is stabilized in its transition to the superconducting state. A small temperature rise results in an increase of resistance that is measured with a SQUID-based readout circuit [4].

A resistive heater structure is attached to each thermometer. This structure is used to provide a constant heating bias as well as for periodic injection of additional heat pulses. Each such injected heat pulse briefly drives the thermometer out of its transition. The amplitude response of the thermometer to the injected heat is controlled to stay at a constant value with small adjustments to the heating bias by means of a feedback system. In this way the operating point of the thermometer is stabilized over several months to within the required µK precision [START_REF] Angloher | The CRESST Collaboration)[END_REF]. This can also be seen from the sharpness of gamma lines discussed here.

Pulse Height Evaluation

Our thermometers yield rather slow pulses with rise times of about 1 ms. These pulses are digitized with a time base of 40 µs. 4096 such samples comprise a record, of which the first 1024 samples are taken from the time before the trigger to define the baseline of the pulse.

A template pulse is fitted to the record in order to evaluate the amplitude of the pulse as a measure of the deposited energy. This template is an average of appropriate pulses, often obtained from a gamma calibration line. Since the SQUIDs only measure relative changes in current, the baseline needs to be defined individually for each recorded pulse. Hence, the free parameters in the fit are the level of the baseline, the onset of the pulse, and its amplitude.

The pulse samples the transition curve of the thermometer from its superconducting to the normal conducting state. For interactions with recoil energies below ∼ 200 keV, the resultant pulse generally samples only the linear region of the thermometer's transition curve between the superconducting and normal conducting state. For interactions with larger energies, the pulse samples a region of the transition curve that is non-linear, distorting the resulting pulse shape. Therefore, in the template pulse fit, only those parts of the pulse are included that show a linear behavior. This linearizes the energy response of the fitted amplitudes over a wide energy range, even up to the MeV range [5].

Energy Calibration

To reach the cryogenic temperatures, the detectors are enclosed by a fivefold thermal shield made from copper. Calibration sources are placed outside this shield, which has a total thickness of 12 mm of copper. This requires gamma energies 100 keV in order to penetrate this barrier. Hence, a method is required to take the energy calibration to lower energies.

To this end, additional heater pulses corresponding to a number of fixed energies are injected into the heater structure every 30 seconds. The pulse height of the response pulse from the thermometer is also evaluated with a template fit. Non-linearities are due to the shape of the transition curve and due to the readout circuit. They manifest themselves in the amplitude of the response pulse, but, for given amplitude, not in its shape. Therefore, the amplitudes of the response pulses give the correct energy assessment even if there are small differences in the pulse shape of heater response pulse and particle pulse.

The heater pulses are calibrated at relatively high energies, typically with the 122 keV or 136 keV line from a 57 Co source. The amplitude response of the thermometer as a function of energy injected into the heater can then be fitted with a low-order polynomial function. Figure 1 depicts this procedure for the calibration of particle pulses in the complete energy range covered by the heater pulses. In the following we identify a variety of lines at the expected energies, thus validating our calibration procedure.

Calibration Spectrum

Figure 2 shows a typical spectrum from a 57 Co calibration. The two strong peaks at 122 keV and 136 keV originate from gamma transitions in the 57 Co source. They show pronounced shoulders towards lower energies, which originate from photons that are Compton scattered in the copper thermal shields and the detector support structures that surround the detectors. Then, for each particle pulse, the particle's energy is found by evaluating the transfer function at the fitted amplitude (dashed line).

β 122-K 1 β 122-K 2 α K 1 α K 1 α 122-K 2 α 122-K 1 β K 2 β 136-K 2 β K 1 β 136-K 1 α 136-K 2 α 136-K 57Co 57Co
A gamma can interact in the crystal via the photoelectric effect, ejecting an electron from its shell. The X-ray produced in the subsequent recombination may escape from the crystal. In this case, we measure the energy of the original gamma minus the energy of the X-ray (escape peak). If the escaping X-ray reaches a neighboring crystal, its energy will be measured there. In figure 2 we observe a set of peaks in the energy range between 50 keV and 80 keV, originating from the dominant K α1 , K α2 , K β1 /K β3 and K β2 escape lines of tungsten X-rays. The origin of the various observed lines are labeled in the figure. In addition, a line is visible at 8 keV from copper fluorescence due to interactions in the copper surrounding the cryostat (see section 6.6).

The calibration spectrum allows confirmation of our energy calibration method. This is shown in table 1, where measured energies are compared to the literature values [START_REF] Firestone | Table of Isotopes (CD ROM Edition)[END_REF] for a few lines. The values agree within the resolution of our detectors. Even the extrapolation of our energy calibration to high energies is validated by the 692 keV line from 57 Co, which is reconstructed only about 1% below the expected value, namely at 682 keV with a resolution of 3 keV (1σ).

Background Spectra

In the following, spectra from exposure to background radiation alone are analyzed, from one de- The energy resolution is seen to be better than 1 keV.

tector during the prototyping phase of the experiment [START_REF] Angloher | The CRESST Collaboration)[END_REF] (figures 3, 4) and two detectors during the commissioning phase [2] (figures 5, 8 and 6). All these crystals were grown by the General Physics Institute, Moscow, Russia, but in different batches.

Strontium and Yttrium

Powders containing calcium are difficult to purify from strontium since both elements belong to the same chemical group. Thus, strontium can easily remain in our CaWO 4 crystals during the Czochralski growth process. Indeed this is the case, as can be seen in figure 3, where we show the spectrum of one crystal up to 1 MeV. We observe a continuous background which can be attributed to the beta decay of 90 Sr (with an endpoint energy of 546 keV) together with the subsequent beta decay of 90 Y (2282 keV endpoint energy). A fit above 80 keV to a parametrization of the expected beta spectrum from 90 Sr in equilibrium with 90 Y [START_REF] Kossert | [END_REF] can be seen to describe the data in this energy range. From the fit we deduce an activity of (4.37 ± 0.15) mBq of 90 Sr for this crystal. In addition to this common and continuous background, the three crystals show quite different features in the energy range below 100 keV.

Lead

Isotopes from natural decay chains are present in the crystals and can be identified by their alpha decays [5]. A step in the natural decay of 238 U is 210 Pb, which beta-decays into 210 Bi with an endpoint energy of 63.5 keV. In 84% of these decays a 46.5 keV gamma is emitted. If the decays happen in the vicinity of our crystals, we observe a line at 46.5 keV (see figure 4), which reflects an external activity of (31 ± 5) µBq. If, on the other hand, 210 Pb is a contamination intrinsic to the crystal, our calorimetric measurement gives the energy of the gamma plus that of the emitted electron. Hence, we observe a beta spectrum starting at 46.5 keV that extends up to 63.5 keV, as can prominently be seen in figure 5. The activity of this decay due to sources both internal and external to the crystal is much higher, namely (672 ± 28) µBq. Figure 6 also shows this lead feature for the third crystal with an activity of (241 ± 25) µBq.

Actinium

A similar case is that of 227 Ac, a step in the natural decay of 235 U. 227 Ac beta-decays with an endpoint energy of 44.8 keV into 227 Th, where two excited levels may lead to the emission of 24.5 keV and 9.3 keV gammas. This results in two beta spectra, one starting at 9.3 keV, and one at 24.5 keV, as can prominently be seen in figure 6 but also in figure 5. This contamination is not observable in the spectrum shown in figure 4, which indicates that this detector contains less impurities from the decay of 235 U.

Activated Tungsten

If the raw material or the crystals themselves are exposed to cosmic radiation, the tungsten in the CaWO 4 crystals can be activated. A possible channel is 182 W(p, α) 179 Ta. 179 Ta then decays via electron capture with a half-life of 1.8 years into 179 Hf. The energy signature of this decay in our crystals is the binding energy of the captured electron, mostly from the K-shell, which is E K,Hf = 65.4 keV. The activation via 183 W(p, t) results in 181 W. This decays with a half-life of 121 days via electron capture into an excited 181 Ta nucleus, which de-excites with emission of a 6.2 keV gamma. Hence, the dominant energy signature of this decay is the binding energy of a tantalum K-shell electron, E K,Ta = 67.4 keV, plus the energy of the gamma, 6.2 keV, adding up to 73.7 keV.

Lines at these energies can be seen in all spectra. For the crystal of which the spectrum is shown in figure 4, a measurement with a similar exposure to background radiation, about three months earlier, is available [START_REF] Majorovits | The Identification of Dark Matter[END_REF]. The crystal was stored underground in between these two measurements, so no re-activation could occur. This allows us to extract the 179 Ta and 181 W activities from two measurements. A comparison to the expected decay times is shown in figure 7 and confirms the origin of these lines.

Activated Calcium

Capture of thermal neutrons from cosmic rays can activate calcium present in the CaWO 4 crystals. 41 Ca has a half-life of 10 5 years and decays via electron-capture to 41 K. The K-shell binding energy of potassium is only 3.61 keV, but our detectors are capable of detecting such low energies. Figure 8 shows again the spectrum of figure 5 but for lower energies and with a finer binning. We observe a prominent line at (3.71 ± 0.02) keV (mean and error from the fit of a Gaussian). This is only 3% above the expected energy from 41 Ca, which we therefore attribute to this line. The width of this line is only (130 ± 19) eV, consistent with the calibration. For illustration, figure 9 shows a typical pulse at this low energy, still clearly visible above baseline noise.

From the observed activity of (26 ± 4) µBq we can readily calculate that there are (1.4 ± 0.2) × 10 8 41 Ca atoms present in the crystal, only a fraction of (2.2 ± 0.3) × 10 -16 of all calcium atoms. The level of activity is consistent with expectations for raw materials mined from a few meters below surface. Our measurement is an order of magnitude more accurate than other measurements of this fraction, which are typically carried out by accelerator mass In addition to the lines from 41 Ca and copper fluorescence there is a line at ∼ 11.5 keV which so far could not be identified in a consistent way (detector Verena, run 30). spectrometry (AMS) [START_REF] Merchel | Nuclear Instruments and Methods B[END_REF][START_REF] Fink | [END_REF]11]. Our present sensitivity is in the 10 -17 range even with CaWO 4 crystals that are not optimized for this purpose. This demonstrates the extreme sensitivity of our detectors for rare decays and opens new possibilities for radioactive dating using 41 Ca [12]. In fact, having an in situ calibration line at these low energies is an advantage for the direct Dark Matter search, since it independently confirms our energy calibration procedure in this otherwise inaccessible region.

Copper Fluorescence

Another line, seen at (8.17±0.04) keV in figure 8, can be attributed to copper fluorescence. This originates from gammas or electrons that eject a photoelectron from the copper of the cryostat that surrounds the detectors, resulting in a copper X-ray absorbed in our crystals. Our energy calibration puts the line only 2% higher than the expected 8.04 keV from the literature. This is consistent with the slight energy overestimate seen when identifying the 41 Ca spectral line, again confirming its origin.

Lutetium

All single crystal spectra above 100 keV are rather featureless. Yet, two crystals were operated simultaneously and in close vicinity during the commissioning of the setup. This allows us to search for coincident events, of which a spectrum is shown in figure 10. At 59 keV, tungsten K α escape events can be seen. In addition, two lines at 202 keV and 307 keV become very prominent, and we attribute those to a contamination with lutetium intrinsic to the crystals. 176 Lu beta-decays to 176 Hf with a total decay energy of 1191.7 keV. Each decay produces a gamma cascade with gamma energies of 306.8 keV, 201.8 keV and 88.3 keV. In particular the higher energy gammas may escape from one crystal and hit the other. Hence, for a coincident 306.8 keV event observed in one detector, the other crystal will mostly absorb the other two gammas, plus the energy of the beta electron. This explains the structure seen in the scatter plot of figure 11, showing the energies of coincident events in each detector. We observe a beta spectrum starting at 201.8 keV + 88.3 keV = 290.1 keV for most events at 306.8 keV in the other detector (lines marked (a) and (b) in figure 11). Figure 12 shows the energy spectrum of these events in one detector, tagged by an energy of 306.8 keV in the other one, in agreement with the expectation from the lutetium beta spectrum.

Conclusions

The excellent energy resolution of less than 1 keV over the full energy range used in the Dark Matter search [START_REF] Angloher | The CRESST Collaboration)[END_REF]2] allows us to identify a variety of different background sources in the CRESST experiment. Since the CRESST experiment can discrimi- nate nuclear recoil signal from electron recoil background based on the light output, the presence of these features does not compromise the search for Dark Matter, yet allows an audit of the performance of the experiment. The energy resolution (1σ) is 130 eV at 3.6 keV, and the observed sharpness of the lines proves the high stability of the detectors during the background measurements running for several months. The observed line from activated calcium provides an independent in situ energy calibration even below the energy range used for the Dark Matter search, and in addition proves the high sensitivity of the experiment to rare decays.

Figure 2 :Figure 1 :

 21 Figure2: Energy spectrum from a high statistics 57 Co calibration (detector Daisy, run 31). The relative strengths of escape peaks to direct gamma lines depend on the position of the crystal, the source, as well as the position of neighbouring crystals and the surrounding cryostat and is in agreement with simulations.

Figure 3 :

 3 Figure 3: Background spectrum observed with one crystal (m = 306.8 g) from an exposure to background radiation of 15.00 kg d early 2004 (detector Daisy, run 28). The black line is a fit above 80 keV of a parametrization of the 90 Sr/ 90 Y beta spectrum to the data.

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Zoom on the energy region below 100 keV of the spectrum shown in figure 3.

Figure 7 :

 7 Figure 7: Two measurements of the 181 W and 179 Ta activities in one crystal during two runs with similar exposure to background radiation. The expectations from decays with T 1/2 = 121 d and T 1/2 = 665 d are drawn as dashed lines and seen to be consistent with the measurement. Error bars are uncertainties of the fit to the spectral lines (each containing tens of counts) and the running time of the respective measurements (detector Daisy, run 27 and run 28).

Figure 8 :

 8 Figure 8: Low-energy part of the spectrum shown in figure 5.In addition to the lines from 41 Ca and copper fluorescence there is a line at ∼ 11.5 keV which so far could not be identified in a consistent way (detector Verena, run 30).

Figure 9 :

 9 Figure 9: A typical pulse with an energy of 3.7 keV (grey), and its 50-sample (2 ms) moving average (black).

Figure 10 :

 10 Figure10: Spectrum of coincident events. To increase the counting statistics, the spectra of the two individual detectors were added to produce this figure. Without the coincidence requirement, the added spectrum would have a count rate of O(500) counts per 2 keV bin (detectors Verena and Zora, run 30).

Figure 11 :

 11 Figure 11: Scatter plot of coincident events in two detectors operated during commissioning. The vertical line (a) going up from 290 keV is due to lutetium decaying in the second detector, in coincidence with the 307 keV gamma escaped to the first detector. The horizontal line (b) is due to the same process with detectors reversed, at 307 keV in the second detector and starting at 290 keV in the first one. At (c) we observe a hint of a horizontal line starting at 395 keV due to lutetium decaying in the first detector with the energy of an escaped 202 keV gamma in the second detector.

Figure 12 :

 12 Figure12: Histogram of all coincident events in one detector that appear at (307 ± 3) keV in the respective other detector: The 176 Lu beta spectrum starting at 290 keV is visible despite the low statistics.

Table 1 :

 1 Literature values, reconstructed energies and resolution of a few lines observed in the cobalt calibration. All values agree and validate our energy calibration procedure.

	line	literature reconstructed resolution
	origin	value	energy	(1σ)
	Cu K α	8.0 keV	8.3 keV	0.1 keV
	122 keV -K β1	54.9 keV	55.0 keV	0.4 keV
	W K α1	59.3 keV	59.6 keV	0.5 keV
	122 keV -K α1	62.8 keV	63.0 keV	0.3 keV
	57 Co	122.1 keV	122.2 keV	0.6 keV
	57 Co	136.5 keV	136.5 keV	0.6 keV

8. Acknowledgments R. F. L. acknowledges useful discussions with G. Korschinek on 41 Ca. This work was partially supported by funds of the DFG (SFB 375 and Transregio 27 "Neutrinos and Beyond"), the Munich Cluster of Excellence ("Origin and Structure of the Universe"), the EU networks for Cryogenic Detectors (ERB-FMRXCT980167) and for Applied Cryogenic Detectors (HPRN-CT2002-00322), and the Maier-Leibnitz-Laboratorium (Garching). Support was provided by the Science and Technology Facilities Council.