Kondo Hloindo

Kossi Pelope

Minimizing Maintenance Cost Involving Flow-time and Tardiness Penalty with Unequal Release Dates

This paper proposes important and useful results relating to the minimization of the sum of the flow-time and the tardiness of tasks or jobs with unequal release dates (occurrence date), with application to maintenance planning and scheduling. Firstly, the policy of real-time maintenance is defined for minimizing the cost of tardiness and critical states. The required local optimality rule (FTR) is proved, in order to minimize the sum or the linear combination of the tasks' flow-time and tardiness costs. This rule has served to design a scheduling algorithm, with O(n 3 ) complexity when it is applied to schedule a set of n tasks on one processor. To evaluate its performance, the results are compared to a lower bound that is provided, in a numerical case study. Using this algorithm in combination with the tasks urgency criterion, a real-time algorithm is developed to schedule the tasks on q parallel processors. This latter algorithm is finally applied to schedule and assign preventive maintenance tasks to processors in the case of a distributed system. Its efficiency enables, as shown in the numerical example, to minimize the cost of preventive maintenance tasks, expressed as the sum of the tasks tardiness and flow-time. This corresponds to costs of critical states and of tardiness of preventive maintenance.

Task scheduling, maintenance planning, flow-time, tardiness, release date (occurrence date), availability, distributed systems.

Notations

N number of sites in the system Nk number of machines on site k Eik machine i of site k Aik(t) availability of machine Eik at time t H analysis period Ctot total cost on period H rm,ik tasks' release dates on machine Eik dm,ik tasks' due dates on machine Eik cm,ik tasks' completion dates on machine Eik fi tasks' strength vector first threshold on machine Eik second threshold on machine Eik duration related to threshold duration related to threshold failure rate of machine Eik repair rate of machine Eik Weibull parameter related to machine Eik Weibull parameter related to machine Eik Weibull parameter related to machine Eik tasks' dominance matrix

INTRODUCTION

The modelling and optimization of maintenance process of interest to academic researchers as well as industry managers. They seek optimal management solutions enabling the machines and devices to reach the highest availability and efficiency, which are necessary conditions to provide high quality products and services to customers. The increasing number of publications in this field during the last decade is a witness of this interest. In the literature, most of maintenance process models are stochastic models often based on state graphs. Most of them relate to the study of systems comprising not more than three components. As usual, the goals are to minimize economic objectives on a time horizon that is sometimes supposed to be boundless. Because of the specificity of the problem and of the various models used to describe the maintenance process, various resolution methods are used.

Whereas works of Calibra et al. [START_REF] Calibra | Optimal reliability and maintenability allocation in manufacturing systems[END_REF] and Fyffe et al. [START_REF] Fyffe | System reliability allocation and a computational algorithm[END_REF] focus on the problem of reliability allocation, in design phase, to components of machines to reduce the maintenance cost, works of Berg [START_REF] Berg | The marginal cost analysis and its application to repair and replacement policies[END_REF], Park [START_REF] Park | Optimal number of minimal repairs before replacement[END_REF] and Stadje and Zuckerman [START_REF] Stadje | Optimal repair policies with general degree of repair in two maintenance models[END_REF] aim at reducing the same cost, through reliability and availability improvement and through optimal replacement of components during exploitation of the machine. Working for the same objectives, Lam [START_REF] Lam | An optimal maintenance model for a combination of second-hand-new or outdated-updated system[END_REF] proposed a solution based on optimal combination of old and new components. Whilst the maintenance cost optimization can also be achieved through optimal resources allocation to tasks, the above mentioned approaches do not consider maintenance job management problems. On the basis of queuing theory and stochastic models, Stadje and Zuckerman [START_REF] Stadje | Optimal repair policies with general degree of repair in two maintenance models[END_REF] proposed a solution of repairmen allocation to failed machines.

A solution based on the queuing theory and a Markov model was also proposed by Derman [START_REF] Derman | On the optimal assignment of servers and a repairman[END_REF],

Frostig [START_REF] Frostig | Optimal policies for machine repairmen problems[END_REF], [START_REF] Frostig | Optimal allocation of machines to distinguishable repairmen in order to maximize some reward functions[END_REF], [START_REF] Frostig | Jointly optimal allocation of a repairman and optimal control of service rate for machine repairman problem[END_REF], Koole [START_REF] Koole | Optimal repairman assignment in two symmetric maintenance models[END_REF] and Seshadri [START_REF] Seshadri | Determination of aggregate preventive maintenance programs using production schedules[END_REF].

Among the works considering the specific problem of maintenance tasks scheduling and resources allocation, one should consider the works of Graves and Lee [START_REF] Graves | Scheduling maintenance and semiresumable jobs on a single machine[END_REF], Lee and Chen [START_REF] Lee | Scheduling jobs and maintenance activities on parallel machines[END_REF], Qi et al. [START_REF] Qi | Scheduling the maintenance on a single machine[END_REF] and Weinstein and Chung [START_REF] Weinstein | Integrating maintenance and production decisions in a hierarchical production planning environment[END_REF]. Their approaches belong to the operations research area. The common aspect of these works is the integration of maintenance policies in production management. Adzakpa et al. [START_REF] Adzakpa | On-line maintenance job scheduling and assignment to resources in distributed systems by heuristic-based optimization[END_REF], [START_REF] Adzakpa | A new effective heuristic for the intelligent management of the preventive maintenance of the distributed systems[END_REF] also proposed in this framework an approach for real-time scheduling and assignment of tasks to processors in a preventive maintenance policy. These algorithms use local optimality rules for the flow-time minimization, and are adapted to the NPhard nature of the problem -it cannot be solved within a deterministic polynomial time, using an optimal algorithm at the sense of the complexity theory -and to the problem of real-time preventive maintenance decision making.

It is well known that the management cost of the maintenance actions depends on the components' reliability and maintainability characteristics, defined in the design phase, on the number of allocated maintenance task processors and their skills, as well as on the adopted policy and maintenance strategies. The preventive maintenance cost is related to costs of the critical states and to tardiness of machines preventive maintenance activities. Most of maintenance models in the literature propose optimal periods for the preventive maintenance tasks, in order to comply with given availability objectives at a minimum cost. In general, these methods make very strong assumptions about the resources or processors availability for the maintenance tasks: resources available as soon as they are needed, instantly repaired machines and so on. Very often, the assumptions are too hard for the realistic conditions in industry. Obviously, it is more difficult in practice to face the failures of a whole system of reparable machines with a limited number of maintenance tasks processors. Indeed, a task processor is not always available at the planned date. Moreover, incidents may delay the planned maintenance actions and contribute to degrading the machines and thus increase the maintenance costs. Compared to the state of the art, this paper proposes a solution to realistic cases of system comprising N machines, where the joint problems of scheduling and real-time assignment to processors of the preventive maintenance tasks are considered. In this case, comparatively to production planning, the maintenance cost minimization implies the flow-time cost minimization with unequal release dates (tasks occurrence dates) under tardiness penalty, while assuming that the processing cost on each machine is constant on the considered time horizon. Parts of these problems have been considered in the production scheduling area ( [START_REF] Chu | A branch-and-bound algorithm to minimize total flow time with unequal release dates[END_REF], [START_REF] Chu | Efficient heuristics to minimize total flow time with release dates[END_REF], [START_REF] Chu | Some new efficient methods to solve the n/1/ri/єTi scheduling problem[END_REF], [START_REF] Chu | A branch-and-bound algorithm to minimize total tardiness with different release dates[END_REF]).

This paper presents an algorithm for real-time maintenance tasks planning in distributed systems. It enables to minimize the preventive maintenance cost, expressed as the sum of the tardiness and flow-time costs. The proposed algorithm is based on a priority rule locally optimal with proof. Section 2 presents the real-time maintenance problem considered in this paper. Section 3 deals with the priority rule in a static scheduling process, a lower bound for performance evaluation and an algorithm of the management decision making. In section 4, the priority rule is applied to the problem of real-time planning of maintenance tasks with a real-time scheduling algorithm, emphasized by experimental numerical results. Lastly, section 5 concludes the paper with some indications on the future extensions.

THE REAL-TIME MAINTENANCE PROBLEM

The aim of this work is to schedule and assign tasks to processors while minimizing the total maintenance actions cost on a finite horizon H, under availability constraint.

Formulation of the problem

One considers a distributed system composed of N sites working in parallel and sharing q processors for the preventive maintenance activities of all the machines of the system. On a given exploitation site, the machines operate in series. A machine may be a production machine or just a component. Obviously, the number of processors is less than the total number of machines in the system and the processors are shared by the whole system, Nk being the number of machines on site k. The set-up times of tasks (task preparation times) are assumed to be included in their processing time.

To model this problem, finite time horizons are considered, and unavailability and maintenance cost, pertaining to each site must be minimized within the finite horizon. The reliability models of machines in the system are assumed to be known. Perfect maintenance is assumed, i.e., at the end of each preventive maintenance action at time T the machine Eik is as new, with an availability reset to its initial value. Moreover, for a machine, one assumes the maintenance task completion time is also the time of restart. It should be recalled that, if, for a machine, the failure and repair rates are constant on the horizon H, then the reliability and availability of this machine depend on the exponential law. In this context, the availability of a machine Eik (machine i of site k) after the completion of a preventive maintenance task at time T is expressed as [START_REF] Villemeur | Reliability, Availability, Maintainability and Safety Assessment[END_REF]: [START_REF] Calibra | Optimal reliability and maintenability allocation in manufacturing systems[END_REF] where the failure and repair rates ( and ) of Eik are constant. The operational availability limit is then defined for each machine, which is converted into a threshold such that . Solving this inequality yields the duration between two consecutive maintenance tasks on Eik

(2)

Under this threshold, when a machine's working time exceeds (i.e.

), one supposes that it operates in a critical state beyond this moment with a high probability of failure. Beyond a second threshold ( ), there are additional costs to the maintenance action. These last costs are expressed as tardiness costs. The time from the completion date of a maintenance task on each machine Eik (supposed to be a restarting date of the machine) to the deadline of the next task is then ( )

å = << N k k N q 1 ( ) ( )( ) [ ] T t t A ik ik ik ik ik ik ik ik ik - + - + + + = µ l µ l l µ l µ exp ik l ik µ 1 ik a 1 ik ik A a ³ 1 ik t ú ú û ù ê ê ë é - ÷ ÷ ø ö ç ç è ae + + - = ik ik ik ik ik ik ik ik l µ l µ a µ l t 1 ln 1 1 1 1 ik t 1 ik ik A a £ 2 ik a 1 2 ik ik a a < (3)
In a more general case where machine Eik's reliability follows the Weibull law, its availability Aik is renewed at the end of a maintenance task completed at time T. This availability is then Aik(t-T) where Aik(t) is expressed with the availability function below with a constant repair rate (which is for machine Eik) [START_REF] Adzakpa | Maintenance of Distributed Systems: Methods for Real-time Decision Making[END_REF] (4)

where the Weibull parameters relative to machine Eik are the time origin parameter (denoted in general) -the scale parameter (denoted in general) -the shape parameter (denoted in general).

In this case and when the shape parameter is greater than or equal to 1, the availability function is a strictly decreasing function which tends to a limit when the time t tends to infinity. Thus it exists a unique duration (resp.

) such that (resp.

). These durations and can be computed from the availability function.

The total cost due to the critical states and to the tardiness of the preventive maintenance tasks on the considered time horizon H is then [START_REF] Stadje | Optimal repair policies with general degree of repair in two maintenance models[END_REF] where rm,ik, cm,ik and dm,ik are respectively the release date (task occurrence date), the completion date and the due date of the m th preventive maintenance task on machine Eik. The processing time of the corresponding task is supposed to be the mean time to preventive maintenance of the machine, which is .

Remark 1:

In the cost , the flow-time cost and the tardiness cost are supposed to be equally weighted. When they are unequally weighted with the weight Wf for the flow-time cost and Wt for the tardiness cost, then is replaced by

A corresponding remark in the sequel allows to still use the approach in its weighted form.

The problem is then (P1): Minimize with the q processors and , and in its weighted form, (P2): Minimize with the q processors and , .

When there is only the flow-time (response-time) cost or critical state cost, and no tardiness cost, and the machines in the system are equally weighted, the problem was dealt with in [START_REF] Adzakpa | On-line maintenance job scheduling and assignment to resources in distributed systems by heuristic-based optimization[END_REF]. The

ú ú û ù ê ê ë é - ÷ ÷ ø ö ç ç è ae + + - = ik ik ik ik ik ik ik ik l µ l µ a µ l t 1 ln 1 2 2 µ ik µ ( ) ( ) ( ) dx x x t t t A t ÷ ÷ ø ö ç ç è ae + + ÷ ÷ ø ö ç ç è ae ú ú û ù ê ê ë é ÷ ø ö ç è ae - + - - = ò - s g b b µs µs s g g µ 0 exp . 1 . exp ik g g 0 > ik s s 0 > ik b b ik b ( ) ¥ ik A 1 ik t 2 ik t ( ) 1 1 ik ik ik A a t = ( ) 2 2 ik ik ik A a t = 1 ik t 2 ik t ( ) ( ) (
)

å å å = = Î ÷ ÷ ø ö ç ç è ae - + - = N k N i H r m ik m ik m ik m ik m tot k ik m d c r c H C 1 1 / , , , , , , 0 max ) ( ik µ / 1 ) (H C tot ) (H C tot ( ) ( ) (
)

å å å = = Î ÷ ÷ ø ö ç ç è ae - + - = N k N i H r m ik m ik m t ik m ik m f w k ik m tot d c W r c W H C 1 1 / , , , , , , 0 max . . ) ( ) (H C tot 1 ik ik A a ³ { } { } k N i N k ,..., 1 , ,..., 1 Î Î ) (H C w tot 1 ik ik A a ³ { } { } k N i N k ,..., 1 , ,..., 1 Î Î
proofs of the following propositions are also based on the same reference. When the different machines in the system are weighted, the problem was considered in [START_REF] Adzakpa | A new effective heuristic for the intelligent management of the preventive maintenance of the distributed systems[END_REF]. Problems (P1) and (P2) have the properties in the propositions below. The first one is straightforward. The second one is proved as a consequence of an analogous proposition in [START_REF] Adzakpa | On-line maintenance job scheduling and assignment to resources in distributed systems by heuristic-based optimization[END_REF].

Proposition 1:

The tasks relative to machine Eik are linked by precedence relationships defined as follow. If cm,ik is the completion time of the m th preventive maintenance task on machine Eik (its m th restarting date), then the release date (resp. the due date) of the (m+1) th task is determined by (resp.

) with defined in equations ( 2) and [START_REF] Berg | The marginal cost analysis and its application to repair and replacement policies[END_REF].

Proposition 2: Problem (P1) is NP-hard.
In the solution, the tasks are processed without pre-emption, i.e. once started, each maintenance tasks is completely processed without interruption. It is based on the development in the following section on an off-line problem and subsequent properties for a real-time use.

THE SCHEDULING PROBLEM AND ITS RESOLUTION

In order to solve the cost minimization problem of real-time preventive maintenance planning and scheduling, one considers, in a first step, the static problem in this section. To this end, one seeks to schedule a set of n tasks having unequal release dates on a single processor or on parallel processors, in order to minimize, through a non-preemtive schedule, an objective function which is a linear combination of the total flow-time and the total tardiness: , where ci is the completion time of task i. Each task i has a release date ri (the time the task is released for processing), a processing time pi, and a due date di; beyond this last date, the task involves a tardiness penalty. Wf and Wt are respectively the weight of the critical states and the weight of the tardiness.

The problem is NP-hard, even in the case considered in this paper where Wf=Wt=1 is assumed. In order to solve it, some tools are provided for local optimality and used with the concepts introduced in the following subsection.

Concepts and definitions

Definition 1: Consider a pair of tasks {i, j} to be scheduled at time t at the end of a partial sequence [ ].

i) It is said that a task i dominates a task j at time t and noted , if scheduling i at the end of [ ] before j, yields a less great value of the total cost (here the total weighted sum of the flow-time and the tardiness) than scheduling j before i.

ii) A matrix denoted related to the tasks i and j is then created as follow: [START_REF] Derman | On the optimal assignment of servers and a repairman[END_REF] Remark 2: In the classical scheduling vocabulary, the term dominance is used for partial schedules. But in the above definition of dominance, it can be noticed that if at time t, [ ] is a partial schedule, saying is equivalent to saying that the sequence [ ] dominates [ ].

Definition 2: At any time t, the strength fi of a task i is defined as the number of tasks that i dominates. In other words,

1 , , 1 ik ik m ik m c r t + = + 2 , , 1 ik ik m ik m c d t + = + 1 ik t ) , 0 max( . ) .( 1 i i t i i n f d c W r c W - + - å = p j i ! p W ( ) î í ì ¹ = W otherwise 0 and if 1 , j i j i j i ! p j i ! ji p ij p (8)
Now one considers a pair of tasks {i, j} to be scheduled at time t at the end of a partial schedule [p].

Remark 3: Ties are broken with the tasks' indices. So, if , one does not have and .

Definition 3:

At any time t, the function FTR (Flow-time and Tardiness Rule) relatively to a pair of tasks {i, j} is defined as

FTR(i, j, t) = max[PRTF(i, t), Q1(i, j, t)]+max[PRTT(i, t), Q2(i, j, t)] (9)
where

PRTF(i, t) = 2.max(ri, t) + pi, PRTT(i, t) = max(ri, t) + max(max(ri, t)+ pi, di),
Q1(i, j, t) = max(ri, t) + max(rj, t) and Q2(i, j, t) = max(max(ri, t), di -pi)+max(max(rj, t), dj -pj).

Then the theorem below can be proved.

Theorem 1: At any time t at the end of a partial schedule [ ] and for any pairs of tasks , i dominates j if and only if (10)

Proof: Task i is defined by the parameters ri, pi and di. Task j is defined by the parameters rj, pj and dj.

By scheduling i and j at time t at the end of the partial sequence [ ], let's denote the additional cost involved by scheduling i directly before j, and the one obtained by scheduling j directly before i. In this case, if one poses and , then and Therefore, Now, with PRTF and Q1 defined as previously, Chu [START_REF] Chu | Efficient heuristics to minimize total flow time with release dates[END_REF] proved that

å W = j i j i f ) , ( j i ¹ j i ! i j ! p { } j i, ) , , ( ) , , ( t i j FTR t j i FTR £ p ij C ji C ) , max( t r R i i = ) , max( t r R j j { } { } )] 0 ( [ ] [ ] 0 [ ] [ , d p ) R , p R max( max r p ) R , p R max( ) , d p R max( r p R C j j j i i j j j i i i i i i i i ij - + + + - + + + - + + - + = { } { } )] 0 ( [ ] ) ( [ )] 0 ( [ ] [ , d p ) R , p R max( max r p R , p R max , d p R max r p R C i i i j j i i i j j j j j j j j ji - + + + - + + + - + + - + = { } { } { } { } { } { } { } { }. )] 0 , ) , max( ( max [ ] ) 0 , max( [ ] ) , max( [ ] [ )] 0 , ) , max( ( max [ ] ) 0 , max( [ ] ) , max( [ ] [ )] 0 , ) , max( ( max [ ] ) , ( max [ )] 0 , ( max [ ] [ )] 0 , ) , max( ( max [ ] ) , max( [ ] ) 0 , max( [ ] [ i i i j j j j j i i i j j j j j j j j i i i i i j j j i i i i i i i i j j i i i j j j j j j j j j j j i i j j j i i i i i i i i ji ij d p R p R d p R r p R p R r p R d p R p R d p R r p R p R r p R d p R p R r p R p R d p R r p R d p R p R r p R p R d p R r p R C C - + + + - + - - + + + - + - - + + + - + + - + + + - + = - + + + - + + - - + + - + - - + + + - + + + - + + - + = -
He also proved in [START_REF] Chu | Some new efficient methods to solve the n/1/ri/єTi scheduling problem[END_REF] that with the above expression of PRTT and Q2, the following equality is obtained where So, By using these different equalities, it is obtained that Hence task i dominates task j if and only if ■

Remark 4:

The FTR function is a local optimality priority rule for the cost expressed as the sum of the total flow-time and the total tardiness penalty.

Remark 5:

It can be noted that if the cost is the weighted sum of the flow-time cost and the tardiness cost, where the flow-time is weighted by Wf and the tardiness is weighted by Wt, then theorem 1 still holds if function FTR is replaced by the following priority function.

FTRW(i, j, t) =Wf .max(PRTF(i, t), Q1(i, j, t))+Wt.max(PRTT(i, t), Q2(i, j, t)

Therefore, all the development that follows can still be used to solve the real-time maintenance tasks scheduling problem. These definitions and properties of the FTR function, combined with the dominance matrix and the task strength, are then used in a scheduling algorithm that minimizes the objective function. This algorithm, designed for off-line implementation, is summarized as follow.

FTR-based algorithm

Algorithm 1: Off-line use of the FTR with one processor

The decision instant t is the time when the processor is available the after the completion of a task k. The decision is made through the following steps.

1) The search set .

2) Compute FTR(i,j,t) for all pairs of tasks where and .

3) Compute the dominance matrix in the set St.

4) Compute the tasks' strength (the vector containing fi for all the tasks not yet scheduled) in the set St .

)

. ( ) ( )) , , ( ), , ( max( ] ) , max( [ ] [ 1 j i j i j j j i i i i i r r p p t j i Q t i PRTF r p R p R r p R + - + + = - + + + - + ) d d ( ) p p ( ) t , i ( PRTT , Q max , d p ) R , p R max(max( ) , d p R max( j i j i ij j j j i i i i i + - + + = - + + + - + ) ( )] 0 [ ] 0 [ ). p p ( ) t , j , i ( Q ) p p ( ) p d , R max( ) p d , R max( ) d , p R max( ) d , p R max( Q j i j i j j j i i i j j j i i i ij + + = + + - + - = + + + = 2 ). d d ( ) p p ( )) t , i ( PRTT ), t , j , i ( Q max( ) , d p ) R , p R max(max( ) , d p R max( j i j i j j j i i i i i + - + + = - + + + - + 2 ] 0 [ ] 0 [ { } ). , , ( ) , , ( )) , ( ), , , ( max( )) , , ( ), , ( max( )
) , ( ), , , ( max( )) , , ( ), , ( max( 2 1 2 1 t i j FTR t j i FTR t j PRTT t i j Q t i j Q t j PRTF t i PRTT t j i Q t j i Q t i PRTF C C ji ij - = + - + = - ) t , i , j ( FTR ) t , j , i ( FTR £ W { } scheduled yet not tasks all S t ¬ { } j i, t S i Î t S j Î W 5
) Evaluate the subset of tasks having the greatest strength.

. Until card(St)=1, go to step 4.

These steps are repeated until all tasks are scheduled.

Remark 6: This algorithm can be adapted to the case of q ( ) parallel processors. In this case, the processor available the earliest is used at each decision date.

Proposition 3:

The algorithm based on the FTR function and the tasks' strength has an O(n 3 ) complexity if n is the number of tasks to be scheduled.

Proof. The proof of this proposition is straightforward.

A lower bound to the cost enables evaluating the performance of this algorithm and the priority rule on one processor. This lower bound is based on the SRPT rule (shortest remaining processing time rule) which gives a lower bound to the mean flow-time cost in presence of the unequal release dates of tasks and a lower bound to the mean tardiness cost in presence of the unequal release dates developed in [START_REF] Chu | A branch-and-bound algorithm to minimize total tardiness with different release dates[END_REF] using modified due dates.

Lower bound

To obtain the lower bound, the non-preemtion assumption is relaxed and so interruption is permitted when the tasks are being processed. This is obtained in the following steps.

1) Determine the series (

) of the tasks' due dates sorted in their non decreasing order (modified due dates).

2) Schedule the preemptive tasks according to the SRPT rule. In the SRPT rule, tasks are scheduled in the non decreasing order of their remaining processing time.

The quantity is then a lower bound to the cost expressed as the sum of the flow-time and the tardiness. In this expression, [i] is the task completed in the i th position and its completion time. This lower bound is due to the fact that the SRPT rule gives a lower bound to the flow-time in presence of tasks with unequal release dates and, as proved in [START_REF] Chu | A branch-and-bound algorithm to minimize total tardiness with different release dates[END_REF], the modified due dates combined with the SRPT rule gives a lower bound to the tardiness in presence of unequal release dates. Before using this rule to solve the real-time maintenance problem, figure 1 shows some results obtained on two sets of data among others. On the curves, MC denotes the average cost (including the flow-time cost and the tardiness cost) obtained by using the FTR function and the dominance matrix. LB_MC is the lower bound to the average cost. Delta_MC is the difference (MC -LB_MC). These results, among others, reveal a very good performance of the rule. Now, one considers the original problem which is the real-time preventive maintenance tasks scheduling and assignment to processors in order to minimize the total cost. It is proposed, in the following section, a solution based on the FTR function and the W matrix.

t S 1 t t S S ¬ 2 ³ q ' n ' ' d ,..., d , d 2 1 ( ) ( ) å å
= = - + - n i ' i ] i [ n i ] i [ i ] i [ , d c max r c 1 1 0 ] i [ c - 0,

REAL-TIME APPLICATION AND EXPERIMENTATION

In addition to the development in the previous section, the properties below are used in the real-time decision-making process.

Additional properties

The following theorem allows the use of the FTR function in accordance with the precedence constraint linking the tasks on a machine.

Theorem 2: By using theorem 1, the precedence relationships linking the tasks on a given machine Eik are also respected. In other words, for any pair of integers m1 and m2 with m1 < m2, and for any machine Eik, theorem 1 allows to schedule the task on Eik before the .

Proof:

The proof is based on the FTR function's definition.

Consider two tasks i and j whose processing times, release dates and due dates are such that pi = pj, ri rj and di dj. The FTR function defined in equation ( 9) is

FTR(i, j, t) = max[PRTF(i, t),Q1(i, j, t)]+max[PRTT(i, t),Q2(i, j, t)]
where PRTF(i, t) = 2.max(ri, t) + pi, PRTT(i, t) = max(ri, t) + max(max(ri, t)+ pi, di), Q1(i, j, t) = max(ri, t) + max(rj, t), and Q2(i, j, t) = max(max(ri, t), di -pi)+max(max(rj, t), dj -pj).

Then, PRTF(i, t) PRTF(j, t) and PRTT(i, t) PRTT( j, t).

As Q1(i, j, t) = Q1(j, i, t) and Q2(i, j, t) = Q2(j, i, t), then the inequality FTR(i, j, t) FTR(j, i, t) is satisfied. Now for two integers m1 and m2 with m1 < m2, the task and the task on a machine Eik are supposed to have the same processing time Mtpik and they satisfy and . So the FTR function gives the priority to the task.

The problem then becomes establishing priorities between the different tasks on machines in the system.

Definition 4: Let t be a decision-making time and denote, at time t, the duration between t and the end of the last preventive task on machine Eik before that time. Let Ut (set of urgent tasks) denote at time t the subset of tasks satisfying condition .

On the basis of these different definitions and properties, one proposes, below, a real-time algorithm for scheduling and assignment of tasks to processors in a distributed system. The algorithm uses the

th m 1 th m 2 £ £ £ £ £ th m 1 th m 2 ik m ik m r r , , 2 1 £ ik m ik m d d , , 2 1 £ th m 1 ) (t ik D t ik ik 1 ) ( t ³ D
FTR function that takes into account the critical state costs and the tardiness costs, the dominance matrix , the strength vector, the urgency of tasks and theorem 2.

Real-time decision making

The real-time decision making consists in establishing priorities between the tasks, as they arrive in the system, and to assigning them to the processors. The algorithm, called OL-MTSA-2T, is summarized as follows.

Algorithm 2: OL-MTSA-2T

At the beginning of the planning horizon H (t = 0): Determine the first critical dates and the due dates of the preventive maintenance tasks on each machine in the system (they are equal to and respectively for machine Eik).

The decisions are then made through the steps below.

1) Consider the time t when a processor is available the earliest.

If two or more processors are available at time t, select the one with the smallest index.

While t<H, go to step 2.

Else, go to step 5.

2) Compute the subset Ut.

3) If then the search set .

Else 4) Apply to the set St the steps described in algorithm 1 based on the FTR function for the tasks scheduling, to select the task that should be scheduled at time t.

Determine its completion time.

Determine the next release date of preventive maintenance task on the corresponding machine by adding the quantity to its completion time.

Determine the next due date by adding its quantity to its completion time.

Go to step 1.

5) End.

The following subsection gives some experimental results highlighting the performances of the approach.

Numerical experimentation

The algorithms based on the FTR rule are programmed in C language and tested on a Compaq AlphaServer ES40 DEC6600 workstation operating under UNIX with 2048 MB of RAM memory. Two versions of the OL-MTSA-2T algorithms have been programmed. The first one uses the urgency criterion while the second one does not. Both programs are implemented on a 365-day horizon.

Experimentations are performed on systems with an overall number of 500 machines. The failure rates data ( ) and the repair rates data ( ) are simulated according to normal laws and the processing time (MTPik) of tasks have been computed on the basis of the repair rates. The first availability thresholds are uniformly generated in the interval ]0, 1[ with respect to the . The durations are then computed on the basis of these thresholds. Instead of the second thresholds ( ), we generated the durations ( ) which are the duration from the end of a maintenance task on machine Eik to the deadline of the next task on the same machine. To obtain , we generated variables uniformly in [0.25, 0.5]. Then is computed as . The total number of processors is varied from 2% to 20% of the total number of machines in the system. It is noted that in the manufacturing industry, the number of processor is scarcely greater than . On graph (a) of figure 2, MC is the average cost of the preventive maintenance tasks undertaken on the 365-day horizon H. It corresponds to the average cost relating only to the machines whose tasks have been processed on the horizon H. MRC is the average cost of all the preventive tasks that have been needed on the horizon H (including machines whose maintenance has not been able to be processed on the horizon). The _U (resp. _NU) extension to the curve names denotes the result obtained with (resp. without) considering the urgency criterion.

Fig. 2. Results on 500 machines

From this figure, one can observe the following. 1) The average cost of the tasks that are processed on the time horizon (MC) is almost the same whatever the tasks' urgency is used or not. But, on the overall set of tasks, including those that have not been able to be processed on the horizon, the average cost (MRC) is far better when the urgency is used. It represents the average cost involved for the whole system on the horizon including the cost of critical states and of tardiness. Moreover, on the enlarged part of the graph (a), one can see that when the number of resources or processors increases, the MC and the MRC (in the case where the urgency is used) converge from about 10%. This means that from a given number of resources, all the preventive tasks on the horizon are processed when the urgency is used (while minimizing their total cost, which is equivalent to minimizing the average cost). This is not the case when the urgency is not used.

2) Graph (b) of figure 2 shows the average time of processors utilisation on the 365-days time horizon. When the urgency criterion is considered, the processors are used at their full capacity as long as all the maintenance tasks are not processed. The average time per processor under the urgency criterion begins to decrease only from 10%, when all the tasks are processed. On the contrary, the mean time per resource continually decreases as the number of resources increases when the urgency is not used, even though there remain tasks that are not processed.

All these remarks emphasize the effectiveness of the priority rule expressed by the FTR function, when it is combined to the urgency criterion of tasks for making the preventive maintenance decisions in real time.

CONCLUSION

In this paper, we proposed a real-time algorithm for the scheduling of the preventive maintenance tasks of a distributed system comprising several sites, and each one containing a set of machines working in series. The main objective is to minimize the maintenance cost on a given horizon, while ensuring the system an availability higher than a given threshold. The maintenance cost comprises the cost of the time spent by the machines operating in critical states (expressed as the flow-time cost or response time cost), the cost due to the maintenance tasks tardiness and tasks processing cost. For this purpose, we propose a local optimality priority rule (FTR) using a priority rule for the flow-time minimization and a priority rule for the tardiness' minimization. This priority rule considers task pairs at decision times. This priority rule has a polynomial complexity of O(n 3 ) when used for off-line scheduling. This polynomial complexity enables to use it as well in production scheduling as in real-time decision making of preventive maintenance tasks on large systems. We also proposed a lower bound which reveals the good performance of the rule with one processor.

Basing on this rule and other concepts introduced in the paper, we proposed an algorithm for the real-time decision making in the preventive maintenance activities management of distributed systems.

Among the different assumptions, we considered exponential availability law for the machines incorporated in the system. This assumption is not restrictive. The method described in this paper can be used in the case of other laws. We also supposed that the processors are identical. In the extension to this work, non-identical processors will be considered.

Fig. 1 .

 1 Fig. 1. The performance of the FTR rule relatively to the lower bound.

1

 1


Lisseth Gavilan

Laurent Remusat

Mathieu Roskosz

Horia Popescu

Nicolas Jaouen

Christophe Sandt

Cornelia Jäger

Thomas Henning

Alexandre Simionovici

Jean Louis Lemaire

Denis Mangin

Nathalie Carrasco

X-Ray-induced Deuterium Enrichment of N-rich Organics in Protoplanetary Disks: An Experimental Investigation Using Synchrotron Light

Keywords: X-rays, stars -Planetary Systems, protoplanetary disks -ISM, dust, extinction -Methods, laboratory, solid state -Techniques, spectroscopic

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The study of organic molecules, dust and ices is essential to understand the origin and evolution of the solar system. The deuterium to hydrogen ratio (D/H) is used as a proxy of the surrounding physico-chemical conditions affecting organics, omnipresent in inter-and circumstellar media (ISM/CSM), and occasionally delivered to Earth and other planets in the form of meteorites. Kerridge (1983) proposed early on that deuterium enrichment in meteorites had an interstellar origin. Yet interstellar clouds can further evolve, reaching a critical mass and density, collapsing, and forming protoplanetary disks (Williams & Cieza 2011) on par with our own protosolar nebula (PSN). Isotopic fractionation can occur at any of these stages via multiple chemical pathways such as isotopic exchange reactions in the ISM (Roueff et al. 2015) and gas-surface reactions in protoplanetary disks (Henning & Semenov 2013).

Deuterated molecules have been detected in all stages of star formation; for a review, see Ceccarelli et al. (2014). In molecular clouds, detected deuterated species include DCO + and DNC (van der Tak et al. 2009), H 2 D + and HD + 2 (Parise et al. 2011) among many others. In proto-stellar cores, detections include DCN and D 2 CO (Roberts et al. 2002;Ceccarelli et al. 2002), CH 3 CCD (Markwick et al. 2005), CH 3 OD (Parise et al. 2006), NH 2 CDO and DNCO (Coutens et al. 2016). In protoplanetary disks such as TW Hydrae, the most wellstudied protoplanetary disk because of its proximity and near face-on orientation, detections include DCO + (van Dishoeck et al. 2003;Qi et al. 2008), H 2 D + (Ceccarelli et al. 2004) and DCN (Qi et al. 2008;Öberg et al. 2012). In the young protoplanetary disk surrounding the solartype protostar DM Tau, detections of H 2 D + (Ceccarelli et al. 2004) and HDO (Ceccarelli et al. 2005) were reported but later contested by Guilloteau et al. (2006), along with the confirmed detection of DCO + . Chemical models including a vast network of deuterated species have pointed at the PSN origin of water found in comets and asteroids (Albertsson et al. 2014), although modeling work by (Cleeves et al. 2014) argued that ion-driven deuterium pathways in the disk are inefficient, implying an interstellar origin for water. Therefore, the dominance of deuterium fractionation at a particular stage in stellar evolution persists as an open question in astrophysics.

In the solar system, deuterated species have been measured in comets and meteorites. Recent in situ measurements in the Jupiter family comet 67P/Churyumov-Gerasimenko by the Rosetta spacecraft found the D/H in water to be approximately three times the terrestrial value (Altwegg et al. 2015), suggesting an asteroidal rather than cometary origin for the Earth oceans and atmosphere. This contrasts with reports on the D/H ratio of the Jupiter-family comet 103P/Hartley 2 (Hartogh et al. 2011), compatible with the value of Earth ocean-like water. These measurements point at mixing processes in the early PSN. Analysis of meteorites have shown deuterium enrichment of their insoluble organic matter (IOM) compared to the solar hydrogen abundance (Robert & Epstein 1982;Yang & Epstein 1983), and suggested an interstellar inheritance (Sandford et al. 2001;Busemann et al. 2006). D-enrichment has also been measured in interplanetary dust particles (IDPs) (Messenger 2000) and more recently in UltraCarbonaceous Antarctic Micrometeorites (UCAMMs) (Duprat et al. 2010). However, isotopic anomalies in meteorites and comets are not necessarily fingerprints of an interstellar origin: organic matter can be further modified in the coldest outer regions of protoplanetary disks (Remusat et al. 2006;Aléon 2010), e.g. via isotopic gasgrain exchange reactions (Remusat et al. 2009).

Laboratory experiments have the potential of unveiling the chemical and physical mechanisms at the origin of isotopic fractionation. Thermal processing experiments can impact the chemical and structural properties of astrophysical dust and ice analogs, e.g. (Allamandola et al. 1988;Ehrenfreund et al. 1999;Collings et al. 2004). However, non-thermal irradiation processes are as important in astrophysical environments where elevated dust or gas densities inhibit thermal annealing. Ionizing radiation includes charged particles (electrons or protons), ions, but also UV photons and Xrays (Feigelson et al. 2002). Recent experiments have investigated the effects of electron irradiation on the structure and the D/H signature of kerogens (natural organic macromolecules) (Le Guillou et al. 2013), synthetic polymeric organics (Laurent et al. 2014(Laurent et al. , 2015) ) and in analogs of protoplanetary hydrous silicate dust (Roskosz et al. 2016). These experiments have revealed that non-thermal irradiation can offer viable pathways to deuterium enrichment in the solid phase.

X-rays from T-Tauri stars are associated to the energy dissipation of magnetic coronal fields (Preibisch et al. 2005). Together with UV photons, these can irradiate the disk surface, mid-plane and outer edges promoting photochemistry and radiolysis (Glassgold et al. 1997;Henning et al. 2010;Aresu et al. 2011). As of now, no laboratory investigations have considered the possible isotopic effects of X-ray irradiation on organic matter. Yet, we know hard X-rays (> 5 keV) can strongly impact the structure of silicates in protoplanetary disks (Gavilan et al. 2016). Thus, our goal is to experimentally unveil the isotopic effects of stellar soft X-rays (0.1 -2 keV) on protoplanetary organic matter.

To this end, we employ nitrogen-rich amorphous organic solids as analogs to refractory organics produced from the thermal and/or photo-processing of N-rich ices in regions beyond the snow-line of protoplanetary disks (Greenberg et al. 1995;Ciesla & Sandford 2012). To simulate the T-Tauri X-ray field, we perform controlled experiments using synchrotron radiation in the soft Xray range. The brilliance of synchrotron X-rays is ideal to perform highly-focused monochromatic studies with remarkable reliability, allowing us to simulate astrophysical fluences. Tholins, complex organic compounds rich in carbon and nitrogen used as analogs to astrophysical solid matter (Sagan & Khare 1979), were produced using the PAMPRE setup, a low pressure (0.95 mbar) radiofrequency (RF) plasma reactor located at LATMOS (Guyancourt, France) (Szopa et al. 2006). A 13.56 MHz RF power source tuned at 30 W generates a capacitively coupled plasma (CCP) fed by a gas mixture of N 2 :CH 4 = 95:5, conditions chosen for optimal tholin production (Sciamma-O' Brien et al. 2012). The plasma was produced within a vertical cylindrical cage and the substrates were placed on the bottom electrode. The produced films, rich in carbon and nitrogen, i.e. N/C = 0.2, H/C ∼0.1 (Carrasco et al. 2016), make them comparable to the organic matter in UCCAMs with N/C = 0.05 -0.12 (Dartois et al. 2013) and to refractory organics that may have formed from the energetic irradiation of nitrogen-rich ices in the outer solar system.

Two types of tholins with different thicknesses were prepared. A thin film (∼300 nm) on an X-ray transparent substrate was prepared for X-ray absorption spectroscopy with the goal of obtaining absolute absorption coefficients in the soft X-ray range. Thick films (∼1200 nm) were prepared for irradiation with the goal of absorbing most of the deposited X-rays within the sample volume. The thin tholin was produced after turning on the plasma for 30 minutes. It is evenly deposited on both the Si 3 N 4 membrane and the supporting frame made of silicon (5 × 5 mm). For this sample, ellipsometry is performed on the reflective Si frame. The Si 3 N 4 itself is 30 nm thick and quasi-transparent to X-rays (Törmä et al. 2014). Thick tholin films were deposited on single-side polished silicon substrates (9 × 4 mm) after leaving the plasma on for 150 minutes.

To determine the film thickness, we used spectroscopic ellipsometry (SE). This is an optical technique where a beam of polarized light is obliquely sent towards a film deposited on a reflective surface. Light is elliptically polarized upon reflection. Both tan(Ψ), the ratio of the incident to the reflected light amplitude, and ∆, the relative phase change for p-and s-polarized light components, are measured, as seen in Fig. 1. These determine the absolute value of the complex reflectance ratio, ρ, defined as,

ρ = r p r s = tan(Ψ)e (i∆) (1)
where r P and r S are the complex Fresnel reflection coefficients of p-and s-polarized light. Using ρ and an appropriate optical model we can find the χ-minimized dielectric properties, layer (or multi-layer) film thickness, and roughness. For this study we used a rotating spectroscopic ellipsometer (M-2000V, J.A. Woollam), which allows all wavelengths to be measured simultaneously. The light source is a tungsten lamp covering the 370 -1000 nm range, incident at an angle of 70 • and with a spot size at the sample of ∼2 × 3 mm. Modeling of the ellipsometric parameters was done with the Complete-EASE T M ellipsometry software. We selected a Tauc-Lorentz model, suitable for amorphous organic films (Jellison & Modine 1996), providing the thickness d, surface roughness r, E g (the Tauc optical gap), and the mean squared error (MSE) between the modeled and the experimental data. The modeled parameters are listed in Table 1, showing the consistent gap value (and thus homogeneity) of films prepared under the same conditions but at different durations, in agreement with previous measurements (Mahjoub et al. 2012).

Table 1. Ellipsometric parameters for two tholin films used for XAS and irradiation respectively. We performed X-ray absorption spectroscopy (XAS) measurements at the COMET experimental station, an ultra-high vacuum chamber placed in the intermediate focal point of the elastic branch of the SEX-TANTS beamline (Sacchi et al. 2013) of the SOLEIL synchrotron. The samples were placed on a motorized sample holder perpendicular to the X-ray beam. Four different gratings were used: 100 -300 eV, 300 -500 eV, 500 -800 eV and 800 -1200 eV with resolution E/∆E >10 4 . For each range two samples were measured: the thin tholin on a Si 3 N 4 membrane and a bare Si 3 N 4 membrane. The optical depth is obtained using the Beer-Lambert law,

τ (ν) = -ln I(ν) I 0 (ν) , (2)
where I(ν) is the intensity of light at a given frequency, ν, transmitted through the tholin film deposited on the membrane and I 0 (ν) is the light intensity through the membrane alone. The imaginary index of refraction can be obtained via,

k = τ (ν) 4πνd (3)
where d, the film thickness, is found by ellipsometry as aforementioned. Similarly, the penetration depth (L P ) of a photon is given by,

L P (µm) = λ(µm) 4πk (4)
The resulting k and L P (µm) in the soft X-ray range for the tholins is presented in Fig. 2. Fluctuations in the k value around 100 eV are due to Si L edges from the Si 3 N 4 substrate. In addition to the expected C Kedge (390 eV) and N K-edge (397 eV) absorptions, we also notice the O K-edge (530 eV) and the Fe L-edges (720 eV). These are attributed to pollution and sputtering of the confining plasma cage where tholins are produced. Both k and L P are the required parameters to optimize the irradiation experiments on thick tholins. While 0.5 keV X-rays penetrate the tholin film up to a thickness of ∼625 nm (depositing most of their energy within the 1.16 µm film), 1.3 keV X-rays have a penetration depth greater than 5 µm (i.e. only a fraction of their energy will be deposited within the film, the other will be scattered). In this paper, we report the incident dose on each experiment and not the adsorbed dose, which varies depending on the chemical nature of the material (reflected in its absorption index) and on its thickness.

X-ray irradiation

We performed the X-ray irradiation experiments at the COMET station using two monochromatic settings, chosen to evaluate the impact of irradiation across the soft X-ray region: 0.5 keV and 1.3 keV. The beam was right-circularly polarized for both experiments with a spot size at the sample of 80 × 50 µm, flux density ∼10 17 photons s -1 cm -2 and resolution of E/∆E > 10 4 up to 1200 eV. For each energy, tholin films were irradiated at spots separated by 300 µm for different durations: 10 minutes, 30 minutes, 1 hour, 4 hours and 10 hours.

Following irradiation, the sample surface was examined with a white-light interferometer (smartWLI ), which allowed us to record 3-D topographies with depth resolution ∼3 nm. The result is shown in Fig. 3, where X-irradiated regions are manifested by volume expansion seen as elevated craters on the tholin surface (with rising edges up to 200 nm from the surface), reticulation and darkening. The irradiated spots extend up to ∼200 × 200 µm 2 , i.e. twice as large as the incident beam size, revealing the diffusion of secondary electrons during irradiation. The central spot, corresponding to the highest X-ray dose, extends even farther out than all others.

NON-DESTRUCTIVE DIAGNOSTICS

We performed ex-situ diagnostics of the irradiated tholins via infrared and Raman microspectrometry on the 0.5 keV irradiated tholin. These diagnostics allow us to further quantify the chemical and/or structural evolution of the irradiated tholin films, revealed by optical interferometry in the previous section. Access to infrared and Raman microscopes was kindly provided by the SMIS beamline (Dumas et al. 2006) at SOLEIL.

Infrared Microscopy

Infrared spectra of irradiated and non-irradiated regions were measured to quantify the chemical evolution of the irradiated films at increasing X-ray doses. We used a iN10 microscope (ThermoFisher Scientific) coupled to a FTIR spectrometer in offline mode. Areas of 100 × 100 µm 2 within the irradiated regions were integrated from 4000 to 1000 cm -1 . A polynomial baseline was subtracted from the infrared spectra to normalize the continuum arising from rough surface scattering and from broad ultraviolet absorptions typical of amorphous carbon materials reaching the near-infrared (Khare et al. 1984;Bernard et al. 2006).

Fig. 4 shows the evolution of the infrared vibrational bands from the initial (non-irradiated) tholin to the increasingly irradiated regions. This is detailed in Fig. 5, where infrared spectra have been split into three main regions: (a) the 3600 -2500 cm -1 region, which includes stretching modes of amine groups (primary/secondary), superimposed by the methyl (CH 3 ) and methylene (CH 2 ) stretching modes, (b) the 2400 -2000 cm -1 region, which includes nitriles and isonitriles (C≡N/N≡C) stretching modes and (c) the 1800 to 1000 cm -1 region which includes several overlapping bands including the sp 2 aromatic C=C and C=N stretching modes, sp 3 -CH 2,3 bending modes, and sp 3 -like C-C or C-N stretching of defects (such as non-hexagonal rings of aromatic compounds (Carpentier et al. 2012)). For a complete attribution to vibrational modes of tholins in the mid-infrared we refer to Imanaka et al. (2004). To analyze the kinetic evolution of the main vibrational modes, we performed the deconvolution of the continuum subtracted 3.4 µm band (3100 -2800 cm -1 ) and on the nitrile bands (2300 -2000 cm -1 ), presented in Fig. 6. Integrated optical depths were fitted as a function of the irradiated dose with the following exponential relation,

y = y o + Ae -kE dose , (5)
where y corresponds to the integrated optical depths , y 0 corresponds to an asymptotic steady-state, A is a constant, k is the rate constant , and E dose is the irradiation dose. Fig. 7 shows the infrared band kinetics during X-ray irradiation. The fitted parameters are reported in Table 2. We note the rapid depletion of the C-H stretching band (2800 -3000 cm -1 ) by the 95% loss of its original intensity and the depletion of the N-H band and the R-N≡C at 2175 cm -1 band both by the ∼67% loss of their original intensity following the maximum X-ray doses. Before the continuum baseline correction, we note that the absorbance at 4000 cm -1 increases from 0.7 to 0.9 at the highest X-ray dose, signaling darkening of the tholin film.

We also note the growth of specific band intensities during irradiation, such as the vibrational mode at 2210 cm -1 , attributed to -CN groups attached to aromatic molecules (Khare et al. 1981;Sagan et al. 1993), and the sp 3 carbon defect band between 1250 -1300 cm -1 , affecting the entire bending mode region between 1800 - 1000 cm -1 . Nitriles are not completely lost but the network bearing them becomes unsaturated and more aromatic. For the longest exposures (4 and 10 hours), a new peak appearing at 2335 cm -1 is attributed to isocyanate R-N=C=O modes, likely the result of photo-oxidation of impurities in the tholin (noticed by the weak O edge in the X-ray absorption spectra, Fig. 2). In Fig 8, we further compare the infrared spectra of our X-irradiated tholins to a thermally annealed tholin (Bonnet et al. 2015) and to a UCAMM (Dartois et al. 2013). While the X-irradiated tholin does not present the high loss of amines as in the other two, the growth of the defect band near 1300 cm -1 occurs in the three cases. Raman spectra of irradiated and non-irradiated regions were measured to quantify the structural evolution (size of aromatic units, graphitization/amorphization) of the irradiated films at increasing X-ray doses. These were measured with a DXR Raman spectrometer (Ther-moFisher Scientific) using a 514 nm excitation laser, with laser power ≤ 0.3 mW focused with a 50× objective and a spot diameter of 2 -3 µm, well within the 200 × 200 µm 2 irradiated regions (Fig. 3).

Raman microspectrometry

The main features in the Raman spectra of carbonrich organic materials are the so-called G and D peaks, around ω G = 1560 and ω D = 1360 cm -1 dominated by sp 2 sites, as visible excitation resonates with π states. We note that the non-irradiated tholin produces an intense fluorescence background in the Raman signal, an issue that could be avoided via photobleaching (which inherently transforms the material) or by employing an ultraviolet excitation source (Bernard et al. 2006). This was not a problem for the X-irradiated, nitrogendepleted, regions. The baseline corrected Raman spectra of the irradiated tholins are shown in Fig. 9, normalized to the intensity of the G band. For the nonirradiated tholin, measurement of the Raman spectra was hindered by strong fluorescence. Raman spectra were deconvolved following the prescription in Ferrari & Robertson (2000), where the G band is fitted with a Breit-Wigner-Fano profile and the D band is fitted with a Lorentzian profile. This procedure provides Raman parameters including the position of each band, their full-width at half maximum (FWHM), and the ratio of their integrated intensities (I D /I G ). These parameters are listed in Table 4.

As proposed in Ferrari (2007), the Raman parameters for completely disordered, amorphous carbons consisting of distorted sixfold rings are related as follows,

I D /I G = C L 2 a , (6)
where C'(514 nm) = 0.0055, and L a is the average correlation length of polyaromatic units in Å. For the X-ray irradiated tholins, the size of polyaromatic units remains almost unchanged, i.e. L a increases very slightly after 30 minutes irradiation and barely decreases from 1.27 ±0.02 to 1.25 ±0.02 nm after 10 hours, implying that no larger polyaromatic structures are created. As the X-ray irradiation dose increases, the intensity of the D band decreases with respect to the G band, while both the FWHM-G and FWHM-D increase. The broadening and position shift are attributed to an increasing number of defects, including bond length and angle disorder at the atomic scale resulting from cross-linking due to new bond formation. This leads to the reorganization in the tholin structure manifested macroscopically as reticulation and volume expansion (Fig. 3). In Fig. 10, we compare the evolution of Raman parameters of tholins undergoing X-ray irradiation (this study) to tholins that are thermally degraded (Bonnet et al. 2015). For the latter, heating leads not only to carbonization (preferential expulsion of heteroatoms) but also to the creation of polyaromatic structures, while in our case polyaromatic structures are almost unaffected by X-ray irradiation. However, X-ray irradiation may favor the creation of isolated over polyaromatic units, as signaled from the infrared analysis of the nitrile modes in Section 3.1.

DESTRUCTIVE DIAGNOSTICS: SIMS

We used Secondary Ion Mass Spectroscopy (SIMS) to characterize the evolution of elemental abundances of the irradiated tholins at 0.5 keV and 1.3 keV and to detect possible isotopic effects. SIMS is a destructive diagnostic, as seen in Fig. 11. We used the SIMS IMS-7f (CAMECA), installed at the Institut Jean Lamour (Nancy, France), a double-focusing magnetic sector SIMS that allows detection of isotopes at high mass resolution at the scale of a few microns. We used a highenergy primary ion Cs + source (impact energy of 15 keV) for efficient sputtering with a beam current of 0.39 nA. The mass spectrometer was tuned to obtain a mass resolving power of m/∆m = 3000. We collected and quantified secondary ions of H -, D -, and 12 C -, from the irradiated and non-irradiated tholin spots.

Before SIMS operation, the irradiated tholin films were gold-coated (thickness ∼15 nm). Samples were inserted in the analysis chamber operating at 10 -9 mbar. The ion beam was rastered over 50 × 50 µm 2 and, by using an aperture diaphragm, only the inner 35 × 35 µm 2 were sampled, to reduce the contribution of surface contamination to the H signal. The analyzed area was small enough to accurately probe the 200 × 200 µm 2 irradiated spots. In between each measurement of irradiated spots, non-irradiated areas on the same tholin film were measured and constitute reference measurements (Fig. 11). Analyses were saved as depth profiles. Sputtering was stopped as soon as the substrate was reached. The sputtering depth is calibrated using the film thickness obtained via ellipsometry, described in Sec. 2.1. The D/H profiles for irradiated spots are reported as isotope fractionation, defined as,

α D/H = D/H i D/H o (7)
where D/H o is the D/H ratio of the reference tholin measured near the irradiated area and D/H i is the ratio after a given X-ray dose. Similarly, the abundance of H/C for an irradiated spot is parametrized as,

α H/C = H/C i H/C o (8)
where H/C o corresponds to the reference tholin and H/C i is the value after a certain irradiation dose. The For 1.3 keV the α H/C profile after 10 minutes irradiation was discarded due to large systematic errors assigned to localized hydrogenated inhomogeneities on the tholin surface (fractures/wedges), as seen in Fig. 3.

irradiation dose is defined as,

E dose = k × F, ( 9
)
where k is the absorption coefficient (Fig. 2) and F is the fluence obtained from the total flux integrated in time and divided by the irradiated area. In Fig. 12 (a) we show the α D/H profiles for the 0.5 keV irradiated tholin, which reach a maximum after 1 hour irradiation, corresponding to α D/H = 1.9. Longer exposures do not induce further D-enrichment. We also note that for shorter irradiation times (10 and 30 minutes) the deuterium enrichment decreases linearly in depth. This is evidence to the cumulative dose effect: the regions where the enrichment is largest, hence more affected by irradiation, are found near the surface, while towards the tholin/substrate interface the enrichment decreases and becomes negligible. We hypothesize that the bulk of the tholin is shielding the lower regions from irradiation for short duration experiments. From Fig. 2, 0.5 keV photons will have a maximum penetration depth in tholins at ∼625 nm. During irradiation, the released secondary electrons (photoelectrons or Auger electrons) will be diffused completely within the tholin volume. On the other hand, in Fig. 12 (b) we note the linearity of the α D/H profiles for the tholin irradiated at 1.3 keV. At this energy the maximum penetration depth is >5 µm, i.e. for our 1.16 µm thick samples only ∼20 percent of the incident energy will be absorbed. The penetration length of X-rays in a material depends on the X-ray energy and the absorption index of the material, both affecting the total energy absorbed. For our 1.1 µm tholin, lower energy X-rays are absorbed more efficiently. This explains the lower fractionation α D/H = 1.7 for 10 h irradiation at 1.3 keV compared to 0.5 keV, and the need for a higher dose to reach the enrichment plateau at 1.3 keV.

Whereas at 0.5 keV the maximum isotope fractionation is reached after 1 h and does not evolve for longer irradiations, at 1.3 keV isotope fractionation continuously increases from 1 h to 10 h irradiation. In addition, the 10 minute experiment at 1.3 keV induces a negligible D-enrichment. Hence, higher energy X-rays appear less efficient drivers of H-isotope fractionation. The α H/C profiles are reported in Fig. 12 (c) and (d), and are related to the evolution of the molecular structure during X-ray irradiation. For the shortest irradiation at 0.5 keV, tholins exhibit a modest Hdepletion at the surface and a negligible effect near the tholin/substrate interface: α H/C = 0.8 at the surface and increases regularly up to 1 near the interface. Hence, trends for α H/C and α D/H as a function of depth are opposite: D-enrichment is associated to H/C decrease (i.e. loss of H). As the dose increases, dehydrogenation of the organic skeleton increases and becomes more homogeneous across the tholin thickness. At 0.5 keV, α H/C reaches a plateau after 1 hour, while at 1.3 keV the plateau is reached after 4 hours. Fig. 13 shows the average α D/H versus the average α H/C parameters, for both the 0.5 and 1.3 keV irradiated tholins. These are found from the average α D/H and α H/C values over the entire film thickness in Fig. 12, and the scatter in this value determines the error bars. The strong correlation between these two parameters confirms that the preferential H over D loss (a reflection of their differences in zero-point energies) during X-ray irradiation is at the origin of the measured deuterium enrichment. Fig. 14 shows the evolution of the deuterium enrichment, α D/H , as a function of the total X-ray dose for both 0.5 and 1.3 keV. We note that the D-enrichment is slightly larger for 0.5 keV X-rays, which depose a larger energy per volume than 1.3 keV X-rays. We compare our results to the evolution observed for 30 keV electron irradiation of organic polymers (Laurent et al. 2015). We first remark that the deposited dose per volume is 1 to 3 orders of magnitude larger in our experiment than for electron irradiation. Secondly, the total maximum D-enrichment (at the plateau) is about 20% larger in our study. X-ray irradiation acts as a bulk process due to the larger penetration depth (compared to electron irradiation) and volume diffusion of secondary electrons, enhancing D-enrichment.

DISCUSSION

The evolution of tholins induced by X-ray irradiation can be chemically quantified by the sequential depletion and growth of infrared bands, structurally by the analysis of Raman bands, and isotopically by the evolution of the SIMS profiles. For each of these processes, the evolution of the main parameters as a function of the dose can be modeled using first-order rate equations (Laurent et al. 2015),

dy dE = Ae -kE dose , ( 10
)
where y is the chemical, structural, or isotopic parameter, A is a constant, k is the rate constant and E dose is the irradiated dose per volume. We normalized all integrated parameters, such as vibrational band in-tensities and isotopic abundances, to the starting value (non-irradiated tholin) when possible. The parameters for these fits are listed in Table 2. We note that the isotopic rate constant is closely matched by the infrared depletion and the amorphization (Raman broadening) rate constants, pointing at coupled structural and chemical effects of X-ray irradiation, which itself drives the isotopic fractionation, i.e. dehydrogenation of H over D bonds is accompanied by a change in the organic structure (creation of defects).

Table 2. Fitted infrared, Raman and isotopic evolution parameters as an exponential function of X-ray dose While thermal degradation (Bonnet et al. 2015) and X-ray irradiation of tholins, both induce a similar loss of amines and alkyls (i.e. dehydrogenation of the tholin), there are important differences in their respective evolutions. Raman analysis shows that thermal degradation leads to an increase of polyaromatic units (Fig. 10), while X-irradiation does not affect the average size of polyaromatic units (even if infrared analysis shows that isolated aromatic units may be created). In both cases, defect formation competes with aromatization. For soft X-ray irradiation, the former is the dominant process.

y 0 A k [ eV -1 cm 3 ] C-
Disks can be further irradiated by cosmic rays (protons and ions), although their effect may be modulated by winds and magnetic fields in T-Tauri systems (Cleeves et al. 2013). The effect of cosmic rays on organic matter has been studied via laboratory experiments, using low energy or swift heavy ions sources, whose energy can be deposited in solids through nuclear or electronic interactions respectively. Brunetto et al. (2009) examined the effects of 200 -400 keV ion irradiation on soots showing the destruction of large polyaromatic units. On the other hand, Costantini et al. (2005) reported the carbonization of Kapton films and growth of polyaromatic structures induced by swift heavy ion irradiation. Experiments by Sabri et al. (2015) similarly showed the graphitization of water-covered amorphous carbon induced by 200 keV protons. Together, these experiments point at the strong influence of the starting material structure on its evolution under ionizing radiation.

Recent disk models by Cleeves et al. ( 2016) including X-rays and cosmic ray ionization investigated the deuterium enrichment of organic molecules. They argued that since D/H ratios corresponding to the most D-enriched organic materials in the solar system were not attained in the disk, there is a need of some interstellar inheritance. However, these models did not include macromolecular organics. Our experiments show a net loss of heteroatoms from the organic matrix, implying that a fraction of organic solids could be transferred into volatiles (hydrocarbons, radicals, ions). We thus propose that radiolysis and fragmentation of organics can eventually transfer the high D/H ratios of enriched solids back into the gas phase. In addition, as seen by the differences from the irradiation tracks for two X-ray energies (0.5 and 1.3 keV) in Fig. 12, extended studies with tholins of varying structure, using both broadband and hard X-rays sources, will be necessary to fully understand the effect of X-rays on organic matter. Stellar X-ray flux density is an important ionization source in protoplanetary disks, impacting its photochemical evolution. In this experimental study, we have shown that soft X-rays can induce a large deuterium enrichment in protoplanetary refractory organics. We compared the non-thermal effects of X-rays to previous experiments of electron irradiation on polymers (Laurent et al. 2015). X-ray irradiation of tholins induces higher D-enrichment than electron irradiation of organic polymers. This may reflect a better efficiency of X-rays to induce H-D fractionation, or a tendency of tholins to accomodate more defects hence record higher D-enrichment. The efficiency of irradiation will also depend on the size of the irradiated grains (Roskosz et al. 2016). Our experiments on µm thick tholins, comparable to the average size of observed dust grains in protoplanetary disks (from µm to cm sizes) (Natta et al. 2007), demonstrate that photon-induced isotopic enrichment of organic solids in disks is possible.

ASTROPHYSICAL IMPLICATIONS

As shown in Fig. 15, protoplanetary disks undergo irradiation from T-Tauri stellar fields having strong X-ray and UV components, although the higher energy X-ray photons are able to penetrate more effectively toward the disk mid-plane beyond 10 AU (Nomura et al. 2007;Walsh et al. 2012). The X-ray spectra of the classical T Tauri star, TW Hydrae, provides an integrated luminosity from 0.2 to 2 keV of L X = 2 ×10 30 erg s -1 (Kastner et al. 2002). At 10 AU, this amounts to a total exposure of 10 26 eV cm -2 in 1 Myr assuming an X-ray optical depth of τ X = 1. Using synchrotron soft X-rays, we have reached total fluences between 10 22 -10 24 eV cm -2 , well within the lifetimes of protoplanetary disks. The experimental plateau for D-enrichment is achieved at fluences of ∼10 24 eV cm -2 providing an upper limit to D-enrichment of protoplanetary organic matter by soft X-rays.

Compared to UV photons and electron irradiation, soft X-rays and their secondary photoelectrons can strongly modify the bulk of µm-sized organic particles. In spite of its non-thermal nature, soft X-rays can affect the chemical, structural and isotopic properties of organic macromolecules and solids via photochemical and radiolytical reactions, impacting the cycle of deuterium from dust to volatiles in cold regions of protoplanetary disks.

Figure 1 .

 1 Figure 1. Ellipsometric parameters Ψ and ∆ measured in the UV-Visible (390 -1000 nm) and modeled to determine the thickness of (a) a thin tholin film, and (b) a thick tholin film, used for X-ray absorption spectroscopy and irradiation.

Figure 2 .

 2 Figure 2. Imaginary index of refraction as a function of energy and calculated penetration length, LP , in µm. The energies selected for irradiation are marked in green.

Figure 3 .

 3 Figure 3. Top: 2D and 3D topographies measured via optical interferometry on the 0.5 keV irradiated tholin showing crater-like features. Bottom: 2D image of irradiated spots following X-ray doses at 0.5 keV.

Figure 4 .

 4 Figure4. Evolution of infrared spectra of irradiated tholins at increasing X-ray doses at 0.5 keV. The peak at 1100 cm -1 is attributed to the Si substrate.

Figure 5 .

 5 Figure 5. Infrared spectra of irradiated tholins at 0.5 keV split into main spectral regions as defined in Section 3.1.

Figure 6 .

 6 Figure 6. Infrared band evolution of the nitrile bands for X-ray doses at 0.5 keV. Deconvolution of the main nitrile modes for (a) the initial tholin, (b) the 10 hour irradiated tholin, showing the loss of the R-N≡C band intensity for the growth of the Ar-C≡N band intensity

Figure 8 .

 8 Figure 8. From top to bottom: Comparison between infrared spectra of an initial tholin to an X-irradiated tholin, thermally annealed tholin (Bonnet et al. 2015), and a UCAMM (Dartois et al. 2013).

Figure 9 .

 9 Figure 9. Raman spectra of X-ray irradiated tholins as a function of X-ray dose at 0.5 keV. The strong fluorescence from the non-irradiated tholin hinders measurement using a 514 nm excitation source.

Figure 10 .

 10 Figure10. Raman parameters of irradiated tholins as a function of the X-ray dose at 0.5 keV compared to the evolution of thermally annealed (300 -1000 • C) tholins(Bonnet et al. 2015). (a) The size of polyaromatic units (∝ ID/IG) is slightly affected by X-ray irradiation, while strongly affected by thermal heating. (b) Defect creation leads to amorphization in X-ray irradiated tholins, while aromatization dominates over defect creation in thermally degraded tholins.

Figure 11 .

 11 Figure 11. X-ray (1.3 keV) irradiated tholin following SIMS analysis of irradiated regions (brighter spots indicated by arrows) and the reference (non-irradiated) regions. Sputtering stops as soon as the Si substrate is reached (blue spots).

  Figure 12. D/H and H/C abundances normalized to the initial tholin (α D/H or α H/C ) as a function of sputtering depth.For 1.3 keV the α H/C profile after 10 minutes irradiation was discarded due to large systematic errors assigned to localized hydrogenated inhomogeneities on the tholin surface (fractures/wedges), as seen in Fig.3.

Figure 13 .

 13 Figure 13. Evolution of average normalized D/H vs H/C ratios for X-ray irradiated tholins at 0.5 keV and 1.3 keV.

Figure 14 .

 14 Figure14. Evolution of the normalized D/H abundances as a function of irradiation doses, for 30 keV electrons on organic polymers: polyethylene (PE), polystyrene (PS) and polyethylene terephthalate (PET)(Laurent et al. 2015) and for 0.5 keV and 1.3 keV X-rays on tholins (this work).

  (±0.05) 0.63(±0.04) 5.0(±1.8)×10 -27 R-N≡C 0.27(±0.04) 0.73(±0.06) 5.0(±1.6)×10 -27 C-defect 1.97(±0.20) -0.77(±0.05) 3.4(±3.0)×10 -27 I D /I G 0.87(±0.08) 0.02(±0.01) 1.6(±1.3)×10 -28 Γ-D 1.07(±0.08) -0.10(±0.01) 1.3(±0.2)×10 -28 Γ-G 1.07(±0.80) -0.09(±0.01) 3.7(±0.2)×10 -28 α D/H,0.5keV 1.86(±0.14) -0.90(±0.02) 8.8(±1.0)×10 -28 α D/H,1.3keV 1.71(±0.21) -0.75(±0.01) 8.7(±0.9)×10 -28 α H/C,0.5keV 0.49(±0.04) 0.52(±0.01) 7.6(±0.7)×10 -28 α H/C,1.3keV 0.66(±0.07) 0.53(±0.02) 7.2(±0.6)×10 -28

Figure 15 .

 15 Figure15. Schema of the thermal and non-thermal processes on an XUV illuminated protoplanetary disk. UV photons and X-rays irradiate the disk surface (ionized and molecular layers), but only X-rays reach the mid-plane (regions of ice formation and planetesimal growth) and outer edges. Ionization by cosmic rays (CRs) is represented in green.

Table 3 .

 3 Normalized integrated infrared modes (αIR) as a function of X-ray dose at 0.5 keV Integrated infrared mode [ cm -1 ] time [ s ] 0.5 keV dose [ eV/cm 3 ]

	C-H	N-H	R-N≡C	C-defect

Table 4 .

 4 Raman parameters as a function of X-ray dose at 0.5 keV. The FWHM of the Raman D and G bands are expressed as Γ-D and Γ-G. time [ s ] 0.5 keV dose [ eV/cm 3] ωD [ cm -1 ] ωG [ cm -1 ] Γ-D [ cm -1 ] Γ-G [ cm -1 ]

	ID/IG

Table 5 .

 5 SIMS elemental abundances as a function of X-ray dose at 0.5 keV and 1.3 keV

			0.5 keV			1.3 keV
	time [ s ] dose [ eV/cm 3 ]	α D/H	α H/C	dose [ eV/cm 3 ]	α D/H	α H/C
	600	3.80 ×10 26	1.17 (±0.12)	0.89 (±0.12)	1.93 ×10 26	1.05 (±0.07)	NA *
	1800	1.14 ×10 27	1.51 (±0.10)	0.74 (±0.05)	5.78 ×10 26	1.17 (±0.10) 1.07 (±0.13)
	3600	2.28 ×10 27	1.83 (±0.06)	0.53 (±0.03)	1.16 ×10 27	1.52 (±0.05) 0.76 (±0.04)
	14400	9.12 ×10 27	1.86 (±0.08)	0.47 (±0.03)	4.62 ×10 27	1.66 (±0.05) 0.69 (±0.02)
	36000	2.28 ×10 28	1.83 (±0.07)	0.53 (±0.02)	1.16 ×10 28	1.72 (±0.03) 0.67 (±0.02)

Acknowledgments. X-ray experiments were performed at the SEXTANTS beam line of the SOLEIL synchrotron (project No. 20150772). We thank Ferenc Borondics for access to the IR and Raman spectromicroscopes of the SMIS beam line and François Nicolas for his help with diagnostics at the Surface Laboratory. N.C. and L.G. thank the European Research Council for funding via the ERC PrimChem project (grant agreement No. 636829.)