Stanislas Chambon
email: stanislas@rythm.co

N Mathieu

Galtier

Alexandre Gramfort
email: alexandre.gramfort@inria.fr

Domain adaptation with optimal transport improves EEG sleep stage classifiers

Keywords: EEG, sleep stage classification, domain adaptation, neural network, optimal transport

Low sample size and the absence of labels on certain data limits the performances of predictive algorithms. To overcome this problem, it is sometimes possible to learn a model on a large labeled auxiliary dataset. Yet, this assumes that the two datasets exhibit similar statistical properties which is rarely the case in practice: there is a discrepancy between the large dataset, called the source, and the dataset of interest, called the target. Improving the prediction performance on the target domain by reducing the distribution discrepancy, between the source and the target domains, is known as Domain Adaptation (DA). Presently, Optimal transport DA (OTDA) methods yield state-of-the-art performances on several DA problems. In this paper, we consider the problem of sleep stage classification, and use OTDA to improve the performances of a convolutional neural network. We use features learnt from the electroencephalogram (EEG) and the electrooculogram (EOG) signals. Our results demonstrate that the method significantly improves the network predictions on the target data.

I. INTRODUCTION

Training a supervised machine learning algorithm requires labeled data which might be limited or even absent in some cases. Using an auxiliary dataset, largely labeled, a.k.a source domain, enables to cope with this issue by training a predictive model that can then be used on a target domain. Yet, seldom are the source and target domains exactly identical in a statistical way: there is a domain discrepancy. Domain Adaptation (DA) methods [START_REF] Pan | A survey on transfer learning[END_REF] have been proposed to cope with this problem. The rationale behind DA is to find/adapt a feature representation that minimizes this discrepancy [START_REF] Ben-David | Analysis of representations for domain adaptation[END_REF].

DA methods can be divided into two categories. Unsupervised DA methods perform adaptation without using labels in the target domain, whereas semi-supervised DA methods use labels from source domain as well as from a few labeled target samples. Furthermore, DA methods can be separated into two types: deep DA methods which aim at relearning the feature representation [START_REF] Tzeng | Adversarial discriminative domain adaptation[END_REF], [START_REF] Ganin | Domain-adversarial training of neural networks[END_REF] while shallow DA methods aim at reducing the domain discrepancy by applying an operation, typically linear, on a fixed feature representation [START_REF] Sun | Return of frustratingly easy domain adaptation[END_REF]- [START_REF] Pan | Domain adaptation via transfer component analysis[END_REF].

Among shallow DA methods, optimal transport methods have exhibited promising performances on a wide range of This work was supported by the French Association Nationale de la Recherche et de la Technologie under Grant 2015/1005 problem [START_REF] Courty | Domain adaptation with regularized optimal transport[END_REF], [START_REF] Courty | Optimal transport for domain adaptation[END_REF], in particular on some EEG signals [START_REF] Gayraud | Optimal Transport Applied to Transfer Learning For P300 Detection[END_REF]. The principle underlying these OTDA methods is finding a mapping that maps source samples into the convex hull of target samples such that the probability distribution of the mapped samples is similar to the target probability distribution. These proposed OTDA methods build on recent progress on computational optimal transport, such as entropic regularization which leads to a fast Sinkhorn-Knopp's algorithm [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF].

Sleep stage classification, a.k.a. sleep scoring, is of great interest for understanding sleep and its disorders. Indeed, the sequence of sleep stages over a night is often used in the diagnosis of sleep disorders such as narcolepsy [START_REF] Stephansen | The use of neural networks in the analysis of sleep stages and the diagnosis of narcolepsy[END_REF]. Traditionally, this exam is performed as follows. First a subject sleeps with a medical device which performs a polysomnography (PSG). It consists of electroencephalography (EEG) signals at different locations over the head, electrooculography (EOG) signals, electromyography (EMG) signals, and eventually more. Second, a human sleep expert looks at the different time series recorded over the night, and assigns to each 30 s time segment a sleep stage following a reference nomenclature such as American Academy of Sleep Medicine (AASM) rules [START_REF] Iber | The AASM Manual for the Scoring of Sleep and Associated Events: Rules, Terminology and Technical Specifications[END_REF]. The AASM rules identify 5 stages: Wake (W), Rapid Eye Movements (REM), Non REM1 (N1), Non REM2 (N2) and Non REM3 (N3) also known as deep sleep. They are characterized by distinct time and frequency patterns. They also differ in proportions over a night. Sleep scoring is a tedious and time consuming task which is furthermore subject to the scorer subjectivity and variability [START_REF] Stephansen | The use of neural networks in the analysis of sleep stages and the diagnosis of narcolepsy[END_REF], yet it gives access to clinically relevant information about a patient.

From a machine learning perspective, the problem of automatic sleep scoring is an imbalanced multi-class classification problem. Several deep learning approaches have been proposed to learn an appropriate feature representation and classify sleep stages based on raw PSG signals or a timefrequency representation of these signals [START_REF] Stephansen | The use of neural networks in the analysis of sleep stages and the diagnosis of narcolepsy[END_REF], [START_REF] Chambon | A deep learning architecture for temporal sleep stage classification using multivariate and multimodal time series[END_REF], [START_REF] Biswal | SLEEPNET: automated sleep staging system via deep learning[END_REF]. Despite leading to state-of-the-art performances on various datasets, training a deep architecture requires a large quantity of data [START_REF] Chambon | A deep learning architecture for temporal sleep stage classification using multivariate and multimodal time series[END_REF] which might not be available, typically for particular clinical populations.

In this context, training a deep network on a large labeled source domain seems to be a good compromise. Yet, applying this method directly to records from a target domain is likely to give lower performances than expected. Indeed cross-dataset generalization is a difficult problem and is referred to as the dataset bias problem [START_REF] Torralba | Unbiased look at dataset bias[END_REF]. Domain adaptation is therefore a crucial step to improve the generalization performance of the model on a target domain.

Our contribution is a benchmark of OTDA approaches versus standard DA methods to enhance generalization performances of "source domain trained" deep sleep stage classifier on a target domain. In this paper, (i) we recall the OTDA principles, (ii) we present the general approach to adapt sleep stage classifier with OTDA, (iii) we perform a benchmark of DA methods on the sleep stage classification task between a large source domain and a target domain. All data are publicly available.

II. METHOD Notation:

With n ∈ N an integer, 1 n (resp. 0 n ) stands for the vector of R n composed of 1 (resp. 0). Let A, B be two matrices of same dimensions, A, B = tr(A T B). Let C ∈ N be the number of possible classes to predict and Y = {1, • • • , C}. X = R d stands for the feature space in dimension d. X is equipped with the Euclidean norm • 2 . For x ∈ X , δ x stands for the Dirac measure at position x. Let n s , n t ∈ N be the number of source training samples and target samples. (X s , Y s ) (resp. (X t , Y t )) stand for the n s source (resp. n t target) samples. They are assumed to be uniformly drawn from a source (resp. target) probability distribution µ s (resp. µ t ). Let P ∈ R ns×nt + be a matrix such that P 1 nt = 1 ns and P T 1 ns = 1 nt . The entropy of P ∈ R ns×nt is defined as: H(P ) = -ns,nt i,j=1 P ij (log P ij -1). I y stands for the set of source samples with label y: I y = 1 ≤ i ≤ n s : Y i s = y . Finally, • p 1 stands for the 1 norm of a vector to the power of p ∈ R + .

A. Domain adaptation with optimal transport

We first recall the principles of OTDA methods as introduced in [START_REF] Courty | Domain adaptation with regularized optimal transport[END_REF], [START_REF] Courty | Optimal transport for domain adaptation[END_REF], [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. It consists in (i) solving a linear program (LP) to obtain a coupling matrix linking the source and target samples, (ii) expressing the transported source samples as a linear combination of the target samples, where the linear combination is given by the coupling matrix.

Consider n s source samples and n t target samples, with X s , X t ∈ X drawn i.i.d. from two probability distributions µ s and µ t . Both X s and X t are associated with two uniform

empirical distributions of probability μs = 1 n s ns i=1 δ X i s and μt = 1 n t nt i=1 δ X i t .
We identify them to the vectors 1/n s 1 ns and 1/n t 1 nt . OTDA methods involve a pairwise distance matrix M ∈ R ns×nt also called cost matrix or ground metric. It is common to define it as

M ij = X i s -X j t 2
2 . Furthermore OTDA methods involve a polytope of possible couplings between the empirical distributions μs and μt : Π(μ s , μt ) def = P ∈ R ns×nt + : P 1 nt = μs , P T 1 ns = μt . The first step consists in finding an optimal coupling P * which minimizes the transported mass P → P, M [START_REF] Courty | Domain adaptation with regularized optimal transport[END_REF], [START_REF] Courty | Optimal transport for domain adaptation[END_REF] P * ∈ arg min P ∈Π(μs,μt) P, M .

(

The second step consists in computing a transportation mapping T that transports each source sample into the convex hull of the target samples [START_REF] Courty | Domain adaptation with regularized optimal transport[END_REF], [START_REF] Courty | Optimal transport for domain adaptation[END_REF]:

T (X i s ) = nt j=1 P * ij X j t nt j=1 P * ij = n s P * i X t (2)
where P * i stands for the i-th row of P * .

B. Entropic regularization

In this paragraph, we recall how the LP part can be solved efficiently using an entropic regularization term. Solving the LP (1) involved in OTDA methods is computationally expensive. To cope with this problem [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] introduced an entropic regularization term to the initial problem. The new optimization problem reads:

P * = arg min P ∈Π(μs,μt) P, M -H(P ), > 0 . ( 3
)
The benefits of such a regularization are threefold: the problem is now strictly convex with a unique solution, the entropic regularization stabilizes the solution and smooths the transportation plan P * , and finally the solution of (3) can be obtained by iterative diagonal scaling of the Gibbs kernel K = exp(-M/ ) associated to M . This step is often referred to as Sinkhorn-Knopp's algorithm [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF].

C. Regularization on source labels

The authors of [START_REF] Courty | Domain adaptation with regularized optimal transport[END_REF], [START_REF] Courty | Optimal transport for domain adaptation[END_REF] introduced an additional regularization to the OT problem (3) to structure the coupling based on the labels of source samples. It is a regularization defined as Ω c (P ) = nt j=1 y∈Y P (I y , j) p 1 with p < 1 (0.5 in practice). It enforces sparsity in the coupling P so that source samples of the same label are transported onto the same target samples, and additionally source samples of opposite labels are transported onto different target samples. With > 0 and η > 0, the optimization problem now reads:

P * = arg min P ∈Π(μs,μt) P, M -H(P ) + ηΩ c (P ) (4)
We will refer in this manuscript to this approach as LpL1.

As the regularization is non-convex, it is solved iteratively by majorization-minimization [START_REF] Courty | Domain adaptation with regularized optimal transport[END_REF], [START_REF] Courty | Optimal transport for domain adaptation[END_REF]

D. General method

To demonstrate the benefits of the OTDA method for a deep learning model we propose the following approach: (i) select samples from the source domains and a record from the target domain (ii) extract the features from the penultimate layer (before the classification layer) for the source data and target data (iii) apply the OTDA method on extracted features to adapt the source features. (iv) retrain the network classification layer on the adapted features (v) predict on the target features.

III. EXPERIMENTS a) Data: We used MESA [START_REF] Dean | Scaling up scientific discovery in sleep medicine: The national sleep research resource[END_REF] as a source domain which contains EEG records with 3 EEG derivations (Fz-Cz, Cz-Oz, C4) et 2 EOG derivations (EOG left -Fpz, EOG right -Fpz). The target domain was MASS-session 3 [START_REF] O'reilly | Montreal Archive of Sleep Studies: an open-access resource for instrument benchmarking and exploratory research[END_REF]. We used 61 records and consider the derivations similar to the source domain derivations: Fz-Cz, Cz-Oz, C4 and EOG left -Fp1, EOG right -Fp1 (since Fpz is not available).

b) Features extraction: We followed the methodology described in [START_REF] Chambon | A deep learning architecture for temporal sleep stage classification using multivariate and multimodal time series[END_REF] for preprocessing. Here, 1200 (resp. 200) randomly chosen source records are used for training (resp. validation i.e. monitoring the progress of training). The deep learning model was the model introduced in [START_REF] Chambon | A deep learning architecture for temporal sleep stage classification using multivariate and multimodal time series[END_REF] with 5 EEG like channels. It was trained by minimizing a cross-entropy loss weighted by the inverse of the stage proportions in order to mimic the balanced sampling used in [START_REF] Chambon | A deep learning architecture for temporal sleep stage classification using multivariate and multimodal time series[END_REF]. The source (resp. target) EEG/EOG 30 s signals were propagated through the network and the activations of the penultimate layer were considered as source (resp. target) features. This leads to a feature space of dimension d = 600. c) Baselines: We compared our approach to various strategies. S stands for not using any adaptation method. CORAL (Correlation Alignment) [START_REF] Sun | Return of frustratingly easy domain adaptation[END_REF] consists of whitening the source data and coloring them with the target covariance matrix. SA (Subspace Alignment) [START_REF] Fernando | Unsupervised visual domain adaptation using subspace alignment[END_REF] consists of projecting the source and target data onto two subspaces of dimension d SA < d, and aligning the source subspace axes onto the target subspace axes. The hyperparameter d SA was optimized in {25, 50, 100} via a grid search. OT (resp. LpL1) has 1 hyperparameter (resp 2 hyperparameters , η) that was also optimized by grid search over 10 -1 , • • • 10 3 (resp.

10 -1 , • • • 10 3 × 10 0 , • • • 10 3 ). d) Classifier: We used the classification layer of [START_REF] Chambon | A deep learning architecture for temporal sleep stage classification using multivariate and multimodal time series[END_REF]'s model as a classifier, a fully connected layer with a softmax activation. The input dimension is equal to the feature representation dimension (either d or d SA for SA) and the output dimension is equal to the number of classes C = 5. For each DA method, the network is trained on adapted source samples by minimizing the categorical cross entropy with stochastic gradient descent and balanced sampling [START_REF] Chambon | A deep learning architecture for temporal sleep stage classification using multivariate and multimodal time series[END_REF]. The training progress is monitored on the validation source data. An early stopping procedure is used to stop the training when no improvement is observed on the validation set after 5 consecutive pass on the data. Once trained, the model predicts the sleep stages on the (adapted) features of a target record.

e) Cross-validation: The classifier is trained on 2500 (adapted) source training samples (500 per class). 2500 additional (adapted) source validation samples are used to monitor training progress and perform early stopping. The classifier then predicts on an unseen target record. The DA methods are fitted on the source training samples and the considered target record. The process is repeated over the 61 target records.

The hyperparameters of SA, OT and LpL1 are optimized using target validation data. This set is composed of 1000 (adapted) target samples from 10 records. For each set of hyperparameters, (1) the DA method is applied to adapt source training / validation samples, (2) the classifier is trained, (3) it is evaluated on (adapted) target validation data. The set of hyperparameters obtaining the highest performance on validation target data is selected. f) Metrics: The balanced accuracy is used as a general metric. Per-class metrics such as normalized (per row) confusion matrix, F1 score, precision, sensitivity and specificity are furthermore used to characterize the effects of DA methods. Each metric is computed per target record. The presented performances are averaged over the target records.

g) Code: We used POT library (https://pot.readthedocs. io/en/latest) for OT algorithms, PyTorch library [START_REF] Paszke | Automatic differentiation in PyTorch[END_REF] for the deep learning part and scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] for general purposes.

h) Benchmark: Besides results with the DA methods, we present the optimal performance that a classifier could reach, should we had access to the labels for 40 records in the target domain. The averaged and per-record balanced accuracies are reported in Figure 1. We observe that OTDA methods outperform SA and CORAL by delivering a gain of performance larger than 5% on the averaged balanced accuracy. Performance per record indicates that most records exhibit higher performances with OTDA methods.

The confusion matrices and the per-class metrics associated to each DA method are reported in Figure 2. They indicate that OTDA methods improve the detection of all stages: the matrices for OT and LpL1 are more diagonal. They also show that OTDA methods offer a good compromise of sensitivity, specificity and precision for each stage. On the other hand, the other DA methods do not exhibit general good performance as they under perform on some particular stages. Indeed SA exhibits low precision and sensitivity on N1 and low precision on REM. CORAL exhibits low sensitivity on N1 and N2 and low precision on N3.

i) Influence of hyperparameters: We focus here on evaluating the influence of the hyperparameters , η for LpL1. We report the performances in Figure 3-A when varying in 10 -1 , • • • , 10 3 and η ∈ 0, 10 -1 , • • • , 10 3 . We also report the performances obtained by setting , η to the maximum value taken by any feature in either source or target domains. We refer to this as the maximum heuristic.

Both and η have influence on the performances of LpL1 but for each η value, there is a large range of leading to similar and good performances, see Fig 3-A. The maximum

Fig. 1 .

 1 Fig. 1. General benchmark of DA methods: OTDA methods address well the DA problem and recover half of the lost performance. A: averaged balanced accuracy. B: balanced accuracy per record (1 dot is 1 record). Dashed-line: optimal performance of the classifier, should we have labels for 40 records of the target domain.

LpL1 precision 0.98 0.38 0.94 0.50 0.52 0.96 0.38 0.94 0.52 0.55 0.97 0.31 0.89 0.23 0.55 0.60 0.28 0.93 0.53 0.68 0.62 0.27 0.94 0.55 0.69 sensitivity 0.51 0.11 0.60 0.98 0.99 0.51 0.16 0.65 0.96 0.98 0.45 0.02 0.24 1.00 0.80 0.83 0.46 0.47 0.94 0.81 0.85 0.46 0.50 0.95 0.83 specificity 0.95 0.93 0.71 1.00 1.00 0.95 0.93 0.73 0.99 0.99 0.95 0.92 0.56 1.00 0.95 0.98 0.95 0.64 0.99 0.95 0.98 0.95 0.66 0. Importantly, the maximum heuristic avoids the need to run a grid search, and does not require labeled target samples at all.

CONCLUSION

In this paper, we focused on improving the performances of an automatic sleep stage classifier when applied to a different target domain. The proposed approach based on a convolutional network, and built on optimal transport, clearly showed improvements on cross-domain generalization, a problem which is relevant for any predictive model applied to new datasets exhibiting distribution shifts.


Hélène Le Bouder
email: lebouder@emse.fr

Gaël Thomas
email: gael.thomas@xlim.fr

Yanis Linge

Assia Tria

On Fault Injections in Generalized Feistel Networks

Keywords: Fault injections attack, differential fault analysis, Generalized Feistel Networks

In this paper, we propose a generic method to assess the vulnerability to Differential Fault Analysis of generalized Feistel networks (GFN). This method is based on an in-depth analysis of the GFN properties. First the diffusion of faults is studied, both at the block level and at the S-box level, in order to have a fault which maximizes the number of S-boxes impacted by a fault. Then the number of faults in an S-box required to find the key is evaluated. By combining these results, a precise assessment of the vulnerability to fault attacks of GFN can be made. This method is then used on several examples of Feistel ciphers.

I. INTRODUCTION

Feistel and generalized Feistel networks are a family of widespread block symmetric ciphers. Fault injection attacks consist in disturbing the circuit behaviour in order to alter the correct processing of the algorithm. Among those, the Differential Fault Analysis (DFA) is a cryptanalytic technique that exploits erroneous results of fault injections. Since injecting faults on a cryptographic circuit may irreversibly damage the device, limiting the number of fault injections is often a critical aspect of DFA. In this paper, a generic method is proposed to evaluate the vulnerabilities of the block registers in generalized Feistel networks. This method is based on the study of the diffusion properties both of the function used inside the Feistel network and of the Feistel network itself.

The fewer single-bit faults are required in a block to find the key, the more vulnerable this block is. This knowledge allows an attacker to choose the best target for its fault injection, but it also enables the designer to focus countermeasures on the vulnerable spots.

No fault attack is presented in this paper, we focus on the diffusion properties to assess a "potential" vulnerability but do not perform a complete fault attack since it is not our objective.

The paper is organised as follows. Section II gives the context of generalized Feistel networks and fault attack injections on them. Our approach is described in Section III. Section IV presents our results on some examples of known generalized Feistel networks. Finally, the conclusion is drawn in Section V.

II. CONTEXT AND STATE OF THE ART

A. Feistel Networks 1) Definition: A Feistel network is an iterative construction allowing to create a bijection from functions that need not be bijective. It was invented at IBM by Horst Feistel in the 1970's, during the development of the Lucifer cipher, later known as the Data Encryption Standard (DES) [18]. A Feistel network divides the input into two blocks of equal size. One of these two blocks goes into a subkey-dependent function, called the Feistel function, and the result is xored to the other block, as depicted on Fig. 1. This construction is then repeated several times, switching the roles of the two blocks each time. The main characteristics of this construction are:

• only half of the state is modified, using the Feistel function.

• the construction is invertible, even if the Feistel function itself is not. • the inverse bijection follows the same construction up to reversing the subkeys order.

L r R r F L r+1 R r+1
Fig. 1. A classical (two blocks) Feistel cipher.

Introduced at CRYPTO'89 by Zheng et al. [START_REF] Zheng | On the Construction of Block Ciphers Provably Secure and Not Relying on Any Unproved Hypotheses[END_REF], Generalized Feistel Networks (GFNs) are a broader class of schemes taking the main ideas of the Feistel network but dividing the internal state into more than two blocks. As done in [START_REF] Berger | Extended Generalized Feistel Networks using Matrix Representation[END_REF], one round of a GFN is divided into two successive transformations: the nonlinear layer and the permutation layer. The non-linear layer is made of one or more key-dependent Feistel functions F. They have some of the blocks of the GFN internal state as input and their output is xored to some other blocks. The permutation layer consists in rearranging the different blocks of the internal state. An example of GFN is given on Fig. 2.

B 0 r B 0 r+1 B 1 r B 1 r+1 B 2 r B 2 r+1 B 3 r B 3 r+1 B 4 r B 4 r+1 B 5 r B 5 r+1 B 6 r B 6 r+1 B 7 r B 7 r+1 F F F F Fig. 2.
Example of a generalized Feistel network with 8 blocks.

The blocks of a generalized Feistel network are denoted B i r where i is the register/block index (ranging from 0 to b -1 for b blocks ) and r is the current round (ranging from 0 to r for r rounds).

2) The Feistel function: We suppose that the Feistel function F is made of three kinds of component.

• The Xor with the subkey. Since in a single round, more than one block can go into an F function, we denote the different subkey blocks K λ r with λ ∈ [[0, Λ -1]] for Λ subkey blocks. For example, the GFN on Fig. 2 has Λ = 4 subkey blocks.

• The S-boxes. They are the non-linear part of the GFN and provide Shannon's confusion. • The "mixing" functions. They are linear and provide diffusion inside a block. It is the most famous Feistel cipher and is depicted on Fig. 3. Prior to the 16-round Feistel network itself, a bit-wise permutation IP is applied to the 64-bit internal state. The inverse permutation IP -1 finishes the encryption process. The Feistel function F is illustrated on Fig. 4. It acts on 32 bits and consists in 4 steps:

B. Examples of

• Expansion E which maps 32 bits into 48 bits by duplicating half of the bits. • Xor with the 48-bit subkey K r , r ∈ [ [START_REF] Beaulieu | The SIMON and SPECK Families of Lightweight Block Ciphers[END_REF][START_REF] Suzaki | TWINE : A Lightweight Block Cipher for Multiple Platforms[END_REF]].

• S-boxes S i , i ∈ [ [START_REF] Beaulieu | The SIMON and SPECK Families of Lightweight Block Ciphers[END_REF][START_REF] Leander | On the Classification of 4 Bit S-Boxes[END_REF]], which substitute a 6-bit input for a 4-bit output. • Bit-wise permutation P of 32 bits.

2) MIBS: MIBS is a cipher presented at CANS'09 [START_REF] Izadi | MIBS: A New Lightweight Block Cipher[END_REF] and aimed at low-resource environments such as RFID tags or sensor networks. It is made of 32 iterations of a classical Feistel networks (see Fig. 1) acting on two blocks of 32 bits each and accepts key of 64 or 80 bits.

As depicted on Fig. 5, MIBS Feistel function operates on 32 bits divided into 8 4-bit words (nibbles) and is made of three consecutives steps:

• Xor with the subkey.

IP Permutation Initiale (64 bits)

L 0 R 0 F L 1 R 1 F L 2 R 2 L 15 R 15 F R 16 L 16
IP -1 Permutation Finale (64 bits) • a layer of 8 parallel 4 × 4 S-boxes.

Fig. 3. Overview of Data Encryption Standard. L r R r E K r S 8 S 7 S 6 S 5 S 4 S 3 S 2 S 1 P R r+1 L r+1
• a linear mixing layer M C acting at nibble level.

3) TWINE: TWINE is a 64-bit block cipher presented at SAC '12 [START_REF] Suzaki | TWINE : A Lightweight Block Cipher for Multiple Platforms[END_REF] aiming at a reduced circuit area while avoiding non-software-friendly operations (such as bit permutations) to maintain decent performances when implemented in software. It is a GFN with 16 blocks, 4 bits each, that can accept either 80-bits or 128-bits keys. The permutation layer that rearranges the blocks at the end of each round was chosen to maximise diffusion between blocks, according to a result of [START_REF] Suzaki | Improving the Generalized Feistel[END_REF]. Compared to the cyclic shift, this permutation requires twice less rounds for an input difference to influence all the blocks. One round of TWINE is depicted on Fig. 6. The whole encryption process consists in 36 rounds of this scheme for both key lengths. The Feistel functions is used 8 times per round and is consecutively made of:

L r R r K r S 7 S 6 S 5 S 4 S 3 S 2 S 1 S 0 M C L r+1
• a 4-bit Xor with a subkey block.

• a single 4 × 4 S-box. 4) CLEFIA: CLEFIA is a 128-bit block cipher presented at FSE '07 [START_REF] Shirai | The 128-Bit Blockcipher CLEFIA (Extended Abstract)[END_REF] with key sizes 128, 192 or 256 bits. It is a GFN with 4 blocks, 32 bits each, as depicted on Fig. 7. At each round, it uses 2 slightly different Feistel functions F 0 and F 1 . Function F 0 is given on Fig. 8. They are both made of three consecutive steps:

• Xor with the subkey, • 4 parallel 8 × 8 S-boxes. Both functions F 0 et F 1 uses the two same S-boxes S 0 and S 1 but in different orders.

B 0 r B 1 r B 2 r B 3 r F 0 F 1 B 0 r B 1 r+1 B 2 r+1 B 3
r+1 Fig. 7. Overview of CLEFIA.

• two linear diffusion layers, M C 0 for F 0 , M C 1 for F 1 based on finite field operations (much like AES Mix-Columns). C. Attacks using Fault Injections 1) Fault Injections: Fault injection attacks consist in disturbing the circuit behaviour in order to alter the correct processing of the algorithm. Unfortunately, in fault injection attacks the circuit can be damaged. Sakiyama et al. in [START_REF] Sakiyama | Information-Theoretic Approach to Optimal Differential Fault Analysis[END_REF] emphasize the importance of limiting the number of fault injections. So, our analysis is conducted with this goal in mind.

B 0 r B 1 r S0 S1 S0 S1 M C 0 B 0 r+1 K 1 r
By design, algorithms are made secure against classical attacks. For example, in the case of the classical Feistel cipher, Luby and Rackoff demonstrated in [START_REF] Luby | How to Construct Pseudo-Random Permutations from Pseudo-Random Functions[END_REF] the following result: no efficient algorithm exists to distinguish a 4-round Feistel cipher with pseudo-random Feistel functions from a random permutation with only the observation of (plaintext, ciphertext) pairs. Despite this security, an algorithm may be vulnerable to physical attacks. There are two families of physical cryptanalysis : SCA (Side Channels Analysis) attacks; and FI (Fault Injection) attacks. In this paper we will only focus on fault injections attacks. One of the first fault attack was proposed against a Feistel cipher (the DES) by Biham and Shamir in [START_REF] Biham | Differential Fault Analysis of Secret Key Cryptosystems[END_REF] in order to recover the cipher-key.

In our approach, the fault attack changes one single-bit of a block B i r . The single-bit fault model is a common one which account for an attack where a low stress injection method is used.

2) DFA on a Generalized Feistel Network: The Differential Fault Analysis (DFA) is a powerful cryptanalytic technique that exploits differences between the correct ciphertext and erroneous results due to fault injections. The attack path for DFA on a Feistel cipher is classically done on the last round as in [START_REF] Biham | Differential Fault Analysis of Secret Key Cryptosystems[END_REF] for DES [18]. The attack path is shown on Fig. 9 and is repeated for each subkey block K λ r . From now on, variables with a superscript * (star) denote the faulted version of a data and the letter ∆ is used to denote the difference between the correct value and the faulted one. E.g. B i * r-1 denotes the faulted version of B i r-1 , and block B i r-1 as input and has its output xored with B j r-1 to produce B o r , we have the following relation:

∆ B i r-1 = B i * r-1 ⊕ B i r-1 . Thus, if F has B j r-1 , B j * r-1 B i r-1 , B i * r-1 F B o r , B o * r K λ r
∆ B o r = B o r ⊕ B o * r = F B i r-1 , K λ r ⊕ F B i * r-1 , K λ r ⊕ ∆ B j r-1 (1)
However it is possible for the subkey to be xored at the end of the Feistel function. In this case, we have:

∆ B o r = F B i r-1 ⊕ K λ r ⊕ F B i * r-1 ⊕ K λ r ⊕ ∆ B j r-1 = F B i r-1 ⊕ F B i * r-1 ⊕ ∆ B j r-1
This equation does not depend on the subkey, so another attack path must be used instead. This case happens for example in the ciphers SIMON [START_REF] Beaulieu | The SIMON and SPECK Families of Lightweight Block Ciphers[END_REF] or Piccolo [START_REF] Shibutani | Piccolo: An Ultra-Lightweight Blockcipher[END_REF]. Another case to deal with is when blocks B i r-1 and B i * r-1 are masked by a whitening key that is xored between the last round of the Feistel construction and the ciphertext as in Piccolo [START_REF] Shibutani | Piccolo: An Ultra-Lightweight Blockcipher[END_REF] for example. These two particular cases are not taken into account in this study but are not mandatory for fault injections on GFNs to work, see e.g. [START_REF] Jeong | Differential Fault Analysis on Block Cipher Piccolo[END_REF] for a successful attack on Piccolo.

Physical attacks are based on a divide-and-conquer strategy: the targeted subkey block K λ r is cut into L pieces. That is, the whole subkey is then cut into Λ blocks K λ r , each block being cut into L pieces K λ,l r . How this second division is done depends on the S-box layer: one piece of subkey per S-box. An hypothesis on K λ,l r is denoted k and k denotes the correct value.

III. DESCRIPTION OF THE GENERIC APPROACH

In this section a generic method is proposed to assess the vulnerability of a particular block in a GFN. In particular the most vulnerable block B e r , the best suited for injecting singlebit faults, is searched (e is the block index and r is the round number). A block is the most vulnerable when injecting singlebit faults at this location minimizes the total number of faults required to find the key. We wish to have a good diffusion of the fault as to attack as many subkey bits as possible for each fault injection. A priori, this depends on two things: the overall structure of the Feistel network but also the Feistel function used.

A. At the Feistel network level

At the Feistel network level, the diffusion of faults in the blocks are studied which allows to predict the propagation of faults in the blocks.

1) Full diffusion delay: We define here the notion of (blockwise) diffusion in a generalized Feistel network, as given in [START_REF] Berger | Extended Generalized Feistel Networks using Matrix Representation[END_REF]. A block

B i 0 influences block B j r if B i 0 effectively appears when expressing B j r as a function of B 0 0 , • • • , B b-1 0 . A block B i 0 has diffused at round r if B i 0 influences all the B j r for j ∈ [[0, b -1]
]. If all the blocks B i 0 have diffused at round r, the Feistel network is said to have reach full diffusion, i.e. each output block B j r depends on all input blocks B i 0 . The minimum number of rounds required to reach full diffusion is called full diffusion delay, and denoted d. By definition of the full diffusion delay, it is unnecessary to inject the fault before the d latest rounds as doing so will not further propagate the fault but may complexify the analysis. We then focus on the blocks

B e r with e ∈ [[0, b -1]], r ∈ [[r -(d + 1), r -1]
] where e indicates the block index where the single-bit fault has been injected.

2) Matrix Representation of Feistel Networks: At SAC '13, Berger et al. [START_REF] Berger | Extended Generalized Feistel Networks using Matrix Representation[END_REF] presented a general framework to represent generalized Feistel networks using matrices. The matrix representation of a GFN is the adjacency matrix of the graph representing influence relations in the GFN: when block B i r+1 is influenced by B j r , a non zero coefficient at position (i, j) in the matrix is written. The non zero coefficients of this matrix can either be 1 or F depending on whether the influence is direct or through a Feistel function F. They then compute powers of that (adjacency) matrix to determine the full diffusion delay. Technically, the computation is done over the coefficient ring of polynomials Z[F]. In this paper, we reuse this representation but in a slightly improved way.

Let D be N ∪ {-∞} endowed with the two internal laws max and +, defined over the integers N and extended to D in the obvious way: ∀a ∈ D, max(a, -∞) = max(-∞, a) = a and (-∞) + a = a + (-∞) = -∞. Then D is an idempotent commutative semiring, also known as commutative dioid. A semiring is an algebraic structure very similar to a ring except that the first law (max in our case) forms only a commutative monoid instead of a commutative group, i.e. the requirement to have a first-law inverse for every element is discarded. A semiring is said idempotent when ∀a ∈ D, max(a, a) = a, and commutative when the second law (+ in our case) is commutative. A dioid is simply an idempotent semiring.

This object represents degrees of polynomials with the convention that the degree of the zero polynomial is -∞. In our case, a polynomial in the Feistel function F is studied. We then define the dioidic matrix representation M of a GFN to be the matrix with coefficients in D such that coefficient at position (i, j) can either be:

• 0, when block B i r+1 is influenced by B j r directly. • 1, when block B i r+1 is influenced by B j r via the Feistel function F.

• -∞, otherwise, i.e. not influenced. For example, the matrix M of the GFN on Fig. 2 is given on Fig. 10. The link with the matrix representation introduced in [2] is as follows. The coefficients in our matrix represent the degree of the coefficient at the same position in the matrix of [START_REF] Berger | Extended Generalized Feistel Networks using Matrix Representation[END_REF]. That is, we are interested in how many Feistel functions are passed through to go from an input block to an output block, which is represented by the degrees in the polynomials in the Feistel function F.

M =             1 0 -∞ -∞ -∞ -∞ -∞ -∞ -∞ -∞ 0 -∞ -∞ -∞ -∞ -∞ -∞ -∞ 1 0 -∞ -∞ -∞ -∞ -∞ -∞ -∞ -∞ 0 -∞ -∞ -∞ -∞ -∞ -∞ -∞ 1 0 -∞ -∞ -∞ -∞ -∞ -∞ -∞ -∞ 0 -∞ -∞ -∞ -∞ -∞ -∞ -∞ 1 0 0 -∞ -∞ -∞ -∞ -∞ -∞ -∞            
3) Using the Matrix Representation for Fault Injections: For a subkey block K λ r to be attacked using fault injections, the block B i r-1 at the input of the λ-th Feistel function F must be faulted. Thus the block where the fault is done must influence that block B i r-1 . We recall that due to the structure of a (generalized) Feistel network, the input of a Feistel function at round r -1 must be outputted as a ciphertext block to have an invertible construction.

Let B e r be the block where the fault is injected and V the e-th standard basis vector (which has then -∞ coefficients on every coordinates except the e-th one which is 0). The vector V F defined by the vector matrix product given on equation (2) allows to know for each block at the output of the penultimate round whether the fault injected in B e r has influenced that particular block and how many Feistel functions F has been crossed doing so.

V F = M r-(r+1) • V (2)
If at the index i, corresponding to block B i r-1 , we observe:

• -∞: the block is not faulted.

• 0 : the block contains the single-bit fault.

• x ∈ N : the block is faulted and the fault has been through

x Feistel functions.

Hence for each fault, the number n λ of subkeys blocks attacked after r rounds can be computed.

B. At the Feistel function level

Recall that a subkey block is split into L pieces (one per S-box). In the Feistel function, a non-linear layer is necessary. This layer is composed of one or several (L) Sboxes. This division has an impact on the propagation of the fault inside the Feistel function. Since an S-box must be faulted for the attacker to be able to gain knowledge on the corresponding subkey piece, knowing into which S-boxes the fault is propagated is of the utmost importance.

Studying the Feistel function allows to identify how many pieces of the subkey blocks are attacked by each fault injection. This number is denoted n l .

One should also be able to find all the possible differentials ∆ B j r-1

. Indeed, as B j r-1 might be faulted, knowing the different values of ∆ B j r-1 is necessary in order to be able to use equation [START_REF] Beaulieu | The SIMON and SPECK Families of Lightweight Block Ciphers[END_REF] to retrieve the key. The rest of this section deals with the diffusion properties of the Feistel function.

1) Number of pieces of subkey attacked: The first thing to determine from the Feistel function is the number of pieces of the subkey which are attacked, i.e. the number of faulted S-boxes. Looking at the Feistel function, it is possible to know what output bits can be faulted with a single-bit fault at the input. And, since we know what S-boxes depend on which input bits of the Feistel function, it is possible to deduce which S-boxes will be faulted in the next Feistel function.

So from the diffusion of the bits in the Feistel function it is possible to compute the number n l of pieces of the subkey block attacked by a fault injection, depending on the number of Feistel functions crossed which is V F (i), the i-coordinate of vector V F . Care should be taken not to forget the case where the fault has not been through any function. For example in DES, a single-bit fault can impact up to 2 S-boxes because of Expansion E.

2) Search for the values of ∆ B j r-1

:

If the fault is not propagated to B j r-1 , the value of ∆ B j r-1
in equation ( 1) is equal to 0. In the other case, the value of ∆ B j r-1 is unknown.

V F (j), the j-th coordinate of vector V F , indicates if B j r-1 has been faulted and how many times the fault passed through the Feistel function.

Simply put, the study of the diffusion of bits in the Feistel function allows to evaluate how many bits are potentially faulted in B j r-1 when the number of times the injected fault crossed the Feistel function is known. For example, if the fault do not cross any Feistel function, the fault being single-bit, there are as many possible differentials as bits in B j r-1 . As a consequence, it is possible to evaluate the information brought by the fault since we can evaluate the incertitude on ∆ B j r-1 .

But this approach can be improved by studying the actual values of the possible differentials. For that purpose, the single faulted bit in B e r is searched. Several hypotheses are often possible.

In order to find it the diffusion result on the next round is used, i.e. the value:

W F = M • V F
The minimal but finite value (i.e. in N) in the vector W F corresponds to the faulted block where the fault crossed the Feistel function a minimum number of times.

If this value is equal to 0, the single-bit fault can be observed in the ciphertext on the corresponding block.

If the value is equal to 1, the differential of the corresponding block ∆ B min(W F ) r is computed. For all single-bit faults, it is often possible to compute exhaustively all possible differentials at the output of a Feistel function F. As a consequence, knowing ∆ B min(W F ) r restricts the set of possible single-bit faults that have been injected.

If the minimal value in W F is greater or equal to 2, the above method might be still possible. However, since for each intermediate difference tuple (i.e difference after the first, second,. . . , penultimate function), all the output differences (after the last function) must be computed, the complexity of such an algorithm grows exponentially fast in the size of one block and in the value of min(W F ) and quickly becomes infeasible. Besides, due to Feistel function being designed for maximazing diffusion inside a block, it is likely that this will not help restricting the possible single-bit faults injected any longer.

If ∆ B min(W F ) r allows to restrict the set of fault candidates, each one of these candidates is propagated to obtain a reduced set of potential differentials ∆ B j r-1 . The union of all possible differentials (for all fault candidates) indicates the number of possible differentials ∆ B j r-1

. Yet this approach is not perfect since two differentials can be xored in the datapath. This problem has no consequences for us since we only try to evaluate the number of faults required to find the key and we do not construct a fault attack. But a cryptanalyst would have to take it into account to conceive a fault attack.

C. At the S-box level

In this section the information on the key brought by a fault crossing an S-box is evaluated. The value n, the number of faults in an S-box to find the corresponding subkey piece, is computed. The attack path is reduced to the S-box level, considering that what happens in the Feistel function F before the Xor with the subkey is known and that the S-boxes are the only non-linear part of the function. This case is summed up by equation [START_REF] Biham | Differential Cryptanalysis of DES-like Cryptosystems[END_REF], where x the input of the S-box, x * the faulted input and ∆ y the output differential are known.

S(x

⊕ k) ⊕ S(x * ⊕ k) = ∆ y (3)
A method inspired by classical differential cryptanalysis is used. Differential cryptanalysis is a chosen plaintext attack introduced by Biham and Shamir in 1991 [START_REF] Biham | Differential Cryptanalysis of DES-like Cryptosystems[END_REF]. It consists in considering plaintexts pairs having a fixed difference and in studying how this difference propagates throughout encryption.

If we write e = x ⊕ x * in equation ( 3), we have:

S(x ⊕ k) ⊕ S(x ⊕ e ⊕ k) = ∆ y (4)
Each time another fault is injected, values of x, e and ∆ y change while k stays the same. We are then interested in the minimal number n of experiences (fault injections) such that being given

x 0 , • • • , x n-1 , e 0 , • • • , e n-1 and ∆ y0 , • • • , ∆ yn-1
, there is only one key for which the n equations (4) simultaneously hold. The S-box S is supposed to be a vectorial boolean function from

F u 2 into F v 2 . Let a ∈ F u 2 and b ∈ F v 2 , S a,b
denotes the set of all inputs z such that if an input difference a is applied, an output difference b is observed. Said otherwise:

S a,b = {z ∈ F u 2 | S(z) ⊕ S(z ⊕ a) = b}.
By symmetry of this equation, #S a,b is even. Let δ denotes the base-2 logarithm of the maximum size of the sets S a,b taken over all a = 0 and b, that is:

2 δ = max a =0, b #S a,b
Since the S a,b 's are of even size and not all empty, it follows that δ ≥ 1. More generally, δ ≥ u -v. Indeed, let a be a nonzero fixed input difference, then b∈F v 2 #S a,b = 2 u . Hence there exists a b so that #S a,b ≥ 2 u /2 v and thus δ ≥ u -v. Finally, if x is drawn uniformely at random in F u 2 then

Pr [x ∈ S a,b ] = #S a,b 2 u ≤ 2 δ-u .
A key k is then a solution of equation ( 4) if and only if x ⊕ k ∈ S e,∆y . So the number of solutions of equation ( 4) is #S e,∆y . The set of z such that z ⊕ x ∈ S e,∆y is denoted by x ⊕ S e,∆y . Then k is solution of equation ( 4) if and only if k ∈ x ⊕ S e,∆y . For n fault injections, we are then interested in the probability

Pr k ∈ n-1 s=0 (x s ⊕ S es,∆y s ) .
The faults e s are independent from each others. The ∆ ys and x s are independent and identically distributed. That is, they are random variables with the same probability law and are mutually independent. Thus we can write:

Pr k ∈ n-1 s=0 x s ⊕ S es,∆y s = n-1 s=0 Pr k ∈ (x s ⊕ S es,∆y s ) ≤ 2 n(δ-u)
Besides, this probability is lower bounded by 2 -u since there exists a solution. Hence the number of faults n searched, the number of faults on the S-box required to find the corresponding subkey piece, is such that

2 n(δ-u) ≤ 2 -u , which leads to n ≥ u u -δ (5)
In practice, there are mainly two cases possible:

• u×u bijective S-boxes, with δ as small as possible to best resist differential attacks and essentially u = 4 or u = 8. We know that δ ≥ 1. However, at present time, the only known example of S-boxes that reach this bound are for odd u or u = 6. Still, there are countless examples of S-boxes with δ = 2. The case u = 4 is very well studied, see e.g. the classifications of 4 × 4 S-boxes of Leander and Poschamann [START_REF] Leander | On the Classification of 4 Bit S-Boxes[END_REF] or of Saarinen [START_REF] Saarinen | Cryptographic Analysis of All 4 x 4 -Bit S-Boxes[END_REF] and it is now very easy to find a 4 × 4 S-box with minimum δ. The other standard case is u = 8. If there is no known global classifications of 8 × 8 S-boxes at present day, there are many examples in the literature, the most famous being the S-box of the AES. In both cases (u = 4 and u = 8), n = 2 is obtained. • 6×4 S-boxes as in DES. In general, a 6×4 S-box verifies δ ≥ 6 -4 = 2 and hence n = 2. However, because of their structure the DES S-boxes do not reach this bound but verify δ = 4 instead, hence n = 3 for DES. Anyway, since the value of δ is also required to assess the security of the cipher against differential cryptanalysis, it is, by design of the cipher, always possible to compute it in reasonnable time. Now what happens when there is incertitude on the ∆ y value is studied. We have the following equation:

S(x ⊕ k) ⊕ S(x ⊕ e ⊕ k) = ∆ y1 or • • • or ∆ y J , (6)
where J denotes the number of possible ∆ y . Since for a fixed input difference a, the sets S a,bt are pair-wise disjoints, we have :

Pr z ∈ J t=1 S a,bt = J t=1 Pr [z ∈ S a,bt ]
As in the case where there is only one ∆ y , k is solution of equation ( 6) if and only if

k ∈ x ⊕ J t=1 S e,∆y t
Notice that the union of subkey hypotheses t (x ⊕ S e,∆yt ) is interesting only when it is different from the set of all the hypotheses.

If one does this experiment n times, we have:

Pr k ∈ n-1 s=0 x s ⊕ J t=1 S es,∆y s,t = n-1 s=1 Pr k ∈ x s ⊕ J t=1
S es,∆y s,t

Therefore the minimum number of faults n J , where J is the incertitude (the number of possible ∆ y ), to retrieve the key is the smallest n such that the term on equation [START_REF] Jeong | Differential Fault Analysis on Block Cipher Piccolo[END_REF] becomes smaller than 2 -u , i.e. small enough to eliminate all the keys but the right one.

IV. ALGORITHM AND RESULTS

A. Algorithm to find vulnerable locations for a fault injection

In this section, a method is given in order to find good locations to inject a single-bit fault into a generalized Feistel network. Algorithm 1 describes this procedure. Its inputs are the Feistel function F and the GFN matrix M. For every block

B e r with e ∈ [[0, b -1]] and r ∈ [[r -(d + 1), r -1]]
where the fault can be injected, the algorithm returns the number of subkey blocks and, for each of them, the number of pieces of subkey and the number of differentials to take into account.

Using the Feistel network, for each block candidate for a fault injection, the number n λ of subkey blocks attacked is first computed. Then for each of the subkey blocks K λ r , the algorithm computes the number of pieces attacked (worst case and best case). Then if V F (j) = -∞, i.e. if the block xored at the Feistel function output is faulted, giving a first estimation of the number of possibles differentials. The algorithm reuses then the Feistel network to find the least influenced ciphertext block (but influenced still) by the fault, to further reduce the number of possible differentials. If that number does not allow to eliminate some subkey candidates, we consider that the subkey piece K λ r cannot be attacked using this fault. Finally, we compute the number n J of faults required to retrieve the remaining candidates from the average number of ∆ y involved.

This algorithm gives a first insight on where to inject the fault. Yet its genericity does not rule out that a better fault attack exists on a particular GFN. Moreover, it is quite possible to imagine an attack that injects faults in more than one block or with a different fault model.

B. Results on some Feistel ciphers

We have tested our approach on several Feistel ciphers. The exposed results focus essentially on three parameters. First, the number n λ of subkey blocks attacked and for each of them the number n l of pieces attacked, which together are a measure of the fault efficiency: the higher these numbers, the better the fault. And second the number of differences ∆ B j r-1 on the output block (for each of the n λ blocks attacked) which is a measure of the fault inefficiency: the more ∆ B j r-1 there are, the more faults will be required to find the subkey, the worst case being that the subkey cannot be found. A good attacker will then have to find where to inject the fault to balance those parameters that tend to grow together but have opposite results.

In this section, the matrices have coefficients in the dioid D, assuming matrix vector product. The value -∞ are replaced by '.' (dot) for readability.

1) DES: Our algorithm is tested on DES since it was the most well known 2-block Feistel cipher.

In the following, details of our analysis are given:

• Number of rounds: r = 16.

Thus even when the fault does not go through a Feistel function, up to 2 S-boxes can still be attacked. If the fault goes through exactly one Feistel function, there are between 2 and 8 faulted output bits. After passing through the Feistel function, 2 to 8 bits are faulted at the output, leading to 32 * (256-1-8) possible differentials. Conversely, assuming a single-bit input differential, observing the output of the Feistel function is characteristic enough to reduce the number possible input differentials from 32 to only 2.

In the case of the DES, the fault injection on L 15 is eliminated since it does not allow to attack the key. Results of our analysis are summed up in Table I where n l is the number of pieces of subkey blocks attacked and ∆ is the number of guesses on ∆ L15 . The number n λ of subkey blocks attacked is of course n λ = 1 since there is only one Feistel function per round.

B V F W F n l ∆ R 15 or L 14 (., 0) (1, 0) 1 ≤ n l ≤ 2 1 R 14 or L 13 (0, 1) (2, 1) 2 ≤ n l ≤ 8 2 R 13 (1, 2) (3, 2) 2 ≤ n l ≤ 8 32 * 247
Our method shows that it is more judicious to inject faults in R 14 for the DES. Yet Biham and Shamir have shown in [START_REF] Biham | Differential Fault Analysis of Secret Key Cryptosystems[END_REF] that, for the DES, a fault in R 13 leads to a more efficient attack. To find this result, they have studied the round function in depth, while our approach is generic.

Moreover, the study by Rivain in [START_REF] Rivain | Differential Fault Analysis on DES Middle Rounds[END_REF] showed that injecting faults before, in rounds 9, 10, 11, 12 (beyond the full diffusion delay), requires a bigger number of faults to find the key.

2) MIBS: This cipher is chosen because, as DES, it is a classical Feistel network and thus it is mainly the Feistel function that matters for our study. Moreover, contrarily to DES, the S-boxes are bijective, which is more common among block ciphers.

In the following, details of our analysis are given:

• Number of rounds: r = 32.

• Number of blocks: b = 2.

• Full diffusion delay: d = 2,

• The GFN M: M = 1 0 0 . . The diffusion function following the S-boxes strictly preserves the partition into 4-bit blocks defined by the S-boxes. Hence the study of MIBS Feistel function is as follows:

• Each input bit faults exactly one S-box. S-boxes are bijective, so modifying one bit in input means modifying at least one bit in output, so 1 to 4 output bits can be faulted. • If the fault has been through exactly one Feistel function, there are 5 or 6 pieces of subkey blocks that can be simultaneously attacked. • If one observe a faulted output (assuming single-bit faulted input), one can deduce which S-box is faulted and thus the input fault is reduces to 4 possibilities. • If the fault has been through more than one Feistel function, then all the S-boxes are faulted. Results of our analysis are summed up in Table II where n l is the number of pieces of subkey blocks attacked and ∆ is the number of guesses on ∆ R31 . As for DES, the number n λ of subkey blocks attacked is n λ = 1.

TABLE II RESULTS OF OUR ANALYSIS ON MIBS Blocks

B V F W F n l ∆ L 31 or R 30 (0, .) (0, 1) 1 1 L 30 or R 29 (1, 0) (1, 2) 5 ≤ n l ≤ 6 4 L 29 (2, 1) (2, 3) 8 112
Results indicate that it seems best to attack on block L 30 . This way though, not all the subkeys pieces can be attacked at the same time. On the contrary, injecting on L 29 allows to attack the whole subkey, but the number of differences on ∆ R31 becomes greater. A more in-depth study of the L 29 case may be possible, but may not be trivial to perform.

3) TWINE: TWINE was chosen because it is the exact opposite of DES and MIBS: the Feistel function is very simple, it is a single S-box. The whole study will be on the Feistel networks itself.

In the following, details of our analysis are given: Since the Feistel function is essentially a single S-box (i.e. n l = 1), if the fault on block B j r-1 , that is xored to the ouput of the targetted function, has been through exactly one Feistel function in the previous rounds, then this fault can take ∆ = 14 of the 2 4 -1 = 15 possible values, restricted to ∆ = 7 if the input fault is known. Hence, there are three favorable cases:

• B j r-1 is not faulted: ∆ = 1. • B j
r-1 contains exactly the value of the fault: ∆ = 1 if the value of the fault is known, ∆ = 4 else. • B j r-1 has been through exactly one function: ∆ = 7 if the value of the fault is known, ∆ = 14 else. The best place for a fault injection on TWINE will then be the register where the fault reaches the most Feistel functions with the condition that their respective outputs must fall in the above three cases. It appears that the best case achievable is to inject a fault at round 31, on any B i 31 at the input of a function. This way it is possible to attack simultaneously n λ = 5 functions on the last round of TWINE, four of them with a non faulted B j r-1 , i.e. with ∆ = 1, and a fifth one with ∆ = 7. Injecting the fault earlier in the algorithm (round 30 and before) only allows to attack at most 4 functions, which is not as interesting as the previous case. Finally, if the fault is injected at round 32 or after, then only up to 3 functions can be attacked.

4) CLEFIA: CLEFIA was chosen as an intermediate case between DES and TWINE: it has more than two blocks and the Feistel function is more than a single S-box. Thus, both the Feistel function and the Feistel network have to be studied and then mixed together.

In the following, details of our analysis are given:

• The study of the CLEFIA Feistel function is as follows:

• Each input bit faults exactly one S-box. S-boxes are bijective, so modifying one bit in input means modifying at least one bit in output, so 1 to 4 output bits can be faulted. • Because of the diffusion properties of matrices M C 0 and M C 1 of CLEFIA, if the fault has been through one Feistel function, then all the S-boxes are faulted. • If one observe a faulted output (assuming single-bit faulted input), one can deduce which S-box is faulted and thus the input fault is reduced to 4 possibilities. • By symmetry in the network, one can restrict the analysis to the case where the fault is injected on block B 0 r . Results of our analysis are summed up in Table III, assuming a 128-bit key, i.e. there are 18 rounds. For each possible location of the fault, the number n λ ∈ [ [START_REF] Beaulieu | The SIMON and SPECK Families of Lightweight Block Ciphers[END_REF][START_REF] Berger | Extended Generalized Feistel Networks using Matrix Representation[END_REF]] of subkey blocks attacked, as well as for each of the two functions, the number n l ∈ [[0, 4]] 2 of pieces of subkey attacked and the corresponding number ∆ of differences on the output block are given. Results indicate that injecting on B 0 17 is not optimal since if the injection is made one round earlier on B 0 16 , it is possible for the same cost, to fault the whole function F 0 instead of a single S-box. The attack by Chen et al. [START_REF] Chen | Differential fault analysis on clefia[END_REF] with 18 faults (for CLEFIA-128) falls in that case.

If one wishes to fault both functions, then the second one will have a non trivial ∆. The case B 0 15 is still manageable, but only one S-box of the second function is attacked. Fault injections on B 15 were used by Takahashi and Fukunaga [START_REF] Takahashi | Improved differential fault analysis on clefia[END_REF] to successfully attack CLEFIA-128 requiring only 2 faults. Even if they make an in-depth study of CLEFIA Feistel functions to reduce the number of faults to its minimum, this result corroborates our analysis: it is best to inject the fault on B 15 .

It might be possible to manage the case where the fault is done on B 0 14 . This would allow to attack all 8 S-boxes at once but it would require very deep study of the Feistel functions to tackle the problem of ∆ on the second function.

V. CONCLUSION

In this paper, it has been shown that some blocks are more vulnerable than others in generalized Feistel networks. Fewer faults are needed to find the key in the most vulnerable blocks, which make them targets of choice for fault attacks. A method has been proposed to identify these blocks allowing designers to deal with them accordingly. The genericness of the proposed method implies that the vulnerability evaluation is not optimal but a method to assess the vulnerabilities automatically has been proposed. Further work will include finding a refined version of this study that can go deeper into the analysis of the Feistel functions while maintaining genericity.

Feistel ciphers 1 )

 1 DES: The Data Encryption Standard (DES) [18] was developed by IBM and standardized by the NIST in 1977.

Fig. 4 .

 4 Fig. 4. Feistel function F of DES.

Fig. 5 .

 5 Fig. 5. Structure of the MIBS Feistel function.

15 r+1Fig. 6 .

 156 Fig. 6. One round of TWINE.

Fig. 8 .

 8 Fig. 8. The Feistel function F 0 of CLEFIA.

Fig. 9 .

 9 Fig. 9. Attack path of a fault injection on a Feistel cipher.

Fig. 10 .

 10 Fig. 10. Matrix of the GFN given on Fig. 2.

• 1 = ∆ R31 . •

 1. Number of subkey blocks: Λ = 1. One has to remark that B i r-1 = L 31 and ∆ B j r-Number of pieces in subkey blocks (number of S-boxes): L = 8. • Number of faults required to retrieve a piece of subkey using equations of the form (3): n = 2.

•

  Number of rounds r = 18, 22 or 26 depending on the key size (resp. 128, 192 or 256 bits). • Number of blocks b = 4. • Full diffusion delay d = 4. • The GFN M: Number of subkey blocks: Λ = 2. • Number of pieces in subkey blocks (number of S-boxes by Feistel function): L = 4. • Number of faults required to retrieve a piece of subkey using equations of the form (3): n = 2.

Acknowledgements. The authors would like to thank Ronan Lashermes and Thierry P. Berger for their valuable contributions to the development and understanding of the issues discussed in the paper.

This work was partially funded by the French DGCIS (Direction Générale de la Compétitivité de l'Industrie et des Services) through the CALISSON 2 project; and partially supported by the French National Agency of Research: ANR-11-INS-011.

Compute the associated number of pieces attacked n l Compute the number J of output differences possible and the number of faults n J required to retrieve a subkey piece end for for

Deduce the number of blocks that can be attacked n λ for All the subkey blocks K λ r attacked do Identify the input block B i r-1 and the block xored at the output of the Feistel function • The GFN M: M = . 0 0 1 .

• Number of subkey blocks: Λ = 1. One has to remark that

• Number of pieces in subkey blocks (number of S-boxes): L = 8. • Number of faults required to retrieve a piece of subkey using equations of the form (3): n = 3.

Then the Feistel function is studied :

• Because of Expansion E, a single bit input difference can impact 2 S-boxes. • A property of DES S-boxes is that a single bit difference at the input of an S-box ensures at least a two-bit difference at the output. • the bit-permutation P is made such that the 4 output bits of an S-box are sent to 4 different S-boxes at the input of the next round.


Dr Chang-Soo R Kim

R Aydé

A.-L Sablong
email: raphael.sablong@univ-lyon1.fr

O Perrier

Beuf

G Gaborit

P Jarrige
email: pierre.jarrige@kapteos.com

L Duvillaret
email: lionel.duvillaret@kapteos.com

Potentialities of an Electro-Optic Crystal Fed by Nuclear Magnetic Resonant Coil for Remote and Low-Invasive Magnetic Field Characterization

Reina Aydé, Gwenaël Gaborit, Pierre Jarrige, Lionel Duvillaret, Raphaël Sablong, Anne-Laure Perrier, and Olivier Beuf

Abstract-In this paper, we demonstrate the use of a LiTaO 3 crystal associated with a typical nuclear magnetic resonant loop coil to perform an optically remote radio frequency magnetic-field characterization. The whole transduction scheme is theoretically and experimentally studied. The measurement dynamics reaches 60 dB. The minimum detectable magnetic field is lower than 1 nT, which corresponds to an induced inner crystal electric field as low as 30 mV/m. To evaluate the spatial potentialities of the sensor, a 1-D mapping of the field along an asymmetric butterfly-shaped loop coil is performed. The result is in good agreement with finite-difference time-domain simulations and demonstrates the vectorial behavior of the sensor device.

Index Terms-Electro-optics, magnetic field sensor, magnetic resonance imaging, nuclear magnetic resonance.

I. INTRODUCTION

M AGNETIC resonance imaging (MRI) is based on a spatial characterization of the nuclear magnetization. Recording the variation of the nuclear magnetization is often realized with a MRI external surface coil. However, spatial resolution and image quality achievable with external surface coils are limited in the case of thin internal organs, e.g. bowel. On the other hand, the use of an endoluminal MRI coil located close to the area of interest could provide both high contrast to noise ratio (CNR) and high spatial resolution images leading, in turn, to a good evaluation of biological diseases [START_REF] Armenean | Optimization of endoluminal loop radiofrequency coils for gastrointestinal wall MR imaging[END_REF], [START_REF] Beuf | In vivo colon wall imaging using endoluminal coils: Feasibility study on rabbits[END_REF].

Since metallic coaxial cables are mainly used with endoluminal coils, also located within the MRI transmit RF coil, heating phenomena may occur inside biological media [START_REF] Armenean | RF-induced temperature elevation along metallic wires in clinical magnetic resonance imaging: Influence of diameter and length[END_REF]. Optical fibre link can be used to overcome heating problem, thus ensuring patients safety [START_REF] Yuan | A 4-channel coil array interconnection by analog direct modulation optical link for 1.5-T MRI[END_REF], [START_REF] Memis | Miniaturized fiberoptic transmission system for MRI signals[END_REF]. Direct modulation of a laser diode remains invasive because its biasing requires a DC power supply [START_REF] Yuan | A 4-channel coil array interconnection by analog direct modulation optical link for 1.5-T MRI[END_REF]. Furthermore, the achieved realizations including active devices may also be disturbed by the static magnetic field [START_REF] Yuan | A 4-channel coil array interconnection by analog direct modulation optical link for 1.5-T MRI[END_REF], [START_REF] Fandrey | Miniature optical signal transmission system for an active intravascular device[END_REF], [START_REF] Du | Comparison of FP, VCEL and DFB laser diode in optical transmission for MR RF coil array[END_REF].

Magnetic field characterization involving passive optical system has already been performed [START_REF] Suzuki | Optical magnetic field probe with a loop antenna element doubly loaded with electrooptic crystals[END_REF]- [START_REF] Mhapatra | Electro-optic electromagnetic field sensor system with optical bias adjustment[END_REF]. This type of magnetic probe is based on an Electro-Optic (EO) material associated to a magnetic loop inducing an electric (E)-field proportional to a magnetic field component. This experimental configuration allows to minimize the invasiveness: the probe can be pigtailed, does not need power supply and includes metallic element much shorter than the wavelength of the field to be measured. This kind of sensor has been studied by Suzuki et al. [START_REF] Suzuki | Optical magnetic field probe with a loop antenna element doubly loaded with electrooptic crystals[END_REF] and exploits a magnetic loop double loaded with LiNbO 3 crystal, acting as the capacitance of the resonant equivalent circuit. The realized set-up is dedicated to measure the magnetic field associated to a microwave electromagnetic field propagated along a stripline. We here propose a similar configuration, dedicated to MRI based RF magnetic field measurement. Nevertheless, as our field of interest concerns proton MRI, the resonant frequency is precisely defined (e.g. 127 MHz at 3 T). Moreover, this resonant frequency is lower in MRI system compared to the microwave frequency probed by Suzuki et al., thus inducing a lower electromotive force ε = dφ dt applied to the crystal. φ represents the magnetic flux. The loop coil is designed to be intrinsically resonant at 127 MHz and the crystal capacitance is kept very low to keep the same frequency response. In the first section of the paper, the crystal orientation and the measurement set-up are studied and optimized in order to maximize the sensitivity, the linearity and the measurement dynamics. Then, the experimental setup used for the characterization of the EO magnetic probe is depicted. The response intensity of the probe is measured versus applied magnetic field amplitude. Finally, last section is dedicated to 1D magnetic field mapping of a butterfly-shaped loop RF coil.

II. THEORETICAL BACKGROUND AND MEASUREMENT SET-UP

A. Endoluminal Coil Design

Based on the previous study by Armenean et al. [START_REF] Armenean | Optimization of endoluminal loop radiofrequency coils for gastrointestinal wall MR imaging[END_REF], we choose the parameter of the coil as those providing the best signal-to-noise ratio (SNR) for MRI together with suitable dimensions for endoluminal investigations. Hence, the coil has a rectangular shape with an external width set to 5.1 mm, respecting medical constraints, and a length set to 47 mm providing good image penetration depth and exploration length. A 35 μm-thick copper coil is engraved on a insulating substrate FR4 (ε r = 4 and thickness T = 0.8 mm). A conductor width of 0.8 mm is defined to minimize ohmic losses then providing the best SNR images. The coil is tuned to a frequency of 127 MHz and is matched at 50 , using capacitors of 282 pF and 39 pF on distal and proximal part of the coil, respectively. A schematic and a picture of the endoluminal coil are presented in figure 1.

The characterization of the coil was performed thanks to a vector network analyzer and the magnitude of S11 reflection coefficient is given in figure 2. The resonant frequency is 126.9 MHz and its bandwidth is 3.5 MHz at -3 dB.

The extracted parameters of the RLC series equivalent circuit are R = 0.85 , L = 46 ×10 -9 H and C = 34.2 ×10 -12 F. These values are very close to the expected ones, obtained thanks to method of moments simulations.

B. Electro-Optic Transducer

The Pockels effect appears as an E-field induced modification of the eigen refractive indices in non-centrosymetric crystals. The developed transducer is based on an x-cut LiTaO 3 crystal. An exhaustive analysis of the properties of EO crystals is described by Duvillaret et al. [START_REF] Duvillaret | Electro-optic sensors for electric field measurements. II. Choice of the crystals and complete optimization of their orientation[END_REF]. The EO tensor, which gives the vectorial dependence of the index ellipsoid with the E-field is the following one:

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 -r 22 r 13 0 r 22 r 13 0 0 r 33 0 r 51 0 r 51 0 0 -r 22 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
. This latter one leads to the refractive index ellipsoid equation (eq. 1):

x 2 1 n x 2 -E y r 22 + E z r 13 + y 2 1 n y 2 + E y r 22 + E z r 13 +z 2 1 n z 2 + E z r 33 + 2E y r 51 yz + 2E x r 51 xz -2E x r 22 x y. ( 1
)
A figure of merit can be calculated for this crystal [START_REF] Duvillaret | Electro-optic sensors for electric field measurements. II. Choice of the crystals and complete optimization of their orientation[END_REF]. As the crystal acts as a polarization state modulator, the relevant property is the E-field induced birefringence δn

= n + ( -→ E ) -n -( -→ E ).
The EO figure of merit is depicted by the sensitivity vector --→ K , linked to the relative variation of refractive indices induced by the applied E field. This sensitivity vector is defined by:

--→ K = -→ ∇ δn( -→ E )|-→ E = 0 . ( 2
)
Its modulus gives the sensitivity of the EO crystal to the applied E field while its direction leads to the unique E-field component that will be probed. These properties are summarized by the following expression:

δn( E) ∼ = δn( E = 0) + --→ K • E. (3)
A three-dimensional representation of --→ K is given in figure 3.

We here consider an x-cut LiTaO 3 , leading to a high sensitivity together with a transverse E-field probing (see Figure 3: = 90°). For an applied E field along z axis, we get:

| --→ K x-cut,E z | = 1 2 n 3 z r 33 -n 3 y r 13 . ( 4
)
The theoretical value of | --→ K x-cut,E z | at the wavelength of 1550 nm is 120 pm/V assuming n z = 2.116, n y = 2.117, r 13 = 7.5 pm/V, and r 33 = 33 pm/V [START_REF] Duvillaret | Electro-optic sensors for electric field measurements. II. Choice of the crystals and complete optimization of their orientation[END_REF].

The E-field induced dephasing between the two allowed polarization state inside the EO crystal writes:

ϕ = 2π λ opt | --→ K x-cut,E z |.E z L x . ( 5
)
λ opt being the wavelength of the probing optical beam and L x , the length of the crystal.

LiTaO 3 :

--→ ΔK

C. Measurement Set-Up

The magnetic field source is constituted of a single wide band loop with a diameter ∅ = 10 cm to provide an homogeneous field onto the endoluminal coil surface. As the dimension of the transmit coil is much shorter than the wavelength of the signal at 127 MHz, B x is thus directly linked to the input power P in,d Bm of the emission coil, via Biot and Savart equation:

B x = μ 0 I c 2 a 2 a 2 + x 2 3/2 (6)
with a being the radius of the emitting loop and where I c = This emission coil is fed by a function generator providing a CW signal at the frequency of 127 MHz. The electrical power P in varies from 14 dBm down to -60 dBm thanks to additional RF attenuators. The induced magnetic field at the center of the loop ranges from 700 nT downto 140 pT.

The endoluminal coil is positioned in front of the emission coil and probe the field B x (see Fig. 4). It induces a rms electromotive force ε onto the z-axis of the crystal thanks to copper deposited parallel electrodes, thus applying an electric field E z = ε h z , h z being the separation distance between electrodes.

The electric field induces a modulation of the laser polarization state. This latter one can be written as a function of the input power P in,d Bm feeding the transmit coil and of the coupling factor κ ant between the transmit loop and the h z [START_REF] Du | Comparison of FP, VCEL and DFB laser diode in optical transmission for MR RF coil array[END_REF] with R load = 50 being the resistive load of the endoluminal coil.

For x = 0 (at the center of the loop), equations 7 and 6 lead to the link between E z and B x :

E z (B x ) = 2a10 κ ant 20 √ R em R load h z μ 0 B x . ( 8
)
The optical source is a low noise DFB laser (λ = 1.55 μm and Relative Intensity Noise R I N = -160 dBc.Hz -1 ) servo-controlled in temperature and in optical power ( Popt = 50 mW). A half waveplate orients the impinging linear polarization with an angle of 45°relatively to the crystal axes to maximize the conversion between the field to be measured and the polarization state modulation of the laser probe beam. The modulated polarization state is then analysed thanks to a quarter waveplate, compensating the intrinsic static birefringence of the crystal. A polarizer converts the relative dephasing ϕ into a modulation of optical power Popt = 1 2 ϕ Popt . This optical modulation is finally converted into an electrical signal using an amplified (G = 37.3 dB) pigtailed high-speed photodiode (sensitivity η = 0.85 A/W, loading resistance R ph = 50 ). The mean optical power received by the photodiode after polarization treatment and fibre coupling is Popt = 14.4 mW. This latter signal carries the modulation and is then recorded using a spectrum analyzer.

III. EXPERIMENTAL CHARACTERIZATION

The analysis of the signals at different transduction stages, from the emission to the optoelectronic output is firstly investigated. Figure 5 illustrates the spectrum of the input signal feeding the transmit coil, the endoluminal coil received signal and the EO converted signal. The resolution bandwidth (RBW) of the measurement is 30 Hz. It can be noticed on figure 5(a), that the noise floor increases by 40 dB after the first stage of transduction. This is due to internal attenuation of the spectrum analyser that automatically sets for high level signals.

Finally, the conversion factors are κ ant = -22 ± 0.3 dB and κ E O = -34.7 ± 0.3 dB for the detection coil and for the EO conversion respectively.

The sensitivity and the linearity of this conversion has been also studied over the 74 dB range of available power. The result is shown on figure 6. The dynamics exceeds 60 dB and is down limited by the noise floor (-112 dBm) of the spectrum analyser, which is much greater than the R I N of the laser and the shot-noise on the photodiode (-140 dBm and -151 dBm in a 1 Hz analysis bandwidth respectively). The measurements are also compared with two curves. The linear fitting curve (gray dashed line on figure 6) is written:

P out,d Bm = A P in,d Bm + κ conv . ( 9
)
A is the linearity factor and is equal to 0.997. κ conv illustrates the conversion factor of the whole transduction scheme. κ conv = -54.1 dB is close to κ ant + κ E O , which is -56.7 ± 0.6 dB. The second equation does not involve any fitting parameter and describes the expected output signal as a function of set-up parameters.

Introducing α P = 10 G+κ ant 10 R ph R load ( π| --→ K x-cut,Ez |L x η Popt λ opt h z ) 2 =
1.2 × 10 -6 , the expression between P in,d Bm and the output signal writes: The calculated value of α P leads to an expected output signal that is fully consistent with the measured one (see inset of figure 6). Using the link between P in,d Bm and B x , given in equation 6, P out,d Bm can also be written as a function of the magnetic field:

P
+ α B B 2 x ( 11
)
with α B = 10 3 R em 2a μ 0 2 α P = 5.9 × 10 7 mW.T -2 . Moreover, the standard deviation between the measurement and the fit given by equation 9 is equal to 2.1 dB for the whole dynamic range and takes a very low value of 0.20 dB for an input signal ranging from -26 dBm to 14 dBm (40 dB dynamic range). Finally, the minimum detectable field, which is defined by the value of the field equalizing the equation 9 and the noise floor value, is B min = 176 pT. This latter value reaches 7 pT considering only the laser RIN contribution to the noise, in a 1 Hz analysis bandwidth.

In order to perform an accurate and reliable measurement of the magnetic field, another relevant parameter is the stability of the experiment set-up. This stability of the output signal is mainly determined by the efficiency of the servo-control system that regulates the optical power and the laser temperature. The output signal is measured along 15 minutes for P in,d Bm = 14 dBm. The result is shown on figure 7(a).

The standard deviation σ E O with respect to the mean signal is deduced for the distribution represented in figure 7(b). We get σ E O < 0.2 dB.

IV. MAGNETIC FIELD MAPPING

In a more global use of the transduction scheme for magnetic field measurement, we here propose to evaluate the ability of the set-up to perform a near field mapping. As a characterization example of inhomogeneous magnetic field, in orientation and magnitude, we choose an asymmetric butterfly-shaped loop coil to be the tested device. Moreover, the typical dimensions of the butterfly-shaped loop coil are similar to the endoluminal coil and the field to be measured will be integrated on its whole radiation pattern. The device under test geometry exhibits a vanishing B x magnetic field at the crossing between the two loops and allows to exploit the available dynamics of the set-up. For this experiment, the referential is the same as indicated on figure 4: the x axes of the two coils are aligned and the endoluminal coil is translated along the y axis (see figure 8(a). The algebraic value of B x static magnetic field is firstly 2D mapped thanks to FDTD simulation (quickfield ® ) for z = 0. From this simulation, we define an interesting 1D path (thick arrows on figure 8(a) and 8(b)) for magnetic field mapping in order to make a comparison between measurement and simulation. Both magnetic field distributions are plotted on figure 8(c).

The experimental result exhibits a dynamics of more than 40 dB and is in good agreement with the simulation, especially for the spatial evolution of B x . Nevertheless, the minimum measured value of -82 dB, corresponding the crossing location, is much higher than the noise floor. This is explained by the residual magnetic field surrounding the crossing point and still probed by the coil (mainly B x for z = 0). This result also demonstrates the immunity of the transducer relatively to the transverse component (B y ) of the magnetic field and illustrate the vectorial behaviour of the sensor.

V. CONCLUSION

A LiTaO 3 crystal is used together with an endoluminal coil to demonstrate the ability of such transducer to perform a remote and reliable characterization of RF magnetic field at 3T NMR frequency. The whole transducing scheme has been theoretically studied. The in-lab implementation demonstrates a very good agreement with the expected results. The experimental minimum detectable field is weaker than 1 nT. The accuracy of the measured signal is better than 0.5 dB along a 50 dB dynamic range with a rather good stability (deviation lower than 0.2 dB for a 15 minutes temporal window). The previous mentioned performances allow to map one component of the magnetic field and a vectorial selectivity greater than 40 dB is demonstrated. Development of such sensor scheme constitutes a first step towards a fully pigtailed probe. This technique would also benefit from the potential miniaturization of the crystal size (h z 10 μm for an EO waveguide), in order to increase the sensitivity (by a factor 400). Furthermore, optical amplifier (+24 dB typically) can be used to increase the optical modulation after polarization treatment to reach lower minimum detectable field, required for actual MRI system, while keeping advantages of passive and remote EO configuration.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Endoluminal RF coil. (a) Schematic. (b) Picture.

Fig. 3 .

 3 Fig. 3. LiTaO 3 figure of merit for --→ K . | --→ K | is determined by the distance to the origin. Color: angle corresponding to the sensitivity axis orientation between --→ K and the optical wave-vector -→ k . Scale at the bottom left. Projections onto the reference planes at the right side. The unit of | -→ K | is pm • V -1 .

10 P

 10 in,d Bm -30 10 R em . R em is the resistive load of the emission loop and takes a value of 8 at 127 MHz.

Fig. 4 .

 4 Fig. 4. Experimental bench. DFB LD: servo-controlled laser diode (λ = 1.55 μm). PM fiber: polarization maintaining fiber. LiTaO 3 : Stoechiometric lithium tantalate crystal. λ/2 and λ/4: half and quarter waveplates.

  endoluminal coil:E z (P in,d Bm ) = R load 10 P in,d Bm +κ ant -3010

Fig. 5 .

 5 Fig. 5. Response at different stage of transduction for three values of power feeding the transmit coil. (a) P in = 14 dBm. (b) P in = -20 dBm. Black dashed line: impinsging signal P in,dBm . Gray line: measured with the endoluminal coil. Black line: measured with the EO crystal coupled to the endoluminal coil.

10 Noise d Bm 10 + α P 10 P in,d Bm 10 .Fig. 6 .

 101010106 Fig.6. EO response as a function of the input power feeding the transmit coil (RBW = 100 Hz). Top axis: input power is linked to the induced magnetic field. Black dots: measurements. Gray dashed line: linear fitting curve. Dashed dotted line: theoretical curve, including noise contributions, obtained without any fitting parameter. Horizontal gray dotted line: noise floor. Inset: zoom of the result.

P

  out,d Bm = 10 log 10 10 Noise d Bm 10

  ween measurement s and mean value dB Dist ribut ion

Fig. 7 .

 7 Fig. 7. (a) Black curve: stability of the EO response. Dashed gray line: average signal is indicated. Gray dotted line: noise floor (-105 dBm). (b) Histogram: distribution of the standard deviation. Black line: associated Gaussian fit.

Fig. 8 .

 8 Fig. 8. (a) Schematic of the asymmetric butterfly-shaped loop coil. Thin arrows: current flow and the thick arrow is the mapping path. (b) FDTD simulation of the butterfly-shaped loop coil. Color scale from blue: lowest negative B x . To red: highest positive B x . Down arrow: measurement location. (c) Spatial distribution of the magnetic field component B x along the y-axis of the asymmetric loop coil. Black dots: measured with EO transducer. Dashed gray line: simulated with FDTD.

part by the Region Rhone Alpes, the DGA (French Military Programs Management and Procurement Agency -Project RAPID SNIFER) and DGCIS.