Joël Goossens
email: joel.goossens@ulb.be

Goossens

Masson

Simulation intervals for uniprocessor real-time schedulers with preemption delay

Keywords: Computer systems organization → Real-time operating systems Scheduling theory, uniprocessor, preemption delay

In the framework of embedded and time-critical systems we consider the scheduling of preemptive real-time periodic tasks upon uniprocessor. We consider the notion of simulation interval, a finite time interval such that the schedule starts to repeat in a cycle. The interest is to design a time interval which includes all possible (reachable) schedule states. Our study focuses on a model where preemption costs are explicitly considered, i.e., the time required by the real-time operating system (rtos) or the hardware to load the context of execution of preempted jobs. We present and prove correct a simulation interval for asynchronous arbitrary deadline tasks with preemption delays and that holds for any deterministic and memoryless scheduler upon uniprocessor. This first contribution is valid for all schedulers one can imagine (but deterministic and memoryless), including non-preemptive schedulers, not work-conserving, not necessarily popular "real-time" ones. We then consider a particular case, regarding the preemption delays, for edf scheduling for which we extend the work of Leung and Merril [25] showing that [0, 𝑂 max + 2 • 𝐻 ) is a simulation interval (significantly shorter than the general result). We also show that [0, 𝑆 𝑛 + 𝐻 ) from [17] remains a simulation interval for fixed task priority (ftp) schedulers. Before concluding we open a discussion on the scope of the results, paths for reducing pessimism as well as other potential results opened by this work.

INTRODUCTION

In this research we study the scheduling of asynchronous arbitrary deadline preemptive hard real-time periodic task sets upon uniprocessor systems. We consider preemptive tasks/jobs in the sense that a higher priority task/job can preempt the current task/job in order to execute the highest priority one. The originality of this work, regarding the state-of-the-art, is the fact that the task model considers preemption delays, i.e., the duration required by the real-time operating system (rtos) or the hardware to resume a task/job after a preemption.

We study the notion of simulation interval with the following definition: Definition 1 (Simulation interval). A finite and safe time interval [0, 𝑏) such that the schedule start to repeat in a cycle.

In this work, we consider the notion of simulation interval for feasible schedule. Note that infeasible schedules are not necessarily periodic, but this question is beyond the scope of this work.

A simulation interval is a window of time during which all reachable states of the system have been observed. A schedule consists of two parts, a first one called the transient phase, the time interval [0, 𝑋 ), during which the succession of states appears only once, and a second one called the steady phase, the time interval [𝑋, 𝑋 + 𝑘𝐻 ), during which the succession of states is periodic. Most of the time 1the period of the schedule is a multiple of the hyperperiod (see Section 3 for a formal definition). This is illustrated in Figure 1, with 𝑋 and 𝑘 as unknown constants and 𝐻 as the hyperperiod. We are interested in finding exact values for 𝑋 and 𝑘 or at least an upper-bound for 𝑋 + 𝑘 • 𝐻 .

The notion of simulation interval should not be confused with the notion of feasibility interval: Definition 2 (feasibility interval). A finite interval [𝑎, 𝑏) such that if all the deadlines of jobs released in the interval are met, then the system is schedulable.

Informally speaking, the feasibility interval includes the worstcase schedule states regarding schedulability like, e.g., the first busy period for synchronous arbitrary deadlines and fixed task priority (ftp) schedulers [START_REF] Lehoczky | Fixed Priority Scheduling of Periodic Task Sets with Arbitrary Deadlines[END_REF] while the simulation interval includes all possible (reachable) schedule states.

We will also use the notion of exact simulation interval:

Definition 3 (Exact simulation interval). An exact simulation interval is a simulation interval [0, 𝑏) such that there is no simulation interval [0, 𝑐) with 𝑐 < 𝑏.

Organization. The remainder of this paper is structured in the following manner. Section 2 presents a motivations of our work. Section 3 presents our task model and assumptions. Section 4 describes related work. Section 5 presents our contributions. Section 6 presents discussions on the scope of the results, paths for reducing pessimism as well as other potential results opened by this work. Lastly, Section 7 concludes our work.

MOTIVATIONS

We motivate here the two axis of our work: first why we are interested in the simulation intervals, then why we consider preemption delays in our model.

Simulation Intervals

The need to have a simulation interval is justified each time we need to study the significant part of a schedule and thus avoid the anecdotal (non representative) part. This could be for performance or benchmark studies. Design tools (such as Cheddar [START_REF] Singhoff | About Real Time Scheduling Analysis of Ada Applications[END_REF][START_REF] Singhoff | Cheddar: a flexible real time scheduling framework[END_REF] or Tina [START_REF] Berthomieu | The tool TINAconstruction of abstract state spaces for Petri nets and time Petri nets[END_REF]) often require to provide the simulation interval as input. When the worst-case scenario is known (or sometimes several cases to cover the worst-case, with a brute-force approaches see e.g. [START_REF] Hladik | A brute-force schedulability analysis for formal model under logical execution time assumption[END_REF]) the notion of simulation interval allows to have an exact schedulability test. For example Nasri et al. [START_REF] Nasri | FIFO with offsets: High schedulability with low overheads[END_REF] use it to have an exact simulation-based schedulability test for the fifo scheduler and tasks with offsets. Unfortunately, regarding the tasks duration, the worst-case scenario does not necessarily corresponds to consider the worst-case execution time (wcet -see Section 3), i.e., schedulers are not necessarily 𝐶-sustainable [START_REF] Burns | Sustainability in real-time scheduling[END_REF]. We will discuss this aspect before our conclusion (see Section 6). In offline approaches such as in the work of Xu and Parnas [START_REF] Xu | Priority scheduling versus pre-run-time scheduling[END_REF], it is fundamental to know a simulation interval in order to build and evaluate the whole schedule (i.e., until the end of the steady state). Regarding our scheduling problem, asynchronous periodic tasks with preemption delays, unfortunately the synchronous case is not the worst-case scenario. This fact motivates specifically the study of asynchronous scenarios.

Preemption Delays

A common assumption in scheduling theory, implicitly or not, since the seminal work of Liu and Layland [START_REF] Laung | Scheduling algorithms for multiprogramming in a hard-real-time environment[END_REF], is that we neglect the time needed for preemption, i.e., the time needed for the system (bare metal or operating system) to switch from one task to another. More generally, all disturbances related to the hardware or the operating system, which could be called scheduling costs, are mostly neglected when studying the schedulability of real-time systems. Precisely, it is often implicitly assumed that these costs are an integral part of the worst-case execution time (wcet) and therefore indirectly included in the analysis [START_REF] Laung | Scheduling algorithms for multiprogramming in a hard-real-time environment[END_REF]. However, in the case of preemption, the cumulative time spent in operating system routines or interrupt handlers depends on the number of preemptions and therefore on the scheduling. One could argue that the scheduling costs can easily be integrated in the wcet with the argument that the price is paid by the task that initiates the preemption and not by the task that suffers it (the price is the cost of the context switch to execute the new task plus the one to resume the preempted task). Then, since with fixed task priority (ftp) schedulers or even fixed job priority schedulers such as edf, a task can only preempt one task per activation (when it is activated), so the cost is limited and can be integrated into the wcet.

A more detailed study of the mechanisms at the origin of what can be called the scheduling cost goes against this argument.

We can distinguish several components to the scheduling cost, for example:

• 𝐶 in the cost of adding a job to the queue, • 𝐶 start the cost of starting a job, • 𝐶 restart the cost of restarting a job after a preemption, • 𝐶 out the cost of terminating a job. Depending on the system, most of these components are constant or can be bounded and occur only once per job. For example, if the system is handling the arrival of tasks with interrupts and has sufficient priority levels for interrupts, 𝐶 in is only paid when the job effectively starts, so it does not interfere with the other jobs (this is called integrated management). 𝐶 start and 𝐶 out depend on the jobs, they correspond respectively to the loading of the execution context of the job and to its release. Depending on the implementation, they can cover different types of overhead for the system, but it is reasonable to think that they can be integrated into the wcet. 𝐶 restart is more problematic, as it depends on each task, and is paid after each preemption. Because of that, it is not necessarily equivalent, in terms of duration, for a task 𝜏 𝑖 to preempt a task 𝜏 𝑗 or another task 𝜏 𝑘 . We so think it is important to add 𝐶 restart to the model and the analysis.

Another issue with these scheduling costs is that they correspond, or may correspond, to the execution of operating system routines, and therefore to non-preemptive or at least "atomic" sections of code, i.e., they can be interrupted, but they will then have to be re-executed in their totality.

The model considered in this paper is a first step towards a model that fully integrates the operating system (or the hardware) interference during task preemption. In order to simplify the problem we have arbitrarily chosen to consider only 𝐶 restart , which we note 𝛼 𝑖 for the task 𝜏 𝑖 (see Section 3). We also chose to consider its execution completely non-preemptive, i.e., once the system starts to reload the context of a previously preempted task 𝜏 𝑖 , no preemption can take place for a duration equal to 𝛼 𝑖 . In future work, we may relax this constraint partially to allow preemption on the condition that we must re-execute the reloading section completely later. We will also have to consider the same constraints for the other components of the scheduling cost.

MODEL AND ASSUMPTIONS

In this work we consider the scheduling of asynchronous periodic task set with preemption delays upon uniprocessor. More precisely, we assume there is a delay

𝛼 𝑖 ∈ {N ∪ 0}

RELATED WORK

We will start by reviewing the works on simulation intervals, and then we will look at the works considering the overheads related to the scheduler and in particular to preemptions. 2 where means "equals by definition".

Simulation intervals

In the uniprocessor context, assuming an utilisation factor less than or equal to 1 and for systems of independent tasks with constrained deadlines and scheduled with edf, the pioneering work [START_REF] Joseph | A note on preemptive scheduling of periodic, real-time tasks[END_REF] shows that the interval [0, 𝑂 max + 2𝐻 ) is a simulation interval. Moreover, the transient phase is included in the window [0, 𝑂 max + 𝐻 ) and the steady phase is contained in the scheduling produced in the interval [𝑂 max + 𝐻, 𝑂 max + 2𝐻 ).

For fixed task priority (ftp) schedulers, it is shown in [START_REF] Leung | On the complexity of fixedpriority scheduling of periodic, real-time tasks[END_REF] that the same interval applies. The result is extended to arbitrary deadlines in [START_REF] Goossens | Feasibility intervals for the deadline driven scheduler with arbitrary deadlines[END_REF]. A better bound is presented in [START_REF] Goossens | The Non-Optimality of the Monotonic Priority Assignments for Hard Real-Time Offset Free Systems[END_REF] and is [0, 𝑆 𝑛 + 𝐻 ) where 𝑆 𝑛 is given by the recursive equation

𝑆 1 𝑂 1 , 𝑆 𝑖 max(𝑂 𝑖 , 𝑂 𝑖 + ⌈ (𝑆 𝑖-1 -𝑂 𝑖 ) 0 𝑇 𝑖 ⌉𝑇 𝑖 )
. A more general result, which holds for resource-sharing tasks with precedence constraints, constrained deadlines, and assumes only a work-conserving scheduler, was first shown in [START_REF] Choquet | Minimal Schedulability Interval for Real Time Systems of Periodic Tasks with Offsets[END_REF] then extended to parallel tasks in [START_REF] Bado | A semi-partitioned approach for parallel real-time scheduling[END_REF]. In that work, it is shown that an exact simulation interval (see Definition 3) is [0, 𝜃 𝑐 + 𝐻 ) where 𝜃 𝑐 is the date of the last idle time of the transient phase. Moreover, it is shown that 𝜃 𝑐 belongs to the interval [0, 𝑂 max + 𝐻 ). It follows that the steady phase corresponds to the first time window of length 𝐻 containing exactly 𝐻 (1-𝑈 ) idle times in the interval [0, 𝑂 max +2𝐻 ).

Several works have also addressed the problem for uniform multiprocessor platforms with global scheduling. In [START_REF] Cucu | Feasibility Intervals for Fixed-Priority Real-Time Scheduling on Uniform Multiprocessors[END_REF], it is shown that for constrained deadlines, in the synchronous case the simulation interval is bounded by [0, 𝐻 ), while in the asynchronous case the interval [0, 𝑆 𝑛 + 𝐻 ) given in [START_REF] Goossens | The Non-Optimality of the Monotonic Priority Assignments for Hard Real-Time Offset Free Systems[END_REF] still held. In [START_REF] Cucu | Feasibility Intervals for Multiprocessor Fixed-Priority Scheduling of Arbitrary Deadline Periodic Systems[END_REF], the authors modify the definition of 𝑆 𝑛 to give a simulation interval for arbitrary deadlines. The bound becomes [0, Ŝ𝑛 + 𝐻 ) where Ŝ𝑛 is given by the

recursive equation Ŝ1 𝑂 1 , Ŝ𝑖 max(𝑂 𝑖 , 𝑂 𝑖 + ⌈ ( Ŝ𝑖-1 -𝑂 𝑖 ) 0 𝑇 𝑖 ⌉𝑇 𝑖 ) + 𝐻 𝑖 .
The result is extended to heterogeneous (unrelated) multiprocessor systems in [START_REF] Cucu-Grosjean | Exact Schedulability Tests for Real-Time Scheduling of Periodic Tasks on Unrelated Multiprocessor Platforms[END_REF]. A result is given in the case of offline schedulers for systems with precedence constraints in [START_REF] Baro | Off-Line (Optimal) Multiprocessor Scheduling of Dependent Periodic Tasks[END_REF]. It is generalized to edf and to independent tasks in [START_REF] Nélis | Feasibility intervals for homogeneous multicores, asynchronous periodic tasks, and FJP schedulers[END_REF]: [0, 𝑂 max + 𝐻 𝑛 𝑖=1 (𝐶 𝑖 + 1)). Finally, a very general result for any deterministic and memoryless scheduler is given in [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF] and is [0, 𝐻 𝑛 𝑖=1 ((𝑂 𝑖 +𝐷 𝑖 -𝑇 𝑖 ) 0 + 1)). This interval works for task systems with multiple dependencies, defined as structural constraints (mutual exclusions, precedence constraints, self-suspension, non-preemptive tasks, etc.).

In a recent paper [START_REF] Lagha | Toward an Exact Simulation Interval for Multiprocessor Real-Time Systems Validation[END_REF], the authors propose a method to reduce the size of the simulation interval. They provide an exact simulation interval. However, the computation of this interval is heavily time consuming, involving a factorial time complexity computation of a set of points defined by an exponential number of linear constraints. The authors show that, despite its intractable time consuming computation, the size of the exact interval grows much slower than the one from [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF] with the number of tasks.

Taking into account the cost of preemptions

Many works have proposed ways to integrate scheduling delays or even other overheads related to the hardware platform or the operating system. Different approaches exist. Some of them consist in bounding this interference in order to add it to the schedulability analysis as part as the wcet or as a blocking term, inputting cost overheads either to the preempted tasks, the preempting ones or 2022-03-08 15:18. Page 3 of 1-10.

both. Some other approaches propose different task models with associated analyses. Finally some works adapt the scheduler algorithm to make it aware of the existence of such delays.

It is interesting to note that a part of the scheduling delays are external to the chosen scheduling policy, i.e., the overheads are imposed by the hardware platform or the operating system. They must nevertheless be taken into consideration in the schedulability analysis or more generally in the scheduler performance evaluation. As we point it out in this paper, they have an impact on the simulation interval and the periodicity of the schedule.

The first works to have been interested in the impact of the scheduler implementation and the kernel on scheduling are those of [START_REF] Daniel | Engineering and analysis of fixed priority schedulers[END_REF]. In that paper, the authors describe four distinct generic methodologies for implementing a scheduler: Integrated Interrupt Event-Driven Scheduling, Nonintegrated Interrupt Event-Driven Scheduling, Timer-Driven Scheduling and Timer-Driven Scheduling with counter. The first one assumes that tasks are activated by external hardware interrupts. The second one is a special case where there is no prioritised interrupts on the platform. The third one corresponds to what is also referred to in the literature as tick scheduling where the scheduler is woken up periodically every quantum of time and processes the events that have occurred since its last wake-up. In addition to the overhead experienced at each tick, due to the granularity of the triggering timer, transient priority inversions can also occur. The last one is a refinement of the previous one where instead of a periodic timer, an updated counter is used allowing a wake-up only when scheduling events occur (next activation, deadline, etc.). For each of them, they detail the cost components inherent to the scheduler, and extend the sufficient popular schedulability condition for ftp and constrained deadlines (∀ 𝑖 , ∃𝑡 ≤ 𝐷 𝑖 s.t. 𝑖 𝑗=1 𝑡 𝑇 𝑗 𝐶 𝑗 + 𝐵 𝑖 ≤ 𝑡) to take them into account either in the 𝐶 𝑗 term or in the 𝐵 𝑖 one. The extended conditions are first given under the assumption of a preemptive kernel, then a non-preemptive kernel, which adds an extra component to the blocking term 𝐵 𝑖 .

The works in [START_REF] Burns | Effective analysis for engineering real-time fixed priority schedulers[END_REF] point out the pessimism of the sufficient conditions developed in the preceding works and gives a finer analysis of the overhead of a scheduler implemented according to the tick scheduling paradigm.

Other works (see, e.g., [5,[START_REF] Burns | Preemptive Priority Based Scheduling: An Appropriate Engineering Approach[END_REF][START_REF] Wang | Scheduling fixed-priority tasks with preemption threshold[END_REF]) have considered the problem of limiting the number of preemptions or the one of computing a bound on the number of possible preemptions, permitting an integration of a maximum blocking time linked to the system interference in the schedulability analysis.

Another source of disturbance related to preemptions analysed in the literature is related to the use of caches. While caches permit a reduction of response times, they also introduce a significant source of variability in wcets due to the possibility for the task to suffer additional delays due to cache misses. The runtime overhead associated with cache misses caused by an early preemption is called Cache Related Preemption Delay (crpd) [START_REF] Busquets-Mataix | Adding instruction cache effect to schedulability analysis of preemptive real-time systems[END_REF]. The strategy for taking crpds into account is to bound them for inclusion in the feasibility analysis of ftp schedulers. To calculate a bound, one can use Evicting Cache Blocks (ecbs, the blocks used by the preempting task and which can replace those of a preempted task [START_REF] Busquets-Mataix | Adding instruction cache effect to schedulability analysis of preemptive real-time systems[END_REF]) or Useful Cache Blocks (ucbs, the blocks used by a preempted task which are reused later, and therefore potentially deleted during a preemption [START_REF] Lee | Analysis of cache-related preemption delay in fixed-priority preemptive scheduling[END_REF]). Some works combine the two approaches [START_REF] Altmeyer | Improved cache related pre-emption delay aware response time analysis for fixed priority preemptive systems[END_REF][START_REF] Staschulat | Scalable precision cache analysis for preemptive scheduling[END_REF][START_REF] Tomiyama | Program path analysis to bound cache-related preemption delay in preemptive real-time systems[END_REF]. The case of edf is studied in [START_REF] Lunniss | Integrating cache related pre-emption delay analysis into edf scheduling[END_REF] by extending the feasibility test based on the Demand Bound Function (dbf). It is shown in [START_REF] Phavorin | Scheduling with preemption delays: anomalies and issues[END_REF] that optimal scheduling with crpds taken into account is an NPhard problem. More recently, a feasibility interval is given for this problem in [START_REF] Nam Tran | Feasibility interval and sustainable scheduling simulation with CRPD on uniprocessor platform[END_REF].

What distinguishes our study from these works is that we consider the overhead of preemption not as a special case where the wcet is increased, but by adding a non-preemptive code block after each preemption.

In his PhD thesis and in [START_REF] Meumeu | Extending rate monotonic analysis with exact cost of preemptions for hard real-time systems[END_REF], Meumeu proposes an approach to include the cost of reloading a task after a preemption in the feasibility analysis of a system scheduled with ftp. In his model, this reloading cost is different for each task, and consists in an "atomic" sequence of operations (it can be preempted but then all the sequence must be re-executed later). He shows that with this model, the critical instant cannot be characterised, that the worst-case response time of a task can occur in the transient phase, and that the Optimal Priorities Assignment algorithm (opa) compatibility conditions [START_REF] Audsley | Optimal Priority Assignment and Feasibility of Static Priority Tasks With Arbitrary Start Times[END_REF] are not respected. He proposes an optimal priority assignment algorithm.

In his PhD Thesis, Bimbard [START_REF] Bimbard | Dimensionnement temporel de systèmes embarqués : application à OSEK[END_REF] considers the kernel overheads of an osek kernel (fixed priority scheduling), taking into account activation, termination and preemption delays. For a given task, activation and termination are extra duration added to the wcet of a task while preemption delays are taken into account at each activation of higher priority task. A worst case response time analysis is proposed with a bound on the number of preemptions that makes the synchronous scenario still the worst case scenario at the price of some pessimism. This Research. In this paper (i) we model preemption delays as task-dependant non preemptive additional code sections executed for a task when it is resumed after a preemption, (ii) we present and prove correct a simulation interval for asynchronous arbitrary deadline tasks with such preemption delays and that holds for any deterministic and memoryless scheduler upon uniprocessor, (iii) for particular case, regarding the preemption delays, and edf scheduling we extend the work of Leung and Merril [START_REF] Joseph | A note on preemptive scheduling of periodic, real-time tasks[END_REF] showing that [0, 𝑂 max + 2 • 𝐻 ) is a simulation interval (significantly shorter than the general result). Lastly, (iv) for ftp schedulers we also extend the work [START_REF] Goossens | The Non-Optimality of the Monotonic Priority Assignments for Hard Real-Time Offset Free Systems[END_REF] and we show show that [0, 𝑆 𝑛 + 𝐻 ) remains a simulation interval.

CONTRIBUTIONS

A first contribution of this paper is to provide a simulation interval (see Definition 1) for the scheduling of periodic tasks with preemption delays. We build upon the work of [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF] where the authors provided a simulation interval and more interestingly a technique to do so.

The former work concerned multiprocessor platforms and dependent tasks. Here the novelty is to consider preemption delays for uniprocessor platform. In Section 5.1 we will introduce the technique of [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF], in Section 5.2 we present and prove correct a simulation interval for tasks with preemption delays and that holds for any deterministic and memoryless scheduler (see Definition 5).

Summarising the technique of [19]

The authors consider multiprocessor scheduling of preemptive periodic real-time tasks without preemption delays. They first formalise the notions of system state, system pre-state and deterministic and memoryless scheduler with Definitions 4 and 5.

Definition 4 (State and pre-state of a system [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF]). The state of a system of 𝑛 tasks can be defined as a (2𝑛)-tuple ). We will see that we need to consider arbitrary deadlines even if the original system under study have constrained deadlines (see Lemma 7 for details) consequently we assume to schedule the various jobs of the same task in the fifo order (the oldest job first).

𝑆
The pre-state of a system of 𝑛 tasks can be defined as a (2𝑛)-tuple Ŝ Ĉrem 1 , . . . , Ĉrem 𝑛 , Ω 1 , . . . , Ω 𝑛 , where:

• Ω 𝑖 is the same local clock as in the state 𝑆 of the system,

• Ĉrem 𝑖 is the remaining work to process for 𝜏 𝑖 not taking the releases at the considered instant into account. We can formalize the remaining work in state and pre-state, for any 𝑡 ≥ 𝑂 𝑖 as follows:

Ĉrem 𝑖 (𝑡) 0, ∀𝑡 ≤ 𝑂 𝑖 𝐶 rem 𝑖 (𝑡) Ĉrem 𝑖 (𝑡) + 𝐶 𝑖 if Ω 𝑖 = 0 Ĉrem 𝑖 (𝑡) otherwise Ĉrem 𝑖 (𝑡 + 1) 𝐶 rem 𝑖 (𝑡) -1 if 𝜏 𝑖 executed on [𝑡, 𝑡 + 1)
𝐶 rem 𝑖 (𝑡) otherwise Definition 5 (Deterministic and memoryless scheduler [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF]).

A scheduler is deterministic and memoryless if, and only if, the scheduling decision 3 at time 𝑡 is unique and univocally defined by the state of the system at time 𝑡.

Consequently, the state of the system must contain all information required by the scheduler to take its decision but only the strictly required and minimum information in order to bound as much as possible the number of distinct states.

Then the authors show that considering only the synchronous scenario of a transformed and arbitrary deadline task set is sufficient to enumerate all feasible schedules of a system, thanks to Lemma 7 (and Definition 6). Definition 6 (Set of feasible schedules [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF]). We define the function F such that F (𝑆) is the set of all feasible schedules obtained by deterministic and memory-less schedulers for task system 𝑆.

Lemma 7 (Lemma 1 in [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF]). Let 𝑆 be a set of independent tasks with ∀𝑖 ∈ 1, . . . , 𝑛, 𝑂 𝑖 ≥ 0. We denote 𝑂 𝑖 the offset of the task 𝜏 𝑖 and 𝐷 𝑖 its relative deadline. Let 𝑆 ′ be the same system, except for the release dates given by 𝑂 ′ 𝑖 = 0 and the relative deadlines 𝐷 ′ 𝑖 = 𝐷 𝑖 + 𝑂 𝑖 .

3 By scheduling decision at time 𝑡 we mean what the scheduler/system decides to do during the slot [𝑡, 𝑡 + 1): idle the processor, execute a task, reload a preempted task.

The set of feasible schedules of 𝑆 is included in the set of feasible schedules of 𝑆 ′ , i.e., F (𝑆) ⊆ F (𝑆 ′ ).

Please notice that, while the original system 𝑆 can be composed of constrained deadlines, the synchronous system 𝑆 ′ could be composed of arbitrary deadlines.

Based on the notions of state and pre-state of a system, the authors show, in Lemma 8, that if two states or pre-states are identical, then the scheduling decision of any deterministic and memoryless schedulers are identical, consequently the schedule repeats. This allows the authors to bound the number of distinct (pre-)states in Lemma 9.

Lemma 8 (Lemma 2 in [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF]). For synchronous task systems, if two pre-states are identical, then the scheduling decision of a deterministic and memoryless scheduler is the same. Lemma 9 (Lemma 3 in [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF]). Any feasible schedule of a synchronous arbitrary deadline independent task system generated by a deterministic and memoryless scheduler reaches a cycle4 at or prior to

𝐻 • 𝑛 𝑖=1 ((𝐷 𝑖 -𝑇 𝑖 ) 0 + 1).
Lastly, the authors combine Lemma 7 and Lemma 9 to provide and prove correct a simulation interval for asynchronous systems with the Theorem 10.

Theorem 10 (Theorem 1 in [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF]). Any feasible schedule of an asynchronous independent tasks system generated by a deterministic and memoryless scheduler reaches a cycle 4 at or prior to

𝐻 • 𝑛 𝑖=1 ((𝑂 𝑖 + 𝐷 𝑖 -𝑇 𝑖 ) 0 + 1).
Recap. For the determination of simulation intervals, the technique of [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF], is based on the formalization of the notion of system state in a mininalist way which combined with a deterministic and memoryless scheduler leads to enumerate and bound the distinct possible (reachable) states.

A first contribution: a simulation interval for deterministic schedulers with preemption delay

In this section we extend the work of Goossens et al. in [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF] in the sense that for each task 𝜏 𝑖 we have an additional 𝛼 𝑖 parameter the reload delay. Following the same arguments as Goossens et al. in [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF] we design and proof correct the following simulation interval: Theorem 11 (Simulation interval). Any feasible schedule of asynchronous arbitrary deadline periodic tasks with preemption delays generated by a deterministic and memoryless scheduler upon uniprocessor platform with a reload delay of 𝛼 𝑖 for the task 𝜏 𝑖 reaches a cycle 4 at or prior to:

𝐻 • (𝑛 + 1) • (𝛼 max + 1) 𝑛 𝑖=1 ((𝑂 𝑖 + 𝐷 𝑖 -𝑇 𝑖 ) 0 + 1). (1)
In order to prove that property we need first to extend Definition 4, the system (pre-)state to consider preemption delay. We have to design it as minimal as possible to limit as much as possible the number of distinct (pre-)states. Our proposition is to consider two additional parameters:

• 𝑟 ∈ {0, 1, . . . ,
𝛼 ℓ if 𝑟 (𝑡 + 1) = ℓ > 0 ∧ 𝑟 (𝑡) ≠ ℓ ∧ (𝐶 rem ℓ (𝑡 + 1)%𝐶 ℓ ) ≠ 0 (3) 0 otherwise (4)
Equation 3 corresponds to the begin of an reload region such that 𝜏 ℓ was preempted in the past. Indeed we consider arbitrary deadlines consequently several jobs of the same task can be active simultaneously, 𝐶 rem ℓ (𝑡 + 1) represents the cumulative remaining computation time. Notice that if 𝐶 rem ℓ (𝑡 + 1) is not a multiple of 𝐶 ℓ a job of 𝜏 ℓ is resumed otherwise this is a new job of 𝜏 ℓ since at time 𝑡 a different task (𝑟 (𝑡) ≠ ℓ) was executed/reloaded5 . Lemma 13. Any feasible schedule of a synchronous arbitrary deadline task system generated by a deterministic and memoryless scheduler with preemption delay reaches a cycle 4 at or prior to

𝐻 (𝑛 + 1)(𝛼 max + 1) 𝑛 𝑖=1 (𝐷 𝑖 -𝑇 𝑖 ) 0 + 1 .
Proof. For parameters in common between Definition 4 and Definition 12 (𝐶 rem 𝑖 's and Ω 𝑖 's) from [START_REF] Goossens | Periodicity of real-time schedules for dependent periodic tasks on identical multiprocessor platforms[END_REF] we know that, regarding the states and pre-states, that in the worst case after 𝑛 𝑖=1 (𝐷 𝑖 -𝑇 𝑖 ) 0 + 1 hyperperiods the state or the pre-state repeats. We must complete the arguments considering the two extra parameters (𝑟 and rem-reload): we need that value 𝑟 repeats as well as the value rem-reload. Since 𝑟 is limited to 𝑛 + 1 values and rem-reload is limited to 𝛼 max + 1 values. We conclude that after a maximum of (𝑛 + 1) • (𝛼 max + 1) 𝑛 𝑖=1 (𝐷 𝑖 -𝑇 𝑖 ) 0 + 1 hyperperiods the schedule (if feasible) repeats. □

In order to have similar result for asynchronous systems we need to extend Lemma 7 for systems with reload delay. Lemma 14. Let 𝑆 be a set of independent tasks with ∀𝑖 ∈ 1, . . . , 𝑛, 𝑂 𝑖 ≥ 0. We denote 𝑂 𝑖 the offset of the task 𝜏 𝑖 and 𝐷 𝑖 its relative deadline. Let 𝑆 ′ be the same system, except for the release dates given by 𝑂 ′ 𝑖 = 0 and the relative deadlines 𝐷 ′ 𝑖 = 𝐷 𝑖 + 𝑂 𝑖 . The set of feasible schedules of 𝑆 is included in the set of feasible schedules of 𝑆 ′ .

Proof. A feasible schedule for 𝑆 is also a feasible schedule for the system 𝑆 ′ , since the absolute deadlines are identical, we just release jobs earlier which doesn't affect the validity of the solution. □ Theorem 11 follows from Lemma 14 and Lemma 13.

Shorter interval for edf and ftp schedulers in a particular case

Introduction.

In this section we will study the particular case of constrained deadlines and binary delays (0 ≤ 𝛼 𝑖 ≤ 1, see Definition 17). We will see that it is not straightforward to define shorter intervals -in comparison with [0, 𝐻 • (𝑛 +1) • (𝛼 max +1) 𝑛 𝑖=1 ((𝑂 𝑖 + 𝐷 𝑖 -𝑇 𝑖 ) 0 + 1)) -with preemption delays. We will first consider the case of edf and see why the Leung and Merril [START_REF] Joseph | A note on preemptive scheduling of periodic, real-time tasks[END_REF] interval, [0, 𝑂 max + 2𝐻 ), for edf scheduling is not a simulation (or feasibility) interval. We will exhibit the phenomenon which occurs because of preemption delays and which invalidates the Leung and Merril's properties. Lastly, we consider the case of ftp schedulers.

Let's first recap the work of Leung and Merril. In [START_REF] Joseph | A note on preemptive scheduling of periodic, real-time tasks[END_REF], the authors consider the edf scheduling of periodic constrained deadline task sets upon uniprocessor. They design and prove correct the [0, 𝑂 max + 2𝐻 ) simulation (and feasibility) interval. They first prove Lemma 1 required to prove Lemma 2 the main property: Lemma 15 (Lemma 1 in [START_REF] Joseph | A note on preemptive scheduling of periodic, real-time tasks[END_REF]). Let 𝑆 be schedule of a task system 𝑅 constructed by the edf algorithm. Then for each tasks 𝜏 𝑖 and for each time instant 𝑡 1 ≥ 𝑂 𝑖 , we have 𝑒 𝑖,𝑡 1 ≥ 𝑒 𝑖,𝑡 2 , where 𝑡 2 𝑡 1 + 𝐻 , and 𝑒 𝑖,𝑡 is the amount of time for which 𝜏 𝑖 has executed since its last request.

Lemma 16 (Lemma 2 in [START_REF] Joseph | A note on preemptive scheduling of periodic, real-time tasks[END_REF]). Let 𝑆 be schedule of a task system 𝑅 constructed by the edf algorithm. If 𝑅 is feasible, then 𝐶 𝑆 (𝑅, 𝑡 1 ) = 𝐶 𝑆 (𝑅, 𝑡 2 ), where 𝑡 1 = 𝑂 max + 𝐻 , 𝑡 2 = 𝑂 max + 2𝐻 and 𝐶 𝑆 (𝑅, 𝑡), to be the n-tuple (𝑒 1,𝑡 , . . . , 𝑒 𝑛,𝑡 ).

In the following, we will show that, unfortunately, Lemma 15 and consequently Lemma 16 are not satisfied with preemption delays.

Consider, for instance, the system described by Table 1. Figure 2 corresponds to the edf schedule, in this case the periodicity starts immediately (at time origin) but the period is 2 • 𝐻 = 24. In Figure 2, we represent a task request by a ↓ (down arrow); we represent the task deadline with a ⃝ (circle); white rectangles represent task execution, grey rectangles represent reload activities (⊲⊳ represents priority inversion a notion that will be introduced later).

This example contradicts Leung and Merril's property (Lemma 15) required to prove that schedule repeats from 𝑂 max +𝐻 with a period 𝐻 , since we have: 0 = 𝑒 4,7 < 𝑒 4,19 = 1.

Informally speaking, the progress of task 𝜏 4 at time instant 7 is strictly smaller than the progress of task 𝜏 4 an hyperperiod later (at time instant 19).

Consequently the system state at time instant 𝑂 max + 𝐻 = 18 is not the same an hyperperiod later (at time instant 30), in particular at time instant 30 task 𝜏 1 must complete its non-preemptive reload Priority inversion. We will now exhibit more precisely the phenomenon which is problematic: the priority inversion. On Figure 2, ⊲⊳ represents priority inversion at time instants 3, 6, 3 + 2𝐻 = 27, 6 + 2𝐻 = 30, 3 + 4𝐻, 6 + 4𝐻, . . . where an higher priority job is delayed by a lower priority one.

p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p

Binary delays. An interesting particular case is to consider binary delays because priority inversion is impossible (we will prove that in this section), and for which shorter intervals can be identified, which is a second contribution of this paper. Definition 17 (Binary delays). In this model the preemption delay are equal to zero or one: ∀𝑖 𝛼 𝑖 ∈ {0, 1}.

Proposition 1 (No priority inversion with binary delays).

For binary delays and discrete time systems priority inversion is impossible, i.e., a lower priority job/task cannot interfere in the execution of higher priority job/task.

Proof. If we resume a task/job during the time-slot [𝑡, 𝑡 + 1) we know there is no higher priority task/job at time 𝑡, moreover the resumption will be completed at time 𝑡 + 1 in order to schedule immediately an higher priority task/job (if any). □

For this particular case we can generalise the properties of Leung and Merril. To do this we will extend the definition of 𝐶 𝑆 (𝑅, 𝑡) (system state). Please notice that we are not interested here in enumerating all distinct states but we need to compare states distant by one hyper-period. Definition 18 ( C𝑆 (𝑅, 𝑡)). Let 𝑆 be a schedule of the task system 𝑅. The state of a system of 𝑛 tasks can be defined as C𝑆 (𝑅, 𝑡) 𝑒 1,𝑡 , . . . , 𝑒 𝑛,𝑡 , 𝑟, and rem-reload where at time instant 𝑡:

• 𝑒 𝑖,𝑡 is the amount of time for which 𝜏 𝑖 has executed since its last request, • 𝑟 ∈ {0, 1, . . . , 𝑛} is the task index of the running/reloaded task (if any) or 0 is the processor is idle, • rem-reload ∈ {0, 1} is the remaining reload time of 𝜏 𝑟 if 𝑟 ≠ 0 and rem-reload > 0. If 𝑟 = 0 rem-reload is undefined.

We have now the material to extend the Leung and Merril lemmata, However, the proofs are much more sophisticated and complex.

Lemma 19. Let 𝑆 be schedule of a task system 𝑅 constructed by the edf algorithm with binary delays. Then for each tasks 𝜏 𝑖 and for each time instant 𝑡 1 ≥ 𝑂 𝑖 , we have 𝑒 𝑖,𝑡 1 ≥ 𝑒 𝑖,𝑡 2 , where 𝑡 2 𝑡 1 + 𝐻 , and 𝑒 𝑖,𝑡 is the amount of time for which 𝜏 𝑖 has executed since its last request.

Proof. We will show the property by contradiction. We assume, for the purpose of the contraction, that ∃ ℓ 1 s.t. 𝑒 ℓ 1 ,𝑡 1 < 𝑒 ℓ 1 ,𝑡 2 . Then there must be some time instant 𝑡 ′ 1 < 𝑡 1 such that task 𝜏 ℓ 1 is active at both instants 𝑡 ′ 1 and 𝑡 ′ 2 = 𝑡 ′ 1 + 𝐻 , and 𝜏 ℓ 1 is executing at time instant 𝑡 ′ 2 but not at time instant 𝑡 ′ 1 . We will consider all cases where this can happen with binary delays. For each case ((a)-(d)) we will show that either the case is not possible or that it produces an infinite progression of requests that precede 𝑡 ′ 1 which is not possible by definition of task offsets and time origin. But first we introduce an additional notation: 𝑆 (𝑡) denotes the schedule decision at time 𝑡, 𝑆 (𝑡) = 𝐸 𝑖 meaning the execution of 𝜏 𝑖 , 𝑆 (𝑡) = 𝑅 𝑖 meaning the reloading of 𝜏 𝑖 , and 𝑆 (𝑡) = 𝐼 meaning that the processor is idle.

Case a. 𝑆 (𝑡 ′

1 ) = 𝐸 ℓ 2 , ℓ 2 ≠ ℓ 1 (and 𝑆 (𝑡 ′ 2 = 𝐸 ℓ 1 )), in this case task 𝜏 ℓ 2 , an edf higher priority task, is active at time instant 𝑡 ′ 1 but not at time instant 𝑡 ′ 2 . But this means that 𝑒 ℓ 2 ,𝑡 ′ 1 < 𝑒 ℓ 2 ,𝑡 ′ 2 , and we may repeat the above argument to produce an infinite progression of task computations 𝜏 ℓ 3 , 𝜏 ℓ 4 , . . . , for which no lower bound will exist for the time at which the task computations in the sequence are Binary delays and ftp schedulers. For fixed task priority (ftp) schedulers such as Rate Monotonic, Deadline Monotonic or opa a much shorter simulation interval than [0, 𝐻 • (𝑛 + 1) • (𝛼 max + 1) 𝑛 𝑖=1 ((𝑂 𝑖 + 𝐷 𝑖 -𝑇 𝑖 ) 0 + 1)) and even than [0, 𝑂 max + 2 • 𝐻 ) could be designed based on the work [START_REF] Goossens | The Non-Optimality of the Monotonic Priority Assignments for Hard Real-Time Offset Free Systems[END_REF] for binary delays.

The benefit of the ftp family is that we can often and easily prove properties by induction on the task set cardinality, but first we recall a definition: Definition 22 (𝑆 𝑛 [START_REF] Goossens | The Non-Optimality of the Monotonic Priority Assignments for Hard Real-Time Offset Free Systems[END_REF]).

𝑆 1 𝑂 1 , 𝑆 𝑖 𝑂 𝑖 + ⌈ (𝑆 𝑖-1 -𝑂 𝑖 ) 0 𝑇 𝑖 ⌉𝑇 𝑖 .
Informally speaking 𝑆 𝑖 (𝑖 > 1) corresponds to the first release of 𝜏 𝑖 not before 𝑆 𝑖-1 .

Theorem 23 (Simulation interval for ftp). Any ftp-schedulable asynchronous constrained deadline periodic tasks with binary delays upon uniprocessor platform reaches a cycle 4 at or prior to:

𝑆 𝑛 + 𝐻 . (6)
Proof. By induction. We first assume that, wlog, the ftp scheduler corresponds to 𝜏 1 > 𝜏 2 > • • • > 𝜏 𝑛 . Let P (𝑖) be the proposition: if the task set {𝜏 1 , 𝜏 2 , . . . , 𝜏 𝑖 } is ftp-schedulable then the schedule is periodic from 𝑆 𝑖 with a period of 𝐻 𝑖 .

Base case: P (1) is true: if feasible the schedule for 𝜏 1 is periodic of period 𝑇 1 = 𝐻 1 from the first release of 𝜏 1 (𝑆 1 = 𝑂 1 ).

Inductive step: For 𝑖 > 1, P (𝑖 -1) ⇒ P (𝑖) is true. Assume P (𝑖 -1) is true for purpose of induction. We need to show that P (𝑖) is true. P (𝑖 -1) means that if the task set {𝜏 1 , 𝜏 2 , . . . , 𝜏 𝑖-1 } is ftp-schedulable then the schedule of {𝜏 1 , . . . , 𝜏 𝑖-1 } is periodic from 𝑆 𝑖-1 with a period of 𝐻 𝑖-1 . Now we consider the scheduling of {𝜏 1 , . . . , 𝜏 𝑖 }, since 𝜏 𝑖 is the lowest priority task and since there is no priority inversion with binary delays (Proposition 1) the schedule and thus periodicity of {𝜏 1 , . . . , 𝜏 𝑖-1 } is unchanged by the requests of task 𝜏 𝑖 . 𝑆 𝑖 is the first release of task 𝜏 𝑖 after (or at) 𝑆 𝑖-1 (𝑆 𝑖 ≥ 𝑆 𝑖-1 ), combining the periodicity of the schedule of {𝜏 1 , . . . , 𝜏 𝑖-1 } and the periodicity of 𝜏 𝑖 (from 𝑆 𝑖 with a period of 𝑇 𝑖 ) we conclude that P (𝑖) is true: if {𝜏 1 , . . . , 𝜏 𝑖 } is schedulable, its ftp schedule is periodic from 𝑆 𝑖 with a period of lcm{𝐻 𝑖-1 ,𝑇 𝑖 } = lcm{𝑇 𝑗 | 𝑗 = 1, . . . , 𝑖} = 𝐻 𝑖 . □

DISCUSSION

In this section, we open the discussion on other properties that could be brought by considering binary preemption delays. We review ways to reduce the pessimism of our simulation interval for any deterministic and memoryless scheduler. Finally, we discuss the future work needed to evolve our model to take into account other interference sources from the platform on the scheduler.

𝐶-sustainability

We believe that the absence of priority inversion will lead to other positive and better results for the particular case of binary delays. We conjecture that work-conserving priority driven scheduler are C-sustainable [START_REF] Burns | Sustainability in real-time scheduling[END_REF] for binary delays. The later property means that considering the worst-case execution requirement (𝐶 𝑖 ) for each task is actually the worst-case regarding schedulability. Consequently, in this case the simulation interval is equivalent to an exact simulation based schedulability test. An open question here is whether the assumption of binary delays is realistic or not? To go further, we could also imagine systems where this hypothesis is forced by design when choosing the clock frequencies on an embedded board, for example, or by the tick value of the scheduler if it is a tick scheduling type.

Reducing the pessimism of the general result

The general result, the simulation interval [0, 𝐻 • (𝑛 + 1) • (𝛼 max + 1) 𝑛 𝑖=1 ((𝑂 𝑖 + 𝐷 𝑖 -𝑇 𝑖 ) 0 + 1)), can be computed very easily (linear complexity), but is large and pessimistic. This is probably the price to pay for such a general result which is applicable for all schedulers one can imagine (but deterministic and memoryless). Please notice that the scheduler is not necessarily preemptive, not necessarily work-conserving, not necessarily "real-time". . .

In order to reduce the size/pessimism of the interval we could reduce the parameters of the system state. In particular, we believe that many real-time schedulers (e.g., edf, rm, dm, ftp) do not need the totality of the parameters. For example edf does not need the 𝑟 parameter because priorities do not change during a reload activity.

Extending the model

To extend this work, we will also need to consider slightly different models where sections corresponding to the return of a preemption remain preemptive, but must be completely re-executed if preempted. We will also need to extend this property to the other components considered in Section 2.2 in order to model not only preemptions but more generally the scheduling cost of the platform/operating system. In such case, we could also consider a task model with deadline prior to completion as proposed in [START_REF] Burns | Fixed priority scheduling with deadlines prior to completion[END_REF]. Finally, a future extension will be to consider multiprocessor platforms.

CONCLUSION

In this paper we have considered the notion of simulation interval for asynchronous arbitrary deadline tasks with preemption delays. We have designed and proved correct such an interval for any deterministic and memoryless scheduler upon uniprocessor. We then studied particular cases where significantly shorter simulation intervals can be designed for popular real-time schedulers. Firstly, for edf we extended the work of Leung and Merril, secondly for ftp we extended the work [START_REF] Goossens | The Non-Optimality of the Monotonic Priority Assignments for Hard Real-Time Offset Free Systems[END_REF], in both cases for binary delays and constrained deadlines. We also discussed the scope of the results, paths for reducing pessimism as well as other potential results and future work opened by this one.

Figure 1 :

 1 Figure 1: A periodic schedule of periodic tasks.

◎

  2 2022-03-08 15:18. Page 2 of 1-10.

  4 2022-03-08 15:18. Page 4 of 1-10. uniprocessor real-time schedulers with preemption delay RTNS '22, June 07-08, 2022, Paris, France

  2022-03-08 15:18. Page 7 of 1-10. uniprocessor real-time schedulers with preemption delay RTNS '22, June 07-08, 2022, Paris, France

  when a job of 𝜏 𝑖 is resumed after a preemption. It represents the time required to load the context of execution of the preempted job. We assume this activity non-preemptive.More formally a task 𝜏 𝑖 is characterised by the tuple (𝑂 𝑖 , 𝐶 𝑖 , 𝑇 𝑖 , 𝐷 𝑖 , 𝛼 𝑖 ). A system Sys = {𝜏 1 , . . . , 𝜏 𝑛 } is a task set, where every task 𝜏 𝑖 is defined by: • 𝑂 𝑖 ∈ {N ∪ 0} the task offset, i.e., the release date of the first job of 𝜏 𝑖 , • 𝐶 𝑖 ∈ N the worst-case execution time (wcet), i.e., the maximum cumulative cpu units required for a job of 𝜏 𝑖 to be executed, • 𝑇 𝑖 ∈ N the task period, the 𝑘 th job 𝜏 𝑖,𝑘 of 𝜏 𝑖 is released at the instants 𝑂 𝑖 + (𝑘 -1)𝑇 𝑖 , 𝑘 ∈ N, • 𝐷 𝑖 ∈ N is the relative deadline and represents the timing constraint of a task: the 𝑘 th , (𝑘 ∈ N) job 𝜏 𝑖,𝑘 of 𝜏 𝑖 must be completely executed in the window [𝑂 𝑖 + (𝑘 -1)𝑇 𝑖 , 𝑂 𝑖 + (𝑘 -1)𝑇 𝑖 + 𝐷 𝑖 ). We consider arbitrary deadline where there is no constraint between task deadline and period. Particular cases include constrained deadlines where ∀𝑖 ∈ {1 . . . 𝑛}, 𝐷 𝑖 ≤ 𝑇 𝑖 , and implicit deadlines where ∀𝑖 ∈ {1 . . . 𝑛}, 𝐷 𝑖 = 𝑇 𝑖 . When considering arbitrary deadline it is admitted to have 𝐷 𝑖 > 𝑇 𝑖 , we assume to schedule the various active jobs of the same task in the fifo order (the oldest job first). • 𝛼 𝑖 ∈ {N ∪ 0} represents the time required by the rtos/hw to load non-preemptively the context of execution of each preempted job of task 𝜏 𝑖 . The following parameters can be deduced 2 : • 𝐻 lcm(𝑇 1 , . . . ,𝑇 𝑛 ) is the hyperperiod of the system, where lcm denotes the least common multiple, • 𝐻 𝑖 lcm{𝑇 𝑗 | 𝑗 = 1, . . . , 𝑖}, • 𝑂 max max 𝑖=1,...,𝑛 (𝑂 𝑖 ) is the largest offset, • 𝛼 max max 𝑖=1,...,𝑛 (𝛼 𝑖 ) is the largest preemption delays,

	• 𝑈	𝑛 𝑖=1 𝐶 𝑖 /𝑇 𝑖 is the processor utilization factor.
	We will use the notations:
	• (𝑎) 0 meaning (𝑎) 0 max(𝑎, 0),
	• [𝑥, 𝑦) denoting an half-open interval: {𝑧 | 𝑥 ≤ 𝑧 < 𝑦}
	Please notice that, in our model, we consider discrete time.

  𝐶 rem 1 , . . . , 𝐶 rem 𝑛 , Ω 1 , . . . , Ω 𝑛 , where • Ω 𝑖 is the local clock of 𝜏 𝑖 , undefined before 𝑂 𝑖 , initialized at 0 at the time 𝑂 𝑖 , being reset at every new request of the task. Formally, at time 𝑡 ≥ 𝑂 𝑖 , Ω 𝑖 (𝑡 -𝑂 𝑖 ) mod 𝑇 𝑖 , • while 𝐶 rem 𝑖 is the cumulative remaining work to process for 𝜏 𝑖 (can be greater than 𝐶 𝑖 since it is allowed to have 𝐷 𝑖 > 𝑇 𝑖

  𝑛} the task index of the running/reloaded task (if any) or 0 if the processor is idle. As we have nonpreemptive reloading phases, the decision at time 𝑡 + 1 depends on the decision taken at time 𝑡.• rem-reload ∈ {0, 1, . . . , 𝛼 max } the remaining reload time (0 means that a task executes). , . . . , 𝐶 rem 𝑛 , Ω 1 , . . . , Ω 𝑛 ,

	More formally:
	Definition 12 (State and pre-state of a system). The state
	of a system of 𝑛 tasks can be defined as a (2𝑛)-tuple
	𝑆	𝐶 rem 1

𝑟 ∈ {0, 1, . . . , 𝑛}, and rem-reload ∈ {0, 1, . . . , 𝛼 max }. See Definition 4 for details (concerning the already defined attributes).

We will now formalize the remaining reload time in state, as follows:

rem-reload(𝑡 + 1) rem-reload(𝑡) -1 if rem-reload(𝑡) > 0

(2)

Table 1 :

 1 Task set characteristics 𝑇 𝑖 𝐷 𝑖 𝐶 𝑖 𝑂 𝑖 𝛼 𝑖

	𝜏 1 12 12 2	0	2
	𝜏 2 6	6	1	1	2
	𝜏 3 12 8	1	3	2
	𝜏 4 12 3	2	6	2
	Figure 2: edf schedule of the system described by Table 1

  time 18 task 𝜏 1 can be preempted (by task 𝜏 4 )immediately. Of course we could extend the Leung and Merril's formal definition of system state (𝐶 𝑆 ) to take account of reloads as we did in the Definition 12 to formally have 𝐶 𝑆 (𝑅, 𝑡 1 ) ≠ 𝐶 𝑆 (𝑅, 𝑡 2 ) but this is not required.

	𝜏 4 𝜏 3 𝜏 2 𝜏 1	? ?	? 3 ⊲⊳	? 6 ⊲⊳ ? 7 8 e 5 6	e	e	? e ? e	?	? 18 19 ? e	e	e	? e ? e	? 27 ⊲⊳	? 30 ⊲⊳ ? 31 32 e	e	p p p p p	p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p e	p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p ? p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p ? p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p e e	p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p ? p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p	p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p	p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p ? p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p ? p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p e	p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p e	p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p	p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p e	p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p ? e
	activity while at

In particular cases, it is possible to define schedulers whose periodicity does not depend on the arrival of tasks, but in practice they make no sense.

i.e., all reachable states are visited. 2022-03-08 15:18. Page

of 1-10.

if the scheduler has decided to idle the processor at time 𝑡 (i.e., 𝑟 (𝑡 ) = 0) we decide arbitrarily to make a reload in this case.

active. But this is impossible since every task 𝜏 𝑖 in the system has an initial request time 𝑂 𝑖 .

Case b. 𝑆 (𝑡 ′

1 ) = 𝐼 (and 𝑆 (𝑡 ′ 2 = 𝐸 ℓ 1 )), this is impossible since 𝜏 ℓ 1 is active at time instant 𝑡 ′ 1 and edf is work-conserving.

Case c. 𝑆 (𝑡 ′ 1 ) = 𝑅 ℓ 2 , ℓ 2 ≠ ℓ 1 (and 𝑆 (𝑡 ′ 2 = 𝐸 ℓ 1 )), task 𝜏 ℓ 2 is an edf higher priority task, active at time instant 𝑡 ′ 1 but not at time instant 𝑡 ′ 2 . Using the very same argument than case a we can produce an infinite progression of task computations which leads to a contradiction.

Case d. 𝑆 (𝑡 ′

1 ) = 𝑅 ℓ 1 (and 𝑆 (𝑡 ′ 2 = 𝐸 ℓ 1 )). We will consider all subcases where this can happen with binary delays.

Case d.1. 𝑆 (𝑡 ′

1 -1) = 𝐼 this impossible since at time 𝑡 ′ 1 -1 there is at least one active task (𝜏 ℓ 1 ).

1 we must have 𝑆 (𝑡 ′ 1 ) = 𝐸 𝑗 (which is not the case), if there is an request of task 𝜏 𝑗 at time instant 𝑡 1 but then we must have 𝑆 (𝑡 ′ 1 ) = 𝐸 𝑗 because it is a new job that begins its execution (which is not the case either).

Case d.3. 𝑆 (𝑡 ′

1 -1) = 𝑅 ℓ 1 this impossible for binary delays.

Case d. [START_REF] Baro | Off-Line (Optimal) Multiprocessor Scheduling of Dependent Periodic Tasks[END_REF].

Once again, using the very same argument than case a we can produce an infinite progression of task computations which leads to a contradiction. □ Lemma 20. Let 𝑆 be schedule of a task system 𝑅 constructed by the edf algorithm with binary delays.

Proof. When 𝑅 is edf-schedulable we will first show (case a) that ∀𝑖, 𝑒 𝑖,𝑡 1 = 𝑒 𝑖,𝑡 2 , and then we will show that (case b) if ∀𝑖, 𝑒 𝑖,𝑡 1 = 𝑒 𝑖,𝑡 2 then that the system states at time instant 𝑡 1 and 𝑡 2 cannot differ on the other two parameters (𝑟 and rem-load).

Case a. We will show the property by contradiction. We assume, for the purpose of the contraction, that ∃ ℓ s.t. 𝑒 ℓ,𝑡 1 ≠ 𝑒 ℓ,𝑡 2 . By Lemma 19, we must have 𝑒 ℓ,𝑡 1 > 𝑒 ℓ,𝑡 2 . We first show that during the time interval [𝑡 1 , 𝑡 2 ) there is no idle unit. We will prove this by contradiction as well. For the purpose of the contradiction we will assume an idle unit during the time-slot [𝑡 1 + Δ, 𝑡 1 + Δ + 1). This imply that all task computations requested prior to 𝑡 1 + Δ have finished their execution. Consequently, ∀𝑖, 𝑒 𝑖,𝑡 1 +Δ = 𝐶 𝑖 . By definition we have ∀𝑖, 𝑒 𝑖,𝑂 max +Δ ≤ 𝐶 𝑖 , by Lemma 19 we have ∀𝑖, 𝑒 𝑖,𝑂 max +Δ ≥ 𝑒 𝑖,𝑡 1 +Δ . We just proved that ∀𝑖, 𝑒 𝑖,𝑡 1 +Δ = 𝐶 𝑖 . Consequently, ∀𝑖, 𝑒 𝑖,𝑂 max +Δ = 𝐶 𝑖 . In other words, all task computations requested prior to 𝑂 max + Δ have finished their execution, as well. Since the task requests in the interval [𝑂 max + Δ, 𝑡 1 + Δ) and [𝑡 1 + Δ, 𝑡 2 + Δ) are the same and since ∀𝑖, 𝑒 𝑖,𝑂 max +Δ = 𝑒 𝑖,𝑡 1 +Δ = 𝐶 𝑖 , the schedule in the interval

and 𝑡 2 -𝑡 1 = 𝐻 ; we conclude that ∀𝑖, 𝑒 𝑖,𝑡 1 = 𝑒 𝑖,𝑡 2 contradicting our assumption that ∃ ℓ s.t. 𝑒 ℓ,𝑡 1 > 𝑒 ℓ,𝑡 2 . We have just proved that there is no idle unit in the time interval [𝑡 1 , 𝑡 2 ). Let us remember that we assume, for the purpose of the contradiction, that ∃ ℓ s.t. 𝑒 ℓ,𝑡 1 > 𝑒 ℓ,𝑡 2 . During a hyperperiod [𝑡 1 , 𝑡 2 ) the processor is fully busy (executing or reloading tasks) but at time 𝑡 2 the task 𝜏 ℓ is delayed in comparison with time instant 𝑡 1 . As edf decisions do not rely on the remaining computation times edf will make the same decisions in the following hyperperiods which will sooner or later lead to a deadline miss. Hence, 𝑅 cannot be edf-schedulable on one processor. In other words, during each hyper-period what the system is able to schedule is less than the processor demand that arrives in each hyper-period. This does not mean that 𝑈 > 1 because the processor spends time on execution and reloads but a backlog will accumulate and sooner or later a deadline miss will occur.

Case b. We now know that ∀𝑖, 𝑒 𝑖,𝑡 1 = 𝑒 𝑖,𝑡 2 we will show that the system states at time instant 𝑡 1 and 𝑡 2 cannot differ on the other two parameters (𝑟 and rem-load). Firstly, notice that if there is no active task at time 𝑡 1 (or equivalently 𝑡 2 ) mandatory edf idles the processor. Secondly, since ∀𝑖, 𝑒 𝑖,𝑡 1 = 𝑒 𝑖,𝑡 2 the active tasks are identical at both time instants 𝑡 1 and 𝑡 2 , edf will choose the same task for execution/reloading (let's say task 𝜏 ℓ ). If the states differ (and consequently the edf decisions) we have to distinguish between two cases:

where 𝑆 (𝑡) denotes the schedule decision at time 𝑡 (𝐸 ℓ meaning the execution of 𝜏 ℓ , and 𝑅 ℓ meaning the reloading of 𝜏 ℓ ).

Proof. This is a direct consequence of Lemma 20: if the system (𝑅) is edf-schedulable on a single processor then C𝑆 (𝑅, 𝑂 max +𝐻 ) = C𝑆 (𝑅, 𝑂 max + 2 • 𝐻 ) and consequently the schedule repeats from

Please notice that, as the original Leung and Merril's result, it is possible to miss a (first) deadline after time instant 𝑂 max + 2 • 𝐻 . To ensure the periodicity/schedulability it is necessary to test the condition: C𝑆 (𝑅, 𝑂 max + 𝐻 ) = C𝑆 (𝑅, 𝑂 max + 2 • 𝐻 ). Unfortunately, the condition 𝑈 ≤ 1 is not sufficient for system with preemption delays.


Brian Stout
email: brian.stout@fresnel.fr

J C Auger
email: augerjc@gmail.com

Alexis Devilez
email: alexis.devilez@fresnel.fr

Recursive T matrix algorithm for resonant multiple scattering : Applications to localized plasmon excitations

Keywords: Multiple scattering, T matrices, morphology dependent resonances, plasmon resonances

A matrix balanced version of the Recursive Centered T Matrix Algorithm (RCTMA) applicable to systems possessing resonant interparticle couplings is presented. Possible domains of application include systems containing interacting localized plasmon resonances, surface resonances, and photonic jet phenomena. This method is of particular interest when considering modifications to complex systems. The numerical accuracy of this technique is demonstrated in a study of particles with strongly interacting localized plasmon resonances.

Introduction

It has been well established that certain kinds of recursive T matrix algorithms (known as RCTMA) [START_REF] Stout | A transfer matrix approach to local field calculations in multiple scattering problems[END_REF][START_REF] Auger | A recursive centered T-Matrix algorithm to solve the multiple scattering equation : numerical validation[END_REF] are numerically stable and can be used to solve the Foldy-Lax multiplescattering equations for particles exhibiting "modest" inter-particle couplings. By "modest couplings", we refer to situations in which the order of orbital number of the Vector Spherical Wave Functions (VSWFs) necessary to describe the field scattered by each particle in an aggregate of particles are not too much larger than that necessary for describing isolated particles. The "modest coupling" criteria apply to a host of multiple scattering situations, including systems of dielectric particles comparable in size to the wavelength and for most packing fractions including dense packing. The modest coupling criteria can also apply to metallic particles under certain conditions.

Like any multiple scattering technique not employing matrix balancing, the RCTMA can encounter numerical difficulties in certain extreme situations of strongly coupled resonant phenomenon. In this work, we present a matrix balanced form of the Recursive Centered T Matrix Algorithm (or RCTMA) that can readily be employed even in the presence of strong (i.e. resonant) inter-particle couplings. The rather extreme situation of "strong couplings" studied here will generally require carefully micro-scaled engineered systems where high Q-factor resonances can occur for particles illuminated in isolation, and in which the particles are sufficiently closely spaced that neighboring particles modify the resonance response properties. Examples of strong inter-particle couplings can be found in particles exhibiting plasmon resonances, surface resonances, or even photonic jet phenomenon.

In section 2, the notation is introduced in a brief review of the relevant multiple-scattering theory. Section 3 describes an analytic matrix balancing procedure used to 'well-condition' the multiple scattering system of equations. A matrix balanced RCTMA is derived in section 4. Essential formulas for applications are summarized in section 5. Their applications are then demonstrated by applying matrix balanced RCTMA calculations to study systems of interacting localized plasmon excitations. Some known and novel aspects of interacting localized plasmon excitations are presented herein.

Multiple-scattering theory -VSWF approach

Let us consider an arbitrary incident electromagnetic field incident on a collection of threedimensional particles (as shown in fig. 1). The particles are considered as 'individual' scatterers if they can be placed in a circumscribing sphere lying entirely within the homogeneous medium (actually this constraint can frequently be relaxed, cf. [START_REF] Doicu | Light Scattering by Systems of Particles[END_REF] ).

The electromagnetic field incident on an N-particle system, E i , is developed in terms of the transverse regular VSWFs developed about some point O arbitrarily chosen as the system 'origin':

x 2 x 3 x N O x y z x 1 R 1 R 2 R N R 3 E inc
E i (r) = E 0 ∞ n=1 n m=-n {Rg [M nm (kr)] a 1,n,m + Rg [N nm (kr)] a 2,n,m } = E 0 2 q=1 ∞ p=1 Rg [Ψ q,p (kr)] a q,p ≡ E 0 Rg Ψ t (kr) a (1)
where E 0 is a real parameter determining the incident field amplitude. Eq.( 1), introduces a condensed notation for the VSWFs, M nm and N nm : Ψ 1,p (kr) ≡ M n,m (kr) and Ψ 2,p (kr) ≡

N n,m (kr). The notation Rg [ ] stands for "the regular part of" and distinguishes these regular VSWFs from the "irregular" scattered VSWFs (cf. appendix A). In the second line of eq.( 1), the two subscripts (n, m) are replaced by a single subscript p defined such that p (n, m) ≡ n(n + 1)m and has the inverse relations [START_REF] Tsang | Theory of Microwave Remote Sensing[END_REF] :

n(p) = Int √ p m(p) = -p + n(n + 1) . (2)
The last line of eq.( 1), adopts the compact matrix notation allowing the suppression of the summation symbols. [START_REF] Chew | Waves and Fields in Inhomogeneous Media[END_REF] The superscripted ( t ) stands for the transpose of a column 'matrix' of composed of VSWFs into a row 'matrix' of these functions.

For points external to all individual circumscribing spheres, the total field, E t (r) can be written as the sum of the incident field, and a set of 'individual' scattered fields, E

s , centered respectively on each of the particle centers:

E t (r) = E i (r) + N j=1 E (j) s (r j ) = E 0 Rg Ψ t (kr) a + E 0 N j=1 Ψ t (kr j ) f (j) N (3)
where each scattered field, E

s , is developed, with coefficients f

(j)
N , on a basis of outgoing VSWFs defined with respect to the associated particle center, denoted x j . The spherical coordinates relative to each scatterer are denoted r j ≡ rx j .

The crucial idea of Foldy-Lax multiple-scattering theory is that there exists an excitation field, E (j) exc (r), associated with each particle which is the superposition of the incident field and the field scattered by all the other particles in the system (excluding the field scattered by the particle itself). [START_REF] Lax | Multiple Scattering of Waves[END_REF] From this definition, the excitation field of the j th particle can be written

E (j) exc (r j ) ≡ E 0 Rg Ψ t (kr j ) e (j) N ≡ E i (r) + N l=1,l =j E (l) s (r l ) = E 0 Rg Ψ t (kr j ) J (j,0) a + N l=1,l =j H (j,l) f (l) N (4)
where e

(j)

N are the coefficients of the excitation field in a regular VSWF basis centered on the j th particle. In the second line of eq.( 4), we have used the translation-addition theorem [START_REF] Stout | A transfer matrix approach to local field calculations in multiple scattering problems[END_REF][START_REF] Tsang | Theory of Microwave Remote Sensing[END_REF][START_REF] Chew | Waves and Fields in Inhomogeneous Media[END_REF] and introduced the notation where J (j,0) ≡ J (kx j ) is a regular translation matrix and

H (j,l) ≡ H [k (x j -x l )
] is an irregular translation matrix. Analytical expressions for the matrix elements of J (kx j ) and H (kx j ) are given in refs. [START_REF] Stout | A transfer matrix approach to local field calculations in multiple scattering problems[END_REF][START_REF] Tsang | Theory of Microwave Remote Sensing[END_REF] .

The other key idea of multiple scattering theory is that the field scattered by the object,

f (j)
N , is obtained from the excitation field e (j) N via the 1-body T matrix, T (j) 1 , derived when one considers the particle to be immersed in an infinite homogeneous medium. This relation is then expressed as

f (j) N = T (j) 1 e (j) N (5)
[The index 1 on the T (j) 1 indicates that this T matrix concerns an isolated particle, henceforth referred to as a '1-body' T matrix.] Employing eq.( 5) in eq.( 4), one obtains a Foldy-Lax set of equations for the excitation field coefficients [START_REF] Stout | A transfer matrix approach to local field calculations in multiple scattering problems[END_REF][START_REF] Chew | Waves and Fields in Inhomogeneous Media[END_REF] :

e (j) N = J (j,0) a + N l=1,l =j H (j,l) T (l) 1 e (l) N j = 1, ..., N (6)
For numerical applications where one is obliged to solve the equations on a truncated VSWF basis, it is advantageous to work with a set of formally equivalent equations involving the scattering coefficients f (j)

N . This set of equations is derived by multiplying each of eqs.(6) from the left by T (j) 1 and using eq.( 5) to obtain

f (j) N = T (j) 1 J (j,0) a + T (j) 1 N l=1,l =j H (j,l) f (l) N j = 1, ..., N (7)
In the RCTMA, one calculates the centered multiple scattering transition matrices, T (j,k) N , which directly yield the scattered field coefficients in terms of the field incident on the system through the expression

f (j) N = N k=1 T (j,k) N a (k) with a (k) ≡ J (k,0) a (8)
In this equation, we have introduced the column matrix a (j) which contains the coefficients of the field incident on the entire system developed on a VSWF basis centered on the j th particle.

Basis set truncation and matrix balancing

Although the multiple scattering formulas of the previous section are expressed as matrix equations on VSWF basis sets of infinite dimension, the finite size of the scatterers naturally restricts the dimension of the dominate VSWF contributions. In order to discuss this phenomenon analytically, we consider the case of spherical scatterers. For non-spherical scatterers, the matrix balancing procedure described below should be applied to the circumscribing spheres of the particles.

The Mie solution for a sphere of radius R j immersed in a homogeneous host medium, can be cast in the form of a 1-body T matrix that is diagonal in a VSWF basis centered on the particle:

T (j) 1 q,p;q ′ ,p ′ = δ q,q ′ δ p,p ′ T 1 (j, n(p), q) (9)
where the T (j)

1 (n(p), q) correspond to the Mie coefficients and depend on q and n (cf. eq.( 1)).

With the objective of matrix balancing, it is helpful to express the Mie coefficients of the scatterers in terms of the Ricatti Bessel and Hankel functions, respectively ψ n (z) ≡ zj n (z) and ξ n (z) ≡ zh n (z), and their logarithmic derivatives

Φ n (z) ≡ ψ ′ n (z) ψ n (z) Ψ n (z) ≡ ξ ′ n (z) ξ n (z) (10)
The T matrix elements of eq.( 9) for a sphere of dielectric contrast ρ j ≡ k j /k can then be cast in the convenient form [START_REF] Stout | Absorption in multiple scattering systems of coated spheres[END_REF] :

T (j, n, 1) = ψ n (kR j ) ξ n (kR j ) µ j µ Φ n (kR j ) -ρ j Φ n (ρ j kR j ) ρ j Φ n (ρ j kR j ) - µ j µ Ψ n (kR j ) ≡ ψ n (kR j ) ξ n (kR j ) T (j, n, 1) T (j, n, 2) = ψ n (kR j ) ξ n (kR j ) µ j µ Φ n ρ j kR j -ρ j Φ n (kR j ) ρ j Ψ n (kR j ) - µ j µ Φ n (ρ j kR j ) ≡ ψ n (kR j ) ξ n (kR j ) T (j, n, 2) ( 11
)
where k is the wavenumber in the external medium. The normalized T matrix coefficients, T (j, n, q), of eq.(MieT) contain a rich resonant structure. The ratios ψ n (kR j ) /ξ n (kR j )

on the other hand have an exponentially decreasing behavior for large, n ≫ kR j as is demonstrated in fig. 2 for kR = 10. One can remark from figure 2 that |ψ n (kR) /ξ n (kR)| become quite small beyond n max = kR + 3 and its value at n = 14 is ∼ 2 10 -4 . Although these factors permit an appropriately truncated VSWF basis set to contain essentially all the physical information necessary for accurate calculations, they also tend to produce illconditioned linear systems when one is obliged to enlarge the VSWF space far beyond ≈ kR + 3 in order to account for strong coupling phenomenon.

A solution to the above problem is to 'balance' the matrix manipulations in section 4 below by defining "normalized" scattering and incident coefficients:

f (j) q,p ≡ ξ n(p) (kR j ) f (j) q,p a (j) q,p ≡ ψ n(p) (kR j ) a (j) q,p (12)
For notational purposes, it is convenient to define diagonal matrices ψ (j) and ξ (j) with Ricatti-Bessel functions along their diagonals, namely ψ (j) q ′ ,p ′ ;q,p ≡ δ q,q ′ δ p,p ′ ψ n(p) (kR j ) and ξ (j) q ′ ,p ′ ;q,p ≡ δ q,q ′ δ p,p ′ ξ n(p) (kR j ). This notation allows normalized or 'balanced' versions of the one-body and many-body T -matrices to be defined respectively as

T (j) 1 ≡ ξ (j) T (j) 1 ψ (j) -1
and In terms of these normalized quantities, eq.( 8) then reads

T (j,k) N ≡ ξ (j) T (j,k) N ψ (k) -1 (13)
f (j) N = N k=1 T (j,k) N a (k) j = 1, ..., N (14)
In the next section, these 'normalized' T

Derivation of a matrix balanced recursive algorithm

In this section, we derive a matrix balanced version of the Recursive Centered T Matrix Algorithm (RCTMA) using purely algebraic manipulations. The recursive algorithm can be invoked once we have a solution for the T (j,k)

N -1 matrices of a N ≥ 1 particle system. If we wish to solely use the recursive algorithm to solve a system, we initiate the recursive process with a single particle solution described by

T (1,1) 1 ≡ T (1) 1 .
One then considers an arbitrarily positioned particle being added to the system. The excitation field on a particle N added to the system can be expressed as the superposition of three fields. The first contribution is simply the field incident on the system, the second contribution results from the scattering of the incident field by the N -1 cluster of particles onto the particle N, and finally the third contribution comes from field scattered by the particle N onto the N -1 cluster and which returns to the N th particle as an excitation field. Invoking the translation-addition theorem and eq.( 8), these three contributions can be expressed in matrix form as 1 e (N )

N = a (N ) + N -1 j,k=1 H (N,j) T (j,k) N -1 a (k) + N -1 j,k=1 H (N,j) T (j,k) N -1 H (k,N ) f (N ) N (15)
Defining now the normalized irregular translation matrices and excitation coefficients respectively as

H (j,k) ≡ ψ (j) H (j,k) ξ (k) -1
and e

(j)

N ≡ ψ (j) e (j) N (16)
the normalized form of eq.( 15) can be written

e (N ) N = a (N ) + N -1 j,k=1 H (N,j) T (j,k) N -1 a (k) + N -1 j,k=1 H (N,j) T (j,k) N -1 H (k,N ) f (N ) N ( 17
)
where we also used the definitions in eqs.( 12) and (13).

Recalling that the excitation field is linked to the scattered field by the 1-body T -matrix via eq.( 5), and invoking the definitions of eq.( 12) and ( 16) we can write e (j)

N = T (j) 1 -1 f (j) N (18)
Employing this relation for particle N on the LHS of eq.( 17) and rearranging we obtain

T (N ) 1 -1 - N -1 j,k=1 H (N,j) T (j,k) N -1 H (k,N ) f (N ) N = a (N ) + N -1 j,k=1 H (N,j) T (j,k) N -1 a (k) . ( 19
)
We now take the normalized T (N,N ) N matrix to be expressed as

T (N,N ) N = T (N ) 1 -1 - N -1 j,k=1 H (N,j) T (j,k) N -1 H (k,N ) -1 (20)
With this assignment, we multiply both sides of eq.( 19) by T (N,N ) N

and obtain an expression consistent with equation ( 14):

f (N ) N = T (N,N ) N a (N ) + T (N,N ) N N -1 j,l=1 H (N,j) T (j,k) N -1 a (k) = T (N,N ) N a (N ) + N -1 l=1 T (N,k) N a (k) = N l=1 T (N,k) N a (k) (21)
where we have assigned the matrix T

(N,k) N , k = N as T (N,k) N = T (N,N ) N N -1 j=1 H (N,j) T (j,k) N -1 (22)
One completes the description of the scattering by the system by remarking that the field scattered by the other particles in the system are the superposition of the field that would be scattered by the N -1 cluster in the absence of the N th particle, plus the field scattered from the N -1 particle originating as a scattered field emanating from the N th particle.

Using again the translation-addition theorem, the field coefficients of f

(j)
N can in turn be expressed in a form consistent with equation ( 14) as

f (j) N = N -1 l=k T (j,k) N -1 a (k) + N -1 k=1 T (j,k) N -1 H (k,N ) f (N ) N = N -1 l=1 T (j,k) N -1 a (k) + N -1 k=1 T (j,k) N -1 H (k,N ) T (N,N ) N a (N ) + N -1 l=1 N -1 k=1 T (j,l) N -1 H (l,N ) T (N,k) N a (k) = T (j,N ) N a (N ) + N -1 k=1 T (j,k) N a (k) = N k=1 T (j,k) N a (k) (23)
where we invoked eq.( 21). In the second and third lines we have defined the T (j,N ) N and T

(j,l) N matrices such that

T (j,N ) N = N -1 k=1 T (j,k) N -1 H (k,N ) T (N,N ) N (24a) T (j,k) N = T (j,k) N -1 + N -1 l=1 T (j,l) N -1 H (l,N ) T (N,k) N (24b)
At this point, all the T (j,k) N matrices have been obtained and the matrix manipulations in eqs.(20), ( 22) and (24) can then be repeated to add as many particles to the system as desired.

A. Relationship with system matrix inversions

Although the recursive algorithm is quite efficient for systems with relatively small numbers of particles, for systems with many particles, one may prefer to try and solve an entire N-particle system directly. A balanced linear system for the entire system corresponding to our recursive algorithm can be obtained by applying the relation of eq.( 5) to the left hand side of eq.( 7), then multiplying both sides of the resulting equations by the ψ (j) matrix and finally rearranging to obtain a system of balanced linear equations for the unknown scattering coefficients:

T (j) 1 -1 f (j) N - N k=1,k =j H (j,k) f (k) N = a (j) j = 1, ..., N (25)
where we used eqs.( 12) and (13). The system of linear equations in eq.( 25) can in principle be directly solved by inverting the balanced system matrix:

        f (1) N f (2) N . . . f (N ) N         =          T (1) 1 -1 -H (1,2) • • • -H (1,N )
-H

(2,1)

T (2) 1 -1 • • • -H (2,N ) . . . . . . . . . . . . -H (N,1) -H (N,2) • • • T (N ) 1 -1          -1         a (1)
a [START_REF] Auger | A recursive centered T-Matrix algorithm to solve the multiple scattering equation : numerical validation[END_REF] . . .

a (N )         (26)
Once we have inverted this system, one can associate each block with a corresponding T (j,k) N matrix as shown below as

        f (1) N f (2) N . . . f (N ) N         =         T (1,1) N T (1,2) N • • • T (1,N ) N T (2,1) N T (2,2) N • • • T (2,N ) N . . . . . . . . . . . . T (N,1) N T (N,2) N • • • T (N,N ) N                 a (1)
a [START_REF] Auger | A recursive centered T-Matrix algorithm to solve the multiple scattering equation : numerical validation[END_REF] . . .

a (N )         (27)
which is the same form as the desired solutions given in eq.( 14).

Summary and applications to localized plasmon excitations

In this section, we will apply the RCTMA to solve systems exhibiting strong interactions between localized plasma resonances. We begin this section by summarizing the balanced recursive algorithm. We then recall some useful formulas for extracting physical quantities from the T matrix. Finally, we carry out some illustrative calculations for strongly interacting systems.

A. Summary of the balanced RCTMA algorithm

In order to implement the RCTMA, one must first solve the 1-body T -matrices, T

1 ,T

1 ,...,

T (Ntot) 1
, for all the particles in the system. Normalized versions of the 1-body T matrices and the irregular translation matrices [START_REF] Stout | A transfer matrix approach to local field calculations in multiple scattering problems[END_REF][START_REF] Tsang | Theory of Microwave Remote Sensing[END_REF] , H (j,k) , are then calculated via

T (j) 1 ≡ ξ (j) T (j) 1 ψ (j) -1 H (j,k) ≡ ψ (j) H (j,k) ξ (k) -1 (28)
where the diagonal matrices ψ (j) q,q ′ ,p,p ′ = δ q,q ′ δ p,p ′ ψ n(p) (kR j ) and ξ (j) q,q ′ ,p,p ′ = δ q,q ′ δ p,p ′ ξ n(p) (kR j ) respectively have Ricatti-Bessel and Ricatti-Hankel functions on the diagonal. (R j being the radius of the circumscribing sphere of the jth scatterer).

The balanced recursive algorithm is that the solution for the T (N,N ) N matrix is obtained from the T matrices of the N -1 system, T (j,k)

N -1 , via the matrix inversion in eq.( 20). All the other matrices T (j,k) N with j = N or k = N are then obtained via matrix multiplications and additions via equations ( 22) and ( 24). This process is then repeated as many times as desired.

B. Physical quantities

When the incident field is a plane wave, it is convenient to express physical quantities in terms of cross sections. Appealing to the far-field approximation of the field, the extinction and scattering cross sections of clusters of N objects can be respectively expressed 8,9

σ ext = - 1 k 2 Re N k=1 a (j), † f (j) N and σ scat = 1 k 2 N j,k=1 f (j), † N J (j,k) f (k) N (29)
It is also possible to produce analytical expression for local field quantities like individual absorption cross sections. For lossy scatterers in a lossless host medium, one can obtain individual particle absorption cross sections by integrating the Poynting vector on a circumscribing sphere surrounding the particle to obtain the formula as

σ (j) a = - 1 k 2 Re f (j), † N e (j) N - 1 k 2 f (j) N 2 (30)
In an analogous fashion, optical forces on the particles can be calculated by integrating the Maxwell tensor on a circumscribing sphere surrounding the particle [START_REF] Moine | Optical force calculations in arbitrary beams by use of the vector addition theorem[END_REF][START_REF] Mishchenko | Scattering, Absorption and Emission of Light by Small Particles[END_REF] . It is frequently convenient to characterize the optical force by vector 'cross sections', -→ σ opt , defined such that the time averaged optical force on particles immersed in a liquid dielectric of refraction index n med can be expressed as

F opt = S inc n med c -→ σ opt ( 31
)
where S inc = 1 2 Re{E * inc × H inc } is the incident irradiance. The binding force and its associated cross section, σ b , between two particles separated by a relative position vector r pos ≡ r 2r 1 , can be defined as

F b ≡ 1 2 (F 2 -F 1 ) • r pos ≡ S inc n med c σ b (32)

C. Interacting localized plasmon excitations

For those conductors, such as the noble metals, that support surface plasmon resonances, one can usually observe localized plasmon resonances in sufficiently small particles. These resonances are typically dominated by absorption if the particles are sufficiently small with respect to the incident wavelength and by scattering for larger particles. We chose to study silver spheres 50 nm in diameter immersed in air (for which both scattering and absorption are non-negligible).

The study is carried out for wavelengths ranging from the near ultra-violet through the visible (300 to 850 nm). We ignore the relatively modest finite size corrections to damping [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF] and simply adopt the bulk dielectric constant of silver from ref. [START_REF] Gray | American Institute of Physics Handbook[END_REF] and extrapolate between the experimental values. The extinction, scattering and absorption cross sections for these particles are readily obtained from Mie theory and are displayed in figure 3 as a function of frequency. These spheres are quite small with respect to visible wavelengths, (size parameters in the 300↔800 nm wavelength range go through kR = 0.52 ↔ 0.20) and the isolated particle cross sections are obtained to high precision with n max = 4. One can also see from figure 3 that the strength of the plasmon resonance for these particles is about half due to absorption and about half due to scattering.

One of the principal sources of interest of the localized plasmon resonances is their capacity to produce large field enhancements in regions much smaller than the incident field wavelength. This property is demonstrated in fig. 4a) with a 2D and 1D plot of the electric field intensity in and near an isolated 50 nm diameter silver sphere illuminated near its resonance peak (λ 0 = 365 nm with N Ag = 0.077 + 1.6i). The plots in Fig. 4 are performed in a plane containing the center of the sphere and perpendicular to k inc (the polarization lies along the horizontal axis). The dimensionless extinction and scattering 'efficiencies',

Q = σ/(πR
λ 0 (nm) σ / ( π R 2 ) Q ext Q abs Q scat
Fig. 3. Total cross section 'efficiencies', Q ≡ σ/(πR 2 ) for an isolated 50 nm diameter sphere.

the T matrix calculated by RCTMA contains information for arbitrary incident fields, we study the physically interesting case of a plane wave perpendicular to the axis separating the particles. As is widely known, the response then depends strongly on the polarization of the incident light. In figures 5a) and 5c), the extinction and scattering cross sections per particle are presented when the polarization is respectively perpendicular and parallel to the symmetry axis. From figure 5c), one sees that the cross section for the parallel to axis polarization presents a two sphere coupled resonance that is strongly red-shifted with respect to the isolated particle resonance. The polarization perpendicular to this axis on the other hand presents only slight modifications with respect to an isolated sphere.

The optical binding force cross sections for these same polarizations are respectively plotted in figures 5b) and 5d). While the binding force for the polarization perpendicular to the particle axis (cf.5c)) is slightly repulsive, the force for polarization parallel to the resonance can be highly attractive with the dimensionless |Q b | attaining amplitudes of three orders of magnitude. There has already been experimental and theoretical evidence supporting the existence of optical force couplings in particles with plasmon excitations 14 although such high precision calculations at such small separations seems not to have been presented before now.

This dimer system dramatically illustrates the 'strong' coupling category since correct calculations require that the VSWF space be enlarged far beyond the predominantly dipolar response characterizing the particles in isolation. The normalized cross sections per particle are given in the table 1 for different values of the VSWF space truncation. Although it was necessary to go to ∼ n max = 30 to achieve 4 digit precision in all the cross sections, the table indicates that results were already quite good at n max = 20.

A base 10 logarithmic intensity field map for a two silver sphere dimer illuminated with light polarized along the symmetry axis (frequency near the coupled sphere resonance maximum (λ 0 = 467 nm and N Ag = 0.048 + 2.827i) is presented in figs.6a) and figs.6b) which are respectively a 2D plot (in the same plane as figure 4) and a 1D plot along the symmetry axis.

The size parameter of the individual spheres is kR = 0.34 and the isolated cross sections at this frequency are Q ext = 0.136 and Q scat = 0.0963. As can be seen in fig. 6, the fact that one had to go so far beyond the dipolar response has a dramatic effect on the field inside and near the the particles. Notably, the fields inside the particles are no longer quasi-constant as was the case for isolated particles.

Q ext = σ ext /(2πR 2 ), Q scat = σ scat /(2πR 2 ) and Q b = σ b /(πR 2 ).
The system is a dimer composed of D = 50 nm diameter silver spheres (1 nm separation) at (λ 0 = 467 nm and N Ag = 0.048 + 2.827i)

An important word of caution should be made at this point. Although 1 nm separation may appear to 'nearly' touching, the coupled resonance is in fact quite sensitive to exact separation details when resonant particles are so closely separated. For example, at a separation distance of 0.5 nm for the silver dimer, the coupled plasmon resonance is displaced to λ 0 ≃ 516 nm as compared with λ 0 ≃ 467 nm for a 1 nm separation, and the multipole order has to be pushed to n max ≃ 50 to achieve four digit accuracy in the cross sections.

Nanometer scale separations are not necessarily theoretical idealizations however as recent experiments with DNA separators have demonstrated [START_REF] Bidault | Plasmon-Based Nanolenses Assembled on a Well-Defined DNA Template[END_REF] . Nevertheless, in applications like DNA separators, one may well have to consider the strong optical forces between these particles on account of the exceptionally strong attractive optical forces efficiencies of these resonances. For instance, the binding force efficiency at 0.5 mathrmnm separation was calculated at Q b = -20174 for λ 0 ≃ 515.6 nm (cf. Q b = -6018 for 1 nm separation at λ 0 ≃ 467 nm). The question of perfect spheres exactly in contact however seems untenable from an experimental standpoint and quite difficult from a theoretical standpoint on account of the singular behavior of the contact point. Theoretical 'separations' of 1 Å for instance require multipole truncations of the order of n max 120 before convergence is achieved, but the idea of 'perfect' spheres separated by atomic scales has clearly gone beyond domain of applicability of our mesoscopique physical model in any case.

It is also important to verify that the recursive algorithm works for more complicated systems. Towards this end, we illustrate in figure 7, the results of calculations for a system composed a line of 5 identical silver spheres spheres separated by 1 nm. For the binding force, we now present Q b,1 which is the binding optical force between outermost spheres and their nearest neighbor and Q b,2 which is the binding force between central sphere and each of its nearest neighbors. It is interesting to remark that addition of other spheres in the chain dramatically lessens the strong binding force interactions between spheres even though the fields between the spheres (cf. fig. 5) can still be almost as high as the dimer case.

We remark that the interactions have continued to red-shift and widen the coupled "chain" resonance. This chain resonance peaks at ≈ 561 nm and N Ag ≃ 0.0564 + 3.685i (Q ext /5 = 14.416 and Q scat /5 = 12.543). It is clear from figure 7c) that the extinction cross section of the chain resonance is increasingly dominated by scattering rather than absorption. The number of VSWF orders necessary for high precision was also seen to decrease slightly for the chain. The calculation of fig. 7 was carried out with n max = 20 since calculations at n max = 24 produced negligible differences on this scale.

Despite the dominance of scattering, a considerable amount of absorption is still present in the 5 sphere chain. Furthermore, from the field maps in figures 6c) and 6d), one can see that the central sphere has the highest field internal field intensities, and one consequently expects increased absorption in the central sphere. This supposition can readily be confirmed quantitatively by using eq.( 30) to calculate the absorption in each individual sphere. The results are given in table 2.

Q a,1 Q a,2 Q a,3 Q a,4 Q a,5
0.8346 2.333 3.030 2.333 0.8346 Table 2. Individual absorption efficiencies Q a,j ≡ σ a,j /(πR 2 ) in a five sphere chain at λ 0 = 561 nm and N Ag = 0.0564 + 3.685i

We conclude this section with some calculations systems concerning larger chains of particles. One can remark that the chain coupled resonance continued to red-shift and widen when passing from the dimer to the five particle chain. Results for the extinction and scattering cross sections chains of 10 and 20 sphere chains are presented in figure 8 for the same polarizations and incident directions as considered previously.

One readily sees that ultraviolet and perpendicular responses per particle seem to have stabilized for large chains. The collective chain response on the other hand continued to broaden and slightly red shift as one passed from 10 to 20 sphere chains and it is an interesting point for future studies to examine the evolution of this phenomenon for even longer chains and to study the impact of defaults in the chains. In a) the polarization is perpendicular to the symmetry axis and in b) it is parallel to the symmetry axis.

Conclusions

The balanced recursive algorithm can give useful and highly accurate information in systems with large numbers of strongly interacting resonances. This has been demonstrated herein for the case of localized plasmon resonances and the studies presented here suggest that chains of closely spaced localized plasmons can have potentially interesting applications with respect to frequency shifting and broadening. Although not demonstrated here, this technique also proves useful for treating closely spaced systems possessing surface resonances of 'whispering gallery' type.

It is worth remarking that matrix balancing seems to be a useful method to employ in almost any Foldy-Lax equation solution scheme, be that for direct system matrix inversion, iterative techniques or linear system solutions. In fact, some modern matrix inversion programs actually integrate numerical matrix balancing into their algorithms. Nevertheless, since the matrix balancing in Foldy-Lax equations can be obtained analytically at relatively

Fig. 1 .

 1 Fig. 1. Schematic of an field incident on a collection of scatterers centered on x 1 , x 2 , ..., x N . The radii of the respective circumscribing spheres are denoted R 1 , R 2 , ..., R N .

Fig. 2 .

 2 Fig. 2. Plot of the spherical Bessel to Hankel function ratio, ψ n (kR) /ξ n (kR) occurring in the Mie coefficients when kR = 10.

  derive a matrix balanced version of the recursive T matrix algorithm.

Fig. 4 .

 4 Fig. 4. Electric field intensity ||E t || 2 /||E inc || 2 in an isolated 50 nm diameter sphere (λ 0 = 365 nm, N Ag = 0.077 + 1.6i). In a) is presented a 2D (hot) plot of the electric field intensity in a plane perpendicular to the wavevector and containing the origin of the sphere (the horizontal axis lies along the polarization direction). Fig b) is a 1D plot of the field intensity along the line in this plane containing the direction of electric field polarization.

Fig. 5 .

 5 Fig.5. Dimensionless cross section 'efficiencies' per particle, Q = σ/(2πR 2 ) and binding force 'efficiencies' for a dimer of 50 nm diameter spheres (1 nm separation). In a) and b) the polarization is perpendicular to the symmetry axis and in c) and d) it is parallel to the symmetry axis.

  n

Fig. 6 .

 6 Fig.6. Logarithmic scale plots of the field intensity for a two sphere dimer with (λ 0 = 467 nm and incident light polarized along the sphere axis (N Ag = 0.048 + 2.827i). In a) is a 2D plot in the plane containing the centers of the spheres and the polarization vector while b) is a 1D Logarithmic plot along the symmetry axis of the spheres. Figures. c) and d) are the same as a) and b) respectively but for a 5 sphere chain of spheres at its resonance maximum (λ 0 = 561 nm and N Ag = 0.0564 + 3.685i). (cf. figure7).

1 Q b, 2 Fig. 7 .

 127 Fig.7. Total cross section and binding 'efficiencies' for a line of 5 'touching' silver spheres 50 nm in diameter (1 nm separation). In a) and b) the polarization is perpendicular to the symmetry axis and in c) and d) it is parallel to the symmetry axis.

Fig. 8 .

 8 Fig. 8. Total extinction and scattering cross section efficiencies per particle in chains of 10 and 20 particles.

Table 1 .

 1 Dimensionless cross section 'efficiencies' per particle in function of the VSWF truncation, n max .

	max	5	10	15	20	25	30	35	40
	Q ext /2	4.60	15.53	17.38	17.20	17.14	17.13	17.13	17.13
	Q scat /2	3.51	10.62	11.30	11.04	10.98	10.97	10.97	10.97
	Q b	-417	-3639	-5530	-5918	-6000	-6015	-6018	-6018

low computational cost, it seems beneficial to carry out this balancing explicitly rather than relying on purely numerical treatments.

The matrix balanced RCTMA has potentially interesting applications for other kinds of resonance phenomenon, notably whispering Gallery modes. Such studies are currently underway. Furthermore, the ability of the matrix balanced RCTMA to study defaults and small modifications in large complicated systems is particularly promising and will be employed in subsequent studies.

Brian Stout and Alexis Devilez would like to thank Ross McPhedran, Evgeny Popov

and Nicolas Bonod for helpful discussions. This work was funded in part by the grant ANR-07-PNANO-006-03 ANTARES of the French National Research Agency.

Appendix A: Vector spherical wave functions

The vector spherical wave functions can be readily written in terms of the Vector spherical harmonics (VSHs) and outgoing spherical Hankel functions :

In the same manner, the regular VSWFs are obtained by replacing the spherical Hankel functions in eq.(A1) by spherical Bessel functions. Our adopted definition of the VSHs is

where the Y nm (θ, φ) are the scalar spherical harmonics.