Bruno Robisson
email: bruno.robisson@cea.fr

Michel Agoyan

Patrick Soquet

Sébastien Le Henaff

Franck Wajsbürt

Pirouz Bazargan-Sabet

Guillaume Phan

Sébastien

SMART SECURITY MANAGEMENT IN SECURE DEVICES

Keywords: Hardware tamper resistance, Formalism to detect and counter fault and side-channel attacks, Architectures for trusted computing, Application of fuzzy logic

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Security is a key component for information technologies and communication. Among the security threats, a very important one is certainly due to vulnerabili-ties of the integrated circuits that implement cryptographic algorithms to ensure confidentiality, authentication or data integrity (such as smartcards). With the access to one of these circuits, the attacker tries either to reconstruct the functionality of the circuit (reverse engineering) or to recover cryptographic materials when the cryptographic algorithm is known (physical or hardware cryptanalysis). Both threats share a set of techniques. The first one, called side channel attacks, consists in observing some physical characteristics which are modified during the circuit's computation [START_REF] Mangard | Power Analysis Attacks Revealing the Secrets of Smart Cards[END_REF][START_REF] Bauer | Horizontal collision correlation attack on elliptic curves --extended version[END_REF]. The second technique, called fault attacks, consists in disrupting the circuit's behavior [START_REF] Ali | Differential fault analysis of aes: towards reaching its limits[END_REF][START_REF] Barenghi | Fault injection attacks on cryptographic devices: Theory, practice, and countermeasures[END_REF]. The third one consists in getting information about the chip design by direct inspection of its structure [START_REF] Torrance | The state-of-the-art in ic reverse engineering[END_REF]. This inspection may be performed by using any kind of imaging techniques or by using destructive means such as abrasion, chemical etching or focused ion beam. Combinaison of those three techniques have also been proposed [START_REF] Roche | Combined fault and side-channel attack on protected implementations of AES[END_REF]. Many protections have therefore been proposed to counter such attacks. Some protections (hereafter referred to as "sensors") give information about the state of the system either by measuring the light, the voltage, the frequency or the temperature of the chip or by detecting errors during computations. This detection is generally based on spatial redundancy (i.e. performing the same computation several times simultaneously), temporal redundancy (i.e. performing the same computation several times) or information redundancy (i.e. performing a computation with more bits than required) [START_REF] Carlet | Higher-order CIS codes[END_REF] or ad-hoc sensors [START_REF] Dutertre | Improving the ability of bulk built-in current sensors to detect single event effects by using triple-well CMOS[END_REF][START_REF] Homma | Em attack is non-invasive? -design methodology and validity verification of em attack sensor[END_REF]. Several mechanisms are also proposed to detect a modification in the software execution flow. In some cases, an additional hardware block is dedicated to this task [START_REF] Mao | Hardware support for secure processing in embedded systems[END_REF]. To render physical attacks more difficult, "noise" has been added, for example, by using an internal clock, by randomizing the order of the instructions, by adding dummy operations or by masking the internal computations that can be predicted by the attacker [START_REF] Rivain | Higher-order masking and shuffling for software implementations of block ciphers[END_REF], [START_REF] Coron | Analysis and improvement of the random delay countermeasure of CHES 2009[END_REF]. Another way to reduce sensitivity to side channel attacks consists in reducing the correlation between physical values (such as power consumption or electromagnetic radiation) and the data processed, for example, by using balanced data encoding and balanced place and route [START_REF] Guilley | Evaluation of power constant dual-rail logics countermeasures against DPA with design time security metrics[END_REF], by using power filters or electromagnetic shields. Finally, some countermeasures modify the functional behavior of the circuit in case of attacks. Such "reactions" may consist, for example, in temporarily stopping the communication with the reader (the card "mutes") and/or resetting parts of the running software. The ultimate reaction consists in permanently destroying (i.e. killing) all the data (including sensitive information) stored in the chip.

In practice, the chip security is achieved by a combination of such countermeasures. Nowadays, state-of-the-art circuits can protect data for weeks, months and even years. But the implementation of these countermeasures not only drastically decreases the performance of the chip in terms of power and speed but also decreases its availability. In order to overcome this limitation, we propose a complementary approach: implementing complex management (through the application of a "strategy of security") of the whole set of countermeasures embedded on a chip. Such a "system-level" approach has already been addressed in [START_REF] Gogniat | Reconfigurable hardware for high-security/ high-performance embedded systems: The SAFES perspective. Very Large Scale Integration[END_REF]. The authors proposed mechanisms which reconfigure the architecture of cryptographic hardware blocks embedded in an FPGA in accordance with performance and energy consumption criteria. Our work addresses this tradeoff but also the trade-off between security and availability. Besides, we consider that these trade-offs have to be adjusted only with technologies which are regularly available in smart cards (i.e. without the use of hardware reconfiguration capabilities provided by FPGAs).

In section 2, a representative case study (a conditional access system for Pay-TV) and the associated expected properties for complex strategy of security are detailed. In section 3, an example of such a strategy, based on fuzzy logic, is proposed. In section 4, an innovative hardware/software prototype enabling the execution of this strategy is described. Finally, in section 5, the results obtained with this prototype are presented.

Case study

In the followin g sections, an implementation of the chosen representative system, the smart card part of a conditional access system (hereafter referred to as the "host system"), is described. Its functional parts and some of its protections are presented. Finally, by expressing the difficulties of designing such a system, the specifications for complex strategies of security are highlighted.

Conditional Access System: Functional Part

The functional part of the conditional access system (or CAS) consists of a JavaCard application which is interpreted by a virtual machine (VM) which runs on a micro-controller.

Application The Conditional Access Systems (CAS) protect content (such as radio, TV, data stream) by requiring certain criteria to be met before granting access to this content. One of the main criteria is to own a smart card that stores sensitive information. This information is used by the receiver to decipher the content. Security of this content is achieved by combining protections at different levels. A management key protects the transfer of access rights (subscriptions, the related period of validity, geographical localization of the receiver, etc.) and the transfert of the exploitation key. This key, which is changed approximately once every month, protects the access criteria associated with the content and the value of the control word. This control word, which is changed approximately once every 5 to 10s, protects the content.

Virtual Machine

The chosen virtual machine (VM) complies with the Java Card 2.2.2 specifications [START_REF]Java card technology website[END_REF], and provides card content management capabilities based on the GlobalPlatform Standard [START_REF]Global platform specifications website[END_REF]. Hence, it is possible to download and install on-card applications (or applets) running on top of the VM and to use standard Java Card and GlobalPlatform Application Programming Interfaces (APIs).

Host's micro-controller The host sub-system is built on a 5-stage pipelined 32bit Harvard RISC microcontroller able to execute one instruction per clock cycle. The instructions are stored in 640kB of ROM and the data in 256kB of RAM and in 128kB of EEPROM. The peripherals include two Universal Asynchronous Receiver Transmitters (or UARTs) compliant with the ISO7816 and RS232 standards and an Advanced Encryption Standard (or AES) crypto-engine (to accelerate the ciphering and deciphering of data) [START_REF]NIST: Announcing the Advanced Encryption Standard (AES)[END_REF].

Protections

The different components of the host described above embed the protections described below. Some of them, such as the security sensors, may have outputs. Other countermeasures may be configured. Configuration parameters are either integer numbers (for redundancy level, insertion of dummy instructions and random power generator) or boolean values (which activate or not a security reaction). As will be explained in section 4.1, even though these protections are physically embedded in the host sub-system, these parameters are not driven by the host.

Security sensors Several sensors have been implemented or emulated. Some of them measure physical characteristics such as voltage, light intensity or chip temperature. Others detect errors during the execution flow of a piece of code, during I/O data transmissions or during computations. One way of detecting errors consists in performing the same computation several times and then comparing the results. If they are identical, the computation is considered error-free. If not, a counter of corrupted executions, noted CE, is increased. The number of times that the same computation is performed is called the "redundancy level" (noted RL). A redundancy level of "i" is noted "×i". It is equal to 1 when the redundancy countermeasure is not activated but when RL ≤ 3, the host has error-correction capabilities.

The various sensors which are taken into account in our study and their different possible outputs are described in Table 1. The hardware (HW) or software components ("VM" denotes the virtual machine and "App" the application) which update these outputs are listed in the last column of this table. For example, the number of times that the check (performed by the application) of a Message Authentication Code (MAC) is false is noted M E. It takes its value between 0 and 10 4 . Another sensor is the number of times that a cryptographic operation is performed with the same key. This number is noted CO.

Insertion of Dummy Instructions (IDI)

Another protection consists in increasing the level of noise by inserting dummy instructions (i.e. that are not useful) during computations. The execution of a program is thus made up of several sequences of valid instructions followed by dummy instructions. In our framework, which has been drawn from [START_REF] Ambrose | RIJID: Random Code Injection to Mask Power Analysis based Side Channel Attacks[END_REF], two parameters D and N configure the countermeasure. D and N denote the maximum number of consecutive useful and dummy instructions respectively in a sequence. Note that N is equal to 0 when the countermeasure is not activated. ). Let us also suppose that the power consumption added by one random generator also follows a Gaussian distribution with a constant mean µ R and with a standard deviation which is equal to σ c . We also suppose that the different RNG are identical and that their power consumptions are statistically independent. The probability density function of the total power consumption, noted x tot (t), of the circuit when R RNG are activated is thus the following:

pdf (x tot (t)) = e - (x tot (t)-µc (t)-R•µ R ) 2 2•σ 2 c •(1+R 2 ) σ c • √ 1 + R 2 • √ 2 • π (1)
Note that R is equal to 0 when the RPG countermeasure is not activated.

Security reactions Two reactions have been taken into account: stopping the communication and resetting the software (hereafter noted "Mute/Reset"), or destroying all the data (hereafter noted "kill"). But several reactions could also be considered at the application level: forbidding the access to one or several TV channel packages, suppressing services such as video on demand, restricting the access to only free TV or program contents, etc.

Other protections In the case of a commercial device, other security features (such as power filters, shields, balanced logic, etc.) which can never be deactivated, are also supposed to protect the circuit, independently of all kinds of strategy of security.

Impact of protections

An estimation of the evaluation of performances in terms of speed, energy consumption and security for different values of the parameters of the countermeasures is proposed in Annexe 1. To estimate the security level, we have only considered the threat related to side channel attacks (SCA) and differential fault attacks (DFA).

-The gain in terms of SCA, noted F SCA, is the ratio between the number of curves needed for the attacker to recover the key when the countermeasure is activated and the number of curves without countermeasure (the greater the better). -The gain in terms of DFA, noted F DF A, is the ratio between the number of experiments needed to recover the key when the countermeasure is activated and the number of experimentations without countermeasure (the greater the better). -The loss in terms of speed, noted F T ime, is the ratio between the duration of a computation with the countermeasure and the duration of the same computation without countermeasure (the smaller the better). -The loss in terms of energy, noted F N RJ, is the ratio between the energy consumption with countermeasure and the energy consumption without countermeasure (the smaller the better).

Figure 1 displays the values of these factors for different values of D and N (with RL = 2 and R = 3). Note that the values do not share the same scale. This figure clearly shows that performance decreases as the security level increases. Note also that the theoretical considerations described above are only approximations and have to be validated with experimental results.

Difficulties

The design of the smart card part of a CAS, as described in section 2.1, is a very challenging task because of the following specifications:

-Its security level has to be very high: the security of a smart card is expected to hold for 2-5 years. But it is generally used for 8 years. In this condition, a card which resists one additional year represents large savings for the content and CAS providers. -Its performance has to be high: For Pay-TV applications, the switch from channel to channel has to be as fast as possible.

-The availability has to be high: A malfunction of the card causes a prejudice for the user (who is no longer able to watch his favorite TV program) but also for the CAS provider who has to deal with expensive field feedback. A malfunction may be due to abnormal system usage (if the card, for example, is put in a very hot place) but also to incorrect system integration (if the card, for example, is inserted in a poor quality reader). Such cases will further be called "anomalies", in the sense that the component operates in abnormal conditions but is not subjected to attacks. -It processes data with various sensitivity levels. For example, the retrieval of the exploitation key (EK) would enable the pirate to decipher the control word (CW) during the period of validity of the key EK. With this CW, he gains access to the content. So, the corruption of the key EK is a very serious security threat. On the contrary, knowing the value of the word CM at a given moment is of little value because it is only valid for a very short period of time. Let us call DS the sensitivity level of data which is handled by the application. -Its power consumption must be low if it is embedded in a mobile phone or in a multimedia tablet. -It has to be inexpensive.

Need for complex strategies of security

Section 2.3 shows that the setup of the countermeasures defines the performances of the circuit in terms of security, speed and power consumption (but also availability) and that those performances are closely linked. For example, the increase in security causes the decrease in speed. Thus, the set of performances which may be reached, despite the wide range of setups of the countermeasures, is restricted. This restriction may be such that none of the performances of this set fulfills the specifications described above. In such a case, the specifications are antagonistic and the designer has to make a trade-off. He has, for example, to choose between availability and security.

We propose to add mechanisms which enable to dynamically modify the setups of the countermeasures, i.e. which enable to switch from a high performance but "poorly" secured state to a low performance but secured one. To draw an analogy, these mechanisms provide a new "degree of freedom" which make it possible to reach new sets of performances. In other words, these mechanisms release the links that exist between the different performances. We consider that a smart strategy of security should be able to control these mechanisms to fulfill all the specifications (or at least as many as possible).

In this section, a complex strategy of security tailored for the application described in section 2.1 is proposed. The chosen strategy is inspired from methods and concepts developped in Intrusion Detection Systems (or IDS) [START_REF] Bace | Intrusion detection[END_REF]. In this framework, two methods are proposed. The first method consists in determining the whole set of states of the system reached by a user with "normal" behavior. Each time the state of the system leaves this set, an "anomaly" is detected. The second method is complementary. It consists in determining the whole set of states which would have to be reached by an evil-minded user to perform an attack. When the system enters this set, misuse (considered an attack) is detected. Contrary to classical IDS, the component is off-line most of the time. So, it has to react autonomously without any outside help. In order to obtain, in this context, both robust and secure behavior, mechanisms must detect both anomalies and misuses.

According to principles developed for IDS, the chosen strategy of security is decomposed into three different processes: The first consists in collecting information about the state of the host system, the second, called analysis, consists in computing (from this information) the anomaly and misuse levels (called resp. M L and AL) and the third in configuring (according to M L and AL) the parameters of the countermeasures.

Information sources

We consider that the sensitivity level of data DS (provided by the application and defined in 2.4) and the outputs {LS, V S, EF E, CE, P E, N E, M E, CO} (provided by the sensors and defined in section 2.2) are the inputs of the analysis algorithm.

Analysis

As the secure circuit designer generally expresses the anomaly and misuse cases with words which are very often vague and inaccurate, fuzzy logic has been used to formalize the analysis mechanism [START_REF] Dubois | Fundamentals of Fuzzy Sets Series: The Handbooks of Fuzzy Sets[END_REF]. The details of such an analysis is out of scope of this article but it is sketched here. Within this framework, the analysis process is mainly described by a set of "IF-THEN" rules and an inference process.

The "IF-THEN" fuzzy rules are expressed with words or "linguistic variables" such as LOW, HIGH and adverbs such that RATHER, VERY whose meanings are precisely defined with "fuzzy sets". For example, the chosen strategy of security is expressed in our system by using a dozen of rules such as :

-R 0 : "IF the number of methods that have processed without error (N E) is VERY HIGH THEN the misuse is LOW " -R 1 : "IF the voltage (V S) is RATHER HIGH and the light (LS) is HIGH THEN the misuse is HIGH " -R 2 : "IF the number of cryptographic errors (CE) is RATHER HIGH THEN the misuse is HIGH " -R 3 : "IF the number of PIN code error (P E) is RATHER HIGH OR the number of triggers of the voltage sensor (V S) is HIGH THEN the misuse is HIGH " -. . . The inference process chosen for computing the misuse and anomaly level from fuzzy rules and inputs was first proposed by Mamdani in [START_REF] Mamdani | Application of fuzzy logic to approximate reasoning using linguistic synthesis[END_REF]. This process consists in fuzzifying the inputs, then computing the degree of truth of each rule, then aggregating the results of all the rules and lastly "defuzzyfying" them to obtain the misuse and anomaly levels (which are scalar values between 0 and 1). This process can be activated, for example, each time an input value is modified.

Configuration of countermeasures

The configuration of the countermeasures has been chosen from the analysis of their impact, described in section 2.3. We have evaluated the impact of the RPG for R ∈ {0; 3; 10} (the value for α is chosen equal to 10%). As the cost of the redundancy technique is very high, we have only considered RL ∈ {1; 2; 3} (the default value for q is 8) and the impact of the IDI for D ∈ {0; 4; 8} and for N ∈ {2; 3; 4} (the value for m is chosen equal to 150). Among all those possibilities, only four configurations of the different countermeasures are taken into account. These configurations are chosen since they cover a large range of security/performance trade-offs. They are further called Saf e, U nsaf e, Critical and F atal and are reported in Table 2. It is important to note that whatever the configurations, the different physical sensors are always activated. It is also important to note that, in practice, the security level is determined as the minimum value between F DF A and F SCA (in bold in the table). The function which assigns the configuration to each value of the anomaly (AL) and misuse (ML) levels is described in Table 3. The different thresholds have been determined according to the parameters chosen for the analysis process.

Prototyping

This section describes a hardware implementation which has been defined for our case study but which could also be used for different strategies of security.

Architecture

We choose an architecture, represented in Figure 4, which separates the calculations related to the execution of the application software and those related to the management of the security. This leads to the addition of a system called "monitor" and strictly dedicated to the management of the security. In the left part of the figure, the host system, described in section 2.1, is represented with its software part (the CAS application running above the JavaCard virtual machine) and its hardware part (the 32-bit RISC microcontroller). Each of these parts implements its own set of countermeasures described in section 2. The software and hardware parts of the monitor and an example of communication between the host and the monitor are described above.

Operating system and strategy of security The monitor uses a specific component-based operating system. It needs a small memory footprint (less than 4kB) and it responds in less than a hundred clock cycles to the host's requests. An application is a set of components linked together through exchange memory zones. A component is a set of ports towards these exchange zones, of private attributes and of methods (which describe its behavior). The strategy of the security formalized by using fuzzy logic and described in section 3 is embedded in such a component. It is important to note that, thanks to a particular choice of fuzzification functions and off-line precomputations, the fuzzy logic mechanisms can easily be implemented in the very light-weight hardware of the monitor. Monitor's micro-controller The micro-controller of the monitor is built around a 5-stage pipelined 32-bit Harvard RISC processor. To increase the security of the system, this processor does not share any piece of hardware with the host system and in particular it has its own addressing space: the instructions are stored in 32kB of ROM and data in 4kB of RAM. The peripherals comprise one UART, two FIFOs and an Interrupt Controller Unit (ICU). These ICU and FIFOs provide two communication channels: The monitor is informed of the events occurring in the host (without having access to the sensitive data themselves) and in return, the monitor is able to change the configuration of the countermeasures embedded into the host.

Communication between the host and the monitor The communication between the host and the monitor is based on a request/acknowledge protocol. The request is always initiated by the host. The host waits after its request until the monitor responds. The following example describes such a communication initiated by the application.

1. The application indicates that the sensitivity of the data that it manipulates (DS) is changing. 2. The virtual machine sends this information to the monitor via the FIFOs. 3. The host stops its current execution. 4. From the fuzzy sets and the "IF-THEN" fuzzy rules defined by the user, the monitor processes the CM's outputs and the DS's value by using the fuzzy reasoning described in section 3.2. The outputs of this reasoning are the misuse and the anomaly levels. These levels are used to select the configuration of the countermeasures thanks to Tables 3 and2 reported in section 3.3. 5. The monitor asks to the host system to configure software (via FIFOs) and hard wired countermeasures (via the ICU). 6. The monitor waits until all the countermeasures are configured and are ready. 7. The monitor asks to the host system to resume its execution.

FPGA synthesis results

The whole system has been implemented in the Xilinx Virtex-5 FPGA of an ML501 evaluation platform. This platform has been extended with an ISO7816 interface adapter (for an easy connection with off-the-shelf smart card readers) and with an additional RS232 serial interface (for trace or debug purposes). The details of the hardware implementation of the host, of the monitor and of the communication channels are described in Figure 9 of Annexe 3. The combinatorial part of the system (i.e. excluded the memory areas) uses 48% (3490 slices out of 28800) of the total FPGA's resources in term of slices. As a reminder, a slice is the smallest hardware element used in a Xilinx FPGA. A slice contains four D flip-flops associated to four look-up tables with six independent inputs. Excluding the monitor, the resources usage falls to 34% (2462 slices out of 7200). So, the combinatorial part of the monitor represents around 40% of the combinatorial part of the whole system. As the resources associated with the memory parts are far larger than the combinatorial ones, the estimated total overhead (i.e. including the memory areas) decreases to about 7%. But it is important to note that this total overhead greatly depends on the size of the host's memories.

Security analysis

Side channel analysis

For this discussion, we distinguish two kinds of information leakage. The first kind is "directly" related to the value of sensitive data (such as cryptographic keys) and is exploited with classic side channel techniques. The second kind of leakage is not directly linked to the sensitive data but only helps the attacker to understand the functioning of the system and so, to focus his attacks. For example, if the circuit blurs the manipulation of sensitive data with the activation of RPG, the attacker could easily detect, by measuring the power consumption, the location in time of this sensitive computation. This information will help him perform his side channel attacks on this location.

The system is protected against direct information leakage by the host's countermeasures. The theoretical efficiency of some countermeasures is estimated in Section 2.3. Several methods have also been proposed to quantify in practice the direct information leakage of circuits in terms of "number of traces to recover the key". One of these methods, described in [START_REF] Standaert | A unified framework for the analysis of sidechannel key recovery attacks[END_REF] is based on a profiled information theoretic analysis. If the security of the host has been first estimated with such a worst-case security evaluation, the user is able to know how many bits per trace the host leaks. From this data, from the number of cryptographic executions which has be done with the same key (sensor CO), he is able to estimate the number of bits of the key which could have been retrieved by the attacker with side channel analysis. The parameters of the strategy of security are to be chosen according to these estimations. Besides, as the host sends to the monitor only the values of the sensors and the sensitivity level of the data handled by the application, there is no additional direct information leakage due to the communication channels or due to the monitor.

The indirect information leakage is linked to the "discretion" of the host's countermeasures. We consider that the computation of the configuration of the countermeasures executed in the monitor is very small compared to the signature of the countermeasure itself and so does not represent an additional indirect leakage of information. Besides, a more complex strategy of security such as the one described in Section 3, could be defined in order to reduce the indirect information leakage (for example, by creating a randomly high level of noise even when no sensitive data is manipulated).

Fault attacks

As for protections against side channel attacks, the system is mainly protected by the host's countermeasures. In particular, the host has to embed highly effective protections in order to be protected against highly effective attacks such as those described in [START_REF] Moradi | A generalized method of differential fault attack against AES cryptosystem[END_REF] (which enable the pirate to retrieve the key with only one faulted execution). In the following discussion, we consider that it is the case, at least when sensitive data is handled by the application. During fault attacks, several cases have been distinguished:

-If the functioning of the monitor and of the communication channels are not corrupted during a fault attack, two sub-cases are possible:

• The functioning of the host is not corrupted during the fault attack. In this case, the system functions normally. • The functioning of the host is corrupted during the fault attack. In this case, as the host is protected by using error detection or correction, the host stops its computation, the counter of corrupted execution is incremented and this information is processed by the monitor. The monitor will decide to increase the redundancy level until providing error correction capabilities (rendering DFA ineffective). The only possibility for the pirate in such a case consists in bypassing the HW/SW error detection or correction capabilities of the host to avoid the increase of the redundancy level.

-If the functioning of the monitor (or of the communication channels) are corrupted during a fault attack, this incorrect functioning does not directly give information about sensitive data because the monitor never handles such data but two sub-cases are possible:

• The monitor is unable to compute a configuration. Two sub-sub-cases have to be taken into account: * The functioning of the host is not corrupted during the fault attack. In this case, as the host waits after any request until the monitor responds, there is no risk of external communications of data. * The functioning of the host is also corrupted during the fault attack.

In this case, as the host is protected by using error detection or correction, the host stops its computation and a request towards the monitor is triggered through the increase of the counter of corrupted execution. As the host waits after this request until the monitor responds, there is no risk of external communications of data. The only possibility for the pirate in such a case consists in bypassing the HW/SW error detection or correction capabilities of the host to avoid any request towards the monitor.

• The monitor computes an incorrect configuration of the countermeasures. In this case, the sensitive data is still protected by the countermeasures that can never be deactivated and by this incorrect set of countermeasures. In the worst case, the attacker will have to face a lower level of security. In the best case, he will have to face a higher level of security. But in both cases, the host has to be attacked after the monitor.

This discussion shows that in order to take advantage of the proposed architecture, the pirate has to attack first the monitor and, in a second step, to attack the host. Such a scenario is unfortunately possible with state-of-the-art attack equipment but we consider that, in practice, it leads to more difficult attacks.

Discussion

As explained in Section 4.2, the proposed architecture leads to a penalty of 5-20% (depending of the size of the host's memories) in terms of area (and to some extent, of power consumption), but we underline that these figures have been obtained with no design optimization. The penalty in terms of performances of the adjunction of the monitor greatly depends on several parameters such as, for example, the rate of communication of data between the host and the monitor, the rate of change of configurations of the countermeasures or the time necessary for the host to reconfigure its countermeasures. But it is important to note that these penalties have to be compared with the penalties due to the countermeasures embedded in the host.

Another difficulty concerns the design of the best strategy of security for a given application. The theoretical considerations of Section 2.3 help the user, first, to evaluate a priori, that is according to a given set of hypotheses, the performances and the security levels of the host for several sets of parameters and, second, to choose an a priori strategy of security for the given application. As the set of hypotheses may be false (and it is unfortunately often the case in the security domain), the performances have to be measured in practice and security evaluations have to be performed with state-of-the-art attack equipment. As most secure designs, these results have to be used to improve the a priori strategy of security. At the end of this refinement process, the strategy should not only increase the availability but also trade-off the performances and the security according to the application's constraints.

Finally, the testing and debugging of the whole prototype is also more complex than for a single component. Each hardware and software component has been debugged one by one. After these unitary tests, the monitor and the host have been tested independently. Finally, the host and the monitor have been tested together.

This article has completed the first step towards "system level" management of the security dedicated to the improvement of the availability and the performance of a circuit without reducing its security. The three main contributions of this article are the following. First, the impact in terms of security and performances of different combinations of well-known countermeasures has been quantified. This study has shown, without surprise, that security and performances are antagonistic. But this study has also shown that, by modifying only a few parameters of these countermeasures, states with very distinct performance and security levels can be reached. Second, we have proposed a strategy of security designed to minimize both the rate of anomalies considered to be attacks (to increase the availability) and the rate of attacks considered to be normal functioning (to increase the security). The strategy is based on mechanisms which dynamically modify the parameters of the countermeasures. Third, we have proposed an HW/SW architecture which implements these mechanisms. The proposed solution is based on a double-processor architecture: one processor embeds a representative set of countermeasures (and mechanisms to define their parameters) and executes the application code. The second processor, on the same chip, applies the strategy, but without requesting sensitive data from the first processor.

Future work will consist in determining the theoretical performances and security levels for combinations of other countermeasures. It will also consist in refining the strategy by performing real attacks on the prototype and in measuring its real performance (preliminary simulation results are described in Annexe 4). We also plan to develop mechanisms to enable over the air upgrades of fuzzy rules according to the evolution of the threat. It could also be interesting to investigate, in parallel, more expressive formalisms such as, for example, those used for the detection of intrusion in complex computer systems.

Annexe 1: Impact of protections

As explained in 2.3, we consider only the threat related to side channel attacks (SCA) and differential fault attacks (DFA). In order to successfully mount such attacks, the attacker proceeds in a divide and conquer manner (i.e. he attacks small pieces of the key one by one). On each iteration of these attacks, he targets the result of one particular piece of computation, hereafter referred to the "targeted result". It is, for example for SCA, the computation of a SBox during the first round. To mount a DFA, the attacker also needs the cipher text resulting from one or several faulty computations.

Impact of Redundancy

FSCA We consider that the redundant computations generate identical power traces which could easily be added in the time domain by the attacker, decreasing the number of curves needed to recover the key by a factor RL.

FDFA As explained in section 2.2, a redundant computation is associated with a comparison of the results in order to increment the counter CE if the results are different. To bypass the redundancy protection, the attacker will have to both avoid the update of the counter in case of error detection and realize several faults of the same value, noted e 0 , during the RL successive computations of the targeted result. The probability of realizing such a set of faults determines the number of realizations which are required to mount the attack. We shall call q the number of bits of the targeted result. If we consider that all the faults on these bits are equally probable, then the probability of realizing the same fault e 0 (whose value does not matter) during the RL successive execution of the targeted instruction is equal to (1/2 q ) RL-1 . In classical DFA schemes, a fault generally has to affect 1 byte. So, the default value is chosen equal to q = 8.

FTime In our framework, the redundant computations are not performed in parallel. We also assume that the comparison of the results is negligible in terms of computation time. Thus, the redundancy countermeasure increases the computation time by a factor RL.

FNRJ We assume that the comparison of the results is negligible in terms of energy consumption. So, the redundancy countermeasure increases the energy consumption by a factor RL.

Impact of Insertion of Dummy Instructions

Let us call D the random variable equal to the number of instructions in the sequences of useful instructions. The domain of D is chosen equal to {1; . . . ; D}. Let us call N the random variable equal to the number of instructions in the sequences of useless instructions. The domain of this random variable is set equal to {0; . . . ; N }. We consider that the random variables N and D follow uniform distributions. Let us also suppose that the m th valid instruction should be the instruction which computes the targeted result. For the sake of simplicity, we consider that each instruction (useless or not) is executed in one clock cycle. The typical value for m is equal to 100-200 which corresponds to a software implementation on a 32-bit processor of a round of an AES [START_REF]NIST: Announcing the Advanced Encryption Standard (AES)[END_REF]. Let us define the random variable X equal to the number of the clock cycle associated with the execution of the instruction m.

F SCARL = 1 RL F DF ARL = (1/2 q ) RL-1 F T imeRL = RL F N RJRL = RL
Each realization x of this random variable, is equal to

x = k i=1 (d i ) + k i=1 (n i ), with k such that k i=1 (d i ) = m,
with i being the i th sequence of useless/useful instructions and with d i and n i the i th realizations respectively of the random variables D and N . We have x = m + k i=1 (n i ). We consider that, because m ≫ D, x could be approximated by x ∼ m + q i=1 (n i ) with q = 2 • m/(D + 1). In these conditions, the density of probability of X follows a normal distribution (µ X ,σ X ) with:

µ X = m + q • µ N = m • (1 + N/2)
(2)

σ 2 X = q • σ 2 N = m • N • (N + 2) 6 • (D + 1) (3)
For the sake of simplicity, we consider, as proposed in [START_REF] Clavier | Differential power analysis in the presence of hardware countermeasures[END_REF], that m is uniformly distributed (with the probability 1 out of 2 • σ X ) between mσ X and m + σ X .

FSCA In these conditions, the SCA peak is reduced by a factor 2 • σ X and the number of curves necessary to retrieve the key increases by a factor 4 • σ 2 X . But by using the sliding window method (with consists in reconstructing the peak by integrating the consumption curves on 2 • σ X samples) also described in [START_REF] Clavier | Differential power analysis in the presence of hardware countermeasures[END_REF], this saving in terms of number of power curves is only equal to 2 • σ X .

FDFA In order to realize a DFA, we suppose that the attacker is able to target clock cycles comprising between mσ X and m + σ X . In these conditions, he has one chance out of 2 • σ X to modify the instruction m. The number of faulty realizations is thus increased by a factor 2 • σ X .

FTime The formula 2 indicates that the computation time is increased by a factor (1 + N/2).

FNRJ Because we consider that useful and dummy instructions consume the same energy, the consumption of the circuit is also increased by a factor of (1 + N/2).

F SCAIDI = 1 if N=0 2 • m • N•(N+2) 6•(D+1) otherwise F DF AIDI = F SCAIDI F T imeIDI = 1 + N/2 F N RJIDI = F T imeIDI

Impact of Random Power Generator

FSCA According to [START_REF] Clavier | Differential power analysis in the presence of hardware countermeasures[END_REF], if we call δ the amplitude of the SCA peak (i.e. the difference of the power consumption or electromagnetic radiation between data) and σ c the standard deviation of the power consumption curve, the number of power curves necessary to recover the key has to be higher than (σ c /δ) 2 . The activation of the R random generators increases this number by a factor (1+R 2 ).

FDFA We suppose that the RPG does not protect against differential fault attacks.

FTime As the RPG processes at the same time as the computation, computation time is not increased by its activation.

FNRJ We consider that the power consumption of an RNG is equal to α times the temporal mean of µ c (t) (with α = 10% because the random number generator is a small piece of hardware). In these conditions, the total energy consumption of the circuit is increased by a factor (1 + α • R).

Combination of countermeasures When the different countermeasures are combined, factors computed above for the different countermeasures are simply multiplied for side-channel attacks, for the time of computation and the energy consumption. For the DFA, the chance of making the same fault (in order to bypass detection of RL) on the targeted result (whose position is blurred with IDI) is equal to the chance obtained in the first occurrence (i.e. 1/(2 q • 2 • σ X )) up to the redundancy level minus one. So we obtain: The process chosen for inferring a decision from fuzzy rules and inputs was first proposed by Mamdani in [START_REF] Mamdani | Application of fuzzy logic to approximate reasoning using linguistic synthesis[END_REF]. As reported in the following pseudo-code, this process consists in fuzzifying the inputs (according to the fuzzy subsets described in section 7.1), then computing (according to 7.2) the degree of truth of the rules described in section 3.2, then aggregating the results of these rules (according to 7.3) and lastly "defuzzyfying" them to obtain the output values (according to 7.4). Only the computation of the misuse level is described in detail below.

F SCARP G = 1 + R 2 F DF ARP G = 1 F T imeRP G = 1 F N RJRP G = (1 + α • R)
F SCA = F SCARP G • F SCAIDI • F SCARL F DF A = F DF ARP G • F DF A RL-1 IDI • F DF ARL F T ime = F SpeedRP G • F SpeedIDI • F SpeedRL F N RJ = F N RJRP G • F N RJIDI • F N RJRL

Fuzzy subsets

Definition The significance of the words (or "linguistic variables") "low" and "high" for the inputs and outputs are defined with "fuzzy subsets". The membership function, denoted µ A , of a fuzzy subset A is a generalization of the characteristic function in classical set theory. It is any function mapping the values s of the input S to the real unit interval [0, 1]. The value µ A (s) is called the membership degree of s in the fuzzy subset A. This degree quantifies the grade of membership of the element s to the fuzzy subset A: the value 0 means that s is not a member of the fuzzy set; the value 1 means that s is a full member of the fuzzy set. The values between 0 and 1 characterize fuzzy members, which belong to the fuzzy set only partially. The degree of membership also quantifies the degree of truth of the assertion "s is A" (which is true with a degree µ A (s)).

Membership functions for inputs Firstly consider the fuzzy sets associated with the value of the input S i . The 8 membership functions (or fuzzy subsets) selected for each input S i are presented in Table 2. For example, let us consider the voltage sensor input whose maximum value is equal to V S max = 10. If this input is equal to 3 (and so is comprised between V S max /5 and 2 • V S max /5) it is considered "rather high" with a grade of truth of 1/4 and "very low" with a grade of truth of 1/2. If this voltage sensor value is equal to 7 (and so is comprised between 3 • V S max /5 and 4 • V S max /5) it is considered "high" with a grade of truth of 2/3 and "very very low" with a grade of 0. Membership functions for outputs The outputs of the analysis mechanisms are the levels of anomaly and misuse. Their values are real values in [0, 1]. The two different membership functions "LOW" and "HIGH", considered for these outputs are described in Table 3 The second analysis consists in defining several attack and normal use scenarios and playing these scenarios on an algorithmic model of the prototype. For example, consider a scenario where the card is connected to a low quality card reader. In such an abnormal case, the electrical connection between the two devices regularly triggers the voltage sensor. In the first part of the scenario (called I), there is no error and the level of security remains low, as represented in Figure 10 (top). In this figure, the x-axis corresponds to the time and the y-axis to the levels of security. In the second part (called II) of the scenario, the voltage sensor alone is triggered. The security level increases slowly (because we consider that these mistakes are not important). In the third part (called III) of the scenario, the MAC error sensor M E is also triggered. In this part, the security level increases quickly. In the last part of the scenario (called IV), the sensors stop being triggered and the security level decreases quickly. On the contrary, let us consider a laser attack scenario. In the first part of this scenario, there are no errors and the level of security remains low, as represented in Figure 10 (bottom). In the second part, the light sensors are triggered because the attacker injects faults in the middle of a long sequence of correct commands. Even with this precaution, the level of security increases quickly (we consider that the triggers of the light sensor are important). We suppose that the attacker is able to detect this increase (the activation of the RPG when the system switches to the "unsafe" configuration is easy to detect) and that, in the third part of the scenario, he stops to inject faults. The level of security tends to decrease. But when the laser attack is re-started in the fourth part of the scenario, the monitor increases the security level very quickly until sensitive data is deleted.

Fig. 1 .

 1 Fig. 1. Decrease of the performances with the increase of the security

Fig. 2 .

 2 Fig. 2. Countermeasure configurations

  2. These countermeasures may be configured via signals noted "CMs config." and may give information about the state of the host via signals noted "CMs outputs". The signal which indicates the sensitivity of the data handled by the application is noted DS. All these signals are exchanged between the host and the monitor through a communication channel. The monitor is represented in the right part of the figure. The software parts of the monitor are the strategy of security and the component-based operating system. Its hardware part is a micro-controller.

Fig. 4 .

 4 Fig. 4. Architecture of the prototype (hardware in dark grey, software in grey)

Fig. 5 .

 5 Fig. 5. Performances for the redundancy level countermeasure

Fig. 6 .

 6 Fig. 6. Performances for the IDI countermeasure

Fig. 7 .

 7 Fig. 7. Performances for the RPG countermeasures

Fig. 8 .

 8 Fig. 8. Performances for combination of the countermeasure

Fig. 9 .

 9 Fig. 9. Hardware part of the prototype

Fig. 10 .

 10 Fig. 10. Response of the system for an anomaly case (top) and for an attack case (bottom)

Table 1 .

 1 Outputs of the CMs

	Name Values	Description	Updated by
	LS	{0, 1, . . . , 5} # of triggers of the light sensor	HW
	V S	{0, 1, . . . , 10} # of triggers of the voltage sensor	HW
	EF E {0, 1, . . . , 10} # of corrupted execution flow	VM
	CE {0, 1, . . . , 10} # of corrupted execution	VM
	P E {0, 1, . . . , 10} # of wrong PIN	App
	N E {0, 1, . . . , 10 3 } # of methods processed without error VM
	M E {0, 1, . . . , 10 4 } # of MAC errors	App
	CO {0, 1, . . . , 10 7 } # of cryptographic execution	App
	Random Power Generator (RPG) In order to blur the power consumption
	of the circuit, several Random Number Generators (RNG) may be activated.
	Assuming that each step t of this power consumption, noted x(t), is drawn from
	a Gaussian shaped probability distribution (with a mean µ c (t) and a constant
	standard deviation σ c

  Algorithm 1 Algorithm for calculating the Misuse Level (ML) Fuzzy Sets for the outputs (see 7.1 ) Require: The set of rules (see 3.2) Fuzzify the values of the inputs Compute the degree of truth of each rule of the rule set (see 7.2) Aggregate all the rules to obtain the membership function of the ML (see 7.3) Defuzzify this membership function to obtain the scalar value of the ML (see 7.4)

	Require: Inputs: Scalar values of the inputs (S)
	Ensure: Output: Misuse Level (ML)
	Require: Fuzzy Sets for inputs (see 7.1)
	Require:

Table 2 .

 2 The 8 membership functions (or fuzzy subsets) for an input S i ( S i max is the maximum value of this input)

					Very	Very
	Name	Rather		Very Very	Very	Very	Rather
		Low Low Low	Low	High High High High
	Acronym s i ∈ [0; S i max /5]	L S i max -1	L S i max --1	L S i max ---L S i max ----H 1 1	S i max ++++ H 0	S i max +++ H 0	S i max ++ 0	H	S i max + 0
	s i ∈]S i max /5; 2 • S i max /5]	3/4	2/3	1/2	0		0		0	0	1/4
	s i ∈]2 • S i max /5; 3 • S i max /5] 1/2	1/3	0	0		0		0	1/3	1/2
	s i ∈]3 • S i max /5; 4 • S i max /5] 1/4	0	0	0		0	1/2	2/3	3/4
	s i ∈]4 • S i max /5; S i max	0	0	0	0		1		1	1	1

DATA CM ctrl reg Timer CODE 32 KB

Name

LOW

HIGH o ∈ [0; 0, 2] 1 0 o ∈]0, 2; 0, 8] -5/3 • o + 2/3 5/3 • o -1/3 o ∈]0, 8; 1] 0 1

Table 3. Membership functions for outputs

Operations on fuzzy subsets The boolean operations such as AND, OR and NOT are also defined in fuzzy logic. In the following, we consider the standard Zadeh operators on fuzzy subsets A 0 and A 1 respectively defined on variable x and y such as:

The extension of Z AN D and Z OR to "n-ary" operators is trivially defined and denoted as:

The degree of truth of the premise is a real number in [0, 1] which depends on the values S of the inputs. This degree of truth of the premise of the rule i is denoted pre i (S).

Compute the degree of truth of the rules

Computing the values of the premises The degree of truth of each premise is computed first by evaluating the degree of membership of the input S i to fuzzy subsets. Then, if necessary, the Zadeh operators are applied according to the premise's formula.

It is important to note that because we have chosen input membership functions that are discontinuous and because we have chosen Zadeh operators, regardless of the values of the entries and regardless of the different formulae used to express the premise i, the result pre i (S) is always an element of the set P = {0; 1/4; 1/3; 1/2; 2/3; 3/4; 1}; Modification of the membership function of the conclusion of a rule In the method proposed by Mamdani, the degree of truth of the premise of a rule modifies the membership function of its conclusion. The modification consists in truncating the membership function A k (y) of the conclusion with the value of the premise, that is:

In fuzzy logic inference mechanisms, all the rules are considered to fire in parallel. So, at first glance, the use of this logic could lead to inconsistency, that is, several rules could lead to different conclusions (for example, the misuse is both HIGH and LOW). The aim of aggregation of rules and defuzzification is to compute a unique value for the decision.

Aggregation of rules

The different rules are considered to be linked together with the operator OR. So, combining the rules consists in taking for all y ∈ [0, 1], the maximum value of the conclusions of the different rules, according to the following formula:

Defuzzification

The operation of defuzzification consists in calculating a scalar value (also called "crisp" value) for the output membership function of the conclusion (µ R (y|S)), obtained given a set of entries and a set of rules. Before explaining this process of defuzzification, this output membership function is rewritten by taking into account the properties of our set of rules. We distinguish the rules which conclude to a LOW value for misuse (noted "LOW-m" rules) and those which conclude to the HIGH values for misuse (noted "HIGH-m" rules). Without loss of generality, we reorder the rules so that those numbered from 0 to q -1 are "LOW-m" rules and those numbered from q to p are "HIGH-m" rules. The considered set of rules is noted R.

Output membership function rewriting

The fuzzy subset of LOWm is defined as: µ LOW -m (y|S) = max(min(pre 0 (S), µ LOW (y)), . . . , min(pre q-1 (S), µ LOW (y))

= min(max(pre 0 (S), . . . , pre q-1 (S)), µ LOW (y))

= min(

In the same way, the fuzzy set associated to the "HIGH-y" rules is defined as:

µ HIGH-m (y|S) = min(

As explained in section ??, the different premises pre k are in the set P . So, let us define p l and p h and the membership functions µ LOW -m (y) p l and µ LOW -m (y) p h such as:

The output membership function for the output y can in these conditions always be written in the following form: µ R (y) p l ,p h = max(µ LOW -m (y) p l , µ HIGH-m (y) p h ) Defuzzification techniques There are different kinds of defuzzification techniques to compute the crisp output from the output membership function. We have considered four of the more popular ones, called "centroid" (CT), "mean of max" (MM), "first of max" (FOM) and "last of max" (LM). For example, the FOM crisp value is computed using the formula:

Then, we have computed their values for the whole set of possible values of premises, that is, for all of the 49 membership functions µ R (y) p l ,p h with P h ∈ P and P l ∈ P . The results for First Of Max ("FOM") is reported in Table 4.


J Boreux
email: jehan.boreux@fundp.ac.betimoteo.carletti@fundp.ac.be

T Carletti

Ch Skokos
email: hskokos@pks.mpg.de

M Vittot
email: vittot@cpt.univ-mrs.fr

Hamiltonian control used to improve the beam stability in particle accelerator models

Keywords: hamiltonian, chaos, control, accelerator mappings, dynamical aperture

We develop a Hamiltonian control theory suitable for a 4D symplectic map that models a ring particle accelerator composed of elements with sextupole nonlinearity. The controlled system is designed to exhibit a more regular orbital behavior than the uncontrolled one. Using the Smaller Alignement Index (SALI) chaos indicator, we are able to show that the controlled system has a dynamical aperture up to 1.7 times larger than the original model.

Introduction

Particle accelerators are technological devices allowing studies at both "infinitely small scale", e.g. particles responsible for elementary forces, and "extremely large scale", e.g. the origin of cosmos. In a simplified approach, such devices are composed of basic elements sequence: focusing magnets, defocusing magnets, accelerating electromagnetic fields and trajectory bending elements in the case of ring accelerators. The resulting total dynamics is thus nonlinear, and can be described, in the absence of strong damping, by a conservative system. This system can be modelled by a symplectic map built from the composition of several elementary maps corresponding to each basic magnetic element. One of the main problems in the dynamics of ring accelerators is to study the stability around the nominal orbit, i.e. the circular orbit passing through the centre of the ring. Each component of the ring can be seen as a nonlinear map that deforms the trajectory at large amplitude. Moreover, such maps possess stochastic layers whose effect is the reduction of the stability domains around the nominal circular orbit (the so-called dynamical aperture -DA) [1]. Such behaviors imply that (chaotic) nearby orbits can drift away after a few ring turns, eventually colliding with accelerator's boundaries, and consequently reducing the beam lifetime and the accelerator performance. The goal of this paper is to present a reliable improvement of beam stability by increasing the DA in a simplified accelerator model, consisting of only one type of element having a sextupole nonlinearity [2,3,4,5]. We work in the framework of the Hamiltonian Control Theory presented in [7,8], where two methods of controlling symplectic maps have been described, namely using Lie transformations and generating functions. In the present paper we use the former method, which allows direct determination of the new controlled map; avoiding the possible problems related to coordinate inversion. The aim of control theory is to improve selected features of a given system by slightly modifying its Hamiltonian with the addition of a small control term, so that the new system and the original one are conjugated namely, they have the same dynamics. This technique is particularly suitable whenever one can directly act on the system and modify it, e.g. in the case of a particle accelerator where the addition of a control term in the Hamiltonian function can be seen as the introduction of a suitable magnet in the accelerator lattice. In our study, we use the Smaller Alignement Index (SALI) [6,[START_REF] Ch | [END_REF]11,12] method, which is an efficient indicator for characterising orbits as chaotic or regular in Hamiltonian flows and symplectic maps. The SALI is computed using the time evolution of two deviation vectors along the studied orbit. The paper is organized as follows: After the introduction of the accelerator model in Section 2, we present a general result for the control of symplectic maps in Section 3 and the application to the studied model in Section 4. Then we briefly recall the SALI chaos indicator in Section 5, while Section 6 presents our results on the chaotic behavior of the model. Finally, in Section 7 we summarize our conclusions. In the appendix, several technical theoretical issues are discussed in detail.

The model

Let us consider a charged particle passing through a simplified accelerator ring with linear frequencies (tunes) q x , q y , having a localized thin sextupole magnet (for more details the interested reader is referred to [2]). The magnetic field of this element induces a modification of the orbit once the particle passes through it, modelled by the following 4D symplectic map

        x 1 x 2 x 3 x 4         =     cos ω 1 -sin ω 1 0 0 sin ω 1 cos ω 1 0 0 0 0 cos ω 2 -sin ω 2 0 0 sin ω 2 cos ω 2             x 1 x 2 + x 2 1 -x 2 3 x 3 x 4 -2x 1 x 3         (2.1) = T         x 1 x 2 x 3 x 4         . (2.2)
where x 1 (x 3 ) denotes the deflection from the ideal circular orbit in the horizontal (vertical) direction before the particle enters the element, and x 2 (x 4 ) is the associated momentum. Primes denotes positions and momenta after the particle passes through the element. The parameters ω 1 and ω 2 are related to the accelerator's tunes 1 q x and q y by ω 1 = 2πq x and ω 2 = 2πq y . The first matrix in (2.1) describes the linear motion of a particle, which corresponds to a simple rotation in the phase space. The nonlinearity induced by the thin sextupole magnet is represented by the 2 nd order polynomial expression in (2.1). The particle dynamics at the n-th turn, can be described by the sequence (x

(n) 1 , x (n) 2 , x (n) 3 , x (n) 4 ) n≥0
, where the (n + 1)-th positions and momenta are defined as a function of the n-th ones by 2.1. Map (2.1) decomposes naturally in an integrable part, the rotation by angles ω 1 , ω 2 , in the planes x 1 , x 2 and x 3 , x 4 respectively, and a quadratic "perturbation". The following observation will be a key point used in the next section: such maps can be obtained as the time-1 flow of the following Hamiltonian systems (see A)

H(x 1 , x 2 , x 3 , x 4 ) = -ω 1 x 2 1 + x 2 2 2 -ω 2 x 2 3 + x 2 4 2 and V (x 1 , x 2 , x 3 , x 4 ) = - x 3 1 3 + x 1 x 2 3 , (2.3)
more precisely using the compact form of the Poisson brackets

2 x = T ( x) = e {H} e {V } x , (2.4)
where x = (x 1 , x 2 , x 3 , x 4 ) T , with T denoting the transpose of a matrix, and by definition, for any function f defined in the phase space, {H}f = {H, f } = (∇H) T J∇f , with J = 0 1 -1 0 , being the symplectic constant matrix,

e {H} f = n≥0 {H} n n! f and {H} n f = {H} n-1 ({H}f ) . (2.5)
The dynamics of map 2.1 has been already studied in [6] where it has been shown that its phase space contains both regular and chaotic orbits. The chaotic orbits reduce the DA to a hypersphere of radius ∼ 0.39 in the 4-dimensional phase space (see Figs. 5 and 6 of [6]). The goal of the present paper is to show that the stability region of the nominal circular orbit can be increased once map 2.1 is controlled by an appropriately designed map.

Control theory for symplectic maps

The general aim of Hamiltonian control theory is to slightly modify a given system, by introducing a small control term, to improve some chosen feature of the system (e.g. to reduce the chaotic regions or to build invariant tori).

In the following we will be interested in controlling a quasi-integrable symplectic map in such a way that it will allow us to obtain a new, controlled map "closer" to the integrable part of the original map, and thus increase the stability region around the nominal circular orbit. Since chaos opens the way to diffusion in the physical space, the controlled map is expected to have a smaller number of escaping orbits and a larger region occupied by invariant curves in a neighbourhood of the origin. For this purpose, we apply the method of [8], except that in the present case the integrable part is not expressed in action-angle variables. In particular, the integrable part is a rotation, so the present theory applies to perturbations of rotations, instead of maps close to identity.

Let us consider an integrable symplectic map defined through its infinitesimal generator

3 H x = e {H} x , (3.6)
and consider the quasi-integrable map perturbation of the former

x = T ( x) = e {H} e {V } x , (3.7)
where x ∈ R 2N and V is a perturbation, namely V = o(H). The aim of Hamiltonian control theory is to construct a third map, the control map, whose generator F is small (it satisfies F = o(V )) and moreover the controlled map

T ctrl = e {H} e {V } e {F } , (3.8)
will be conjugated to a map T * , closer to e {H} than T (see 3.12 below). We note that the use of the exponential of a Poisson bracket, ensures that such maps are symplectic by construction. More precisely, let us define the unperturbed map

A -1 = e -{H} , (3.9)
and observe that (1 -A -1 ) is not invertible, since its kernel contains any smooth function of H. Thus we assume the existence of a "pseudo-inverse"operator, G, that should verify (see [8] for details)

G 1 -A -1 G = G . (3.10)
At this point we can define the non-resonant and the resonant operators

N := 1 -A -1 G and R := 1 -N , (3.11)
which are projectors, i.e. N 2 = N and R 2 = R.

Our main theoretical result is the following theorem:

Theorem 1 Under the above hypotheses and defining S = GV we have

e {S} T ctrl e -{S} = e {H} e {RV } := T * , (3.12)
where

T ctrl = e {H} e {V } e {F } , (3.13)
with a control term given by e {F } = e -{V } e {(N -G)V } e {RV } e {GV } .

(3.14)

Remark 2 (Warped addition) Let us define as in [8] the warped addition, {A} ⊕ {B}, of two operators by e {A} e {B} := e {A}⊕{B} .

(3.15)

An explicit formula can be obtained using the Baker-Campbell-Hausdorff formula [START_REF] Bourbaki | Eléments de Mathématiques: Groupes et Algèbres de Lie[END_REF], where {A} ⊕ {B} is a series whose first terms are

{A} ⊕ {B} = {A} + {B} + 1 2 ({A}{B} -{B}{A}) + . . . , (3.16)
hence the warped addition is a deformation of the usual addition between operators.

Proof: Using this warped addition, we can rewrite the controlled map

T ctrl = e {H}⊕{V }⊕{F } , (3.17)
where the control term 3.14 can also be rewritten as

{F } = -{V } ⊕ {(N -G) V } ⊕ {RV } ⊕ {GV } . (3.18)
From 3.11 we have:

N -G = -A -1 G . (3.19)
One can easily prove (see appendix A of [8]) that

(-{H}) ⊕ (-{S}) ⊕ {H} = -{A -1 S} , (3.20)
hence, recalling the definition of S and 3.19, we can rewrite 3.18 as

{F } = -{V } ⊕ (-{H}) ⊕ (-{S}) ⊕ {H} ⊕ {RV } ⊕ {S} . (3.21)
By rearranging the terms we can easily get

{S} ⊕ {H} ⊕ {V } ⊕ {F } ⊕ (-{S}) = {H} ⊕ {RV } , (3.22)
which is nothing but 3.12 rewritten using the warped addition.Q.E.D.

Remark 3 Let us observe that the control term is, as required, small compared to V . In fact from 3.18 and by using the approximated formula for the warped addition 3.16, we obtain

{F } = -{V } ⊕ {(N -G) V } ⊕ {RV } ⊕ {GV } (3.23) = -{V } + {N V } -{GV } + {RV } + {GV } + o(V ) = o(V ) , (3.24)
where we use the relation

N + R = 1.
Under the assumption of absence of resonances, i.e. RV = 0, the o(V ) term in the control map can be explicitly computed to give (see Appendix A)

F = 1 2 {V }GV + o(V 2 ) . (3.25)
4 The control term for the non-resonant map (2.1)

Let us now calculate the control map for system (2.1). The variables (x 1 , x 2 , x 3 , x 4 ) in which the generators 2.3 are given are not suitable to develop the control term, roughly because the operator {H} is not diagonal in these variables. So, we introduce complex variables

ζ 1 = x 2 + ix 1 and ζ 2 = x 4 + ix 3 , (4.26)
and rewrite H, still denoting it with the same letter, as

H(ζ 1 , ζ 2 ) = - ω 1 2 |ζ 1 | 2 - ω 2 2 |ζ 2 | 2 . (4.27)
Using these complex variables the Poisson bracket with H takes the form

{H} = iω 1 ζ1 ∂ζ 1 -ζ 1 ∂ ζ 1 + iω 2 ζ2 ∂ζ 2 -ζ 2 ∂ ζ 2 . (4.28) Hence, for any n = (n 1 , n 2 ) ∈ N 2 and m = (m 1 , m 2 ) ∈ N 2 we obtain {H}ζ n ζ m = i (ω 1 m 1 -ω 1 n 1 + ω 2 m 2 -ω 2 n 2 ) ζ n ζ m = i ω • ( m -n)ζ n ζ m , (4.29)
where we introduce the vector ω = (ω 1 , ω 2 ) and use the compact notation

ζ n = ζ n 1 1 ζ n 2 2 for the complex vector ζ = (ζ 1 , ζ 2 )
. The operator {H} is diagonal in these variables and thus map 3.9 is straightforwardly obtained as

A -1 ζ n ζ m = e -{H} ζ n ζ m = e -i ω•( m-n) ζ n ζ m . (4.30)
Once we have this map, we can compute the operators G, N and R. For all n and m ∈ N 2 \ {0},

such that n = m and ω • ( m -n) = 2kπ ∀k ∈ Z (4.31)
(which defines the non-resonance condition), we get

Gζ n ζ m = 1 1 -e -i ω•( m-n) N ζ n ζ m , (4.32) with N ζ n ζ m = ζ n ζ m (4.33) if ω • ( m -n) = 2kπ ∀k ∈ Z (4.34)
and is 0 otherwise In the rest of the paper, we will assume the above non-resonant condition 4.31 to hold for the considered values of q x and q y . Let us remark that this is not a limitation of the actual theory, but just a working assumption. We could equivalently have chosen to work in the resonant regime, using a different control term suitable for resonant dynamics.

We are now able to explicitly compute the control term 3.14. Even under the non-resonance assumption, RV = 0, the computation is cumbersome. In fact, we need to express V in the complex variables, then compute S = GV , go back to the original variables, and finally (what is really the difficult task) compute the exponential e {GV } . Actually, even if GV is a polynomial of degree three in the variables x 1 , x 2 , x 3 , x 4 the map e {GV } is given by an infinite series. The terms of this series can be sequentially computed, but the degree of complexity (i.e. the number of involved terms) increase very fast. That's why we decided to use the approximated generator of the control map 3.25 truncated at order 2

F 2 = 1 2 {V }GV . (4.35)
We note that F 2 is composed by about twenty terms. A detailed discussion of the whole procedure needed to obtain F 2 , as well as its explicit formula, are presented in A.

We now face another difficulty, namely the computation of the control map from the generator F 2 . This is equivalent to performing the sum of

e {F 2 } = 1 + {F 2 } + 1 2 {F 2 } 2 + . . . , (4.36)
whose complexity once again grows very fast. We thus introduce a second approximation to our construction, by computing only a finite number of terms in the above sum. So, we define a truncated control map of order k

C k (F 2 ) = k l=0 {F 2 } l l! , (4.37)
and a truncated controlled map of order k

T k (F 2 ) = e {H} e {V } C k (F 2 ) = T C k (F 2 ) . (4.38)
Let us now discuss the properties of the truncated controlled map. Since the exact control map e {F 2 } given by 4.36 is symplectic, the controlled map 3.8 will also be symplectic. On the other hand, we cannot expect the k-th order control map C k (F 2 ) to be symplectic. We know that a map is symplectic if its Jacobian matrix A verifies (in its definition domain) the equality

A T JA -J = 0 , (4.39)
Thus, in order to check the symplecticity defect of T k (F 2 ) we compute the norm

D k of matrix A T k JA k -J, where A k is the Jacobian of T k (F 2 )
given by 4.38. The results are presented in Fig.

1 for orders k = 1 up to k = 5, in the region (x 1 , x 3 ) ∈ [-1, 1] × [-1, 1], x 2 = x 4 = 0. These results show that T k (F 2
) is a good approximation of a symplectic map for k ≥ 4, because we get D k 10 -4 for a large portion ( 53 %) of variables values. We note that in the central region of the truncated controlled map, where the actual physical process of beam's evolution occurs, the symplectic character of the map is established quite well since there D k 10 -8 . As expected, the larger the order k, the closer to symplecticity the approximation is. The main objective of the addition of a control term is to increase the size of the stability region around the central periodic orbit. This increase leads to decreasing the number of escaping orbits 4 , as we can see from the results presented in Fig. 2, where we plot in black the initial conditions on the square (x 1 , x 3 ) ∈ [-1, 1] × [-1, 1], x 2 = x 4 = 0, giving rise to orbits that do not escape up to 105 iterations of the map. In particular, we consider in Fig. 2(a) the original uncontrolled map 2.1, and in Figs. 2(b) to 2(d), the k order controlled map T k (F 2 ) for k = 1 to k = 5 respectively. One can easily see that the region of non-escaping orbits for the original map is smaller than the one of the controlled maps. This observation can be quantified by considering initial conditions inside a circle centered at the origin of each panel of Fig. 2 (which represent the actual physical plane since the initial momenta are x 2 (0) = x 4 (0) = 0) with radius r 2 = x 1 (0) 2 + x 3 (0) 2 , and evaluate the number of escaping and non-escaping orbits as a function of the circle radius for T k (F 2 ) with k = 1 up to k = 5. Results reported in Fig. 3 support the previous claim, by clearly showing that controlled maps of orders 3, 4 and 5 behave very similarly and lead to an increase of the nonescaping region. Let us note that the behavior of the controlled maps of orders k = 1 (Fig. 2(b)) and k = 2 (Fig. 2(c)) is somewhat misleading if it is not analyzed together with the information from the symplecticity defect (see Figs. 1(a) and (b) respectively). In fact, these maps are strongly dissipative and produce a strong shift of orbits towards the origin, preventing them from escaping. This dissipation effect is not physical, as it is not observed in real accelerators, and therefore we do not discuss further the k = 1 and k = 2 controlled maps. From the results of Figs. 2 and3 we see that the addition of even the lower order (k = 3) control term, having an acceptable symplecticity defect, increases drastically the size of the region of nonescaping orbits around the central periodic orbit. A further increase of the order of the control term results to less significant increment of this region, while the computational effort for constructing the controlled map increases considerably. In fact, T 1 (F 2 ) contains around 100 elementary terms, i.e. monomials in x 1 . . . x 4 , while this number is almost doubled for each order, so that T 5 (F 2 ) contains around 2000 terms. Also the CPU time needed to evolve the orbits increases with the order. For example, while the integration of one orbit using T 1 (F 2 ) takes about 1.4 times the CPU time needed to integrate the original map (2.1), the use of T 5 (F 2 ) needs almost 21.5 times more. Thus we conclude that the T 4 (F 2 ) controlled map, which can be considered quite accurately to be symplectic, is sufficient to get significant increment of the percentage of non-escaping orbits, without paying an extreme computational cost.

The SALI method

The Smaller Alignment Index (SALI) [START_REF] Ch | [END_REF] has been proved to be an efficiently simple method to determine the regular or chaotic nature of orbits in conservative dynamical systems. Thanks to its properties, it has been already successfully distinguished between regular and chaotic motion both, in symplectic maps and Hamiltonian flows [11,12,13,14,[START_REF] Manos | Proceedings of Semaine de l' Astrophysique Française Journées de la SF2A[END_REF]. For the sake of completeness, let us briefly recall the definition of the SALI and its behavior for regular and chaotic orbits, restricting our attention to 2N -dimensional symplectic maps. The interested reader can consult [START_REF] Ch | [END_REF] for a more detailed description. To compute the SALI of a given orbit of such maps, one has to follow the time evolution of the orbit itself and also of two linearly independent unitary deviation vectors v1 (0), v2 (0). The evolution of an orbit of a map T Figure 3: Percentages of non-escaping orbits for the controlled map T k (F 2 ) 4.38 as a function of the distance from the origin in the physical space (x 1 , x 3 ). We iterate initial conditions in a circle of radius r centred at the origin of plane (x 1 , x 3 ), with x 2 (0) = x 4 (0) = 0, and compute the percentages of non-escaping orbits during n = 10 5 iterations, for the controlled map T k (F 2 ) with k = 1, 2, 3, 4, 5, as a function of r.

is described by the discrete-time equations of the map

x (n+1) = T ( x (n) ) , (5.40)
where

x (n) = (x (n) 1 , x (n) 2 , ..., x (n)
2N ) T , represents the orbit's coordinates at the n-th iteration. The deviation vectors v

(n) 1 , v (n) 2
at time n are given by the tangent map

v (n+1) i = A( x (n) ) • v (n) i i = 1, 2 , (5.41)
where A denotes the Jacobian matrix of map 5.40, evaluated at the points of the orbit under study. Then, according to [START_REF] Ch | [END_REF] the SALI for the given orbit is defined as

SALI(n) = min v(n) 1 + v(n) 2 , v(n) 1 - v(n) 2 , (5.42)
where • denotes the usual Euclidean norm and vi = v i v i , i = 1, 2 are unitary normalised vectors. In the case of chaotic orbits, the deviation vectors v1 , v2 eventually become aligned in the direction defined by the maximal Lyapunov characteristic exponent (LCE), and SALI(n) falls exponentially to zero. An analytical study of SALI's behavior for chaotic orbits was carried out in [12], where it was shown that SALI(n) ∝ e -(σ 1 -σ 2 )n , (5.43) with σ 1 , σ 2 being the two largest LCEs. On the other hand, in the case of regular motion the orbit lies on a torus and the vectors v1 , v2 eventually fall on its tangent space, following a n -1 time evolution, having in general different directions. This behavior is due to the fact that for regular orbits the norm of a deviation vector increases linearly in time. Thus, the normalization procedure brings about a decrease of the magnitude of the coordinates perpendicular to the torus, at a rate proportional to n -1 , and so v1 , v2 eventually fall on the tangent space of the torus. In this case, the SALI oscillates about non-zero values (for more details see [11]).

The symplecticity of SALI's definition, its completely different behavior for regular and chaotic orbits, and its rapid convergence to zero in the case of chaotic motion are the main advantages that make SALI an ideal chaos detection tool. Recently a generalization of the SALI, the socalled Generalized Alignment Index (GALI) has been introduced [16,17], which uses information of more than two deviation vectors from the reference orbit. Since the advantages of GALI over SALI become relevant in the case of multi-dimensional systems, in the present paper we apply the SALI method for the dynamical study of the 4D map (2.1).

Dynamics of the controlled map

As already mentioned, the goal of constructing the controlled map T ctrl = T e {F } is to increase the percentage of regular orbits up to a given (large) number of iterations, or equivalently increase the size of the stability region around the nominal circular trajectory (i.e. the DA). Because the presence of chaotic regions can induce a large drift in the phase space, that eventually could lead to the escape of orbits, the achievement of a larger DA can be qualitatively inspected by checking via the SALI method the regular or chaotic nature of orbits in a neighborhood of the origin (see Fig. 4). We note that we define an orbit to be chaotic whenever SALI(t) < 10 -8 , and regular for the contrary. In Fig. 4 (which should be compared with Fig. 5 of [6]) we observe a strong enlargement of the region of regular orbits. This region is characterized by large SALI values. In particular, for 10 4 iterations of the T 4 (F 2 ) map, 54% of the considered orbits are regular, while for the uncontrolled map this percentage reduces to 33%. This improvement can be also confirmed by visual inspection of Fig. 2, where the regions of non-escaping orbits are shown for different orders of the controlled map (4.38). In Figs. 4(a) and (b), we see that there exist orbits of the T 3 (F 2 ) map, which are characterized as regular up to n = 10 4 iterations, while they show their chaotic character once they are iterated up to n = 10 5 . Such orbits correspond to the dark regions marked by a black circle in Fig. 4(b) (for comparison this circle is also plotted in all panels of Fig. 4). This discrepancy is absent for the T 4 (F 2 ) map, which shows almost the same geometrical shape for the non-escaping region when we pass from 10 4 to 10 5 iterations. In order to better understand this behavior we followed the evolution of a single orbit with initial condition x(0) = (-0.50, 0, -0.65, 0) T -inside the black circle in Fig. 4-for both the T 3 (F 2 ) and the T 4 (F 2 ) controlled maps, computing the corresponding SALI values up to 2 × 10 5 iterations. The results are reported in Fig. 5 and clearly show that the orbit behaves regularly up to n ≈ 10 5 iterations of the T 3 (F 2 ) map since its SALI values are different from zero, but later on a sudden decrease of SALI to zero denotes the chaotic character of the orbit. This behavior clearly implies this is a slightly chaotic, sticky orbit, which remains close to a torus for long time intervals (n ≈ 10 5 ), while later on it enters a chaotic region of the phase space. It is interesting to note that iterating the same initial condition by the T 4 (F 2 ) map we get a regular behavior at least up to n = 2 × 10 5 . In order to provide additional numerical evidence of the effectiveness of the controlled map (4.38) in increasing the DA, we consider initial conditions inside a 4D sphere centered at the origin

x 1 = x 2 = x 3 = x 4 = 0 of the map, with radius, r 2 = x 1 (0) 2 + x 2 (0) 2 + x 3 (0) 2 + x 4 (0) 2 .
We compute the number of regular, escaping and chaotic orbits as a function of the sphere radius. The corresponding results are reported in Fig. 6(b), while in Fig. 6(a) we reproduce Fig. 6 of [6] for comparison. From this figure we observe a strong increase of the DA, since the largest sphere containing 100% regular orbits has a radius r ≈ 0.66, while this radius was r ≈ 0.39 for the original uncontrolled map. We also observe that increasing the total number of iterations from but the radius of the 4D sphere containing only regular orbits does not change significantly. and chaotic [(r) red curves] orbits after n = 10 4 (dashed curves) and n = 10 5 iterations (solid curves) for initial conditions in a 4D sphere centred at the origin x 1 = x 2 = x 3 = x 4 = 0, as a function of the sphere radius r. Each point corresponds the average value over 5000 initial conditions. The largest radius at which the percentage of regular orbits is still 100%, is marked by an arrow in each panel.

Conclusions

In this paper we considered a simple model of a ring particle accelerator with sextupole nonlinearity that can be described by a symplectic map. In the framework of Hamiltonian control theory, we were able to control the dynamics of the original system by providing a suitable control map, resulting in a small "perturbation"of the initial map. This control map has been constructed with the aim of DA enlargement of the particle accelerator, and thus improving the beam's life-time and the accelerator's performance.

In particular, the theoretical framework we developed allows a 1-parameter family of approximated controlled maps. We performed several numerical simulations in order to choose "the best"approximated controlled map T k (F 2 ) (4.38), taking into account the complexity of the map, i.e. the number of terms by which it is composed, the CPU time needed to perform the numerical iteration of orbits, and the accuracy of the results in terms of the symplectic character of the map. We find that the 4 th order controlled map T 4 (F 2 ) is an optimal choice for the controlled system. Using this controlled map we succeeded in achieving our initially set goal, since the T 4 (F 2 ) map exhibits a DA with a radius more than 1.7 times larger than the one for the original map (see Fig. 6).

A Computation of the control term.

The aim of this section is to introduce further details for the construction of the control term and of the controlled map, and to provide explicit formulas for the interested reader.

Notations

Let us first introduce some notations and recall some useful relations.

• Lie brackets and operators. Let X be the vector space of C ∞ real or complex functions of 2N variables (p, q). For any F, G ∈ X , the Lie bracket is given by

{F, G} := N i=1 [∂ p i F ∂ q i G -∂ q i F ∂ p i G] , (1.44)
where ∂ x i f ≡ ∂f ∂x i denotes the partial derivative with respect to the variable x i . Using the above definition, we can define a linear operator, induced by an element F of X , acting on X {F } : X → X G → {F }G := {F, G}

This operator is linear, antisymmetric and verifies the Jacobi identity

∀F, G ∈ X {{F }G} = {F }{G} -{G}{F } (1.46)
• Exponential. We define the exponential of such an operator {F }, by

e {F } := ∞ k=0 {F } k k! , (1.47)
which is also an operator acting on X . The power of an operator is the composition :

{F } k G = {F } k-1 ({F }G).
We observe that in the case of the Hamiltonian function H, the exponential provides the flow, namely e {H} x 0 = x(t), of the Hamilton equations ṗ = -∂ q H (1.48) q = ∂ p H .

(1.49)

• Vector Field. The action of the above defined operators, can be extended to vector fields "component by component"

∀ F, G, H ∈ X {F } G H := {F }G {F }H . (1.50)
1.2 Mappings as time-1 flows.

We show now that map (2.1) can be seen as the time-1 flow of a given Hamiltonian system. More precisely we show that

T     x 1 x 2 x 3 x 4     =     cos(ω 1 ) -sin(ω 1 ) 0 0 sin(ω 1 ) cos(ω 1 ) 0 0 0 0 cos(ω 2 ) -sin(ω 2 ) 0 0 sin(ω 2 ) cos(ω 2 )         x 1 x 2 + x 2 1 -x 2 3 x 3 x 4 -2x 1 x 3     = e {H} e {V }     x 1 x 2 x 3 x 4     (1.51)
where

H(x 1 , x 2 , x 3 , x 4 ) = -ω 1 x 2 1 + x 2 2 2 -ω 2 x 2 3 + x 2 4 2 (1.52) and V (x 1 , x 2 , x 3 , x 4 ) = - x 3 1 3 + x 1 x 2 3 . (1.53)
Let us observe that H is the sum of two non-interacting harmonic oscillators with frequencies ω 1 and ω 2 , hence its dynamics is explicitly given by

x 1 (t) = A cos(ω 1 t) -B sin(ω 1 t) (1.54) x 2 (t) = B cos(ω 1 t) + A sin(ω 1 t) (1.55) x 3 (t) = C cos(ω 2 t) -D sin(ω 2 t) (1.56) x 4 (t) = D cos(ω 2 t) + C sin(ω 2 t) .
(1.57)

By definition y = e {H} x is the solution at time 1 with initial condition x = (x 1 , x 2 , x 3 , x 4 ) T , hence we obtain

y 1 = cos(ω 1 )x 1 -sin(ω 1 )x 2 (1.58) y 2 = sin(ω 1 )x 1 + cos(ω 1 )x 2 (1.59) y 3 = cos(ω 2 )x 3 -sin(ω 2 )x 4 (1.60) y 4 = sin(ω 2 )x 3 + cos(ω 2 )x 4 , (1.61)
that is

e {H}     x 1 x 2 x 3 x 4     =     cos(ω 1 ) -sin(ω 1 ) 0 0 sin(ω 1 ) cos(ω 1 ) 0 0 0 0 cos(ω 2 ) -sin(ω 2 ) 0 0 sin(ω 2 ) cos(ω 2 )         x 1 x 2 x 3 x 4     .
(1.62) From 1.53 and the definition 1.45 we easily get

{V } := ∂ x 2 V ∂ x 1 -∂ x 1 V ∂ x 2 + ∂ x 4 V ∂ x 3 -∂ x 3 V ∂ x 4 (1.63) = (x 2 1 -x 2 3 )∂ x 2 -2x 1 x 3 ∂ x 4 .
This means that once applied to a vector x only the second and fourth components of {V } x will be non-zero and moreover they only depend on the first and third components of x, hence {V } 2 x = 0. We can thus conclude that ∀ k ≥ 2 and ∀ x ∈ R 4 , we get {V } k x = 0. Finally using the definition (1.47) we obtain Using the warped addition with equation 3.14 of Theorem 1, we obtain e {F } = e -{V } e {(N -G)V } e {RV } e {GV} = e -{V }⊕{(N -G)V }⊕{RV }⊕{GV} , (1.67) and thus

e {V } x = ∞ k=0 {V } k k! x = Ix + {V } x =     x 1 x 2 + x 2 1 -x 2 3 x 3 x 4 -2x 1 x 3     (1.64)
{F } = -{V } ⊕ {(N -G)V } ⊕ {RV } ⊕ {GV} = -{V } ⊕ {(1 -G)V } ⊕ {GV } + o(V 2 ) =     1 2 ({V }{GV } -{GV }{V }) (1.46) = {{V }GV } -{GV }     ⊕ {GV } + o(V 2 ) = 1 2 {{V }GV } + 1 4 {{{V }GV }GV } =o(V 2 ) +o(V 2 ) = 1 2 {{V }GV } + o(V 2 ) , (1.68)
where we explicitly used the assumption RV = 0 to remove the third term on the the right hand side on the first equation and hence to write N V = V . We are thus able to define the non-resonant control term, up to order V 2 , to be

F 2 = 1 2 {V }GV = 1 2
{V, GV } .

(1.69)

The operator G

To get the explicit formula for To construct it, it will be more convenient to use complex variables

ζ 1 = x 2 + ix 1 and ζ 2 = x 4 + ix 3 . (1.71)
Then the function H becomes

H(ζ 1 , ζ 2 ) = - ω 1 2 ζ 1 ζ1 - ω 2 2 ζ 2 ζ2 , (1.72)
and using

∂ ∂x 1 = ∂ζ 1 ∂x 1 ∂ ∂ζ 1 + ∂ ζ1 ∂x 1 ∂ ∂ ζ1 = i ∂ ∂ζ 1 -i ∂ ∂ ζ1 (1.73) ∂ ∂x 2 = ∂ζ 1 ∂x 2 ∂ ∂ζ 1 + ∂ ζ1 ∂x 2 ∂ ∂ ζ1 = ∂ ∂ζ 1 + ∂ ∂ ζ1 (1.74) for (x 1 , x 2 ), the operator {H} becomes ∂ x 2 H∂ x 1 -∂ x 1 H∂ x 2 = 2i ∂ ζ 1 H∂ζ 1 -∂ζ 1 H∂ ζ 1 = iω 1 ζ1 ∂ζ 1 -ζ 1 ∂ ζ 1 , (1.75)
with a similar expression holding for (x 3 , x 4 ). Hence for any n = (n 1 , n 2 ) ∈ N 2 and m = (m 1 , m 2 ) ∈ N 2 we obtain

{H}ζ n ζ m = i (ω 1 m 1 -ω 1 n 1 + ω 2 m 2 -ω 2 n 2 ) ζ n ζ m = i ω • ( m -n)ζ n ζ m , (1.76)
where we introduced the vector ω = (ω 1 , ω 2 ) and we used the compact notation

ζ n = ζ n 1 1 ζ n 2 2 , for the complex vector ζ = (ζ 1 , ζ 2 ).
We note that from the knowledge of the operators' action on such monomials ζ n ζ m , we can reconstruct the operator action on any regular function. The operators are linear and they will be applied on polynomials in the x variable, which are nothing more than polynomials in the complex variables. Let us now compute the time-1 flow of {H} by using complex variables. Starting from 1.76 and then proceeding by induction, we can easily prove that for all k ∈ N The function F 2 is defined by 1.69, where V is a known function. The term GV will be computed starting from the previously obtained expression of G. To construct F 2 we first have to express V in the complex variables (1.71):

{H} k ζ n ζ m = (i ω • ( m -n)) k ζ n ζ m , ( 1
V (ζ 1 , ζ 2 ) = - 1 24 iζ 3 1 + 1 8 i ζ1 ζ 2 1 - 1 8 i ζ2 1 ζ 1 + 1 24 i ζ3 1 + 1 8 iζ 1 ζ 2 2 - 1 4 iζ 1 ζ 2 ζ2 + 1 8 iζ 1 ζ2 2 - 1 8 i ζ1 ζ 2 2 + 1 4 i ζ1 ζ2 ζ 2 - 1 8 i ζ1 ζ2 2 .
(1.89)

By the linearity of the operator, and by using (1.80), we easily compute GV . In particular we apply G to each term of (1.89). Then using the inverse change of coordinates ) cos(ω 1 ) -x 1 x 2 sin(ω 1 )) -1/6 csc(3/2ω 1 ) x 2 cos(1/2ω 1 ) + x 1 sin(1/2ω 1 ) (2x 2 + 2x 2 cos(ω 1 ) -x 1 sin(ω 1 )) -1/4 (x 2 3 -x 2 4 ) sin(ω 2 ) cos(ω 2 ) -cos(ω 1 + ω 2 ) + 1/4

x
x 1 x 3 (2x 1 x 3 -2x 2 x 4 ) sin(ω 2 ) cos(ω 2 ) -cos(ω 1 + ω 2 ) .

(1.92)

Figure 1 :

 1 Figure 1: The simplicity defect of the controlled map T k (F 2 ) 4.38. Plot of log 10 D k , where D k = ||A T k JA k -J|| and A k is the Jacobian of the k-order controlled map T k (F 2 ) given by 4.38, for 16000 uniformly distributed values in the square (x 1 , x 3 ) ∈ [-1, 1] × [-1, 1], x 2 = x 4 = 0, for (a) k = 1, (b) k = 2, (c) k = 3, (d) k = 4 and (e) k = 5. The percentage of orbits with log 10 D k < -4 is 0.5%, 11%, 28%, 53% and 73% for k = 1, 2, 3, 4 and 5 respectively. The gray scale corresponds to the value of log 10 D k : the darker the color, the smaller the value of log 10 D k is, and hence the closer the map is to a symplectic one.

Figure 2 :

 2 Figure 2: Non-escaping regions of controlled map T k (F 2 ) as a function of the truncation order k. 16000 uniformly distributed initial conditions in the square (x 1 , x 3 ) ∈ [-1, 1] × [-1, 1], x 2 (0) = x 4 (0) = 0 are iterated up to n = 10 5 using (a) the uncontrolled map 2.1 and (b)-(f) the k = 1 to k = 5 order controlled map T k (F 2 ) 4.38, respectively. Initial conditions corresponding to non-escaping orbits up to n = 10 5 are coloured in black, while escaping orbits are color ed in white.

Figure 4 :

 4 Figure 4: Stability analysis in the (x 1 , x 3 ) plane. 16000 uniformly distributed initial conditions in the square (x 1 , x 3 ) ∈ [-1, 1] × [-1, 1], x 2 (0) = x 4 (0) = 0 are integrated using the T 3 (F 2 ) ((a) and (b)), and the T 4 (F 2 ) controlled map ((c) and (d)), up to n = 10 4 ((a) and (c)) and n = 10 5 iterations ((b) and (d)). The gray scale represents the value of log 10 SALI for each orbit at the end of the integration time. The lighter the color the more stable is the orbit, while white color denotes that an orbit escaped before the total number of iterations was reached. The black circle indicate the initial condition of the orbit studied in Fig. 5.

Figure 5 :

 5 Figure 5: (Color online) Dynamics of two orbits with the same initial conditions for the 3 rd and 4 th order controlled maps. Time evolution of the SALI for the orbit with initial conditions x(0) = (-0.50, 0, -0.65, 0) T (see Fig. 4), using the T 3 (F 2 ) [(r) red curve] and the T 4 (F 2 ) controlled map [(b) blue curve].

Figure 6 :

 6 Figure 6: (Color online) Dynamical aperture of the (a) original map (2.1) and (b) the T 4 (F 2 ) controlled map (4.38). The percentages of regular [(b) blue curves], escaping [(g) green curves]and chaotic [(r) red curves] orbits after n = 10 4 (dashed curves) and n = 10 5 iterations (solid curves) for initial conditions in a 4D sphere centred at the origin x 1 = x 2 = x 3 = x 4 = 0, as a function of the sphere radius r. Each point corresponds the average value over 5000 initial conditions. The largest radius at which the percentage of regular orbits is still 100%, is marked by an arrow in each panel.

1. 3 2 (

 32 Computation of the generator F under the assumption RV ≡ 0. Let us recall that the composition of maps expressed by exponential defines the warped addition e {A} e {B} := e {A}⊕{B} , (1.65) whose first terms are {A} ⊕ {B} = {A} + {B} + 1 {A}{B} -{B}{A}) + . . . (1.66)

F 2

 2 we need to compute the expression of G. From definition 3.10 the operator G should satisfy G 1 -e -{H} G = G . (1.70)

e 1 - 1 1 . 5

 1115 {H} ζ n ζ m = ∞ k=0 (i ω • ( m -n)) k k! ζ n ζ m = e i ω•( m-n) ζ n ζ m . (1.78) Similarly e -{H} ζ n ζ m = e -i ω•( m-n) ζ n ζ m . Assuming a non-resonance condition ω • ( m -n) = 2kπ ∀ n = m ∈ N 2 \ {0} and ∀k ∈ Z ,(1.79)a possible choice for the operator G is the followingGζ n ζ m := 1 1 -e -i ω•( m-n) N ζ n ζ m (1.80) with N ζ n ζ m = ζ n ζ m (1.81) if ω • ( m -n) = 2πk (1.82)and is 0 otherwise It is easy to check that the operator G defined by (1.80) verifies (1.70). In order to do so we introduce the compact notationZ n,m := ζ n ζ m and n,m := e -i ω•( m-n) .(1.83)Developing the left hand side of 1.70 and using the linearity of all operators, we getG 1 -e -{H} GZ n,m = G -Ge -{H} 1 n,m ) 2 Z n,m -n,m(1n,m ) 2 Z n,m (Expression of the control term F 2

  (similar expressions hold for (x 3 , x 4 ) and (ζ 2 , ζ2 )), we can go back to the original variables x. The obtained expression after some algebraic simplifications isGV = -1/6 csc (3/2 ω 1 ) [x 2 cos (1/2 ω 1 ) + x 1 sin (1/2 ω 1 )] x 1 2 -3 x 3 2 + x 2 -ω 2 ) -cos (ω 2 ) -1/4 (-x 2 x 4 2 + x 2 x 3 2 + 2 x 1 x 4 x 3 ) sin (ω 2 ) cos (ω 2 ) -cos (ω 1 + ω 2 ) (1.91)Then the explicit expression of the control term F 2 is

		1 =	1 2	i( ζ1 -ζ 1 ) and x 2 =	1 2	(ζ 1 + ζ1 )	(1.90)
							2
	+ 2 x 1	2 -6 x 3	2 + x 2	2 cos (ω 1 ) -x 1 x 2 sin (ω 1 )
	sin (ω 2 ) 3 ) -1/6 csc(3/2ω 1 ) cos(1/2ω 1 )(x 2 +1/4 cos (ω 1 F 2 = 1/2(x 2 1 -x 2 1 -3x 2 3 + x 2 2
	+(2x 2 1 -6x 2 3 + x 2 2

Such parameters have been fixed throughout this work to the values q x = 0.61803 and q y = 0.4152, corresponding to a non-resonant condition (see[5]).

In the literature one can sometimes find the alternative equivalent notation {H} = L H .

In[8] a similar theory has been developed for a general symplectic map devoid of an infinitesimal generator.

an orbit (x 1 (k), x 2 (k), x 3 (k), x 4 (k)) 0≤k≤N is defined as non-escaping if for all k ≤ N x 2 1 (k) + x 2 2 (k) + x 2 3 (k) + x 2 4 (k) ≤ R 2for a certain R (in the simulations we used R 2 = 10 and N = 10 4 , 10

) and escaping otherwise.

Acknowledgements

Numerical simulations were made on the local computing resources (Cluster URBM-SYSDYN) at the University of Namur (FUNDP, Belgium).