Franc ¸ois Chériaux

Laurence Picci

Julien Provost

Jean-Marc Faure

Conformance test of logic controllers of critical systems from industrial specifications

This paper presents a synthesis of works performed in the frame of an industry/academia cooperative research. The overall objective of this research is to automate the construction of test sequences for conformance test of industrial logic controllers when the expected behavior is described in industrial specification languages. The first contribution is aiming at preventing from combinatorial explosion by preliminary verifications on the implementation so as to check whether it satisfies some structural properties. The second contribution is a method to translate an industrial specification into a formal model to take benefit from theoretical results on conformance test of formal models for discrete event systems description. These two contributions are exemplified on case studies from the domain of energy production.

INTRODUCTION

Conformance test is aimed at checking whether an implementation, seen as a black-box with inputs/outputs, behaves correctly with respect to its specification. This kind of test is mandatory for the components of the control system of critical processes, like power production and distribution or transport, and particularly for the software-based logic controllers that are more and more often used in these systems, even to implement safety-related functions.

Conformance test of a logic controller is a noninvasive test that is performed (Figure 1) before operation and consists in sending to the controller an input sequence and comparing the observed output sequence, controller's response to the input sequence, to the expected output sequence so as to build a test verdict (the implemented controller conforms to its specification or not). The set of the input sequence and expected output sequence is termed test sequence.

Numerous theoretical results have been published in the domain of conformance test of logic systems, assuming that the specification is formally described, for instance in the form of a finite state machine [START_REF] Lee | Principles and methods of testing finite state machines -a survey[END_REF], a transition system [START_REF] Tretmans | Model based testing with labelled transition systems[END_REF] or, more recently, a particular class of Petri net [START_REF] Von Bochmann | Testing k-safe petri nets[END_REF]. Generally speaking, these results provide a way to build automatically Unfortunately, these results cannot be used directly for conformance test of industrial logic controllers, because, in industrial practice, the specifications of the behavior of the controllers are not given in formal languages but in tailor-made, often standardized, industrial specification languages, like Logic Functional Diagrams (LFD), Grafcet or state-charts. Test sequences are then built manually, what is a very tedious, time-consuming and error-prone task. Moreover, applying theoretical methods for conformance test sequences generation to industrial-sized examples may easily lead to combinatorial explosion. This synthesis paper presents the results of a work performed in the frame of an industry/academia cooperative research 1 , to contribute to solve these two issues. Section 2 describes shortly the main features of the two industrial specification languages that were selected in this study and underlines their complementary. Section 3 shows how the size of the test sequence can be reduced by verifying whether the implementation satisfies some structural properties, while section 4 presents a method to translate an industrial specification into a formal model, without semantics loss. These two contributions are exemplified on industrial case studies coming from the sector of energy production. Perspectives for further works are given in the last section.

INDUSTRIAL SPECIFICATION LANGUAGES

FOR LOGIC CONTROLLERS Specifications of the logic control of industrial systems are not expressed in formal languages but in languages that are tailored-made for automation and often (officially or de facto) standardized. If it is indeed possible to describe the behavior of small-sized control systems with languages issued from the theory of discrete event systems, like those above-mentioned, this is no more the case when industrial-sized systems are considered. Modeling the expected behavior of the logic control of a real industrial system with these languages leads easily to a huge, complex and not readable and understandable model that cannot be used as a specification. On the opposite, industrial specification languages allow to capture the know-how of engineers and to express clearly what they expect. In the case of energy production, two graphical specification languages: LFD and Grafcet, are widely used and have been selected for this study. It matters to note that these languages are also used in other industrial domains that deal with critical processes, for instance railway transport, chemistry, environment. Hence, the results of this study are not limited to a particular industrial domain.

These two languages are complementary because they are used to specify different kinds of functions, as pointed out in sub-section 2.1. Last, it must be underlined that LFD and Grafcet are functional diagrams that do not depend on implementation technology, i.e. a given diagram may be implemented in different ways, for instance in different programming languages or even using hardware solutions.

Control system architecture

The control system of a power plant is a hierarchical system that comprises four levels (Figure 2): planning and maintenance management (L3), monitoring 1 This research was funded by the French Research Agency (TESTEC project). and operation (L2), automation and control (L1), data acquisition and actuators (L0).

Figure 2: Architecture of the control system Within each level, analogical and logic control functions are implemented. Experience has shown that the analogical functions are less critical than the logic ones. Hence, the parts of the control architecture that implement analogical functions are validated only by simulation, using a plant simulator, while conformance test of parts which perform logic control functions is compulsory. This explains why this work focuses only on logic controllers.

Moreover, only level L1 is considered in this study. It is closer to the plant, owns shorter response time; hence it impacts safety and dependability more strongly than the upper two levels. The expected behavior of the logic control functions at this level is specified either in LFD or in Grafcet, two languages that are briefly presented below, according to the following rules:

• basic protection functions are described in LFD,

• tasks control functions are described in Grafcet.

2.2 Brief description of the LFD language LFD is a graphical specification language that is well suited to specify the behavior of elementary reactive functions, such as actuators control functions in accordance with sensors values. This explains why this language is used to specify basic protection functions.

A LFD (Figure 3) is a network of interconnected logic functions that may include combinatorial (NOT, AND, OR), sequential (set-and reset-dominant (SR and RS) memories) and timed functions (on-delay (TON) and off-delay (TOFF) timers). A LFD specifies the expected relations between its inputs (sensors values or commands from the upper levels) and outputs (commands to actuators).

The example of figure 3 will be used in the remainder of this paper to illustrate the contribution on test sequence size reduction. Grafcet is a standardized graphical specification language (IEC 60848 2002) to describe the behavior of logic sequential systems. Grafcet was developed from the results of the Petri Nets community, and particularly from those on Interpreted Petri Nets. This language was first standardized in France at the beginning of the 80s, and at the international level in 1988. Since this date, several extensions have been proposed to enhance the modeling possibilities; they are included in the last version of the standard (IEC 60848 2002). A good scientific presentation of the main features of the current version of the Grafcet standard can be found in [START_REF] Guéguen | Extensions of grafcet to structure behavioural specifications[END_REF].

This language is widely used in several industrial domains because it allows engineers to describe easily tasks sequences, tasks selection, parallelism and synchronization. A Grafcet is composed of steps, represented by squares, and transitions, represented by horizontal lines; a step can be linked only to transitions and a transition linked only to steps. Actions may be associated to a step; an action associated to a step is performed only when this step is active and then acts upon an output variable. A transition condi-tion must be associated to each transition; this condition is a Boolean expression which may include input variables and steps activity variables.

The example of figure 4 will be used in the remainder of this paper to illustrate the contribution on translation of an industrial specification into a formal model.

BUILDING TEST SEQUENCES FROM LFD

The overall objective of the work described in this section is to avoid (or limit) combinatorial explosion during test sequence construction and to develop a simple, while efficient, method that allows conformance test of a logic controller with a size-reduced test sequence, while providing trustworthy test results. The principle of this method is to verify structural properties of the implementation before the test. This work focuses on conformance test of an implementation that performs highly critical, but simple, logic control functions, e.g. actuators start and stop, protection of actuators, alarms generation, and is specified in LFD.

Method overview

According to the standard on safety-related instrumentation and control systems for nuclear power plants (IEC 60880 2006), system validation is the confirmation by examination and provision of relevant evidences that the system fulfills its specification. Conformance test is clearly a validation technique for control systems and it could be thought that an exhaustive test sequence is to be constructed to provide relevant evidences, when critical systems are considered. Nevertheless, attempting to build an exhaustive test sequence may easily lead to combinatorial explosion; for a LFD with n i logic inputs and n m internal memories, this sequence can include up to 2 n test steps, with n = n i + n m , for instance. This explains why reduction of the test sequence size must be investigated.

The method that was developed to construct automatically a size-reduced conformance test sequence from a LFD comprises two main phases:

1. Verification of the following six structural properties:

a) The specification and implementation own the same memories and timers.

b) Each output of a memory is also an output of the system2 . c) There is no 'backwards loop', i.e. an output of a function must not be connected to an in- put of an upstream function or be input of this same function. For instance, the output of the SR function in figure 3 cannot be connected to an input of the combinatorial functions that are located upstream of this sequential function. d) The dependency relations that link the outputs of the system to its inputs and the stored variables (outputs of memories) in the implementation are identical to those of the specification. e) The dependency relations that link the inputs of the timers and memories to the inputs of the system in the implementation are identical to those of the specification. f) The dependency relations that link the outputs of the system to its inputs and the timers outputs in the implementation are identical to those of the specification.

2. Construction of the test sequence by using the results of the first phase. The core idea of this phase is to take benefit from the knowledge on the structure of the implementation to generate the only test steps that are necessary to test its behavior. The dependency relations that have been stated point out indeed how to control and observe the memories and timers states; among the whole set of test steps that can be obtained by simulation of the specification, the relevant test steps are those which satisfy these two conditions.

It must be noted that the first phase permits to avoid a conformance test be performed on an implementation whose structure does not satisfy some structural properties, what reduces implementation validation duration. This method has been implemented in a prototype tool [START_REF] Salaün | Prospects for model-based testing of discrete safety systems[END_REF][START_REF] Chériaux | Functional test of control systems ensuring a high coverage rate[END_REF]. The results obtained with this tool will be illustrated in the next sub-section on the basis of the example of figure 3. Although limited to logic systems, this approach corresponds to an innovative prospect for model-based conformance test, taking benefit from structural properties verification to reduce test complexity.

3.2 Verification of the structural properties of the implementation The LFD specification of figure 3 includes one SR memory, 2 TON (set to 5 and 10 seconds) and one 3 seconds TOFF timers; the output of the memory is connected to output AAD and there is no backwards loop. The following relations between the three outputs and the nine inputs can be easily obtained by dependency analysis of the diagram:

• AAD = SR(AND(OR(TON 5S(VIBR),TEMP), NOT(VER)),AND(ROP,VAL))

• HOR = TOFF 3S(SR(AND(OR( TON 5S(VIBR),TEMP),NOT(VER)), AND(ROP,VAL)))

• S0 = TON 10S(AND(OR(A1,B1,C1),NOT(D1)))

From these global relations, class d, e and f properties can easily be obtained. For instance, the following properties can be stated from the first relation:

• AAD depends on SR Q,

• TON 5S I depends only on VIBR,

• AAD depends on TON 5S O, TEMP, VER, ROP and VAL, where SR Q, TON 5S I and TON 5S O represent respectively the output of the memory and the input and output of the 5 s TON timer. Once the structural properties stated, they are checked on the implementation. Then, the prototype tool developed during this work includes a configurable parser to check structural properties for logic controllers developed in the following languages: two languages of the IEC 61131-3 standard for programming of programmable logic controllers (Ladder Diagram (LD) and Function Block Diagram (FBD)), C and VHDL). If verification shows that the implementation satisfies the structural properties that have been stated from the specification, conformance test can be performed with the test sequence whose construction is explained in the following sub-section. This approach has been assessed successfully with several examples that describe real industrial LFD specifications. The main limitation is that the LFD must not include 'backwards loop'. However, this limitation is not really significant because industrial experience has shown that basic protection systems behavior can always be specified with LFD without backwards loop.

Test sequences construction

This construction takes benefit from the results of I/O dependency analysis and the properties that have been stated in the first phase. For room reasons, it is not possible to present all the details of the construction; only the two main principles will be given.

Dependency relations between outputs and inputs permit to pinpoint the inputs that are used to compute a given output; then, it is possible to construct a test sequence that contains only values of these inputs and expected values of the considered output. For the example, dependency relations show that the first two outputs only depend on the same five inputs and the third output on the remaining four inputs. Hence, the test sequence will comprise two parts; the first one will be composed of values of the only first five inputs (VIBR to VAL) and the corresponding expected values of the first two outputs, while the second part will contain values of the only last four inputs (A1 to D1) and the corresponding expected values of the last output. This splitting principle may shorten strongly test sequences.

Moreover timers behavior can be tested more easily when the relations between their inputs-outputs and the inputs-outputs of the system are known; this For this specification (nine inputs, one memory and three timers), a test sequence constructed without taking benefit from the structural properties analysis includes more than thousand test steps and is executed in about 5 hours. When the results of this analysis are introduced, less than hundred test steps are necessary and the test is performed in about 1 minute; thus, in this case, the test sequence length is reduced by one order of magnitude and the test duration by more than two orders.

Discussion

It has been shown in this section that a size-reduced test sequence can be used for conformance test from a LFD specification that does not contain any backwards loop, provided that some structural properties of the implementation have been checked previously; these properties are obtained by analysis of the specification model.

This approach combines verification techniques and conformance test. Moreover, when the verification results are negative, non-conformance of the implementation with regard to the specification can be stated without test, what is time-saving. These results have been implemented in a prototype software tool; experiments with several diagrams have shown that the overall time to build and verify the structural properties complies with the requirements of industrial processes of design and implementation of control systems.

It has been shown in this section that a size-reduced test sequence can be used for conformance test from a LFD specification that does not contain any back-wards loop, provided that it has been checked previously that the I/O dependency relations of the implementation are identical to those of the specification. This approach combines verification of I/O dependency relations and conformance test. Moreover, when the result of this verification is negative, nonconformance of the implementation can be stated without test, what is time-saving.

BUILDING TEST SEQUENCES FROM

GRAFCET SPECIFICATIONS The aim of this section is to show how a test sequence can be automatically built from a Grafcet specification by translating this specification into a formal model, without semantics loss. This work considers only non timed models because the first concern of engineers during conformance test is functional correctness; conformance test for time correctness is a second concern, once functional correctness is ensured, and motivates further research. Mealy machines have been chosen as the target formalism because it is suited to non timed logic systems with inputs/outputs and that numerous results have been previously obtained in the field of conformance test of this kind of finite state machine. Hence, it will be possible to take benefit from these results, once the Mealy machine that is equivalent to the initial Grafcet obtained,

Translation of a Grafcet model into an equivalent Mealy machine

Method overview

The translation method of a Grafcet into an equivalent Mealy machine, without semantics loss, comprises two phases (Figure 5) that are detailed the next two sub-sections:

1. Construction of the automaton, termed Stable Location Automaton (SLA), that represents formally all stable states of the logic system described by the Grafcet as well as all evolutions between these states;

2. Translation of this automaton into an equivalent Mealy machine.

Construction of the SLA

The reader is reminded that several steps may be simultaneously active in a Grafcet, for instance because they are located in two parallel sequences, like steps 1 and 4, or 2 and 7 in figure 4. Moreover, according to the evolution rule #4 of the Grafcet standard, A location is characterized by a set of simultaneously active steps, termed situation, and the set of outputs that are emitted when this situation is active. For the example of figure 4, location ({A3, 2, 5}, {IL -0, Act1, P rod1}) means that steps A3, 2 and 5, and only these steps, are active and outputs IL -0, Act1 and P rod1, and only these outputs, emitted. It might be argued that it is not necessary, for this example, to give the list of emitted outputs, once the list of active steps given. However, different sets of outputs may be emitted for the same situation, when the Grafcet contains stored actions or continuous actions that depend on inputs values (conditional actions). Moreover, only the outputs of the logic controller whose behaviour is specified by the Grafcet are observable during conformance test. These two reasons justify the location definition.

Regarding the evolutions in response to inputs events, the Grafcet standard introduces the concepts of transient and non-transient evolutions. A transient evolution is characterized by the firing of several successive sets of transitions on the occurrence of a single input event; on the opposite, a non-transient evolution is one and only one firing of a set of simultaneous transitions when an input event occurs. For the example of figure 4, the evolution from (A3, 2, 5) to (A3, 2, 7) is non-transient (only the transitions t8 is fired) while the evolution from (A3, 2, 5) to (A1, 0) is transient (the two transitions t6 and t8 are simultaneously fired, then transition t11 is fired, then transition t12 is fired, and finally transition t3 is fired. During a transient evolution, several situations are crossed from the source to the target situation, ((A3, 3, 7), (A3, 8), (A3, 0)) in the above example). The standard pinpoints that the continuous actions associated to these 'transient' situations are not executed and consequently the associated outputs are not emitted. As only the outputs of the implementation are observable during conformance test, the 'transient' situations are useless in this work and only the source and target situation of any evolution, named stable situations, must be kept.

Hence, the SLA obtained from a Grafcet is a finite state machine that contains:

• All stable locations, i.e. the locations defined from the stable situations;

• All evolutions from stable location to stable location. Each evolution is caused by a change of the inputs values and is labeled with an evolution condition, Boolean expression on the inputs values.

A software tool was developed during this work to construct automatically a test sequence from a Grafcet model; the first treatment performed by this tool is the automatic construction of the SLA. For the example of figure 4, this automaton comprises 16 stable locations and 102 evolutions.

Construction of the equivalent Mealy machine

A Mealy machine is a finite state machine that is defined on input and output events alphabets. Every transition of this machine is labeled with a couple of input and output event. Then, the issue to solve, to translate a SLA into a Mealy machine, is to translate a state machine whose evolution conditions are defined by Boolean expressions into an event-based state machine, without any semantics loss. This can be achieved by first defining the elements of the input and output alphabets as the different combinations of the logic inputs and outputs of the Grafcet3 . For the example of figure 4, the input alphabet of the equivalent Mealy machine contains 2 8 = 256 elements, because there are 8 logic inputs, and the output alphabet 2 10 = 1024 elements, because the Grafcet model owns 10 outputs. As a SLA is already a state machine, each location of the SLA gives rise to a state of the Mealy machine. Transition and output functions, which define respectively the transitions between the states and the elements of the output alphabet associated to the transitions, can then be computed from the SLA evolution conditions. For room reasons, it is not possible to give a detailed definition of these two functions; the interested reader is referred to [START_REF] Provost | Test sequence construction from SFC specification[END_REF]. The size of the equivalent Mealy machine is easily computable; if n states and n transitions represent respectively the number of states and transitions of this machine, and n locations and n input-var the number of locations and input variables of the SLA, it comes:

n states = n locations n transitions = n locations • 2 n input-var (1)
It must be noted that the number of evolutions of the SLA has no influence on the size of the Mealy machine.

The Grafcet model of figure 4 can then be translated into an equivalent Mealy machine that comprises 16 states and 4, 096 = 16 • 2 8 transitions.

Test sequence construction

Once the equivalent Mealy machine built, a test sequence for conformance test of a logic controller that is supposed to implement the initial Grafcet can then be obtained by using one of the methods surveyed in [START_REF] Lee | Principles and methods of testing finite state machines -a survey[END_REF]. This method will have to provide an exhaustive and minimum-length test sequence. When control of critical systems is considered, as this is the case in this work, the test sequence must be indeed exhaustive, i.e. all evolutions, for all combinations of the logic inputs, of the specification must be tested to avoid some expected behaviors be missing in the implementation. If the specification is formally represented with a Mealy machine, an exhaustive test sequence crosses every transition at least once. Last, a minimum-length test sequence is searched to lessen conformance test duration.

The method to obtain the minimum-length, while exhaustive, test sequence of a Mealy machine is called transition tour method [START_REF] Naito | Fault detection for sequential machines by transitions tours[END_REF] and is a particular solution, for a graph that represents the structure of a Mealy machine (states and transitions between states), of a well-known problem in graph theory: the Chinese Postman Problem [START_REF] Mei-Ko | Graphic programming using odd or even points[END_REF]. The general formulation of this problem is the following: "Find a minimum-length closed path that visits each edge in the graph at least once". As the graph which describes the structure of a Mealy machine is directed, but not weighted, the problem is simplified.

Using the software tool developed during this work and that integrates the transition-tour method, a 6, 436 steps test sequence was generated from the example of figure 4. If a step of this sequence is done every 20 ms during conformance test, the whole test execution will last approximately 2 minutes and 10 seconds, what is quite reasonable.

Discussion

This section has shown that an exhaustive conformance test sequence can be automatically generated from a non-timed Grafcet specification. The main contribution of this work is a method to translate automatically the industrial specification model into a formal model, without semantics loss. The principle of this method is to generate the whole state space of the specification and can be applied to other industrial languages for behavior specification.

All results were implemented in a prototype software tool. Experiments with several Grafcet models shown that the overall time to construct the test sequence from the initial Grafcet, including the translation from Grafcet to the equivalent Mealy machine is very short (some seconds at most), what guarantees control engineers acceptance.

CONCLUSION AND PERSPECTIVES

Safety and dependability of critical systems depend more and more on the correct behavior of logic controllers. This explains the strong demand from industry for conformance test techniques that can be used during the development of industrial control systems; to meet this objective, the test sequences must be built automatically from specifications in industrial specification languages. The two contributions that are presented in this paper are aiming to provide solutions to this issue. The first contribution focuses on specifications in LFD, a language for basic protection functions specification, and on reduction of the size of the test sequence by preliminary verifications of structural properties of the implementation. The second one considers specifications in Grafcet, a standardized language for more complex control functions for which the previous contribution cannot be applied, and on translation of the specification into a formal model so as to build an exhaustive and minimumlength test sequence. The test sequences that were obtained in both cases have been used successfully with real test-benches that integrated commercially available PLCs (Programmable Logic Controllers).

The perspectives of this research are both technical and scientific. First, integration of these results within an existing environment for controllers' specification, design and implementation is planned in a near future. In a longer term, the current limitations (no backwards loop in the LFD, only non-timed Grafcet) should be removed to increase the scope of the results and their interest for safety improvement.

Figure 1 :

 1 Figure 1: Conformance test of a logic controller

Figure 3 :

 3 Figure 3: Example of LFD specification

Figure 4 :

 4 Figure 4: Example of Grafcet specification

Figure 5 :

 5 Figure 5: From Grafcet to Mealy machine

Table 1 :

 1 Testing timers behavior is then controllable and observable and test patterns can be defined, as exemplified in table 1. The second test sequence must be performed after the first one because it assumes that the timer TON 5S, observed from the output AAD, behaves correctly.

	Test of TON 5S by observing the output AAD
	TEMP ROP VAL	VER VIBR	AAD	Remarks
	0	1	1	0	0	0	SR, then
							AAD reset
	0	0	0	0	0	0	TON 5S out-
							put remains
							at 0
	0	0	0	0	1	0 during	Timer
						5,000 ms,	behavior
						then 1	observation
	Test of TOF 3S by observing the output HOR
	TEMP ROP VAL	VER VIBR	HOR	Remarks
	1	0	0	0	0	1	SR set
	0	1	1	0	1	1 during	Timer
						3,000 ms,	behavior
						then 0 during	observation
						2,000 ms,
						then 1
	0	0	0	0	1	1	SR remains
							set
	0	1	1	0	0	1 during	Timer
						3,000 ms,	behavior
						then 0	observation

This restriction is not too strong because this configuration is very common for real LFDs; however, if this is not the case, an additional output must be added in a 'design to test' approach.

The inputs and outputs of the Grafcet and the corresponding SLA are identical.


Sophie Sakka Irccyn
email: sophie.sakka@irccyn.ec-nantes.fr

Chris Hayot
email: hayot@univ-poitiers.fr

Patrick Lacouture
email: lacouture@univ-poitiers.fr

A generalized 3D inverted pendulum model to represent human normal walking

This paper compares different inverted pendulum models to represent the stance phase of human normal walking. We have developed a model which takes into account the mechanism of the foot during the single support phase, by defining a pivot point under the ground level. Similarly to other models, the pivot point as well as the rod length remain constant during the complete single support phase. Lowering the position of the pivot point allows reducing the vertical amplitude of the center of mass (CoM) trajectory and therefore approaching the real CoM trajectory. Another big advantage of this model is that it can be easily adapted to pathological walking as the walking symmetry hypothesis -needed to build classical inverted pendulum models, is not used. We have measured the whole body kinematics of 5 healthy male subjects and set a reference CoM trajectory based on multibody modeling of the human body (16 segments). Then, we have determined a common mathematical definition of three inverted pendulum models: the classical IP-3D, the EIP-3D and our GIP-3D model. To insure continuity of the approached trajectories, a double support phase is represented by a simple pendulum which is a homothetic transformation of the inverted pendulum defining the previous single support phase, except for our GIP-3D model which may include asymmetric walking parameters.

I. INTRODUCTION

The calculation of the mechanical energy used for normal walking is based on the 3D trajectory and velocity characteristics of the center of mass (CoM) displacements. It looked then natural to propose evaluation models of normal walking focusing on the estimation of these trajectory and velocity characteristics, avoiding the use of heavy experimental setting needed by a multibody modeling of the human body. Two main approaches emerged: the inverted pendulum models and the major determinants of walking.

Indeed, locomotion can be approximate from the external forces point of view as the displacement of the center of mass (CoM) in 3D space. During normal walking, the CoM trajectory follows a cyclic sinusoidal pattern, and it seemed natural representing it with pendulum models: an inverted pendulum for the single support phases of the walking cycle, and a simple pendulum for the double support ones [START_REF] Hayot | Eip-3d: An extended 3d inverted pendulum model to represent normal human walking cycle[END_REF], [START_REF] Zijlstra | Displacement of the pelvis during human walking: experimental data and model predictions[END_REF]. When not pathological, walking is assumed perfectly regular and symmetric: one is able to predict the CoM trajectory knowing the characteristics of the subject anatomy and gait comfortable speed. Then the approximation of the CoM trajectory is available for healthy subjects walking at different speeds using a pendulum model from the spatial and temporal gait parameters [START_REF] Jordan | Walking speed influences on gait cycle variability[END_REF], [START_REF] Lelas | Predicting peak kinematic and kinetic parameters from gait speed[END_REF], [START_REF] Stoquart | Effect of speed on kinematic, kinetic, electromyographic and energetic reference values during treadmill walking[END_REF].

In this paper, we focused on the pendulum representation of normal walking. Our main motivation was to improve the existing models by getting closer to a more realistic 3D trajectory of the CoM while reducing the heaviness of the experimental setting needed by multibody models. Studying the dynamics experimental data: the reference CoM (calculated using multibody modeling) and the CoP (measured using force plates) trajectories during the swing phase of normal walking, we have observed a convergence of the straight lines joining these two points. We could then extract a pivot point for the pendulum representation based on a real mechanism of normal walking. Still, its position is not on the ground level as assumed by all the existing pendulum approaches, but about 38 cm below. This situation was already suggested by some authors, assuming that if the motion of the center of pressure is taken into account, it is equivalent to considering the pivot point of walking under the ground level [START_REF] Bertos | The vertical mechanical impedance of the locomotor system during human walking with applications in rehabilitation[END_REF], [START_REF] Gard | What determines the vertical displacement of the body during normal walking[END_REF], [START_REF] Mcgeer | Passive dynamic walking[END_REF]. To our knowledge, no further studies were realized to actually put this assessment in evidence. We have set a three dimensional pendulum model based on the observed pivot point and compared our results to the reference CoM trajectory and two approximate CoM trajectories using Zijlstra's IP-3D and Hayot's EIP-3D models. We first define the three pendulum models in section 2 using a common mathematical representation based on homogeneous transformation matrices. Then we expose our results in section 3 and discuss them in section 4.

II. PROTOCOL AND POST-TREATMENT OF THE DATA

A. Subjects, materials and experimental protocol

The experimentations were realized with 5 healthy male subjects aged XX±XX who gave their oral agreement to participate in this study. Their instructions consisted in walking at their comfortable speed (1,30 ± 0,22 m/s) following a 10-meter walk way. Average step length was 1,322 ± 0,041 m.

Each subject was equipped with 58 markers and performed three walking trials. 3d trajectories of the markers were captured at the frequency of 100 Hz using an optical motion capture system (Vicon, Oxford, UK). The ground reaction force and momentum were registered using two 3d force plates (Kistler 5233A and 9287B, Winterthur, Switzerland) at the frequency 900 Hz. The efforts data were synchronized with the kinematics data and the measurements were realized for a complete walking cycle.

B. Use and post treatment of the data

The human body was modeled by 16 rigid segments representing the feet, legs, thighs, hands, forearms, arms, pelvis, abdomen, thorax and head (including neck). Local frames were defined for each segment in agreement with the recommendations of the International Society of Biomechanics [START_REF] Wu | Isb recommendation on definitions of joint coordinate system of various joints for the reporting of joint motion -part i: ankle, hip, and spine[END_REF], [START_REF] Wu | Isb recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motionpart ii: shoulder, elbow, wrist and hand[END_REF], and bodies masses, CoM trajectory and moment of inertia were calculated using predictive equations and data from Zatsiorsky [START_REF] Zatsiorsky | Kinetics of Human Motion[END_REF]. the reference kinematics of the movement could be calculated based on the homogeneous matrix formulation of the equations of motion [START_REF] Legnani | A homogeneous matrix approach to 3d kinematics and dynamics. part 1: theory[END_REF]. The reference trajectory of the CoM was obtained from the weighted mean of the bodies CoM trajectories and the respective masses:

OC 1 R0 = 1 M 16 i=1 m i T i 0 O i C i 1 Ri (1)
OC denotes the absolute position of the CoM, m i is the mass of the i th body and M = 16 i=1 m i is the total mass of the subject. T k j denotes the homogeneous transformation matrix from R j to R k , R 0 being the reference frame.

Force and torque data were used to determine the CoP trajectory as well as to detect temporal characteristics of the walking cycle: contact instants at heal strike of the leading (LHS) and trailing (THS) limb, and at toe off (LTO and LTO, respectively).

The data extracted from all the trials were averaged after normalization with time and the reference static CoM height.

III. PENDULUM MODELS OF HUMAN WALKING

A. Pendulum models

Inverted pendulum models usually describe the single support phase of normal walking under the assumption that

• The gait is in steady state: there is no variation of CoM walking speed or trajectory amplitude during observation; • The walking is symmetric: similar motion is performed by the left and the right leg; • The system is conservative; Painlevé first integral is valid (conservation of the mechanical energy during walking) The double support phase is rarely represented, or more precisely is widely approached either by a straight line linking two single support phases in the case of the inverted pendulum representation [START_REF] Kuo | The six determinants of gait and the inverted pendulum analogy: A dynamic walking perspective[END_REF], [START_REF] Zijlstra | Displacement of the pelvis during human walking: experimental data and model predictions[END_REF], or just ignored if using the compass gait representation [START_REF] Saunders | The major determinants in normal and pathological gait[END_REF]. A direct consequence is that the CoM velocity is not continuous at the transition instants between single and double support phases. In 1997, Zijlstra et al. [START_REF] Zijlstra | Displacement of the pelvis during human walking: experimental data and model predictions[END_REF] introduced a continuous representation

P ′ 0 E ′ 0 E0 P(THS) E(THS) E(LTO) C ′ (TTO) C ′ (THS) P0
Fig. 1. Pendulum representation of normal walking, projection in the sagittal plane. The amplitude of the CoM trajectory is overestimated using [START_REF] Zijlstra | Displacement of the pelvis during human walking: experimental data and model predictions[END_REF]'s approach, underestimated using [START_REF] Hayot | Eip-3d: An extended 3d inverted pendulum model to represent normal human walking cycle[END_REF]'s approach.

of the CoM trajectory assimilating its oscillations in steady state to the ones of a vibratory system. In this model, the amplitude of the vibrations relies on the walking subject's morphological characteristics and the frequency is set according to the walking speed and the length of the steps. The pivot point of the pendulum, P 0 , is set in the middle of the trajectory of the center of pressure (CoP) under the foot as illustrated by the blue lines in Fig. 1. This model insures the continuity of the CoM trajectory over the gait cycle, but shows the main disadvantage to widely overestimate the amplitude of its trajectory. Dealing with this major drawback and under the assumption that the determinants of walking had to be taken into account in a realistic pendulum representation of normal walking, Hayot et al. [START_REF] Hayot | Eip-3d: An extended 3d inverted pendulum model to represent normal human walking cycle[END_REF] extended the model by introducing the mechanism of the foot during single support phase. The resulting model is represented in black lines in Fig. 1: similarly to Zijlstra, the length of the pendulum rod is measured at TTO (trailing toes off) instant and remains constant all over the single support phase, but the pivot point on the ground moves during the phase and is superposed with the CoP at each instant. Hayot et al. showed that the amplitude of the CoM trajectory got closer to the actual one by comparing their model with a reference trajectory extracted from multibody modeling of the tested subjects. Nevertheless, the approached trajectory underestimates, this time, the trajectory amplitude. Another difficulty of using this model is the motion of the pivot point on the ground, which requires the computation of the position of the CoP under the foot at each instant of the single support phase. Both Zijlstra's and Hayot's models rebuilt the double support phase using an homothetic transformation of the single support phase. The center of the homothety is the position of the approached CoM (C ′ ) at the transition instants between a single support phase and the next double support phase (THS: Trailing Heal Strike). The homothety ratio k 0 is the quotient of the distance performed by the CoM during double support phase (S d ) over the distance performed during the single support phase (S s ) (Fig. 1):

k 0 = S d S s . ( 2
)
Therefore knowing the characteristics of the inverted pendulum at TTO (the positions of the CoM and CoP) allows completely rebuilding the normal gait cycle of the subject.

B. Mathematical formulation

Let R 0 = (O, x 0 , y 0 , z 0 ) be the reference frame, R P k = (P k , x P k , y P k , z P k ) the frame attached to the considered inverted pendulum model during single support and (respectively E k = E ′ 0 ) according to the notations described in Fig. 1. C denotes the actual CoM of the observed subject, and C ′ is the approached CoM calculated by one of the pendulum models. At t =TTO, C ′ and C are considered identical which allows setting the length of the pendulum rod. Then the trajectory of C ′ in the reference frame is defined during the complete walking cycle as follows.

R E k = (E k , x E k , y E k , z E k )
OC ′ 1 R0 =          T P k 0 P k C ′ 1 RP k t ∈ [T T O, T HS] T E k 0 E k C ′ 1 RE k t ∈ [T HS, LT O]
(3) The matrices T i 0 are the homogeneous transformation matrices from the frame R i to the reference frame R 0 .

T i 0 =    
cos α j sin α j sin β j sin α j cos β j x i 0 cos β j -sin β j y i -sin α j cos α j sin β j cos α j cos β j z i 0 0 0 1 are defined in their respective frames as follows.

    (
P k C ′ = ℓ k s z P k E k C ′ = ℓ k d z E k (5)
where ℓ k j (j = {s, d}) defines the rod lengths which remain constant over their respective support phases and are such that 1) ℓ k s denotes the length between the pivot point and the CoM at t =TTO; 2) at t =THS, the three points P k , C ′ and E k are aligned and ℓ k d = k 0 ℓ k s .

TABLE I 3D PENDULUM MODELS CHARACTERISTICS DEFINING THE PIVOT POINTS AND THE ROD LENGTH FOR ZIJLSTRA'S CONTINUOUS MODEL (IP-3D), HAYOT'S EXTENDED MODEL TAKING THE MECHANISM OF THE FEET INTO ACCOUNT (EIP-3D) OR OUR GENERALIZED MODEL

(GIP-3D).

model k P k E k ℓ k s ℓ k d IP-3D P 0 E 0 ℓ s ℓ d EIP-3D P (t) E(t) ls ld GIP-3D P ′ 0 E ′ 0 ls ld
Concerning the EIP-3D model, the horizontal variations of the pivot point E during the double support phase are relative to the horizontal variations of the CoP during the single support phase:

δ(OE) = k 0 δ(OP) (6)
The three models we have described are 3D models and may represent the complete walking cycle of a subject by only measuring kinematics and dynamics of the normal walking during a single support phase. Still , the fact that the length of the pendulum rods remain constant all over the phases does not guaranty the continuity of the approached CoM trajectories at Leading Toe Off (LTO, Fig. 1) instant unless the measured support distance at TTO and LTO are exactly equal. In what follows, we will consider this hypothesis true.

IV. EXPERIMENTAL RESULTS

We have compared the different inverted pendulum models to the reference in terms of CoM approached trajectories. For instance, Fig. 3 shows the average CoM trajectories in the sagittal, frontal and horizontal planes for the IP-3D, EIP-3D and GIP described in Fig. 1, as well as the reference CoM trajectory obtained using multibody representation of the human body (the black line). The main differences between these four cases are described in the following section.

To introduce the GIP approximation, we have used the observation that the support lines joining the CoM and the CoP were converging to a small area located 38 cm under the ground level, as shown in Fig. 2. We have considered that the pivot point of the GIP could be defined using the boundary support lines at TTO and THS and detecting the middle point of the shortest distance in space between these two lines. For this affirmation, we have closely observed the relative movements of the support lines with time during the complete single support phase, extracting the instantaneous rotation vector from one time step to the next one, which is illustrated in Fig. 6, as well as the minimal distances between the support lines (Fig. 4 and5). We have calculated the weighted mean of the instantaneous pivot points which represent the crossing points between two support lines considered at two adjacent time steps. We have observed that it is very close to the GIP pivot point P ′ 0 presented in section II-A (Fig. 1). Fig. 4. The minimal distance between two adjacent support lines does not exceed 0.3 mm. Fig. 5. Position of the pivot points with time. The black lines represent the minimal distance between two adjacent support lines at each time step; the red line is the weighted mean and the blue line is the minimal distance between the first (TTO) and last (THS) support line of the single support phase. The distance between the pivot points calculated using the weighted mean and the first and last support mean is represented by the green line.

V. DISCUSSION

We have presented different inverted pendulum models adapted to normal walking under the hypotheses that walking is symmetrical and regular on flat floor. The classical continuous 3D inverted pendulum model (IP-3D), such as the one proposed by Zijlstra et al. [START_REF] Zijlstra | Displacement of the pelvis during human walking: experimental data and model predictions[END_REF], presents the main drawback to overestimate the amplitude of the CoM trajectory by representing it as an arc of circle which center is fixed on the ground level. An intermediary model between the inverted pendulum and the six determinants of gait representations of human walking was proposed by Hayot et al. [START_REF] Hayot | Eip-3d: An extended 3d inverted pendulum model to represent normal human walking cycle[END_REF]. The approached CoM trajectory gets closer to the real one as it is flattened compared to the IP-3D with the known consequences on energy and work estimations [START_REF] Kuo | The six determinants of gait and the inverted pendulum analogy: A dynamic walking perspective[END_REF]. Still, this model slightly underestimate the CoM amplitude over the gait cycle [START_REF] Hayot | Eip-3d: An extended 3d inverted pendulum model to represent normal human walking cycle[END_REF]. Moreover, it is based on a moving pivot point which requires the measurement of the CoP position at each time step. We have then proposed a generalized pendulum model after having observed that the support lines joining the CoM and the CoP would converge on a small area located 38 cm underground (Fig. 2). From these models, we have rebuilt the CoM trajectory using the boundary conditions of the single support phase, at TTO and THS.

The resulting trajectories are given in Fig. 3 and compared to the reference trajectory in the frontal, sagittal and horizontal planes, respectively. In Fig. 3, while the left figures isolate the CoM trajectories for a more precise observation, the right ones show a general view including the positions of the pivot points and the boundary rods (TTO and THS) for each model. One can notice wide differences around the vertical axis, as show the sagittal and the frontal views, while the respective positions of the pivot points allow the horizontal trajectories (in the X-Y plane) being very similar. Another observation shows that while the reference trajectory is not as circular than the approached ones in the sagittal plane, oppositely it is much more circular in the frontal plane. This is directly related to the fixed length of the pendulum rods in our representations. Indeed, implementing the reference support distance measured between the CoP and the CoM in each models would lead to much more precise approximations: the EIP model would be naturally identical to the reference by definition, while the IP and GIP models would get much closer than when using fixed-length pendulum rods.

Let us focus on the GIP model, more particularly on the convergence area of the support lines joining the CoP and the CoM. These two points may be considered as representative parameters of the motion dynamics of the subjects. We have observed that the support lines do not cross in the three dimensional space. We have then calculated the minimal distances between two support lines observed at two consecutive time steps, all over the single support phase. Let us define the instantaneous pivot point as the middle of these calculated minimal distances, which are represented in Fig. 4 for the complete support phase. In this figure, we can see that the values of the minimal distances are bounded between -0.026 and +0.27 mm, oriented positively in the direction obtained by the cross product between the two consecutive support lines. We have represented the distribution in space of these minimal distances and their respective orientation in Fig. 5. We have observed for all the walking subjects that the distances are almost null up to mid stance, and become larger after mid stance, following a linear length growth as shown in Fig. 4.

The minimal distances are mainly supported along the lateral axis, and their distribution (length and position in space) with time over the single support phase is illustrated in Fig. 5 using the black lines. Using these data, we have located a mean pivot point of the stance phase by calculating the weighted mean of the minimal distances which is represented in the same Fig. 5 using the red line. The mean distance in very small: XX mm. The blue line represents the minimal distance between the first and the last support lines of the stance phase (measured at TTO and THS, respectively). We call it the minimal boundary distance as it is based on the boundary support conditions. The distance between the mean pivot point (middle of the mean minimal distance) and the middle of the minimal boundary distance is represented by the green line: this distance remains smaller than 2 mm (maximum distance observed on the five subjects), which allowed us to consider that taking the boundary conditions of the support lines leaded to a good approximation of the position of the pivot point P ′ 0 . Moreover, comparing the minimal boundary distance and the sum of all the instantaneous minimal distances led to a difference of 0.012 mm which is negligible.

The evolution of the instantaneous rotation vector between two support lines is shown in Fig. 6. It is mainly oriented along the sagittal axis, its orientation remaining within a conic shape which center is the lateral Y -axis. The norm of the rotation vector linearly decreases from 4.85 10 -3 to 4 10 -3 at 60% of the single support phase, stabilizes until 90% and suddenly increases up to 4.4 10 -3 . As mentioned previously, its main orientation is along the sagittal axis. Nevertheless, as can be observed in Fig. 6, its projection on the sagittal plane follows a triangular pattern first oriented in the frontal direction, then slowly axing toward the vertical direction and back to the frontal one, then suddenly moving toward an axis of equation xz. As for the location of the pivot point, we have compared the weighted mean rotation vector to the average one between TTO and THS. Once again, using an average rotation vector between the bounding support conditions is very similar to the mean average of the rotation vector along the comlplete single support phase.

Finally, the pivot point P ′ 0 of the GIP-3D model introduced in this paper can be completely defined using the boundary support lines of the single support phase, measured at TTO and LTO respectively. The position of the pivot point is defined by the middle of the smallest distance between these two lines, and the rotation vector results from the cross product between P ′ 0 C T T O and P ′ 0 C T HS .

VI. CONCLUSION

We have proposed an inverted pendulum model, denoted GIP-3D, to represent human normal walking. Similarly to other existing models, it is based on the cyclic sinusoidal motion of the center of mass during normal walking. Our model also takes the mechanism of the foot into account by considering that the pivot point of the inverted pendulum is located under the ground level. We have observed more precisely the behavior of the convergence area of the support lines (i.e. the straight lines passing by the CoM and the CoP at each time step) and could determine the position of the pivot point only by considering the support lines at TTO and at LTO. The double support phase is then modeled using an homothetic transformation of the single support phase, which center is the approached CoM C ′ at LTO and ratio k 0 defined in Eq. This model presents several advantages: similarly to [START_REF] Hayot | Assessing the respective effects of the gait determinants on the vertical trajectory of the center of mass using the compass gait representation[END_REF], the mechanism of the foot is considered and allows reducing the vertical amplitude of the CoM motion, in comparison to classical approaches which fix the pivot point at ground level and located at the middle of the trajectory described by the CoP. But, in [START_REF] Hayot | Assessing the respective effects of the gait determinants on the vertical trajectory of the center of mass using the compass gait representation[END_REF]'s approach, the pivot point is not fixed and require precise measurement of the CoP position during the whole stance phase. In our approach, the pivot point position only requires knowing the state of the support lines at the boundary instants of the stance phase, LTO and THS, and is fixed for the whole single support phase. Another important advantage is that this model does not rely on the symmetric characteristics of normal walking, moreover it may be applied to pathological walking as long as the boundary states of each stance phase are known. To realize such applications, the length of the pendulum rod should also be considered as variable with time. The next steps in this study will be to extend and validate the GIP-3D model to pathologic walking using variable-length pendulum rods and continuity equations at the transition instants between the single and the double support phases.

Fig. 2 .

 2 Fig. 2. Views of the support lines (joining the CoM and the CoP) at each time step of the single support phase. The black line represents the trajectory of the CoM and the red line the trajectory of the CoP.

Figure 2

 2 Figure 2 shows the three dimensional evolution with time of the segment lines [P C] during single support, and the three associated views in the X-Z (sagittal view), Y-Z (frontal view) and X-Y (top or horizontal view) planes. In this figure, the reference CoM and CoP trajectories are also represented. These segments are called support lines in what follows, as they join the two points characterizing the dynamics of the walking motion (the CoM and the CoP) from a global point of view. The extension of these support lines in space show that they are converging in an area located around mid stance and 34 centimeters under the ground level. Using this observation, we have defined a generalized pendulum model represented in the sagittal plane in Fig. 1 using red lines. Observing the three projections in the X-Y, Y-Z and X-Z planes, it is very difficult to extract a regular pattern of the pivot area which would allow determining a single pivot point for our inverted pendulum representation. We will then focus on this particular area in order to precisely define the pivot point P ′ 0 represented in Fig 1.

  the frame attached to the linear pendulum during the double support phase. Here, P k (respectively E k ) may be replaced by the pivot point matching the considered model as summarized in Tab. I: either Zijlstra's continuous model (IP-3D): P k = P 0 (respectively E k = E 0 ), Hayot's extended model (EIP-3D): P k = P (t) (respectively E k = E(t)) or our generalized model (GIP-3D): P k = P ′ 0

  [START_REF] Hayot | Eip-3d: An extended 3d inverted pendulum model to represent normal human walking cycle[END_REF] with i = {P k , E k } and j = {s, d}. The indexes s and d refer to the single and double support phases, respectively. α(t) and β(t) are the sweep angles, around the lateral and the frontal axis, respectively. The two vectors P k C ′ and E k C ′

Fig. 3 .

 3 Fig. 3. Reference trajectory of the CoM (black line) and its approximation using IP-3D (red line), EIP-3D (blue line) and GIP-3D (green line). The straight lines represent the respective pendulum rods.

Fig. 6 .

 6 Fig. 6. Evolution of the instantaneous rotation vector of the support lines from a time step to the next one.


Nicolas Bonod

Alexis Devilez

Brice Rolly

Sebastien Bidault

Brian Stout

Ultracompact and unidirectional metallic antennas

Keywords: numbers: 78.20.Bh, 73.20.Mf, 78.67.-n, 32.50.+d

We investigate the angular redistribution of light radiated by a single emitter located in the vicinity of dipolar silver nanoparticles. We point out the fundamental role of the phase differences introduced by the optical path difference between the emitter and the particle and demonstrate that the polarizability of the metallic nanoparticle alone cannot predict the emission directionality. In particular, we show that collective or reflective properties of single nanoparticles can be controlled by tuning the distance of a single emitter at a λ/30 scale. These results enable us to design unidirectional and ultracompact nanoantennas composed of just two coupled nanoparticles separated by a distance achievable with biological linkers.

I. INTRODUCTION

Metallic nanostructures are key elements in the control of light interaction with quantum emitters. They can both focus light in tiny volumes and enhance the radiative decay rates of nearby emitters. [1][2][3][4][5][6][7][8][9] The latter property has been thoroughly studied in the case of a single emitter coupled to a single nanoparticle. [10][11][12][13][14][15][16][17][18][19][20] It has been shown that the radiative decay rates depend strongly on the distance between the emitter and the nanoparticle. Furthermore, at very short distances from the surfaces (a few nanometers), non radiative decay channels dominate and the quantum efficiency drops. The coupling efficiency between emitters and optical antennas also depends on the orientation of the dipolar source with respect to the dipolar modes supported by the metallic particles. A longitudinal coupling geometry significantly enhances the emission decay rates while a transverse interaction leads to moderate enhancements. [19][20][21][22] More recently, the ability of nanoantennas to control the angular emission of single molecules has been investigated. [23][24][25] This possibility is particularly important since a high directivity facilitates both the excitation of a quantum emitter by a collimated beam as well as the collection of the radiated light. 26 Yagi-Uda antennas have been successfully introduced at optical frequencies 27,28 and their high angular directivity has recently been confirmed experimentally. 29 This antenna geometry combines a director element generally consisting of a finite chain of identical particles and a reflector element typically based on a slightly larger particle. 30,31 A dipolar emitter can be coupled longitudinally to the antenna by utilizing a nanoparticle located near the reflector element. 32 In this configuration, the emitter is off-axis. For an on-axis emitter, the weak transverse coupling with the chain of particles can be reinforced by employing the so-called 'super emitter'consisting of a dimer of nanoparticles perpendicular to chain axis. 12,31 Pakizeh et al. have recently proposed an ultracompact antenna made of two identical metallic particles. 32 The dipolar emitter is then coupled to a dark mode characterized by opposite phase dipolar modes induced in the two neighboring particles. In that case, it has been shown that the emitter radiates predominantly in only one half space and can thus be characterized as unidirectional.

This paper is dedicated to studying how nearby spherical nanoparticles modify the angular distribution of light emitted by an oscillating dipole. In particular, we provide a thorough study of the phase differences between the dipolar source and the dipolar mode induced in the particle. We emphasize that these phase differences must take into account the optical path between the emitter and the particle in addition to the polarizability of the particle.

We show that by tuning the position of a single emitter from a single nanoparticle by a few tens of nanometers, we can sufficiently modify the dephasing to control the reflective or collective properties of the particle at a given frequency. When the phase difference between the exciting and induced dipoles is strictly equal to π, an equal part of the energy is radiated into each of the half spaces surrounding the emitter (the separation plane being perpendicular to the axis containing the dipoles). We then apply these results to the design of highly unidirectional antennas composed of two nanospheres separated by a mere 50 or 60 nm.

The basic concept is similar to the idea underlying the design of Yagi-Uda optical antennas, which associate the collective and reflective properties of nanoparticles, 27,28 but at much smaller interparticle distances and in simpler geometries. Taking into account the phase lag induced by the distance between the emitter and the nanoparticle, highly directional antennas can be designed with only two nanoparticles. We will first emphasize the role of the distance between the emitter and the nanoparticles by designing a directional antenna composed of identical particles, i.e. with strictly identical polarizabilities. Furthermore, by tuning the relative size of the two particles, it is possible to design an antenna smaller than λ/2 that channels light radiated by a single emitter in the angular aperture of commercial microscope objectives. 33

II. RADIATION PROPERTIES OF A DIPOLAR EMITTER COUPLED TO A SINGLE NANOPARTICLE

We first investigate the radiation pattern of a single emitter located near a single 90 nm silver particle. The dipolar emitter is polarized along the z-axis in order to provide a transverse coupling with the nanoparticle as sketched in Fig. 1. The emission properties of the dipolar emitter are calculated in the framework of rigorous Lorentz-Mie theory, and combined with multiple scattering theory in configurations where more than one silver particle is present. [34][35][36] This analytical method is particularly well suited to tackle light scattering by an ensemble of nanospheres. In order to insure an accurate modeling of the short range couplings, the calculations presented in this study are carried out with 30 multipole orders.

Nevertheless, let us emphasize that the electromagnetic response of the metallic particles under consideration is predominately dipolar in nature, resulting in a rapid multipole convergence. Consequently, a dipole approximation (1st order) would qualitatively exhibit all the underlying physics observed in this work.

The radiation patterns are obtained from the radial component of the Poynting vector in the far field. In order to estimate how much light is collected or reflected by the metallic particle, we define the reflection coefficient R as the ratio of the power emitted in the x ≥ 0 hemisphere with respect to the total radiated power. In this work, we pay particular attention to the dephasing between the emitting and induced dipoles. When the dipole approximation dominates (as it does here), it is sufficient to calculate the electric field at the center of the metallic particle. The induced dipole moment of the nanoparticle is then obtained by multiplying the total electric field by

V s ϵ 0 (ε s -ε b )
where V s is the volume of the sphere, ϵ 0 the permittivity of vacuum and ε s and ε b are the relative permittivities of the metal and the background media respectively. For small particles, the quasi-static approximation applies and we can express the resulting phase differences as the sum of the phase differences due to the optical path difference field from the emitter and the polarizability of the particle. It must be stressed that in previous works, attention was focused on the phase difference of a nanoparticle polarizability with respect to its local excitation fields, while in this work we emphasize importance of taking into account the additional phase difference induced by the (small but non-negligible) optical path between the emitter and the nanoparticle. We point out that we chose the common convention of the phase differences defined to lie in the range from -π to π. The phase differences ∆ϕ between the emitter and dipole moment of a (D = 90 nm) silver nanoparticle are displayed in Fig. 2 (circles) as a function of the emission wavelength.

The emitting and induced dipoles are precisely in opposing phase at λ r = 600 nm, ∆ϕ = π, and remain roughly in opposite phase for longer wavelengths. For wavelengths smaller than λ r = 600 nm, ∆ϕ varies as a function of λ, and the induced dipole is generally out-ofphase with respect to the emitter. The phase of the polarizability of the nanoparticle is also displayed (squares) and the phase value of π/2 at λ = 500 nm indicates the plasmon resonance. Let us now investigate the radiation properties of the coupled system in terms of the reflection efficiency of the nanosphere. We present in Fig. 2 the reflection efficiency, R, of the nanoparticle (i.e. defined as the power emitted in the x ≥ 0 half space over the total emitted power) as a function of the emission wavelength λ. When the reflection efficiency is lower than 0.5, the dipole preferentially radiates towards the x ≤ 0 half space and the nanoparticle behaves as a collector (c.f. Fig. 1: the nanoparticle is located on the negative

x-axis). The full line in Fig. 2 clearly indicates that depending on the phase differences between the emitting and induced dipoles, a metallic nanoparticle can either collect or reflect light radiated by a single emitter. Let us mention that the dipolar approximation (triangles) exhibits almost all the features of the radiation properties, meaning that the coupling between the single emitter and the metallic particle is almost perfectly dipolar for λ > 450 nm. For wavelengths λ > λ r , the emitting and induced dipoles are nearly opposite in phase and the emitter radiates preferentially towards the metallic nanoparticle with R dropping down to 0.4 at λ = 665 nm. For λ < λ r , the emitting and induced dipoles are out-of-phase and the nanoparticle reflects with a rather high efficiency the radiated light since R reaches 90% at λ = 510 nm.

The emission pattern of the light radiated by the exciting dipole at λ = 510 nm is reconstructed in Fig. We consider that in the far-field limit (r ≫ λ), |r -

x j | -r = √ (x -x j ) 2 + y 2 + z 2 -r ≈ r( √ 1 -2xx j /r 2 -1) ≈ -x j (x/r) ≈ -x j sin(θ) cos(φ).
The electric and magnetic fields produced in the far field by p 1 (j = 1) and p 2 (j = 2) then write :

E j (r, θ, φ) = ( ω c ) 2 1 4πrϵ 0 e ik|r-x j | [(e r ×p j ) × e r ] (1) = ( ω c ) 2 1
4πrϵ 0 e ikr e -ikx j sin(θ) cos(φ) p j sin(θ)(-e θ )

(2)

H j (r, θ, φ) = kω 4πr
e ikr e -ikx j sin(θ) cos(φ) p j sin(θ)(-e φ )

with e r , e θ , and e φ the unit vectors of the spherical basis (Fig. 1(b)). The resulting far-field, time-averaged Poynting vector of the sum of these fields writes :

P(r, θ, φ) = 1 2 Re[(E 1 + E 2 ) * × (H 1 + H 2 )] (4)
=

ω 3 k 32π 2 ϵ 0 c 2 r 2 (p *
1 e -ikr e ikx 1 sin(θ) cos(φ) + p * 2 e -ikr e ikx 2 sin(θ) cos(φ) )

(p 1 e ikr e -ikz 1 sin(θ) cos(φ) + p 2 e ikr e -ikz 2 sin(θ) cos(φ) ) sin

2 (θ) e r = ω 3 k 32π 2 ϵ 0 c 2 r 2 [ |p 1 | 2 + |p 2 | 2 + 2 Re(p 1 p * 2 e ik(d+a) sin(θ) cos(φ) ) ] sin 2 (θ) e r ( 6
)
To study the symmetry of the radiation pattern, we compute the sum of the Poynting vectors in one direction and in the opposite direction :

∆P(r, θ, φ) = P(r, θ, φ) + P(r, π -θ, φ + π) (7) = ω 3 k 16π 2 ϵ 0 c 2 r 2 Re{2ip 1 p * 2 sin[k(d + a) sin(θ) cos(φ)]} sin 2 (θ) e r (8)
In our case, we are interested by the evolution of ∆P with respect to the relative phase (9) Hence when ϕ = kπ with k ∈ N, ∆P = 0 for any value of the dipolar amplitudes |p 1 | and |p 2 |. This calculation demonstrates that when the emitted and induced dipoles are in phase or in opposite phase, the emission from the two dipoles is perfectly symmetric with respect to the origin, while the electromagnetic environment of the emitter can be highly asymmetric. This confirms that the wavelength of the plasmonic resonance (λ = 500 nm) taken alone cannot predict the directionality of the emission and that the distance between the emitter and the nanoparticle plays a crucial role. Moreover, this model exhibits a very interesting property: in a given direction defined by θ and ϕ, the change of the sign of sin(ϕ) will change the sign of ∆P, for every distance d. In order to confirm this assumption, phase differences λ (nm) (Fig. 3(a)) and reflection efficiencies (Fig. 3(b)) are now displayed as a function of distance d and λ. These graphs confirm the clear correlation between reflection efficiency and the phase differences of the emitting and induced dipoles. The isoefficiency line is plotted in Fig. 3(b) when the directionality is null and it matches the isodephasing line (see Fig. 3(a)) plotted for ∆ϕ = π. Calculations performed for a 60 nm silver particle (Fig. 3(c)-(d)) show that a similar behavior is obtained, but that the opposite-phase wavelength, λ r , is shifted towards shorter wavelengths. These graphs evidence that in the 500 nm -600 nm range (with d = 30nm), the smaller particles (60 nm) mostly collect electromagnetic radiation while larger particles (90 nm) act as reflectors. This property allows the design of Yagi Uda antennas with a reflector made of slightly larger particles and a collector made of an array of identical smaller particles. [27][28][29][30][31] More importantly, these calculations show that at a given frequency, a single sphere can act as a reflector or a collector depending on its distance from the emitter, and that this behavior can be controlled inside a very small range of distances (a few tens of nm). For example, we can observe in Fig. 3(d) that at λ = 550 nm, the nanoparticle behaves as a collector when it is at a distance of 40 nm from the emitter while it behaves as a reflector when this distance is reduced to 10 nm. Consequently, it is possible to design unidirectional antennas by assembling two identical particles, i.e. with identical polarizabilities and by tuning the distances between the emitter and both particles. However, Fig. 3 suggests that a stronger directionality can be achieved by assembling two particles of different diameters with minimum and maximum reflection efficiencies at the emission wavelength, property that cannot be achieved with equal diameters. We can thus design an ultracompact antenna made of 90 nm and 60 nm particles at an equal 30 nm distance from the emitter, geometry chosen to optimize the dephasing between the dipolar moments of the nanospheres and the emitter.

ϕ = ϕ 1 -ϕ 2 between p 1 = |p 1 |e iϕ 1 and p 2 = |p 2 |e iϕ 2 : ∆P(r, θ, φ) = ω 3 k 16π 2 ϵ 0 c 2 r 2 {-2|p 1 ||p 2 | sin(ϕ) sin[k(d + a) sin(θ) cos(φ)]} sin 2 (θ) e r
d (nm) (a) (b) (c) (d)

III. ULTRACOMPACT AND UNIDIRECTIONAL NANOANTENNA

Before discussing antennas with non equal diameters, let us begin this section by designing an ultracompact antenna composed of two identical nanoparticles of diameters 60 nm closely separated by a distance of only 50 nm. The emitter is located at 10 nm from the first particle (which thus acts as the reflector) and 40 nm from the second particle (the collector) (see Fig. 4a). Fig. 4b shows the reflection efficiency as a function of the wavelength of emission.

A good directivity can thus be achieved with strictly identical particles since 70% of the emitted energy is radiated in the positive x half space. Let us note that both distances, respectively 10 nm and 40 nm are much smaller than the emitting wavelength and achievable with biological linkers. 7,18 Let us now design an asymmetric antenna made of two silver particles with different diameters to optimize their reflective and collective properties. The antenna, consisting of 60 nm and 90 nm diameter silver spheres, was optimized with respect to the particle sizes while keeping d equal to 30 nm (see Fig. 5(a)). In a first step, we compute the phase differences between the emitter and the two induced dipolar modes supported by the particles. The wavelength range in which one particle acts as a collector while the other reflects radiation exceeds 100 nm, which is larger than the width of a typical fluorescent emitter. Fig. 5(b) shows that this antenna geometry fulfills these conditions for wavelengths ranging from 475 nm to 600 nm. As expected, the ratio of the radiated power towards the x ≥0 half space is maximum when the phases are of opposite sign, and it can reach more than 97% at λ = 520 nm. The emission pattern is reconstructed in Fig. 5(c) at λ = 520 nm as a function of the polar angle (in the plane xOz). It confirms the high unidirectionality of this antenna since the emitted power towards the left half plane is unobservable at this scale. Moreover, this antenna narrows the angular redistribution of emitted light compared to a single nanoparticle since the angular width of the emission cone is less than ±30 degrees. Such angular openings are easily achievable with commercial microscope objectives 33 . In practice, this asymmetric nanoantenna, although much smaller than the vacuum emission wavelength of the oscillating dipole, is almost perfectly unidirectional. For the sake of completeness, the evolution of the radiative and total decay rates and the quantum efficiency of the emitter in the vicinity of the optical antenna are displayed in Fig. 5(d). For that purpose, the total emitted power, P tot , and the radiated power, P , are calculated by integrating the radial component of the Poynting vector over a spherical surface surrounding the source, at respective distances of 1 nm and 50 µm. The total (Γ tot ) and radiative (Γ rad ) decay rate enhancements are then obtained by normalizing the emitted power in the presence of the antenna by the emitted power (P 0 ) in the homogeneous background medium: Γ tot = P tot /P 0 and Γ rad = P rad /P 0 .

The quantum efficiency is then defined as η = Γ rad /(Γ tot +(1-η i )/η i ) where η i is the intrinsic quantum efficiency. We consider in this work a perfect emitter (η i = 1). Fig. 5(d) shows that the high directionality achieved at λ = 520 nm is not associated with a drop of the radiative decay rates which confirms that it relies on the association of the reflective and collective features of nanoparticles rather than on opposing phases between the induced dipoles of the nanoparticles. 32 The radiative decay rates obtained are comparable with those observed when dealing with Yagi-Uda antennas 30,31 and they are significantly enhanced by coupling the source dipole to a "super-emitter" as shown in Fig. 6(a). In this last case, we consider a more complex antenna geometry to combine unidirectionality and strong radiative rate enhancements by introducing two dipolar particles coupled longitudinally to the emitter.

Fig. 6(b) shows that a super emitter can strongly enhance the radiative decay rates by up to 3 orders of magnitude. In particular, at λ = 610 nm, the quantum yield is maximum, the radiative decay rate is enhanced by more than 500 and the unidirectionality is preserved as shown on Fig. 6(d).

IV. CONCLUSION

The reflection or collection behavior of the nanoparticle depends on the total phase difference between the emitting and induced dipoles which includes both the polarizability of the metallic particle and the optical path between the emitter and the nanoparticle. We showed the importance of the role of the optical path between the emitter and a metallic particle on the redistribution of light for distances smaller than λ/30. One consequences of this observation was to remark that when the emitting and induced dipoles are exactly in opposing phase, the ratio between the radiated powers in the backward and forward directions is precisely equal to unity. We unveiled the importance of the optical path by designing a directional antenna composed of two identical nanoparticles. The reflective and collective properties were tuned by controlling the distance between the emitter and the nanoparticles at a scale of λ/30. We also presented a means to design highly directive and ultracompact nanoantenna by tuning the relative sizes of the silver particles (while still keeping the overall size much smaller than the vacuum emission wavelength). While angular openings obtained with single plarticles are around 40 degrees, the dimer nanoantennas narrow the angular opening of emitted radiation to around 30 degrees, rendering the radiation readily collectible by commercial microscope objectives. Finally, we showed that the radiative decay rate of an emitter can be increased by three orders of magnitude by introducing a dimer antenna longitudinally coupled to the emitter while preserving a high directivity. * nicolas.bonod@fresnel.fr

FIG. 1 .

 1 FIG. 1. (a) Sketch of a dipolar emitter oriented along the z-axis and located at a distance d from the surface of a silver nanosphere. The refractive index of silver is taken from Ref. 37. Silver nanospheres are embedded in a polymer of refractive index 1.5. (b) The spherical coordinates used in the analytical expressions.

FIG. 2 .

 2 FIG. 2. (a) A dipolar emitter oriented along the z-axis is located at a distance d = 30 nm from a silver sphere of diameter 90 nm. The phase difference between the emitting and induced dipoles (circles, left scale) and reflection efficiency (right scale. Full line: complete calculation, triangles: dipolar approximation) as a function of the emission wavelength. Emission patterns of the oscillating dipole at (b) λ = 510 nm and (c) λ = 600 nm in the xOz plane, (d = 30 nm) calculated by plotting the radial component of the Poynting vector normalized by the forward emitted value as a function of the polar angle.

  2(b) by plotting the radial component of the Poynting vector normalized by the forward emission value as a function of the polar angle (in the plane xOz). This result demonstrates the good unidirectionality offered by a single nanoparticle. The angular width of the emission pattern defined as the angle between the on-axis maximum value and the direction of half maximum value is of the order of ±40 degrees. Finally, it is interesting to note that for λ = λ r , the ratio of the radiated energy in both half spaces is precisely equal to 1. It may seem surprising that a highly asymmetric environment (a single particle at the left of the emitter) results in perfectly symmetric radiation (see the radiation pattern in Fig.2(c)). To fully understand this counter-intuitive result, let us consider the analytic expression of the emission of two dipoles with moments denoted p 1 and p 2 . They are placed along the x-axis with separation d + a, and oriented along the z-axis (see Fig.1(a)). Let us note x 1 = -(d + a)/2 and x 2 = (d + a)/2 the positions of the dipoles along the x-axis.

FIG. 3 .

 3 FIG. 3. (a) Phase differences normalized by π between the dipolar emitter and the induced dipolar moment of the nanoparticle and(b) reflection efficiency as a function of the wavelength of emission and the distance d between the emitter and the metallic surface of a 90 nm silver sphere. Figs. (c) and (d) show similar results as (a) and (b) respectively but for a 60 nm sphere.

FIG. 4 .

 4 FIG. 4. (a) The nanoantenna is composed of two identical silver nanoparticles of diameter 60 nm. The emitter is located at 10 nm from the left particle, and at 40 nm from the right particle. (b) Reflection efficiency as a function of the emitted wavelength. (c) Emission pattern of the oscillating dipole at λ = 550 nm in the xOz plane.

FIG. 5 .

 5 FIG. 5. (a) The nanoantenna is composed of silver particles of diameter 60 nm and 90 nm. The emitter is located at 30 nm from both particles. (b) Full line, right scale: Reflection efficiency as a function of the emitted wavelength; Dotted and dashed lines, left scale: dephasing of the emitting and induced dipoles supported by spheres of diameter 90 nm (dashed line) and 60 nm (dotted line). (c) Emission pattern of the oscillating dipole at λ =550 nm in the xOz plane. (d) Full and dashed lines, left scale: radiative (dashed line) and total (full line) decay rates enhancements as a function of the wavelength emission. Full line+squares, right scale: quantum efficiency (the intrinsic quantum efficiency is equal to 1).

FIG. 6 .

 6 FIG. 6. (a) Sketch of the ultracompact nanoantenna + super emitter: the dipolar source is longitudinally coupled with two 60 nm silver particles on the z axis, with an emitter-particle distance of 8 nm. The emitter-particle distances on the x axis are equal to Fig. 5(a) at 30 nm. (b) Radiative and total decay rate enhancements as a function of the emission wavelength. (c) Quantum efficiency as a function of the emission wavelength. (d) Emission pattern of the oscillating dipole at λ =610 nm in the xOz plane.

1

  H. Metiu, Progress In Surface Science, 17, 153 (1984).