Victor Lomné

Amine Dehbaoui

Philippe Maurine

Michel Robert

Lionel Torres
email: lionel.torres@lirmm.fr

A Dehaboui

Side Channel Attacks

à la diffusion de documents scientifiques de 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

In the past 100 years, we have seen the emergence of modern cryptography, along with many cryptographic primitives and protocols. The development of new theoretical cryptanalytic techniques to try to defeat the main cryptographic algorithms has increased knowledge of how to design cryptographic primitives and schemes.

In theoretical cryptanalysis today, a cryptographic algorithm is considered as a black box. Even when attackers know the cryptographic algorithm, they only have access to pairs of plaintexts/ciphertexts, and their goal is to guess the cryptographic key.

The robustness of modern cryptographic algorithms is based on these assumptions, and a cryptanalytic attack is a method that allows an assailant to guess the key with a complexity (in time and/or in memory) lower than a brute force attack.

But if a cryptographic algorithm is modeled as a gray box, i.e. in such a way that an attacker can obtain intermediate information during the cryptographic operation requiring the key, these assumptions no longer hold true.

In this context, and with the increasing use of embedded cryptographic devices with embedded cryptographic secrets, several tamper attacks have appeared since the beginning of the 1990s. These enable an attacker to obtain intermediate information during the cryptographic operation, and then to deduce the secret key with a complexity much lower than classical theoretic cryptanalytic attacks.

Figure 3.1 shows different models from theoretical and hardware cryptanalytic points of view. In the gray box model, physical leakages from the device can be identified or the computation running on the device can be disturbed. As we explain in this chapter, these different approaches enable cryptographic secrets to be extracted. Cryptanalytic techniques for hardware are classified as invasive attacks, semi-invasive attacks, or non-invasive attacks. We describe these attacks in detail in the following sections.

Invasive Attacks

Invasive attacks are a tamper attack in which the device is completely destroyed in order to extract secret information. The best known invasive attack is hardware reverse engineering, as described in [START_REF] Torrance | The state-of-the-art in IC reverse engineering[END_REF]. The goal is to retrieve the layout of the circuit using chemistry techniques and/or high resolution microscopy. This technique is often used for cloning or anti-cloning purposes.

The first step consists in decapsulating the chip. A high resolution picture is then taken of each metal layer, and chemistry techniques are used to remove each succeeding metal layer to penetrate deeper into the chip right down to the transistor level. When attackers have a picture of each metal layer, they can use dedicated tools to go back to the netlist of the circuit. Finally, they obtain the behavioral description of the chip.

This attack has been conducted in the case of the MyFare stream encryption algorithm, which was kept secret by NXP. An invasive attack revealed the algorithm's functionality [START_REF] Nohl | Reverse-engineering a cryptographic RFID tag[END_REF]. The method is shown in Fig. 3.2.

Although this technique is very powerful, but it has two drawbacks. Firstly, it requires very expensive equipment and highly qualified engineers. Secondly, thanks to technology shrinking techniques, it requires more and more sophisticated microscopy. In the case of SRAM and Flash based FPGAs, the configuration file, also called the bitstream, is stored in a dedicated off or on chip Flash memory. Thus, considering that the layout of an FPGA is roughly an array of reconfigurable cells, invasive attacks that allow attackers to obtain the hardware description of the attacked chip still provide no information about the configuration file, because the reconfigurable cells lose their configuration data as soon as the FPGA is switched off.

For instance, Fig. 3.3 shows the layout of an FPGA; as can be seen, optical observation of the floorplan provides no clues about the FPGA's configuration.

Semi-invasive Attacks

Semi-invasive attacks consist in retrieving physical information or disturbing the cryptographic operation computed by the device after decapsulation. These kinds of attacks are called semi-invasive because they require physical modification of the device, but do not destroy it. Indeed, attackers usually uses chemistry techniques to remove the package of the chip. In some cases, other sophisticated techniques may be required for example thinning down the substrate. For this purpose, different approaches can be used.

The first approach belongs to the fault attacks category. One or several faults are generated during the cryptographic operation using, for instance, light pulses. A concrete light based fault attack was first reported in [START_REF] Skorobogatov | Optical fault induction attacks[END_REF]. Since these first concrete results, more sophisticated techniques have been elaborated using a laser (for instance red or infra-red diode laser), allowing the fault to be generated on either the front or back of the chip, and/or the fault to be focused on one or a few bits [START_REF]Diode Laser Station[END_REF]. Figure 3.4 shows a commercial laser attack setup from Riscure (on the left), and zoomed view of memory with the tiny laser pulse in pink (on the right).

Another technique, called probing attack [START_REF] Handschuh | Probing attacks on tamper-resistant devices[END_REF], consists in placing a small needle on a bus line of the microcontroller to retrieve information. Different authors have shown that even if attackers can only spy on one random bus line, they can nevertheless obtain enough information to reveal a cryptographic key. As an example, the authors of [START_REF] Schmidt | A probing attack on AES[END_REF] describe a probing attack applied on AES. Recently, a general side-channel method based on single bit probing [START_REF] Dinur | Generic analysis of small cryptographic leaks[END_REF] was introduced as a way to optimize cube attacks [START_REF] Dinur | Side channel cube attacks on block ciphers[END_REF]. Such attacks use a learning step, and can consequently even work on a proprietary algorithm.

A recent kind of semi-invasive attack consists in measuring, on the back side, photons emitted by transistors when they switch. This attack [START_REF] Di-Battista | When failure analysis meets side-channel attacks[END_REF] requires thinning down the substrate and a high performance camera. Then, after the acquisition of photons emitted by the chip during the cryptographic computation of different data and the same unknown key, a statistical treatment, similar to that used in SCA reveals the cryptographic key. Moreover, this powerful technique follow the displacement of the current to be monitored by computing a movie disclosing the switching of transistors. However, this attack requires very expensive equipment that is available only in a few advanced laboratories.

Non-invasive Attacks

Non-invasive attacks are powerful techniques to extract secret information from cryptographic devices without physical modification of the device. As explained above, we distinguish two approaches, one consists in identifying physical leakages during the cryptographic operation (called Side Channel Attacks (SCA)), the other in disturbing computation (called Fault Attacks (FA)).

Side Channel Attacks consist in measuring physical leakages from the device during a cryptographic operation. The main idea behind this approach is that CMOS technology has an interesting property from a cryptanalytic point of view: it leaks physical information correlated with processed data. This class of attacks is very powerful, and requires only affordable measuring and testing equipment. Different physical leakages can be measured. Typical leaks are the computational time of the cryptographic operation [START_REF] Kocher | Timing attacks on implementations of Diffie-Hellman, RSA, DSS, and other systems[END_REF], power consumption by the device [START_REF] Kocher | Differential power analysis[END_REF], or the electromagnetic radiations emitted by the chip [START_REF] Gandolfi | Electromagnetic analysis: concrete results[END_REF].

In some cases secret information can be extracted via a physical leakage from only one cryptographic operation, in others the attacker needs to record physical leakages from several cryptographic operations and to apply a statistical treatment on these records. On a classical smart card without dedicated countermeasures, several thousand records of cryptographic operations enable both symmetric ciphers (like DES or AES) and asymmetric ciphers (like RSA or ECC) to be broken in a very short computational time (a few hours using a desktop computer).

The first type of Side-Channel Attacks is the Timing Attack, which exploits differences in the computational time of a cryptographic operation. But from a practical point of view, methods have been proposed to implement a cryptographic algorithm in both hardware or software, ensuring that its computational time is constant.

However, power and electromagnetic side channels are physical leakages that are tightly correlated with the processed data due to CMOS properties, and these kinds of attacks are considered to be very powerful by manufacturers and government agencies. Moreover, using a tiny electric or magnetic sensor that is smaller than the chip, a 3D electromagnetic radiation map can be computed. This is achieved by running the same set of instructions several times, while placing the sensor in different positions.

The second approach belongs to the fault attacks category. The initial idea, from [START_REF] Boneh | On the importance of checking cryptographic protocols for faults (extended abstract)[END_REF], is to run the cryptographic operation twice, one safe and another faulty. Dedicated cryptanalytic techniques enable one pair of safe/faulty cryptographic computations to be exploited to extract the key, with a low complexity in comparison to classical cryptanalytic techniques. In other cases a differential approach is required to retrieve the key [START_REF] Biham | Differential fault analysis of secret key cryptosystems[END_REF], in which case the attacker needs several pairs of safe/faulty cryptographic computations to retrieve the cryptographic key.

The first concrete way proposed in the literature to generate faults during a computation is the glitch attack [START_REF] Anderson | Low cost attacks on tamper resistant devices[END_REF]. The idea is to modify the supplying signal of the device during a sharp time shorter than a clock period. The device can be overpowered or underpowered. As a result, a setup time or hold time violation will occur, and the attacker generates a logical error in the result of the cryptographic operation. A variant consists in suddenly modifying the clock period.

There is another way to generate faults in cryptographic devices without physically modifying the circuit. The idea is to place an electromagnetic sensor above the device and to generate an electromagnetic field through the chip. Although few works have reported concrete results of this kind of attack, concrete experimental results can be found in [START_REF] Schmidt | Optical and EM fault-attacks on CRT-based RSA: concrete results[END_REF].

Some circuits can even be disturbed by a steady stress. For instance, paper [START_REF] Khelil | Fault analysis attack on an FPGA AES implementation[END_REF] shows how to obtain incorrect computations from an AES implemented in an FPGA.

The method consists in underfeeding the device, so as to generate setup time violations. A comprehensive study of this way of injecting such faults in an FPGA is available in [START_REF] Bhasin | Security evaluation of different AES implementations against practical setup time violation attacks in FPGAs[END_REF] for various FPGA families. The outcome of these experiments is that the description of the algorithm can greatly improve the resistance of the implementations against these "global faults".

In the rest of this chapter, we only deal with Power and Electromagnetic based Attacks will be detailed [START_REF] Standaert | Report on DPA and EMA Attacks on FPGAs[END_REF]. For more information on Timing Attacks, the reader can consult [START_REF] Kocher | Timing attacks on implementations of Diffie-Hellman, RSA, DSS, and other systems[END_REF] or [START_REF] Brumley | Remote timing attacks are practical[END_REF]. For information on Fault Attacks, [START_REF] Biham | Differential fault analysis of secret key cryptosystems[END_REF][START_REF] Boneh | On the importance of checking cryptographic protocols for faults (extended abstract)[END_REF][START_REF] Giraud | A survey on fault attacks[END_REF] and [START_REF] Bar-El | The sorcerer's apprentice guide to fault attacks[END_REF] are good starting points.

The following explanations are applicable to for all kind of ICs used in embedded devices, and especially for FPGAs that use hardware or software cryptographic primitives.

Power and Electromagnetic Measurement Platform

In this section, we describe the platform used to perform power and/or electromagnetic measurements of an FPGA running cryptographic operations.

Equipment

One example of an overall platform is given in Fig. 3.5, it is made up of several parts described in the following sections. Other boards are described in [START_REF] Standaert | Report on DPA and EMA Attacks on FPGAs[END_REF]. Also of particular interest are the so-called "Side Channel Analysis Standard Evaluation BOards" (SASEBO) [START_REF] Satoh | Side-channel Attack Standard Evaluation Board[END_REF].

Attacked Device

For the experiments described in this chapter, the attacked device is an Xilinx FPGA Spartan3-1000 embedded in a Digilent board. The Spartan die is encapsulated in a Cavity Up Ball Grid Array (BGA) package, and is roughly 7 mm by 7 mm in size, whereas the whole package is roughly 17 mm by 17 mm.

Oscilloscope

The main component of the equipment required for a Side Channel Attack, (and the most expensive device in a Side Channel Attack platform) is an oscilloscope. A Lecroy 735Zi scope was used for the experiments described in this chapter. It belongs to the latest generation of Lecroy scopes, and has impressive properties: Moreover, like other high-performance scopes, it allows remote control via Ethernet connections. Thus, once the scope is on, it can be fully controlled remotely. Furthermore, specific advanced modes for the trigger can be very useful for SCA purposes, but these functionalities were not exploited in our experiments.

Current Probes

To accurately measure the current consumed by the attacked device, current probes have interesting properties and are affordable. Tektronix current probes CT1 and CT6 [1] were used in the platform described here. First, the voltage regulator on the Digilent board supplying the power to the FPGA core was removed. Power was supplied by a battery, and the current probe was placed on the VDD or GND wire (between the battery and the pins of the core).

The current probe CT1 has a frequency bandwidth of 1 GHz and the CT6 of 2 GHz, and an accuracy of 5 mV/mA. Low-noise amplifier (described below) was plugged at the output of the current probe to increase the amplitude of the signal enabling accurate measurement of power consumption.

Magnetic Sensors

The size and the shape of a magnetic sensor are the most important characteristics that determine the accuracy of the sensor. Hand-made magnetic sensors can be built. Attempts have often been made to build the smallest coiled wire, made of copper, and with several turns. Although these sensors produce good results (for instance some DEMA attacks succeeded with several hundred measurements), commercial magnetic sensors are often preferred [START_REF] Sauvage | Cartographie électromagnétique pour la cryptanalyse physique[END_REF].

Consequently, we used the Rhode and Schwarz HZ-15 probe set which contains different electric and magnetic sensors. Among them, the H field probe RSH2.5-2 gives accurate measurements, with good precision (its diameter is around 0.5 mm, thus when placed close to the IC, one is about 1 mm when working on encapsulated chips). Moreover, it allows precise measurements with a huge frequency bandwidth (spanning from 1 MHz to 3 GHz) even through the package of the FPGA. It consequently appears to be a good compromise to measure the magnetic activity of packaged chips.

Low-Noise Amplifier

To ensure sufficient power and accurate electromagnetic measurements, a low noise amplifier can be used. We used one with a frequency bandwidth of 1 GHz, and a gain of 63 dB.

Motorized X-Y -Z Stage

Using a tiny magnetic sensor to measure a particular area of the attacked chip requires positioning the sensor at the exact location required. For this task, we used a motorized X-Y -Z stage. The three axes have a precision of 50 µm, and the stage can be fully controlled remotely by serial connection.

Computer

A personal computer was used to control the whole platform. It was connected to the board embedding the attacked device via one serial port, to the scope via Ethernet, and to the X-Y -Z motorized stage via another serial port to drive it.

Acquisition Software

To automate the acquisition of measurements, we developed a software in Matlab. It runs on the computer controlling the whole platform. Typically, to perform an EM cartography, i.e. to obtain EM measurements of the same calculation(s) with the magnetic sensor placed at different positions above the chip to draw a map, the following steps are recommended:

• the initial position (X init , Y init , Z init ), the number of positions (nb X , nb Y ), the displacement step (dis), the key (K) and the number of plaintexts (nb PTIs ) are chosen by the user; • once the script is launched, PTIs are chosen randomly;

• the motorized stage is positioned at the first position (X init , Y init , Z init );

• the computer sends to the scope and adjusts it so that the entire computation is correctly recorded; it also sends the key K (optional) and the PTI 1 to the chip; the chip sends to the scope a trigger signal on the channel 1 to launch the recording of the magnetic signal on channel 2 (the magnetic sensor is connected to the low noise amplifier, which is connected to channel 2 of the scope); • once the computation is finished, the scope sends the acquired EM trace to computer; • the two previous steps are repeated nb PTIs times;

• the motorized stage is positioned at the second position (X init + dis, Y init , Z init ), and the three previous steps are repeated; • the same steps are repeated until all the positions have been covered.

The localization of the cryptographic modules is an advanced topic that is beyond the scope of this introduction concerned with side channel attacks. Interested readers should refer to [START_REF] Dehbaoui | Magnitude squared incoherence EM analysis for integrated cryptographic module localisation[END_REF][START_REF] Réal | Enhancing correlation electromagnetic attack using planar near-field cartography[END_REF][START_REF] Sauvage | Electromagnetic radiations of FPGAs: high spatial resolution cartography and attack of a cryptographic module[END_REF][START_REF] Sauvage | Cross-correlation cartography[END_REF].

Leakage Models

Power and Electromagnetic (EM) based Side Channel Attacks is a hardware cryptanalytic technique that exploits physical leakages emitted by a cryptographic device. In this section, we explain how physical measurements can be realized and how leakage model can be used.

Power Consumption Leakage

Complementary Metal-Oxide-Semiconductor (CMOS) is the most used technology in microprocessors, microcontrollers, static RAM, and other digital logic circuits. Two important characteristics of CMOS devices are high noise immunity and low static power consumption. Significant power is only drawn while the transistors in the CMOS device are switching between on and off states. Consequently, CMOS devices do not produce as much waste heat as other forms of logic, like NMOS. CMOS also enables a high density of logic functions on a chip. It was primarily for these reasons that CMOS won the race in the 1980s and became the most widely used technology for VLSI chips.

The power consumption of a CMOS circuit is the sum of the power consumption of the logic cells that make up the circuit [START_REF] Guilley | Differential power analysis model and some results[END_REF]. Hence, the total power consumption mainly depends on the number of logic cells in a circuit, the connections between them, and how the cells are built. The power consumption of a CMOS gate is generally considered in terms of two components [START_REF] Coron | Statistics and secret leakage[END_REF]:

• The dynamic power component: this is mainly related to the charging and discharging of the load capacitance at the gate output, but also to the short circuit current that flows during the transition of the input from one voltage level to another. Indeed there is a short period during which both PMOS and NMOS transistors are on simultaneously, thus creating an electrical path between the VSS and VDD rails. • The static power component: this is due to leakage in the substrate that flows even when the gate is not switching. In turn, this leakage is made up of several components including gate to source leakage, which flows directly through the gate insulator, mostly by tunneling, and source-drain leakage attributed to both tunneling and sub-threshold conduction. However, so far no research paper has reported an attack using this static leakage.

While the transistors that comprise the CMOS gate are switching between on and off states, a significant amount of power is consumed, due to the dynamic power consumption. This characteristic of the CMOS technology is interesting from a cryptanalytic point of view. Indeed, an attacker can use the fact that the dynamic power consumption is tightly correlated with the number of switching bits to guess a secret value in a cryptographic operation.

Moreover, electrical simulations in recent works [START_REF] Lin | Analysis and mitigation of process variation impacts on power-attack tolerance[END_REF] showed that, in deepsubmicron CMOS technologies, the percentage of static power consumption among the global power consumption of a circuit increases with a decrease in the technology shrinking techniques used. As the static power consumption of a gate mainly depends on the biasing of its inputs and its internal nodes, static power consumption is correlated with the last value computed by this gate. An attacker can thus exploit the relation between static power consumption and a computed value to guess a secret value in a cryptographic computation [START_REF] Quisquater | Physically secure cryptographic computations: from micro to nano electronic devices[END_REF].

Different practical techniques can be used to measure the power consumption of a chip, depending mainly on how the chip is encapsulated, and how it is supplied:

• if a small resistor is placed between the VDD (GND) of the circuit and the VDD (GND) of the power supply, the consumption of the circuit can be measured at the pins of the resistor; • another way is to use a dedicated current probe [1].

Electromagnetic Leakage

EM radiations emitted by an integrated circuit are mainly due to the displacement of current through the rails of the metal layers. This phenomenon has been formalized by Maxwell's equations.

Practical experiments in [START_REF] Ordas | Near-field mapping system to scan in time domain the magnetic emissions of integrated circuits[END_REF] showed that the power and ground (P/G) network represent the majority of EM radiations of a CMOS device. Thus, when transistors are switching between on and off states, current crosses the P/G network to supply the pins of the switching transistors, and this displacement of current creates variations in the EM field induced by the electrical behavior of the circuit.

Practically speaking, either the electric or the magnetic field can be measured. But for EM based SCA, measuring the magnetic field generally gives better results. The magnetic field emitted by a chip is measured using a near field magnetic sensor [START_REF] Peeters | Power and electromagnetic analysis: improved model, consequences and comparisons[END_REF].

Hamming Weight vs. Hamming Distance Models

In Kocher's original paper [START_REF] Kocher | Differential power analysis[END_REF], the leakage model used is based on the Hamming Weight of the bit that the attacker tries to guess. This model, called the Hamming Weight leakage model, considers that a 0 does not lead to excess of power consumption (or EM radiations), whereas a 1 involves a significant amount of power consumption. Thus:

• transitions 0 → 0 and 1 → 0 are considered as not leading to excess power consumption; • transitions 0 → 1 and 1 → 1 are assumed involve an excess of power consumption.

Actually this model does not exactly match with the reality, except in precharged logics (where each register is precharged at 0 before being updated, in this case the Hamming Distance of a such word is equal to its Hamming Weight).

So the Hamming Weight (HW) leakage model could be improved by considering the switching state rather than the output state of the word concerned [START_REF] Brier | Correlation power analysis with a leakage model[END_REF] (the Hamming Distance between the previous and the new state):

• transitions 0 → 0 and 1 → 1 do not lead to excess power consumption; • transitions 0 → 1 and 1 → 0 involve excess power consumption. This leakage model is called the Hamming Distance (HD) leakage model. Figure 3.6 summarizes how these two leakage models rank transitions of one bit, considering its four possible transitions.

Other leakage models are possible, for instance a leakage model that distinguishes rising transitions and falling transitions has been studied in [START_REF] Natale | An integrated validation environment for differential power analysis[END_REF][START_REF] Peeters | Power and electromagnetic analysis: improved model, consequences and comparisons[END_REF], but the results obtained were similar to those obtained with the HD model.

Side-Channel Attacks

As shown previously, both the power consumption and the EM radiations of a circuit are tightly correlated with the data it processes. In this context, different cryptanalytic methods have been proposed to exploit this behavior, and to guess secrets involved in a cryptographic operation. The four main methods that exploit power or EM leakage presented here are respectively SPA/SEMA, DPA/DEMA, Template attacks, and Information theoretic attacks.

Note that, in the rest of the chapter, we do not distinguish between power consumption and EM radiations of cryptographic operations. Our explanations focus on the algorithmic treatment of physical leakages, mainly considering power based attacks. Indeed, up to now most published works are related to the power side channel but the algorithmic treatment is the same as for the EM side channel.

SPA/SEMA

The Simple Power Analysis (SPA), originally described in [START_REF] Kocher | Differential power analysis[END_REF], exploits a single power consumption trace of a cryptographic operation to guess the secret used in the computation (SEMA stands for Simple ElectroMagnetic Analysis, in the case of an EM trace of a cryptographic operation).

Actually, this method works directly when applied on some asymmetric cryptography operations. For instance, in the case of an RSA decryption, the computation of the modular exponentiation:

m = c d (mod n) (3.1)
with m being the plaintext, c the ciphertext, d the private key exponent and n the modulus, can be naively calculated using the right-to-left binary method (d i is the ith bit of d), as described in Algorithm 1.

Depending the value of each bit of the private key exponent d, one computes only a squaring operation or both a multiplication and a squaring operation. Hence the power consumption (or EM radiations) of each elements of the for loop will be different depending on the value of the bit of the private key exponent d, and by simply analyzing one trace, an attacker can easily guess bit values of the private key d. Figure 3.7 shows a typical trace (after demodulation at the FPGA clock frequency [START_REF] Meynard | Enhancement of simple electromagnetic attacks by pre-characterization in frequency domain and demodulation techniques[END_REF]).

A similar attack is possible on ECC, but rather than attacking the square and multiply algorithm, one can attack the double-and-add algorithm involved in the point multiplication. Concerning symmetric key cryptography, SPA attacks focusing on specific implementations of the key schedule have been published. [START_REF] Mangard | A simple power-analysis (SPA) attack on implementations of the AES key expansion[END_REF] reports an SPA against the AES Key Expansion, and [START_REF]A simple power analysis attack against the key schedule of the camellia block cipher[END_REF] describes an SPA against the key schedule of the Camellia block cipher.

DPA/DEMA

The Differential Power Analysis (DPA), originally described in the Kocher's pioneering paper [START_REF] Kocher | Differential power analysis[END_REF], is a powerful technique that allows the attacker to guess secret keys used in a lot of cryptographic primitives. A lot of authors have proposed improvements of this attack along with countermeasures. When the attacker uses the EM radiations of the chip rather than its power consumption, DPA becomes DEMA, for Differential Electromagnetic Analysis.

Original Algorithm

DPA is a known plaintext or known ciphertext attack. The adversary enciphers (resp. deciphers) N PlainText Inputs (PTI) (resp. N CipherText Outputs, CTO) with an unknown key stored in the device, and monitors the power consumption of the device Note that these traces have to be well aligned; this means that the time index of the beginning of ciphering has to be the same for all measurements.

If the measurements are not well aligned (due for instance to countermeasures like random clock frequency or dummy instructions), different preprocessing techniques enables traces to be re-synchronized [START_REF]A method for resynchronizing a random clock on smartcards[END_REF][START_REF]Improving the DPA attack using wavelet transform[END_REF][START_REF]High-resolution side-channel attack using phase-based waveform matching[END_REF][START_REF] Clavier | Differential power analysis in the presence of hardware countermeasures[END_REF][START_REF] Kafi | Deconvolving protected signals[END_REF][START_REF] Réal | Defeating classical hardware countermeasures: a new processing for side channel analysis[END_REF].

Figure 3.8 shows several EM radiation traces corresponding to the DES encryptions of different PTIs with the same key monitored on an FPGA.

The second stage is statistical processing of the N PTIs (resp. N CTOs) with the N traces.

In the rest of the chapter, the DES [2] is used as example, because it is the best most known block cipher and the principles that apply to the DPA also apply to other cryptographic ciphers, such as AES [START_REF]Advanced Encryption Standard[END_REF]. For the sake of convenience, we consider that the adversary is using a known plaintext attack and is trying to guess the roundkey 1 of the DES (the remaining 8 bits could then be discovered through a brute force attack). A similar algorithm allows the round-key 16 in a known ciphertext attack to be discovered.

Since the set of all possible values for the round-key 1 is too large to test all of them, the adversary usually divides the round-key 1 into 8 parts of 6 bits each (here called sub-keys) and attacks each sub-key independently and sequentially. Thus, for each sub-key, there are 64 possible values.

Moreover, attacking each sub-key independently allows all kind of implementations to be targeted. For instance, in software implementations, depending on the size of the data bus (generally 8, 16 or 32), sbox processing can be computed se-quentially, and so sbox 1 is not processed at the same time index as sbox 8. Unlike in some hardware implementations, all the sboxes are processed at the same time, meaning other approaches are possible, as explained below.

The adversary makes hypothesizes on the 6 bits of the attacked sub-key, and for each PTI, he computes the output (4 bits) of the corresponding sbox. This value is called the Intermediate Value (IV). In the state-of-the-art mono-bit DPA [START_REF] Kocher | Differential power analysis[END_REF], the adversary targets out of the four, for instance the Less Significant Bit (LSB).

If the LSB of IV 1 (corresponding to the first plaintext, PTI 1 ) is equal to 0, the associated trace (T 1 ) is ranked in set A. If the LSB of IV 1 is equal to 1, T 1 is ranked in set B. As explained above, the adversary ranks all the traces, in sets A or B, and then computes the difference in the means of sets A and B. The resulting curve is called a differential curve, and corresponds to a sub-key hypothesis. The adversary computes the 64 differential curves corresponding to the 64 possible values for a given sub-key hypothesis.

As explained in [START_REF] Kocher | Differential power analysis[END_REF], the differential curve, denoted , for a sub-key hypothesis K s , is calculated as follows:

K s [j ] = N i=1 D(PTI i , K s )T i [j ] N i=1 D(PTI i , K s ) - N i=1 (1 -D(PTI i , K s ))T i [j ] N i=1 (1 -D(PTI i , K s )) , ( 3.2)
where K s [j ] is the j th sample of the differential curve, N is the number of traces used, PTI i is the ith plaintext, T i [j ] is the j th sample of the trace and D the decision function than ranks traces in set A or B, also called selection function.

If the sub-key hypothesis is wrong, all the computed intermediate values will be wrong with respect to the data really processed on the chip. In this case, the traces will be randomly classified in sets A and B. The mean curves of sets A and B will be similar, and the differential curve will look like a thick horizontal line (mainly composed of noise).

On the other hand, if the sub-key hypothesis is correct, all the computed intermediate values will match the data really processed on the chip, and traces in set A will have the same characteristic: at the time index where the intermediate value is computed, the LSB of IV equal to 0 will not lead to excess power consumption. Conversely, when the LSB of IV is equal to 1, a bit more energy will be consumed at the same time index and spikes corresponding to the clock cycle where IV is computed will be greater on traces in set B than on traces in set A. When the difference in the means of the 2 sets is being computed, a spike will appear at the time index of the differential curve concerned indicating that the sub-key hypothesis is correct.

Figure 3.9 shows the 64 differential curves computed following guesses of subkey 1 of round-key 1 of the DES, using 500 EM traces (each one averaged 20 times). Differential curves corresponding to wrong guesses of the sub-key are in cyan, whereas the curve corresponding to the correct guess of the sub-key is in black, and has the greatest peak. To guess the eight parts of round-key 1, processing is applied sequentially on each sub-key. Moreover, an iterative implementation is used in most hardware implementations of block ciphers. Thus, only one generic round is implemented, and at each clock cycle the intermediate result of the previous round is used as input with the current round-key generated by the key schedule.

Hence, rather than considering an output bit of an sbox, we can consider the bit linked to the previous one in R1 [2]. Thus, as the adversary knows the value of this bit in R0 (because he/she knows the PTI), and assuming that R0 and R1 are stored in the same register which is updated at each round, we can compute the Hamming Distance between the value of the targeted bit in R0 and R1. In this case, it is generally not the power consumption of an inverter that is estimated, but the power consumption of a flip-flop.

Multi-bit DPA Another way to improve DPA attacks consists in considering the four bits linked to the output of the sbox rather than only one. As proposed in [START_REF] Bevan | Ways to enhance differential power analysis[END_REF],

we can compute the differential curve for each bit out of the four from the output of the sbox, and sum the four differential curves.

We can also rank traces using another criterion [START_REF] Messerges | Investigations of power analysis attacks on smartcards[END_REF]: by summing the number of switching bits out of the four considered, and if the sum is smaller than 2, we rank the associated trace in set A, whereas if the sum is greater than 2, we rank the trace in set B.

Distinguishers

In the original DPA algorithm, the distinguisher used to correlate predictions with measurements is called Difference of Mean (DoM). Different works have proposed other methods, based on different mathematical tools. Here, we describe the best known.

Partition Distinguishers

As explained previously, the Difference of Means (DoM) distinguisher consists in ranking, according to a key hypothesis, all the traces in two sets, following a criterion. Then the difference in the means of the two sets of traces is computed. Applying this method to the different key hypotheses reveals the correct key hypothesis.

In the previous paragraph, we described different improvements of Kocher's attack, like the multi-bit DPA proposed by Messerges [START_REF] Messerges | Investigations of power analysis attacks on smartcards[END_REF], or the method suggested by Bevan [START_REF] Bevan | Ways to enhance differential power analysis[END_REF].

Another proposal, presented in [START_REF] Le | A proposition for correlation power analysis enhancement[END_REF], consists in ranking traces in more than two sets. In the case of attacking a DES, if the attacker focuses on the four output bits of an sbox, he/she can rank traces in five sets, one per Hamming Weight. Unfortunately, no efficient method has been proposed to determine these coefficients in the general context. Indeed, due to the process variations, foundries cannot ensure that all the bus lines have exactly the same geometrical dimensions. For example, bit b 0 will not consume exactly the same amount of power (or radiate exactly the same field) as bit b 1 . This is all the more true in the presence of countermeasures. For instance, in a circuit where the nets are balanced (for instance thanks to a dualrail logic [START_REF] Tiri | A logic level design methodology for a secure DPA resistant ASIC or FPGA implementation[END_REF]), the coefficients a i are given in Fig. 3.10 [START_REF] Sauvage | Exploiting dualoutput programmable blocks to balance secure dual-rail logics[END_REF].

Blind techniques enable identification of the coefficients. The optimal coefficients in the context of masking are derived in [START_REF] Prouff | Statistical analysis of second order differential power analysis[END_REF]. In a general context, the stochastic approach, discussed in Sect. 3.4.3.2, always applies.

Comparison Distinguishers

In a differential SCA, the attacker wants to estimate the relation between predictions, depending on a key hypothesis, and on measurements. In [START_REF] Brier | Correlation power analysis with a leakage model[END_REF], the authors suggest using a well-known statistical tool, Pearson's correlation coefficient. Let X and Y be two random variables, cov(X, Y ) be the covariance between X and Y , and σ X and σ Y be the standard deviations of X and Y , Pearson's correlation coefficient is computed as follows: This mathematical tool measures the dependence between two quantities, and for SCA, its application is straightforward. Let N be the number of traces, K s the key hypothesis, PTI i the ith PlainText Input, H (K s , PTI i ) the number of switching bits according to the key hypothesis K s and the plaintext PTI i , T i (j ) the j th sample of the ith trace, Eq. (3.3) becomes:

ρ(X, Y ) = cov(X, Y ) σ X × σ Y . ( 3.3)
ρ K s [j ] = N. N i=1 H (K s , PTI i ).T i (j )- N. N i=1 H (K s , PTI i ) 2 -( N i=1 H (K s , PTI i )) 2 × N i=1 H (K s , PTI i ). N i=1 T i (j ) N. N i=1 T i (j ) 2 -( N i=1 T i (j )) 2 . (3.4)
The variant of the DPA using the Pearson's correlation has been called Correlation Power Analysis (CPA), and was first introduced in [START_REF] Brier | Correlation power analysis with a leakage model[END_REF].

However, Pearson's correlation measures only linear dependencies between two quantities. Other tools, like Spearman's rank correlation, or Kendall's tau rank correlation, can measure linear dependencies between two random variables, but practical experiments have shown that these tools give very similar results [9].

Another author [START_REF] Gierlichs | Mutual information analysis[END_REF] introduced a different statistical tool, mutual information. Mutual information allows both linear and non-linear dependencies between two random variables to be measured. This is called Mutual Information Analysis (MIA). This technique is discussed in Sect. 3.4.4.1.

Key Search

Generally, when the attacker has computed the 64 differential curves corresponding to guesses of the sub-key S k using N traces, he/she searches for the maximum (in absolute) of each differential curve. Among the 64 maxima, the biggest one reveals the guessed during the attack using N traces. A classical optimization for the key search step consists in reducing the proportion of the differential curve inspected to detect the maximum. The key search may concentrate on a part of the curve corresponding to the attacked round.

Template Attack

The Template Attack (TA) is considered to be the most powerful type of SCA. The scenario of such an attack is slightly different than that required for other kinds of SCA. First, to build templates for the different possible values of the secret, attackers have to be in possession of a circuit identical to the circuit under attack. Then, they measure one trace on the attacked device, and using appropriate methods, compare this trace with the different templates they built to guess the value of the secret.

Original Template Attack

This attack, firstly described in [START_REF] Chari | Template attacks[END_REF], happens in two stages, the template building stage, and the template matching stage.

Template Building Phase In the first stage, attackers use the circuit under their full control, which is identical to the circuit under attack, to build templates for each pair of data and key.

Power traces can be characterized by a multivariate normal distribution, which is fully defined by a mean vector M and a covariance matrix C. In the rest of the paragraph, the pair (M, C) is referred to as a template. The attacker builds a template (M, C) d i ,k j for each pair of data and key (d i , k j ) following these steps:

• for each pair of data and key (d i , k j ), acquire p traces T d i ,k j 1 , . . . , T d i ,k j p ; • compute a mean vector M d i ,k j for each pair (d i , k j ) from the p traces as follows:

M d i ,k j = 1 p p l=1 T d i ,k j l ; (3.5)
• (optional) compute pairwise differences between the mean vectors M d i ,k j in order to identify and select only points P 1 , . . . , P n at which large differences show up. This optional step significantly reduces the processing overhead with only a small loss of accuracy; • for each pair d i , k j and for each acquired trace T d i ,k j l , l = 1, . . . , p, compute the noise vector N d i ,k j l , l = 1, . . . , p as follows:

N d i ,k j l = (T d i ,k j l (P 1 ) -M d i ,k j (P 1 ) . . . T d i ,k j l (P n ) -M d i ,k j (P n )); (3.6)
• for each pair d i , k j , compute the noise covariance matrix C d i ,k j between all pairs of components of the noise vectors N d i ,k j l s, l = 1, . . . , p, for all the p traces, as follows:

C d i ,k j [u, v] = cov(N l (P u ), N l (P v )).
(3.7)

Thus, for each data and key pair (d i , k j ), the attacker has built a template, that characterizes the deterministic part of the power trace, but also its noise part, modeling it using a multivariate Gaussian model.

Template Matching Phase

In the second stage, the attackers use a power trace from the device under attack to determine the key. To this end, they evaluate the probability density function of the multivariate normal distribution with (M, C) d i ,k j and the trace of the device under attack.

In other words, given a power trace t of the chip under attack, and a template (M, C) d i ,k j , the attacker computes the probability:

p(t; (M, C) d i ,k j ) = 1 √ (2π) n • det(C) e -1 2 (t-M) T •C -1 •(t-M) . ( 3.8)
The attacker does this for every template. As a result, he/she obtains the probabilities p(t; (M, C) d 1 ,k 1 ), . . . , p(t; (M, C) d D ,k K ). The probabilities measure how well the templates fit to a given trace. Intuitively, the highest probability should indicate the correct template. Because each template is associated with a key, the attacker also gets a clue to the correct key. This intuition is also supported by the maximumlikelihood decision rule.

The points of interest can be either selected heuristically [START_REF] Chari | Template attacks[END_REF], by ad hoc techniques [START_REF] Gierlichs | Templates vs. stochastic methods[END_REF], or by dimensionality reduction tools, such as principal components analysis [START_REF] Archambeau | Template attacks in principal subspaces[END_REF].

Stochastic Method

An improvement of TA is suggested in [START_REF] Schindler | Advanced stochastic methods in side channel analysis on block ciphers in the presence of masking[END_REF][START_REF] Schindler | A stochastic model for differential side channel cryptanalysis[END_REF]. It consists in using a stochastic method to approximate the real leakage function within a suitable vector subspace. The attack requires some engineering skills to introduce a relevant parametric model. Then, the on-line attack basically consists in simultaneously:

• estimating the best parameters, using a linear regression, and • deciding which model is the closest to the side channel measurements, using the minimal Euclidean distance as a distinguisher.

An empirical comparison of template and stochastic attacks is provided in [START_REF] Standaert | How to compare profiled side-channel attacks?[END_REF].

Information Theoretic Attacks

MIA-Mutual Information Analysis

At CHES 2008, a generic side channel distinguisher, MIA, was proposed in [START_REF] Batina | Mutual information analysis: a comprehesive study[END_REF]. It is an attractive alternative to the above-mentioned attacks since some assumptions about the adversary can be disregarded. In particular, it does not require a linear dependency between the leakage and the predicted data, as is the case for DPA and CPA, and so it is not only able to exploit any kind of dependency but does not need to profile the leakage as it is the case for template attacks [START_REF] Chari | Template attacks[END_REF].

The rationale of MIA is to compare distributions of observations with random distributions rather observations with a model. This is why the MIA measures the Kullback-Leibler divergence between observations knowing the correct sub-key and observations not knowing anything about the internal values [START_REF] Veyrat-Charvillon | Mutual information analysis: how, when and why?[END_REF].

EPA-Entropy Power Analysis

Mutual information analysis has been tested in noisy real world designs. It indeed appears to be a powerful approach to break unprotected implementations. However, the MIA fails when applied on a DES cryptoprocessor with masked substitution boxes (sboxes) in ROM [START_REF] Standaert | FPGA implementations of the DES and triple-DES masked against power analysis attacks[END_REF]. However, this masking implementation remains sensitive to Higher Order Differential Power Analysis (HO-DPA). For instance, an attack based on variance analysis clearly reveals the vulnerabilities of a first-order masking countermeasure. A novel approach to information theoretic HO attacks, called Entropy-based Power Analysis (EPA), has therefore been proposed. This new attack gives greater importance to highly informative partitions and distinguishes the key hypotheses better [START_REF] Maghrebi | Entropy-based power attack[END_REF].

VPA-Variance Power Analysis

Information theoretic attacks are extremely powerful, as they exploit any deviation from a random probability density function (PDF). However, estimating PDFs is a hard task [START_REF] Prouff | Theoretical and practical aspects of mutual information based side channel analysis[END_REF][START_REF] Veyrat-Charvillon | Mutual information analysis: how, when and why?[END_REF], since for the PDFs to be accurate, it ideally requires a lot of measurements. To simplify the attack, articles [START_REF] Li | Power variance analysis breaks a masked ASIC implementation of AES[END_REF][START_REF] Maghrebi | Evaluation of countermeasures implementation based on Boolean masking to thwart first and second order side-channel attacks[END_REF][START_REF] Standaert | Partition vs. comparison side-channel distinguishers: an empirical evaluation of statistical tests for univariate side-channel attacks against two unprotected CMOS devices[END_REF] suggest limiting the analysis of the PDFs to their second cumulant. The estimation is much faster and the distinguisher thus gains strength. This is referred to as variance power analysis (VPA). This metric to compute distances between PDFs is extremely useful against implementations protected by first-order masking.

Conclusions

More than ten years after the first publication of Kocher's attack, a lot of improvements have been proposed for power and EM based Side Channel Attacks. Firstly, it has been shown that both power consumption and EM radiations of a circuit can leak sensitive information. Secondly, different methods have been proposed to extract a secret key used in a cryptographic operation from physical leakages emitted by the device. In this chapter we described SPA/SEMA, DPA/DEMA and Template Attacks. We also described different improvements of differential SCA that give better results than the original attack, especially when the attacker has some knowledge of the implementation. To thwart these attacks, a lot of methods have been proposed. Due to the specificities of FPGAs, some of these countermeasures are not applicable, while others are relatively easy to integrate. The next chapter describes the main methods used to protect cryptographic algorithms implemented in FPGAs against Side Channel Attacks.

Fig. 3 . 1

 31 Fig. 3.1 Black-box and gray-box models from a cryptanalytic point of view

Fig. 3 . 2 Fig. 3 . 3

 3233 Fig. 3.2 Reverse-engineering of the MyFare CRYPTO1 algorithm. On the left, pictures of the circuit. On the right, netlist reconstruction

Fig. 3 . 4

 34 Fig. 3.4 Commercial laser attack setup from Riscure (on the left)-zoomed view of memory with the laser pulse in pink (on the right) (color online)

Fig. 3 . 5

 35 Fig. 3.5 Power and Electromagnetic Measurement Platform

Fig. 3 . 6

 36 Fig. 3.6 Ranking of a transition of one bit following the Hamming Weight and Hamming Distance leakage models

Algorithm 1

 1 Right-to-left binary method for modular exponentiation Inputs: c, d = (d l-1 d l-2 . . . d 0 ) 2 , n Output: m m = 1 for i = 0 to l -1 do if d i = 1 then m = m * c (mod n) end if c = (c * c) (mod n) end for Result: m

Fig. 3 . 7

 37 Fig. 3.7 Example of an SEMA on RSA

Fig. 3 . 8

 38 Fig. 3.8 Electromagnetic radiation traces of three different DES encryptions monitored on an FPGA

Fig. 3 . 9

 39 Fig.3.9 Example of a successful DEMA: 64 differential curves computed following hypotheses of sub-key 1 of round-key 1 of the DES, using 500 EM traces (each one averaged 20 times). Differential curves corresponding to wrong hypotheses of the sub-key are in cyan, the curve corresponding to the correct hypothesis of the sub-key is in black (color online)

  Considering four bits b 0 , b 1 , b 2 and b 3 , the Hamming Weight of the word b 3 b 2 b 1 b 0 is between 0 and 4, and leads to five possible values. Then the attacker has to choose coefficients a i , i = 0, . . . , 3, assigned to each set. This method has been called Partitioning Power Analysis (PPA).

Fig. 3 . 10

 310 Fig. 3.10 Optimal weights found for the sbox output of DES in dual rail with precharge logic


Marc Alías-Rodríguez
email: marc.alias-rodriguez@univ-amu.fr

Miquel Huix-Rotllant

Coen De Graaf

Coen De

HAL is

First row transition metal complexes with d 4 to d 7 electronic configuration exhibit spin-crossover (SCO), which can be induced by external stimulus such as temperature, pressure, light. The low-spin to high-spin transition has been largely studied, but very little is known about the reverse process. Here, we present a theoretical study of thermal and light-induced high-to-low spin crossover in prototypical Fe(II) complexes. The lifetime of the high-spin state in the thermal process is determined using Fermi's golden rule. With this methodology, we have accurately computed the transfer rate of the HS state thermal relaxation at several time scales (from sub-nanosecond to a few seconds) in two different iron complexes. The use of quasi-degenerate perturbation theory (QDPT2) in the analysis of the LS-HS spin-orbit coupling have allowed to identify 3 T 1 as the main intermediate state coupling LS and HS states. The light-induced process has been studied using wavepacket quantum dynamics along the main vibrational coordinates (one symmetric and two asymmetric Fe-N stretchings). The study suggests that after the initial excitation from 5 T 2g to 5 E g state, the population is transferred back to a vibrationally hot 5 T 2g from which a small amount of population is transferred to the 1 A 1g via the intermediate 3 T 1g . Most of the population remains trapped in the HS state at the time scale of the simulation.

Spin-crossover (SCO) is a molecular phenomena wherein the spin state of the system changes because of an external stimulus such as pressure, temperature, light irradiation or the influence of a magnetic field. This process commonly takes place from the diamagnetic thermally stable low-spin state to the paramagnetic meta-stable high-spin state. Despite SCO can be observed in some molecules, 1,2 the phenomena is typically appreciated in d 4 to d 7 first row transition metal complexes. [3][4][5] Among the different transition metal complexes exhibiting spin-crossover, Fe(II) complexes have probably been the most widely studied ones. Fe(II) has a d 6 valence configuration and is usually presenting a (quasi-)octahedral coordination. [6][7][8][9] In ferrous complexes, the magnetization varies between the diamagnetic singlet electronic state (S=0) and the paramagnetic quintet electronic state (S=2).

The magnetic bistability gives to these complexes a wide range of possible applications such as thermochromic paints, molecular switches or nanophotonic devices to name a few. 10 For technological uses, low and high spin states must show a reasonable lifetime in the working temperatures, ideally room temperature. The lifetime of meta-stable high-spin states strongly depends on their adiabatic energies which makes it vary from a few picoseconds to several seconds or even days.

Hauser et al. described the thermal relaxation of the high-spin state as a tunneling process, which does not need to overcome the classical barrier. 11 They approximated the quantum mechanical expression for the rate constant in a tunneling process, assuming that the LS-HS spin-orbit coupling is a fixed value of 150 cm -1 , the vibrational frequency along the Fe-N symmetric stretching is for both the HS and LS states 250 cm -1 and these states are only displaced along this mode at a constant distance ∆r HL ≈ 0.2Å. The so-obtained analytical expression to determine the rate of the high-spin to low-spin transition depends only on the adiabatic energy difference between these states, and is known as the energy gap rule. This approach was validated for isotropic complexes, for which the reaction coordinate is mainly a totally symmetric Fe-N breathing mode such as [Fe(bpy) 3 ] 2+ .

The light-driven mechanism for the formation of the HS state in SCO complexes, (commonly known as light-induced excited spin state trapping, LIESST) has been largely studied in Fe(II) complexes using [Fe(bpy) 3 ] 2+ as prototype complex. Auböck et al. determined a direct transfer from the metal to ligand charge transfer ( 1,3 MLCT) band to the 5 T state in the sub-50-fs time scale through ultrafast pump-probe spectroscopy. 12 Zhang et al. based on X-ray emission spectroscopy concluded that the population of the 5 T state was in two sequential steps, 1,3 MLCT-→ 3 T in 150 fs and 3 T-→ 5 T in 70 fs. 13 From a theoretical perspective, one of us performed several studies based on Fermi's golden rule. The deactivation path via the 3 MC states was confirmed in the symmetric structures, 14 but the presence of geometrical distortions reopened the direct transfer from the 1,3 MLCT band, becoming both paths competitive. 15 Recently, Pápai published a study about the LIESST in [Fe(NCH) 6 ] 2+ , where only MC states play a role, based on quantum dynamic simulations. 16 He constructed a vibronic Hamiltonian along the three main modes (Fe-N stretchings) and determined an ultrafast first ISC 1 T 1g -→ 3 T 2g in 100 fs and a second slower ISC to the 5 T.

Despite the important amount of LIESST studies, of which only a few are mentioned in the above paragraph, the number of studies on the reverse process is scarce. Marino et al. studied the reverse-LIESST after the excitation to the 5 E g in [Zn 1-x Fe x (ptz) 6 ](BF 4 ) 2 , where ptz=1-propyltetrazole, using ultrafast pump-probe spectroscopy. 17 Following the time-resolved excited-state absorption (ESA), they established that two different transitions are involved in the process, whose rates were determined by a double exponential fitting. The first one was attributed to the 5 E g -→ 3 T 1g in 1.7 ps ISC and the second to the 3 T 1g -→ 1 A 1g in 39 ps. The quantum yield was determined to be 0.1 for the reverse-LIESST and 0.8 for the LIESST in [Fe(ptz) 6 ](BF 4 ) 2 at 10K, with a branching ratio 3 T 1 -→ 5 T 2 / 3 T 1 -→ 1 A 1 of 4:1. 18 In the current study, we simulated the thermal and light-induced conversion from the high-spin to low-spin in [Fe(bpy) 3 ] 2+ , where bpy=2,2'-bipyridine and [Fe(mtz) 6 ] 2+ , mtz = 1-methyl-1H-tetrazole. Two different processes have been studied, the HS-→LS thermal relaxation and the photophysical process upon the absorption of a photon from the HS ( 5 T 2g -→E g ), i.e. the reverse-LIESST. The study of the thermal process is based on Fermi's golden rule and the photo process has been performed using wave-packet quantum dynamics with a vibronic model Hamiltonian.

Thermal relaxation. The HS-→LS direct transfer has been studied in two Fe(II) complexes with significantly different HS-LS relative energies. [Fe(bpy) 3 ] 2+ is a commonly studied complex with short-living HS state and [Fe(mtz) 6 ] 2+ is either a paramagnetic or diamagnetic complex depending on the temperature. The intersystem crossing rate has been calculated using the Fermi's golden rule within the Condon approximation, which allows to express the rate as a direct product between the vibrational overlap and the (effective) spin-orbit coupling.

k a→b ISC = 2π Ω a ∑ α ∑ γ ψ α a | ĤSO |ψ γ b 2 ∑ j=0,k | {ν a j }| |{ν bk } | 2 δ (E a j -E bk )
(1) Herein ψ α a and ψ γ b are the spin-free electronic wave functions where a and b label electronic states of different multiplicity and α and γ are, respectively, the MS-components of the initial and final states. j and k are the vibronic levels, E a j and E bk are, respectively, the energies of the electronic state a in vibrational level j and electronic state b in vibrational level k, E bk =E b +(k+1)hω k /2. Ω a is the degeneracy of the initial state and {ν a j }| |{ν bk } is the two-center Franck-Condon integral. Note that this is the vibrationally cold approximation and the vibrational level of the initial state is fixed, j=0.

The effect of temperature is taken into account by the introduction of the Boltzmann distribution in the initial vibronic states.

k a→b ISC = 2π ZΩ a ∑ α ∑ γ ψ α a | ĤSO |ψ γ b 2 × ∑ j,k e -β E j | {ν a j }| |{ν bk } | 2 δ (E a j -E bk ) (2)
where Z = ∑ j e -β E j is the canonical partition function for vibrational term in the initial state, E j is the vibrational energy of the level j and β is 1/(k B T), in which k B is the Boltzmann constant. We have mainly used the time-dependent approach, obtained after Fourier transform the δ function in Eq. 1, implemented in the VIBES program. [START_REF] Etinski | [END_REF] Then, the rate is determined analytically and the number of grid points and the spacing between them to perform the discrete Fourier transform were selected upon convergence of the correlation and cumulant functions. In one case, we have used the time-independent approach for comparison. η, the threshold that determines the number of vibronic states included, was chosen as 1000 cm -1 after some testing of its influence. The basic ingredients for these calculations of the rates are: i) the vibrational normal modes and frequencies of the states to account the vibrational overlap; ii) the adiabatic electronic energy difference between the states, and iii) the modulus of (effective) spin-orbit among the different spin components of each state.

The vibrational modes and frequencies are computed after a geometry optimization for the lowest singlet (LS) and quintet (HS) states using the TurboMole package. [START_REF]TURBOMOLE V6.4.2 2012, a development of University of Karlsruhe and[END_REF] The Perdew-Becke-Ernzerhof (PBE0) hybrid functional is used in combination with the triple-ζ valence polarized Gaussian-type basis set, def2-TZVP. Two-electron integrals are approximated by the resolution of identity (RI) using the auxiliary basis def2-universal-JFIT. The hessian matrix is calculated through analytical second derivatives and no symmetry restrictions are imposed during the optimization procedure.

The B3LYP adiabatic energies are unfortunately not accurate enough in Fe(II) complexes to be used for further analysis. Therefore, we determine the adiabatic energies with CASPT2 and shift the DFT potentials accordingly. Unfortunately, the DFT minimum considerably differs from the CASPT2 one and we carry out a relaxed scan along the main reaction coordinate, the Fe-N symmmetric stretching. The geometries at each point are optimised through a constrained optimization at B3LYP/def2-SVP level. C 2 symmetry is imposed in the case of [Fe(bpy) 3 ] 2+ and no symmetry restrictions are applied in [Fe(mtz) 6 ] 2+ . The Fe-N distances for the LS and HS states in the relaxed scan minima are in good agreement with the (XMS-)CASPT2 fully optimized geometries by Finney et al. [START_REF] Finney | [END_REF] Further information about the relaxed scan is found in SI.

CASPT2 calculations are performed using (Open)Molcas. 22 The basis sets used are ano-rcc type with the following contractions: (7s6p5d4f3g2h) for the Fe, (4s3p1d) for N bonded with Fe, (3s2p) for C and the other N and (2s) for H atoms. Cholesky decomposition for the two-electron integrals is used to speed up the calculations. The reference wave function is a CASSCF stateaveraged (SA) wave function, including the four lowest states for the singlets, the six lowest for the triplets and the five lowest for the quintets. The active space is formed by 10 electrons distributed over 12 orbitals, three more orbitals are added to those calculations that include MLCT states. The molecular orbitals for these systems are well-established 23,24 and include two σ orbitals mainly corresponding to a small bonding interaction between ligands and Fe, the three 3d (t 2g ) orbitals of Fe, the two 3d (e g ) orbitals of Fe with a small anti-bonding contribution respect to the ligands, an additional set of five 4d orbitals on Fe to introduce radial electron correlation 25 . The influence of the MLCT states has also been studied expanding the active space by including three π* orbitals corresponding to the aromatic system of the ligands. The perturbational treatment is carried out in all the electrons except the deep core ones Fe (1s2s2p), C and N (1s) and an imaginary shift of 0.15 E h is used.

The spin-orbit coupling is calculated using Molcas 22 in the singlet-quintet minimum energy crossing point geometry and in several conformations of the metadynamics for the LS and HS in [Fe(bpy) 3 ] 2+ Ṫhe settings for CASPT2 calculations are the same as for the adiabatic energy differences. The direct coupling is determined using the effective one-electron Douglas-Kroll-Hess (DKH)

Hamiltonian. The off-diagonal spin-orbit matrix elements are calculated by state interaction 26 using the CASSCF wave functions and the diagonal matrix elements are replaced by the CASPT2 energies. The direct coupling between states with ∆S=2, such as singlet and quintet states, is strictly zero using the standard effective one-electron spin-orbit operators such as the DKH Hamiltonian. However, these states may have a non-vanishing effective coupling through intermediate states, which are states with non-zero direct coupling simultaneously with both states. In this case, triplet states can couple singlet and quintet states simultaneously. The effective Hamiltonian is constructed using the orthonormalised projections of SO wave functions Ψ⊥ k and their respective eigenvalues E k using Bloch's formula.

H SO,e f f i j = ψ i | ĤSO,e f f |ψ j = ψ i | ∑ k∈S 0 | ψ⊥ k E k ψ⊥ k | |ψ j (3)
Herein, Ψ⊥ k is the des Cloizeaux orthonormalised spin-orbit wave function expressed in the basis of the spin-free states {ψ} spanning a smaller subspace S 0 and E k are the spin-orbit energies.

The effective Hamiltonian theory is an accurate an elegant way to determine effective couplings but it is difficult to extract information about the intermediate states that allow the coupling. For this reason, the weight of the intermediate triplet states in the singlet-quintet coupling was studied using second order Quasi-Degenerate Perturbation Theory (QDPT2) in different conformations. QDPT2 allows to decompose the spin-orbit coupling as a sum over the different intermediate states.

H SO,QDPT 2 i j = ψ i | ĤSO |ψ j + ∑ µ =i, j ψ i | ĤSO |ψ µ ψ µ | ĤSO |ψ j E µ -E j (4)
Here, {ψ} are the spin-free states and {E} the spin-free energies.

ψ i and ψ j are, respectively, the singlet and quintet states and ψ µ the intermediate states.

The results in Table 1 show that 3 T 1g , the lowest set of 3 MC states, have the largest influence in the coupling of the singlet and quintet states. This can be easily rationalised because it is the intermediate state closest in energy to the singlet and quintet and has a strong direct coupling with both states. In addition, the decomposition also showed an almost negligible influence of the 3 MLCT states in the coupling of the 1 A 1g and 5 T 2g states.

Table 1 also show discrepancies between the effective Hamiltonian theory and the quasi-degenerate perturbation theory. These differences are mainly associated to two intrinsic characteristics of the perturbative method used. The first is that QDPT2 loses accuracy when the initial and final states have larger energy differences. Eq. 4 shows that the larger the difference between E i and E j , the larger the difference becomes between H i j and H ji . The relative similarity of the couplings with both methods in the ISC geometry where singlet and quintet are quasi degenerate in comparison with the disagreement in most of the conformations where the energy difference is larger helps to corroborate this observation. The other shortcoming is that in QDPT2 the interactions are expanded up to second order and only the interactions through one intermediate state are accounted for, while all higher order terms are taken into account in the effective Hamiltonian theory.

The influence of the thermal disorder in the spin-orbit coupling, and therefore, in the intersystem-crossing rates was also studied. Previously, it was demonstrated by one of us that the geometrical distortions have a strong effect on the size of the coupling, especially when MLCT states are involved. These changes in the couplings changed the rates of some of the intersystem-crossings by a few orders of magnitude in Fe(II) polypyridinal complexes. 15 However, in the case of the 1 MC-5 MC coupling, the effect of geometrical distortions is almost negligible, giving in all the conformations the same order of magnitude for the couplings and, therefore, the rates (or lifetimes). The labels a and b refers to conformations from the singlet and quintet simulations, respectively.

The results show that the spin-orbit coupling strength is almost identical for both complexes. However, the HS-LS vibrational overlap is completely different in the two complexes, what is the responsible for the different HS lifetimes in these complexes.

Temperature is also important an important factor; the rates at 300 K are one order of magnitude faster than those at 0K using the vibrationally cold approximation. The HS state lifetime in [Fe(bpy) 3 ] 2+ computed with the methodology described here is in good agreement with the 700 ps determined by Consani et al. for the same complex in aqueous solution 27 , while the very slow ISC rate for [Fe(mtz) 6 ] 2+ is indicative of a stable HS spin, as observed in experiment. Fermi's golden rule methodology is presented as an appropriate method to determine the order of magnitude of the HS state lifetime in Fe(II) complexes.

Photophysical relaxation.

The reverse-LIESST is studied from the 5 E g state, after an excitation from the high spin state ( 5 T 2g -→ 5 E g ). The process is studied in [Fe(mtz) 6 ] 2+ because of the large HS state lifetime in this complex and the presence of experimental studies in an analogue complex, [Fe(ptz) 6 ] 2+ . 17 In the high-spin minimum geometry, the 5 MLCT states are clearly higher in energy than the 5 MC states, therefore the former are highly unlikely to play a role in the deactivation and only the latter are included in the study. The reverse-LIESST was in a first attempt studied using Fermi's golden rule, as a follow-up of the satisfactory results reported in the literature for the singlet-quintet conversion. The only energetically favourable intersystem crossing transition from the 5 E g state to the triplets is the 5 E g -→ 3 T 1g . However, the rate for this transition using the Fermi's golden rule is about a few miliseconds, several orders of magnitudes slower than the experimental reported value of 1.7 ps. For this reason, the process is studied using wave-packet quantum dynamic simulations with a vibronic model Hamiltonian.

The simulations are carried out with the Multi-Configuration Time-Dependent Hartree (MCTDH) method (Heidelberg package, version 8.5.13). 28 The basis of the method is to solve the timedependent Schrödinger equation using the following ansatz for the nuclear wave function:

Ψ(q 1 , q 2 , . . . , q p ,t) = n 1 ∑ j 1 =1 ... n p ∑ j p =1 A j 1 ... j p (t) p ∏ k ϕ (k) j k (q k ,t) (5)
Herein, A j are the expansion coefficients and ϕ j are the singleparticle functions (SPFs), which are linear combinations of primitive functions. In this case, the primitive functions take the shape of the n j 's lowest solutions for an harmonic oscillator. The multiset formalism is employed, where a different wave function is written for each state

Ψ m (q 1 , q 2 , . . . , q p ,t) = σ ∑ α=1 Ψ (α) |α ( 6
)
where σ is the number of electronic states. We have included 5 SPFs for each state along each mode and 151 primitive functions for each mode.

The vibronic-coupling model Hamiltonian takes the following shape:

Ĥ({q}) = 2 ∑ S=0 N S ∑ n=1 |S| ∑ α=-|S| T N ({q}) +V S n ({q}) |n S,α n S,α | + N S ∑ n =n W NA n,n ,S ({q})|n S,α n S,α | + 2 ∑ S=0 N S ∑ n=1 |S| ∑ α=-|S| 2 ∑ S =1 N S ∑ n =1 |S | ∑ α =-|S | W SO nSα,n S α |n S,α n S ,α | +W SO, * n S α ,nSα |n S ,α n S,α | (7)
Herein, S is the spin quantum number, N S is the number of states with S spin quantum number and α is the M S component.

T N ({q}) is the nuclear kinetic energy, V S n ({q}) is the diabatic potential, W NA n,n ,S ({q}) is the non-adiabatic coupling and W SO nSα,n S α and its hermitian conjugate W SO, * nSα,n S α are the spin-orbit coupling terms. Note that in our model spin-orbit is coordinate independent while the other terms are coordinate dependent.

The nuclear kinetic energy is defined (in atomic units) as

T N = - N vib ∑ i ω i 2 ∂ 2 ∂ q 2 i ( 8
)
where ω i is the frequency of the HS state in mode i and N vib is the number of vibrations included in the model Hamiltonian.. The diabatic potential is defined using a quartic polynomial expression

V S n = E S n + N vib ∑ i 4 ∑ j=1 1 j! k n,S j,i q j i (9)
Here E S n is the diabatic energy for the state n with spin S at the Franck-Condon region and k n,S j,i is the j th expansion term along vibrational mode i for this state.

The non-adiabatic coupling term is described as

W NA n,n ,S = N vib ∑ i λ n,n ,S i q i ( 10
)
where λ n,n ,S i is a linear-expansion coefficient for the nonadiabatic coupling between states n and n' with spin S along mode i.

The spin-orbit coupling term is defined as follows

W SO nSα,n S α = n Sα | ĤSO |n S α ( 11
)
where ĤSO is the DKH Hamiltonian, n Sα | is the state n with spin S and M S α and |n Sα is the state n' with spin S' and M S α'.

The model is constructed by the expansion along the most representative modes of the HS state optimized at PBE0/def2-TZVP level and includes the electronic states 1 A 1g , 1 T 1g , 3 T 1g , 3 T 2g , 5 T 2g and 5 E g , which generates a model of 47 states taking into account the spatial and spin symmetries. The terms for the Hamiltonian are obtained through a fitting procedure of the XMS-CASPT2 energies. These calculations are carried out using an active space of 10 electrons distributed over 12 molecular orbitals. The per-turbative treatment of the electron correlation is applied to all electrons except the deep core ones Fe (1s2s2p) and C, N (1s) and the 30 highest orbitals. An imaginary level shift of 0.25 E h is applied to avoid intruder state problems. From the adiabatic energies along the main modes, we obtain the diabatic energies as a quartic polynomia and the non-adiabatic couplings between states of the same M S value within the linear-expansion approximation. The off-diagonal terms between states of different M S are included with the spin-orbit coupling at the DFT optimised geometry. The nuclear kinetic energy operators in mass-frequencyweighted coordinates are added using the frequencies of the HS state at DFT level. The simulations are initialised after a relaxation of the wave-packet in the 5 T 2g state, which starts the dynamics at the minimum energy geometry of the HS state and the wave-packet is distributed over the different components of the 5 E g state according to the oscillator strengths. Further (numerical) information about the couplings is given in the SI.

The selected normal modes are chosen based on those modes with the largest reorganization energy connecting the minima of the different states involved in the process. The difference between the optimized geometries for each state at DFT level is projected into the HS state normal modes, following the expression

q A k = ω 1/2 k ∑ l L HS kl m 1/2 l (R A,l -R HS,l ) (12)
Herein, q A is the normal mass-frequency-weighted Cartesian coordinate expanding the geometrical differences between state A and the HS state projected in the basis of HS state reference normal modes. m l the mass of atom l, R A and R HS are the Cartesian geometries for states A and the HS state, respectively and L HS is the matrix that contains the eigenvectors of the HS state normal modes.

Table 3 Coordinates q A of the main electronic states involved in reversed LIESST projected on the HS state normal modes. In parenthesis the value of q A , those modes with values larger than 1 are reported here. The modes in the reaction coordinate are in bold. According to this analysis, there are three main modes representing the reverse-LIESST. These modes have Fe-N stretching character, q 42 is the symmetric motion and q 31 and q 33 correspond to Fe-N asymmetric stretching vibrations, which are specially important for Jahn-Teller active states. This model with three Fe-N modes coincides with the one proposed by Pápai in his study of the LIESST for [Fe(NHC) 6 ] 2+ . In our model, we include the two asymmetric modes q 31 , q 33 and the reaction coordinate mode (q rc ) that connects HS and LS geometries which is a com-Fig. 2 Schematic motions and diabatic potentials along the modes q 31 , q 33 and q rc . 1 A 1g in black, 1 T 1g in cyan, 3 T 1g in green, 3 T 2g in gold, 5 T 2g in red and 5 E g in violet. The different components are represented with solid, dashed and dotted lines.

Included modes Excluded modes

bined mode described as: q rc = 0.95q 42 + 0.21q 30 -0.18q 10 -0.17q 39 (13) This reaction coordinate mode describes the Fe-N symmetric breathing (q rc ) but slightly assisted for N-Fe-N bending modes (q 30 , q 10 , q 39 ) which aids to adapt the N-Fe-N angles to each Fe-N distance.

The reaction coordinate mode does not break the pseudo-O h symmetry. In this mode, the triple degeneracy of 1 T 1g , 3 T 1g , 3 T 2g , 5 T 2g and the double degeneracy of 5 E g is maintained. The Fe-N asymmetric stretching (q 31 ) is described as a compression of the four Fe-N distance for the N in the equatorial plane with an elongation of the two Fe-N distances for the axial N. This breaks the pseudo-O h symmetry to a pseudo-D 4h , which breaks the triple degenerate states into doubly degenerate states and remove the degeneracy for the 5 E g states. In the other asymmetric stretching (q 33 ), the axial N remain in place while the two equatorial N come closer and recede in pairs. After this motion, the pseudo-O h symmetry is broken to a pseudo-D 2h that removes the degeneracy for all the states.

The efficiency of an intersystem-crossing depends on the energy difference between a pair of states and the size of their coupling. In this complex, the energetically accessible intersystem-crossing from the 5 E g is the 5 E g -→ 3 T 1g . Near the absorption geometry, that is the 5 T 2g minimum, 3 T 1g is slightly lower in energy than 5 E g , therefore, a crossing between these states is found after the absorption. However, the gradients to the 5 E g minimum are high and the system will evolve to this region where the energy gap between the quintet and the triplet is large and the efficiency of the ISC small. The transition to the 3 T 2g is discarded because these states are above the 5 E g states in the absorption region.

Despite the energy gap does not seem to favour the ISC transfer, this could still be efficient when there exists a large coupling between the states. This is not the case for the 5 E g -3 T 1g coupling, which is about 10 cm -1 in average at the 5 T 2g DFT minimum. The small coupling can be easily rationalized observing the two electrons of difference in the main configurations of the 5 E g (t 3 2g e 3 g ) and 3 T 1g (t 5 2g e 1 g ), which are strictly uncoupled in the usual effective one-electron spin-orbit operator.

The minimum of the 5 E g is energetically far away from the other states, 5 E g is lower in energy than 3 T 1g and above the 5 T 2g , well separated from both. Despite the large energy gap with these other states, the nuclear kinetic energy may bring the wavepacket to regions with nearly degenerate 5 T 2g -5 E g , in a sloped conical intersection. Unfortunately, this intersection could not be determined because of the lack of an accurate computational method to determine precisely at the same time geometries (for the minima of states and minimum energy crossings) and energies at an affordable computational cost. The plausibility of the crossing in the quintets is also not favoured by the energy difference between the states in the regions around the 5 E g minimum but at least there is a strong coupling with values larger than 250 cm -1 . This coupling is strong, as expected for the one-electron of difference between the states in their dominant configurations (t 4 2g e 2 g and t 3 2g e 3 g ). These ingredients induce to think that the most favourable transition is the 5 E g -→ 5 T 2g .

The evolution of the diabatic population (Fig. 3) shows that almost half of the population of the 5 E g is rapidly transferred to the 5 T 2g state. Then, part of the remaining population in the 5 E g evolves to the 5 T 2g when the wave-packet visits regions with small energy gap between both states. From, this state part of the population is gradually transferred to the 3 T 1g , above 10% in 2 ps, from which the LS may be populated. This last event is not appreciated in our dynamics because it happens in a larger time scale than our simulation time.

Our results indicates a faster 5 E g -→ 5 T 2g transition than experimentally observed. This unexpectedly fast rate is probably because of the limitation of the model Hamiltonian. The complete vibrational space is not represented but it is reduced to the three main coordinates. Moreover, the diabatic coupling terms between modes are not included in the model. These terms may be especially important in the description of Jahn-Teller modes. The lack of this terms gives stronger gradients in the model and the wavepacket may visit rapidly regions with small 5 T 2g -5 E g energy gaps with high efficient population transfer. Unfortunately, the inclusion of these terms dramatically increase the computational cost of these calculations. In addition, the spin-orbit coupling is fixed at the HS geometry and the variation along coordinates, vibronic spin-orbit coupling, is not included in our model.

The quintet-quintet transition observed in the wave-packet quantum dynamics is also determined by Fermi's golden rule. The calculated lifetime for the 5 E g at 0K is 2.3 ps, in good agreement with the first transition determined by Marino et al. in Fe(ptz) 6 . 17 Then, calculations at different level of theory point in the same direction. The first transition into the reverse-LIESST process is the relaxation to the 5 T 2g . After that, in a longer time scale, part of the population is transferred from the vibrationally hot HS state to the intermediate 3 T 1g and finally the 1 A 1g .

To conclude, we have applied two different methodologies to study the thermal and photophysical conversions from the highspin to the low-spin state. The order of magnitude of the HS state lifetime in the thermal relaxation process has been determined satisfactorily using Fermi's golden rule. The influence of the 3 T 1g as the main intermediate state coupling HS and LS has been demonstrated, while the role of the 3 MLCT states in this coupling can be neglected. This rate theory has been used to determine the HS state lifetime in Fe(II) complexes with nearly O h symmetry. This may be also tested to determine the HS state lifetime in complexes breaking the symmetry such as [Fe(terpy) 2 ] 2+ , where the empirical energy gap rule is not valid. 11 The determination of the adiabatic energies is the key step in this process and would be done using an appropriate DFT functionals or with the recent algorithms to account for the (XMS-)CASPT2 29,30 optimization implemented in quantum codes such as BAGEL.

The calculations indicate that the intersystem-crossing from the 5 E g to the triplets in the reverse-LIESST is not favourable. The mechanism that we propose according to the quantum dynamics is 5 E g -→ 5 T * 2g -→ 3 T 1g -→ 1 A 1g . The intersystem-crossing is produced from the vibrationally hot, 5 T * 2g which is formally as a thermal relaxation at an effective high temperature. This path could explain the limitation on the results observed in Fermi's golden rule when it should work fine because it happens from the minimum of the excited state due to the rate of the transition. This also would explain the two different rates reported in the reverse-LIESST and the quantum yield reported by Marino et al. 17 According to the LIESST, the rate for the 3 T 1g -→ 5 T 2g / 3 T 1g -→ 1 A 1g should be about 1.5 ps and the branching ratio 4:1. This does not compare with the 37 ps and the quantum yield of 0.1 associated to these transitions. For this reason and based on our results, we propose an alternative mechanism from the vibrationally hot 5 T * 2g state. Our results do not match completely with the experimental reported rates, this may probably be associated with a limitation in our model Hamiltonian. New improvements in the generation model Hamiltonian such as the inclusion of diabatic coupling terms or vibronic spin-orbit may be explored.

  J o u r n a l Na me Quantum dynamics simulations of the thermal and lightinduced high-spin to low-spin relaxation in Fe(bpy) 3 and Fe(mtz) 6 † Marc Alías-Rodríguez, * a,b Miquel Huix-Rotllant b and Coen de Graaf a,c

Fig. 1

 1 Fig. 1 Schematic diagram of the main electronic states involved in thermal HS-→LS transition in [Fe(bpy) 3 ] 2+ (top) and in the thermal and photoinduced relaxation in [Fe(mtz) 6 ] 2+ (bottom).

( a )

 a Average of the ten structures extracted from the singlet metadynamics. (b) Average of the ten structures extracted from the quintet metadynamics.

Fig. 3

 3 Fig. 3 Time-evolution population of the diabatic states along the dynamics. The states are represented by the colours: 1 A 1g black, 1 T 1g cyan, 3 T 1g green, 1 A 1g gold, 5 T 1g red and 5 E g violet.

Table 1

 1 Influence of the intermediate states in the singlet-quintet spinorbit coupling for [Fe(bpy) 3 ] 2+ . All the couplings and energies are expressed in cm -1 .∆E=E(HS)-E(LS) in cm-1 ISC*: CAS(10,15); active space formed by 10 electrons distributed over 15 orbitals. γ SµQ is the 2 nd order perturbative correction from state µ. H SO Sµ and H SO µQ are the average matrix elements ψ S | ĤSO |ψ µ and ψ µ | ĤSO |ψ Q , respectively. The labels a and b refers to conformations from the singlet and quintet metadynamics simulations, respectively.

	Conf.	ψ µ	H SO Sµ	H SO µQ	γ SµQ H SO QDPT	H SO eff	∆E
	ISC	3 T 1g 3 T 2g	529 284 20.8 55 232 1.0	21.8	32.7	-186
	ISC*	3 T 1g	548 449 27.9
		3 T 2g	56	248	1.0	29.0	29.9	-317
		3 MLCT	25	58	0.03
	1a	3 T 1g 3 T 2g	526 252 40.3 92 291 1.9	42.3	30.1 2615
	2a	3 T 1g 3 T 2g	516 464 32.2 123 275 2.2	34.4	28.2	696
	1b	3 T 1g 3 T 2g	524 465 19.4 70 294 1.1	20.6	43.0 -5623
	2b	3 T 1g 3 T 2g	519 452 23.8 99 273 1.6	25.4	37.0 -3168

Table 2

 2 Singlet quintet effective spin-orbit coupling in cm -1 at different conformations and its corresponding lifetime using different approaches in [Fe(bpy) 3 ] 2+ .

				T = 0K T=300K	T=0K
	Sys.	Conf.	H SO,eff	τ /ns	τ/ns	τ FC
		ISC	32.7	77.9	2.14	-
		ISC(10,15)	29.9	92.8	2.54	-
		1a	30.1	91.7	2.51	-
	bpy	2a Avg. (a)	28.2 28.1	104.7 106.1	2.87 2.9	--
		1b	43.0	45.0	1.24	-
		2b	37.0	60.7	1.66	-
		Avg. (b)	37.5	61.1	1.67	-
	mtz	ISC	31.4	8.9 s	-	∼10 years

| 1-8

J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Acknowledgements MHR acknowledges financial support by the Agence Nationale pour la Recherche through the project MULTICROSS (ANR-19-CE29-0018). MAR and CdG acknowledge the Spanish ministry of Science and Innovation (Project PID2020-113187GB-I00). The project leading to this publication has received funding from the Excellence Initiative of Aix-Marseille Université -A*Midex, a French "Investissement d'Avenir" programme, AMX-21-ERC-01, and from the funding "ERC Booster" from Région Sud Provence-Alpes-Côte d'Azur.

Notes and references

Conflicts of interest

The authors declare no conflicts of interest.