Amina Guermouche
email: amina.guermouche@telecom-sudparis.eu

Anne-Cécile Orgerie
email: anne-cecile.orgerie@irisa.fr

Experimental analysis of vectorized instructions impact on energy and power consumption under thermal design power constraints

Keywords: Power consumption, energy efficiency, SIMD instructions, TDP, memory

Vectorized instructions were introduced to improve the performance of applications. However, they come with an increase in the power consumption cost. As a consequence, processors are designed to limit the frequency of the processors when such instructions are used in order to maintain the thermal design power.

In this paper, we study and compare the impact of thermal design power and SIMD instructions on performance, power and energy consumption of processors and memory. The study is performed on three different architectures providing different characteristics and four applications with different profiles (including one application with different phases, each phase having a different profile).

The study shows that, because of processor frequency, performance and power consumption are strongly related under thermal design power. It also shows that AVX512 has unexpected behavior regarding processor power consumption, while DRAM power consumption is impacted by SIMD instructions because of the generated memory throughput.

I. INTRODUCTION

The race for computing performance has led the computer vendors to introduce many features and new techniques to run computations as fast as possible. Such hardware improvements allow the supercomputers in the TOP500 to gain 6 orders of magnitude in terms of performance in the last 25 years [START_REF]TOP500 Supercomputer Site[END_REF].

Turboboost allows the processor to run at higher frequencies than the base one in order to increase performance. Simple Instruction Multiple Data (SIMD) model is another of these hardware techniques. In the SIMD model, the same operation is executed simultaneously on different elements of a vector or different data points. For instance, several iterations of the same loop for a vector/vector addition can be processed at once. Figure 1 shows the difference between an SIMD processor and a scalar processor.

This performance race comes with side effects in terms of power consumption and heat dissipation. Indeed, the power efficiency of supercomputers in the TOP500 has gained 3 orders of magnitude in the last 15 years -this metric being collected since fewer years than performance [START_REF]TOP500 Supercomputer Site[END_REF]. It means that the power consumption of supercomputers keeps increasing resulting in larger and larger heat dissipation at the processor b [START_REF]TOP500 Supercomputer Site[END_REF] a [START_REF]TOP500 Supercomputer Site[END_REF] b [START_REF] Cebrian | Improving Energy Efficiency through Parallelization and Vectorization on Intel Core i5 and i7 Processors[END_REF] a [START_REF] Cebrian | Improving Energy Efficiency through Parallelization and Vectorization on Intel Core i5 and i7 Processors[END_REF] b [START_REF]Intel ® xeon ® processor e5 v4 product family. thermal mechanical specification and design guide[END_REF] a [START_REF]Intel ® xeon ® processor e5 v4 product family. thermal mechanical specification and design guide[END_REF] b [START_REF]Intel ® xeon ® processor e5 v3 product family. specification update[END_REF] a [START_REF]Intel ® xeon ® processor e5 v3 product family. specification update[END_REF] b [START_REF] Lento | Optimizing performance with intel ® advanced vector extensions[END_REF] a [START_REF] Lento | Optimizing performance with intel ® advanced vector extensions[END_REF] + (b) SIMD processor: several simultaneous operations Fig. 1: Difference between SIMD and scalar operations. Example: addition of two vectors level, and consequently, because of thermal limits, to a growing fraction of dark silicon [START_REF] Cebrian | Improving Energy Efficiency through Parallelization and Vectorization on Intel Core i5 and i7 Processors[END_REF]. To prevent physical damages due to heat, processors are designed to respect a thermal design power (TDP). TDP is the average power that the processor dissipates when operating at the base frequency with all cores active [START_REF]Intel ® xeon ® processor e5 v4 product family. thermal mechanical specification and design guide[END_REF]. The base frequency is defined as the frequency when turboboost is disabled.

The performance of hardware improvements can be limited by the TDP. Indeed, when the processor detects SIMD instructions, additional voltage is applied to the core. With the additional voltage applied, the processor could run hotter, requiring the operating frequency to be reduced to maintain operations within the TDP limits. This is also the case for turboboost frequencies [START_REF]Intel ® xeon ® processor e5 v3 product family. specification update[END_REF].

The TDP limit enforcement impacts the execution of applications in a non trivial manner. As an example, we run HPL, a CPU-intensive HPC benchmark, with and without turboboost on server chifflet. This experiment will be described in details later in the paper. Here, we only provide the execution time in Table I as an example. It illustrates that, unexpectedly, turboboost does not significantly increase the performance (0.26% difference) on this server for this application, which is optimized for vectorized instructions. The goal of this paper is to experimentally study such behaviors and provide scenarios where unexpected results may be obtained. The idea consists in identifying the impact of thermal design power on SIMD instructions. In order to do so, we will examine the behavior of different application profiles, when using different SIMD instructions on different hardware platforms. Note that we will only focus on Intel architectures in the remainder of the paper. This study aims at analyzing and comparing the power and energy consumption of each type of instruction and presenting cases of abnormal and unexpected behaviors.

The paper is organized as follows. Section II describes the experimental testbed. The used applications are detailed in Section III and experimental results are analyzed in Section IV. Related work is presented in Section V. Section VI concludes this experimental analysis.

II. ARCHITECTURE AND EXPERIMENTAL TESTBED

This section briefly presents SIMD instructions characteristics. Then it describes the processors used for the experiments and the power measurement methodology.

A. SIMD instructions

As stated before, SIMD instructions allow the same operation to be executed simultaneously on different elements of a vector or different data points. The number of simultaneous operations depends on the registers' size provided by the processors. Intel implements floating point SIMD instructions since the late 90's with the introduction of Streaming SIMD Extensions (SSE). Registers of 128 bits (16 Bytes) were used to hold 2 double-precision floats or 4 single-precision floats. Advanced Vector Extensions (AVX) appeared in 2010. The registers' size was doubled (256 bits) for floating point operations. However, the 128-bits integer SIMD instructions were not expanded. Finally, since the Haswell processor (2013), AVX2 extensions were introduced. They expand most 128bits SIMD integer instructions to 256 bits. Moreover, AVX2 extension adds fused multiply-add instructions (FMA) [START_REF] Lento | Optimizing performance with intel ® advanced vector extensions[END_REF]. AVX-512 is a 512-bits extension of the 256 bits operations (for floating point and integer operations). They are implemented in the Intel Xeon Phi and Skylake CPUs (2015).

B. Target platform

For the experiments, we used three servers: two from the Grid'5000 [START_REF] Balouek | Adding virtualization capabilities to the Grid'5000 testbed[END_REF] platform (nova and chifflet) and a Skylake processor located in our lab (skylake), at Telecom SudParis. We chose these nodes because they do not provide the same characteristics: nova runs at the same frequency regardless of the SIMD instructions being used. chifflet and skylake frequencies are impacted by SIMD instructions. Besides, skylake provides AVX512 instruction set. The nodes are described below and their characteristics in terms of TDP and frequencies are summarized in II provides the idle, base frequency, and characteristics for each server. It also shows the turboboost frequency depending on the SIMD instruction. Note that the turboboost frequency is not the frequency of the processor during the whole execution, but rather the frequency when SIMD instructions are used. Thus, higher frequencies may be observed during the execution time. On nova, the turboboost frequency is independent from the SIMD instructions being used while skylake shows the same base and turboboost frequencies when using AVX512. Moreover, when reaching TDP, processor may run at lower frequencies on chifflet and skylake as we will observe in Section IV. skylake data were extracted from [START_REF]Intel ® xeon ® processor scalable family[END_REF] while nova and chifflet data can be found in [START_REF]Intel ® xeon ® processor e5 v4 product family[END_REF]. One should note that while TDP constitutes a physical hard limit, a processor can briefly goes beyond it due to thermal inertia.

C. Power measurements methodology

In our experiments, we measure the power consumption and the execution time. Execution time is provided by the applications themselves.

In our study, power measurements rely on LIKWID [9] (version 4.3.0 1 ). LIKWID is a set of command line tools that are able to read performance counters, pin threads, . . . . In order to measure the power consumed during the execution of an application, we use likwid-perfctr which reads the corresponding hardware counters. We measure both sockets and DRAM power consumption of the target platforms. Note that the words package, socket or processor will be used in the document. likwid-perfctr relies on Running Average Power Limit (RAPL) counters. RAPL was introduced by Intel in order to stay within the power limit. It uses Dynamic Voltage and Frequency Scaling (DVFS) to guarantee the desired power limit. RAPL interface describes Model Specific Registers (MSRs). These registers provide energy consumption information for components such as the processor and the DRAM. On processors like Intel Sandy Bridge, RAPL was based on a modeling approach. Since the Intel Haswell generation, the processors have fully integrated voltage regulators (FIVR) that provide actual measurements and allow for more accuracy [START_REF] Hackenberg | An Energy Efficiency Feature Survey of the Intel Haswell Processor[END_REF]. Note that the literature provides many studies on RAPL accuracy for both DRAM [START_REF] Hackenberg | An Energy Efficiency Feature Survey of the Intel Haswell Processor[END_REF], [START_REF] Desrochers | A Validation of DRAM RAPL Power Measurements[END_REF] and processor [START_REF] Khan | RAPL in Action: Experiences in Using RAPL for Power Measurements[END_REF].

In all our experiments, the measurements are performed every second, with no overhead for all the applications. The mean of the power consumption is computed over the whole execution. The results presented in Section IV represent the average over 5 runs. Regarding measurements errors, all configurations show a small standard deviation. The maximum package power difference observed was on skylake with SVD Bulge using SSE (4W over 121W). The variation in DRAM measurements are very low (< 1%).

III. APPLICATIONS

This section describes the applications and the configuration parameters that we used. It also provides a characterization of these applications.

A. Applications description and configuration

In order to study AVX impact on power and energy consumption, we target 4 different HPC applications, which have automatic vectorization (by setting a compilation flag or an environment variable). We used applications with different CPU behavior and/or available options: HPL [START_REF] Petitet | Hpl -a portable implementation of the high-performance linpack benchmark for distributed-memory computers[END_REF] and Plasma svd [START_REF] Haidar | An improved parallel singular value algorithm and its implementation for multicore hardware[END_REF] which use the Math Kernel Library (MKL), AFF3CT [START_REF] Cassagne | Beyond Gbps Turbo decoder on multi-core CPUs[END_REF] and SToRM [START_REF] Noé | Seed design framework for mapping solid reads[END_REF]. Note that AFF3CT and SToRM only use integer. The following paragraphs present a short description of the applications and their configurations.

• High Performance Linpack (HPL) [START_REF] Petitet | Hpl -a portable implementation of the high-performance linpack benchmark for distributed-memory computers[END_REF]: HPL is a software package that solves dense linear algebra systems. It is used as a reference benchmark to compute the performance of the supercomputers in the TOP500 [START_REF]TOP500 Supercomputer Site[END_REF]. All the parameters are presented in • PLASMA Singular Value Decomposition (SVD) [START_REF] Haidar | An improved parallel singular value algorithm and its implementation for multicore hardware[END_REF]: The SVD decomposition computes the singular values of a matrix. It is performed in three steps. Readers can refer to [START_REF] Haidar | An improved parallel singular value algorithm and its implementation for multicore hardware[END_REF] for more details on the algorithm. The first step is referred to as SVD Band in the remainder of the paper while the second is referred to as SVD Bulge. In our configuration, we used the SVD version implemented in PLASMA [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF]. PLASMA is a software package for solving problems in dense linear algebra using multicore processors and Xeon Phi coprocessors. We used PLASMA version 2.8.0 compiled with MKL library version l mkl 2018.3.22. Thus, setting the desired SIMD instruction is done the same way as HPL. Table IIIb details the parameters values used for the target nodes. We fixed the size such that the first two phases last long enough to have several power measurements. The first phase is CPU-intensive while the second phase is memory-intensive. Note that the third phase is very short (few seconds) for these configurations. For this reason, we will not present an analysis for this phase. • Seed-based Read Mapping tool for SOLiD or Illumina sequencing data (SToRM) [START_REF] Noé | Seed design framework for mapping solid reads[END_REF]: SToRM is a read mapping tool based on mapping data between reads and a reference genome. It runs several phases. We present the results for the search algorithm used in the application. As SToRM uses integers in the SIMD parts, only SSE, AVX2 and AVX512 results will be presented. Note that SToRM is a CPUintensive application. We will no further detail the software.

The user can refer to [START_REF] Noé | Seed design framework for mapping solid reads[END_REF] for more details. Section ?? is an appendix detailing how we generated the input reads for SToRM. Note that the SIMD instructions are handled within the code. Setting the desired SIMD instruction set is done using the right compilation flag. • A Fast Forward Error Correction Toolbox (AFF3CT) [START_REF] Cassagne | Beyond Gbps Turbo decoder on multi-core CPUs[END_REF]:

AFF3CT is a library used for forward error correction. Forward Error Correction (FEC) is used to control errors during data transmission in order to enable efficient communications over noisy channels. It is done through encoding (by the sender) the data frame and decoding (by the receiver). We will no further detail Forward Error Correction, but we will present how the authors used vectorized instructions in their decoding solution. The decoder takes a set of frames as input. The decoding of the frames is vectorized. In order to do so, the frames are first buffered and then the vectorized algorithm is applied on the frames. Thus, depending on the SIMD instruction being used, the number of loaded frames differs. For instance, when using AVX2 instructions, twice the number of frames are loaded simultaneously compared to using SSE instructions. Just like SToRM, the vectorization is handled within AFF3CT code and a compilation flags allows setting the desired SIMD instruction. Note that this application is CPU-intensive and uses integers which are not supported by AVX instructions. As a consequence, there will be not results for AVX with AFF3CT. On nova and chifflet, all applications were compiled against gcc 4.9.2 and the machines were running a Linux Debian 3.16.43-2+deb8u5 (2017-09-19) x86 64 GNU/Linux with a kernel version 3.16.0-4-amd64. On skylake, applications were compiled against gcc 8.2.1 and the machine is running an Arch Linux 4.18.9-arch1-1-ARCH x86 64 GNU/Linux. On all architectures, hyperthreading was disabled.

B. Applications characterization

In order to better understand these applications, we start by characterizing them as CPU or memory-intensive. Note that since AFF3CT and SToRM use only integer operations, we cannot rely on the usual FLOPS and the roofline model to characterize the applications. For this reason, in order to characterize the applications, we vary the CPU frequency and observe the impact on the execution duration of the different applications. This does not determine if one application is more CPU-intensive than another, but rather whether it is CPU-intensive or memory-intensive. Note that HPL mainly computes matrix-matrix multiplications which are CPU-intensive. SVD Band and SVD Bulge have been characterized in the literature [START_REF] Haidar | An improved parallel singular value algorithm and its implementation for multicore hardware[END_REF].

Figure 2 shows the execution time according to the CPU frequency on chifflet. We used the same configurations as the ones described in Section III-A for all applications. The impact on SVD Bulge is as expected. Since it is a memoryintensive phase, CPU frequency should have little impact on its performance. SVD Band and HPL behaviors perform also as expected: the lower the frequency, the slower the application. AFF3CT and SToRM seem to have a similar behavior. This indicates that SToRM and AFF3CT are also CPU-intensive. This section presents SIMD instructions behavior when the power consumption reaches TDP and when it does not. To do so, for each application and each platform, we use the available SIMD instructions on each machine. For each configuration, we will explain the observed behavior of SIMD instructions and if TDP is impacting this behavior. This study is on applications socket power (Section IV-A) consumption, performance (Section IV-B), energy consumption (Section IV-C) and DRAM power (Section IV-D) and energy (Section IV-E) consumption. The goal here is not to compare the different platforms, but rather to study the behavior of SIMD instructions under different configurations. Figure 3 shows the execution time ratio (Figure 3a), socket power ratio (Figure 3b), socket energy consumption ratio (Figure 3c), DRAM power ratio (Figure 3e) and DRAM energy consumption ratio (Figure 3f) of the different applications on nova while using the different SIMD instructions when turboboost is activated and not activated. Figures 4 and5 show the same ratios for chifflet and skylake. For each application, the ratio is computed over the default SIMD instruction which is AVX2 for nova and chifflet and AVX512 on skylake.

A. Impact on applications socket power consumption

In this section, we present how SIMD instructions impact power consumption. We describe the results for each platform separately since they are not similar.

1) nova: SIMD at the same frequency: On nova (Figure 3b), for all CPU-intensive applications except AFF3CT, the larger the SIMD instruction vector size, the larger the power consumption. This is the typical behavior described in the literature [START_REF] Lento | Optimizing performance with intel ® advanced vector extensions[END_REF]. Regarding AFF3CT, SSE and AVX2 have the same power consumption. This is because even if AVX2 is known to consume more power (as shown for the other applications), the performance of SSE are better (as shown in Figure 3a and will be discussed in Section IV-B) and this balances the power consumption with AVX2. We will observe a similar behavior on chifflet.

Regarding SVD Bulge, which is memory-intensive, one can see that SIMD instructions have a very small impact on power consumption (an increase of 2.8% from SSE to AVX2 and 2.07% from AVX to AVX2). This slight impact is explained by the small impact on the performance as shown in Figure 3a which shows the use of SIMD instructions in SVD Bulge.

Finally, disabling turboboost shows the same behavior. As the frequency is lower, the power consumption when disabling turboboost is lower for all applications which is the normal behavior. As expected, for all applications, the power ratios of all instructions over AVX2 when enabling turboboost are the same as the ratios when turboboost is disabled. This indicates that for nova, since all instructions run at the same frequency, the type of used SIMD instructions is the major factor impacting power consumption and shows that, as expected, as the buffer size increases, the power consumption increases.

2) chifflet: SIMD at different frequencies: On chifflet (Figure 4b), the behavior is similar to nova: the larger the SIMD buffer size, the larger the power consumption for all CPU-intensive applications except AFF3CT. However, the difference between the power consumption of the different instructions is less important on chifflet compared to nova for HPL and SVD Band. As a matter of fact, on nova using AVX2 consumes 28% more than using SSE for HPL and 21% more for SVD BAND. On the other hand, on chifflet, using AVX2 only consumes 5% more than using SSE for HPL and 7% more for SVD BAND.

For these two applications, AVX is almost at the thermal design power while AVX2 reaches it. Because of that, for HPL and SVD Band, the frequency is reduced for AVX and AVX2. This explains why they have the same power consumption. Table IV shows the observed frequencies for all applications on chifflet and skylake. For HPL, AVX2 frequency is equal to the base frequency (2.4GHz) while AVX frequency is 2.6GHz. This explains why disabling turboboost shows the same values for performance and power consumption for AVX2 (as hinted in Table I). For SVD Band, the frequency is 2.6GHz for AVX2 and 2.8GHz for AVX. Thus, unlike HPL, AVX2 does not run at the base frequency, which explains the difference in performance and power consumption between AVX and AVX2 for SVD Band compared to HPL. This explains why disabling turboboost shows a greater impact on power consumption for SVD Band compared to HPL (9.76% of power consumption increase when using AVX2 compared to AVX for HPL 16.24% for SVD Band).

Regarding SVD Bulge, SIMD instructions have no impact on the power consumption. Figure 4b shows that AVX consumes slightly less power than to AVX2 or SSE, but this difference is 1.6%, which is in the error measurement range. When setting all processors to the same frequency (No-Turboboost plot), one can see a slight difference between the instructions (6.5W for SSE and 4W for AVX compared to AVX2). This is due to the use of SIMD instructions in SVD Bulge (as will be shown in section IV-B). However, the variation is small (at most 3.5%).

SToRM shows an interesting behavior. According to Table IV, AVX2 reaches the same frequency as SSE. Thus, the ratio SSE/AVX2 is the same regardless of turboboost. Moreover, just like AFF3CT, it uses only integers, which means that for all instructions, only half of the buffer is used. This impacts power consumption. As a matter of fact, SToRM consumes 101W per socket while HPL consumes 124W for AVX2. Note that we do not compare the applications since they may not have the same computation intensity, we only state that using half of the buffer most likely leads to less power consumption and thus higher frequencies can be used.

Finally, regarding AFF3CT, unlike the other applications, the larger the SIMD buffer size, the lower the power consumption. For instance, using SSE consumes 4.3% more power than using AVX2. This is due to the fact that SSE outperforms AVX2 for this application. Moreover, when we look at the power consumption when turboboost is disabled, it is the same. This means that for AFF3CT on chifflet, the frequency is the parameter influencing the power consumption because SSE outperforms AVX2.

3) skylake AVX512: Applications on skylake exhibit a completely different behavior compared to the other platforms: larger SIMD buffer does not mean more power consumption. This means that AVX512 shows better power consumption compared to other instructions. This is different from the behavior observed on nova and chifflet.

On Figure 5b, for all applications, the power consumed when using AVX512 (thus SIMD instruction with the largest buffer) is lower or equal (equal for HPL) compared to using other buffers. The reason why SSE, AVX and AVX2 consume more power than AVX512 is that they reach higher frequencies as shown in Table IV. This can be already be observed for AFF3CT and SToRM when disabling turboboost. We will start by describing the behavior for these two applications before moving to HPL and SVD.

AFF3CT and SToRM have a specific behavior on skylake: AVX2 and AVX512 do not run at the same frequency as specified in Table II. As a matter of fact, AVX2 frequency is the same as SSE whereas AVX512 frequency is 2.4GHz when turboboost is enabled and is the base frequency when turboboost is disabled. Thus, when turboboost is disabled, SSE, AVX2 and AVX512 run at the same frequency.

AFF3CT shows a similar behavior to its behavior on chifflet: SSE has a larger power consumption compared to AVX2 despite running at the same frequency (an increase of 3.8%). When running at the same frequency (when turboboost is disabled), we observe a similar behavior where SIMD instruction (SSE) with the best performance consumes more than the others (using SSE consumes 6.4%more than using AVX2 and 11.6% more than using AVX512).

Regarding SToRM, AVX2 reaches the thermal design power and consumes more energy than SSE (by 3%). Besides, SSE consumes 8.65% more than AVX512 because of AVX512 frequency. This can be verified when observing the results when turboboost is disabled. Since all instructions run at the same frequency, AVX512 consumes more power than SSE. AVX2 and AVX512 have roughly the same power consumption when disabling turboboost.

For HPL, all SIMD instructions reach the TDP which is why they have the same power consumption. Even when turboboost is disabled, AVX2 is almost at TDP.

Regarding SVD Band, SSE, AVX and AVX2 are at the thermal design power. One interesting observation is that, as shown in Table IV, on average, SSE and AVX run at lower frequencies compared to the values presented in Table II for HPL and SVD Band. This is the behavior described by Intel when processors reach the thermal design power [START_REF] Lento | Optimizing performance with intel ® advanced vector extensions[END_REF].

Finally SVD Bulge also shows a similar behavior where SSE is the most power consuming compared to AVX and AVX2 since it consumes 7.45% more power compared to AVX and 6.5% compared to AVX2. This is due to frequency which is higher for SSE as shown in Table IV. AVX512 consumes less than the other instructions. The frequency impact can also be verified when observing the plots where turboboost is disabled. In this case, the power consumption of SSE, AVX and AVX2 is the same, while AVX512 is still below because its frequency is lower when turboboost is disabled.

For HPL, SVD Band and SVD Bulge, most of the time AVX2 consumes more than AVX512 even when turboboost is disabled. We believe that this is because of the frequency. Recall that AVX2 frequency is 2.4GHz while AVX512 frequency is 1.9GHz. In order to validate our assumption, we ran these three applications at 1.9GHz. Note that we focus only on these three applications since on SToRM and AFF3CT, all SIMD instructions already run at the same frequency when turboboost is disabled. Note also that differences will be observed since running applications with AVX512 instructions will have the frequency vary depending on the CPU load and performed instructions, whereas we fixed the frequency to 1.9GHz so it cannot vary. Figure 6 shows the power consumption of HPL and SVD applications when running on skylake at 1.9GHz. Regarding power consumption, HPL and SVD Bulge have the expected behavior observed on nova. For SVD Band however, AVX2 and AVX512 shows the same power consumption. We believe that this is because they have the same performance. As a conclusion, for CPU-intensive applications, both SIMD instructions and frequency seem to impact power consumption. On the other hand, for Memory-intensive applications, frequency seems to be the most impacting factor.

B. Impact on performance

In this section, we study the performance behavior of applications. Note that in this section, we will provide explanations on why applications do not have the expected behavior. We will however not explain why an application running using AVX2 is not at least half as slow as running AVX512 for instance. This is because sometimes this problem is due to the algorithm and how the application is designed, which is out of our expertise since we did not design any of these applications.

As expected, on all platforms, the behavior of CPUintensive applications is similar: the larger the buffer size, the better the performance, except for AFF3CT. We will explain AFF3CT's behavior in the end of this section.

On nova and chifflet, all CPU-intensive applications show great performance gain when increasing the vector size. For instance, HPL shows a ×3.54 improvement on nova and a ×2.89 on chifflet when using AVX2 compared to SSE. On chifflet, HPL is the only application which power consumption reaches the thermal design power. As such, its frequency is lowered as shown in Table IV. The average frequency observed with HPL is the same as the base frequency. This is why the performance when disabling turboboost for AVX2 are so close to the performance with turboboost. Using AVX2 over SSE when disabling turboboost reduces HPL execution time by 71%. This is mainly due to SSE performance which are reduced by a factor of ×1.19 when disabling turboboost.

For SVD Band and SToRM, disabling turboboost shows an even larger impact when increasing the SIMD buffer size. As a matter of fact, using AVX2 over SSE reduces SVD Band execution time by 67% and SToRM execution time by 48%.

On skylake however, using AVX512 does not seem to improve performance as much as using AVX2 compared to SSE for all CPU-intensive applications (except AFF3CT). For HPL, one of reasons lies behind the memory bandwidth. As a matter of fact, the configuration HPL + AVX512 reaches the maximum sustainable memory bandwidth observed using STREAM [START_REF] Mccalpin | Memory bandwidth and machine balance in current high performance computers[END_REF] (approximately 30GB/s) which explains why its performance are limited compared to AVX2. Moreover, because of TDP, HPL frequency when using AVX512 is lower than when using AVX2 as stated in Section IV-A. Note that we tried running HPL on skylake while forcing the frequency to 1.9GHz. This showed a small improvement of the performance ratios (1.16 with the normal behavior and 1.21 when forcing the frequency to 1.9GHz), but not as great as expected. The frequency has also an impact on SToRM performance since disabling turboboost shows a better improvement of AVX512 over AVX2 (9.25% using the default configuration and 30% when disabling turboboost). This is because for SToRM, as shown in Table IV, all instructions run at the same frequency when disabling turboboost on skylake. Finally, SVD Band shows similar performance results with the configuration that we used. Even when forcing AVX2 to run at 1.9GHz, it still provides the same performance as AVX512. We could not find a reason behind this behavior. At 1.9GHz, SVD Band does not reach the thermal design power with AVX512, the frequency does not have an impact and neither does the memory bandwidth or the cache misses. We will no further investigate this behavior.

For SVD Bulge, AVX has little to no impact on the performance on nova and chifflet. On skylake however, one can see that there is an impact. We studied the performance (in flop/s) of SVD Bulge using LIKWID and FLOPS DP group. The group provides double precision floating point performance in addition to vectorization performance and ratios. Vectorization shows an impact of few GFlops/s on SVD Bulges while the total performance does not exceed 44 GFlops/s. This means that despite being memory intensive, this application still computes floating point operations and still uses vectorization which explains why SIMD instructions have an impact on its performance.

AFF3CT exhibits a behavior that is completely opposite to the other applications: the smaller the SIMD instruction buffer size, the better the performance. In other words, for this application, with the parameters that we used, it is better to use SSE rather than AVX2 or AVX512. As stated in Section III-A, AFF3CT handles frames which are loaded in the memory. As a consequence, the larger the buffer (AVX2 or AVX512), the greater the number of loaded frames. Therefore, for the configuration that we used for AFF3CT, for AVX2 and AVX512, the load exceeds the cache. We compared the ratio of cache misses of the different applications when using SSE, AVX2 and AVX512 on skylake (recall that AVX is not supported by AFF3CT). These results are shown in Table V and represent an average over all sockets. Note also that the results show only CPU-intensive applications (thus we exclude SVD Bulge). Although all applications show a difference when comparing their cache miss ratios between AVX512 and SSE or AVX2, AFF3CT shows the highest difference, especially between SSE and AVX512. As a matter of fact, with AFF3CT, using AVX512 generates more than 16.5 times more L3 cache misses than using SSE while this ratio is at most 1.80 for the other applications. Further details are presented in [START_REF] Cassagne | Beyond Gbps Turbo decoder on multi-core CPUs[END_REF].

C. Impact on applications socket energy consumption

Regarding energy consumption, on all platforms, the performance seems to be the major factor impacting energy consumption. This is because SIMD instructions have a large impact on performance compared to power consumption. Thus, except for AFF3CT, for which SSE is better, it is better to use the SIMD instruction with the largest buffer.

Regarding frequency, it has more impact on chifflet and skylake since each SIMD instruction runs at a different frequency. For these two platforms, disabling turboboost shows little to no impact on energy consumption.

D. Impact on applications DRAM power consumption

Figures 3e, 4e and 5e present the DRAM power on nova, chifflet and skylake.

For all applications, on all platforms (except SToRM on skylake), the larger the buffer size, the higher the DRAM power consumption. In order to explain this behavior, we studied the impact of vectorized instructions on the applications memory throughput. Table VI shows the mean memory throughput of all applications on all platforms over all sockets. Note that not all applications show the same memory throughput over all socket. For AFF3CT, the sockets may have different memory throughput but they remain in the same order of magnitude (a difference of 4000 MByte/s over 30000 MByte/s between the sockets on skylake). For SToRM however, the difference between the sockets is noticeable. It is especially the case when using SSE since it varies from 330 to 1900 MByte/s between socket 0 and socket 4 on skylake. For SVD Bulge, only socket 0 has a large memory throughput compared to the other sockets. For SVD Bulge, we present only the memory throughput for socket 0.

The results show a strong correlation between the DRAM power consumption and the memory throughput. First of all, one can notice that, for CPU-intensive applications, as the size of the instruction buffer increases, the memory throughput increases as well [START_REF] Inoue | How SIMD width affects energy efficiency: A case study on sorting[END_REF] which also leads to larger power consumption. SToRM seems to exhibit different behavior on skylake where SSE has a larger DRAM power consumption. This behavior is still correlated to memory throughput. Note that for SToRM, the power consumption is low compared to the other applications: it is roughly half the power consumed by the other CPU-intensive applications for AVX2 and AVX512. SToRM throughput (Table VI) also shows the small memory movement of SToRM compared to the other applications on all platforms. For instance, on skylake using AVX512, SToRM memory throughput is ×0.022 HPL throughput and ×0.037 SVD Band throughput. Finally, as stated before, SToRM memory throughput is different between the sockets. However, there is still a correlation between each socket DRAM power consumption and memory throughput. In other words, for each socket, the larger the memory throughput the larger the DRAM power consumption.

Memory-intensive application SVD Bugle shows almost no impact of SIMD instructions on DRAM power consumption (with a the highest ratio is of 1.01 between SSE and AVX512). Moreover, only the first socket has an impact on power consumption. The other sockets have a very low power consumption (15.9W for socket 0 and between 4.3W and 4.9W for the other three sockets).

Finally, turboboost has also an impact on the DRAM power consumption. This is also due to memory throughput which is lower when turboboost is disabled (as there are less requests per unit of time since the processor is slower).

E. Impact on applications DRAM energy consumption

Figures 3f, 4f and 5f present the DRAM energy consumption on nova, chifflet and skylake.

In order to measure the DRAM energy consumption, we just multiplied the total execution time by the DRAM power consumption. The results show that on all platforms and for all applications, the most performing application is the least energy consuming. Thus, for all applications except AFF3CT, using the SIMD instruction with the largest buffer provides the best energy consumption.

F. Key findings

This study provided an insight on how thermal design power and SIMD instructions impact performance, power and energy consumption for both processor and memory. From our observations, one can conclude that :

• Performance and power consumption are more and more related. Because of TDP and SIMD instructions, core frequency may be lowered which directly impacts performance. • For most CPU-intensive applications (except applications with special design like AFF3CT), the larger the buffer size, the better the performance and the energy consumption. Moreover, turboboost has no impact on energy consumption when using the instruction with the largest buffer. Thus, if power is more important than performance, one can start with disabling turboboost. • Frequency has more impact on power and performance of memory-intensive applications than SIMD instructions. • The larger the SIMD buffer size, the larger is the power consumption, except for AVX512 since it runs at a lower frequency.

• DRAM power consumption is strongly correlated to memory throughput. Thus, larger SIMD buffer sizes lead to larger DRAM power consumption.

V. RELATED WORK

Many studies focused on the energy consumption of vectorization. [START_REF] Hackenberg | An Energy Efficiency Feature Survey of the Intel Haswell Processor[END_REF], presents recent changes in Intel Haswell processors such as AVX frequencies and voltage regulators. In [START_REF] Chen | An empirical study of performance, power consumption, and energy cost of erasure code computing for HPC cloud storage systems[END_REF], the authors compare the time, power consumption and energy consumption of an erasure code running on traditional x86 and more recent SIMD platforms. They showed that using SIMD instructions on such applications reduces execution time and energy consumption. In [START_REF] Schuchart | The Shift from Processor Power Consumption to Performance Variations: Fundamental Implications at Scale[END_REF], the authors present a comparison of Sandy Bridge Processors and Haswell-EP processors in order to show how the new changes, like the TDP, challenges the performance analysis. In [START_REF] Cebrian | Improving Energy Efficiency through Parallelization and Vectorization on Intel Core i5 and i7 Processors[END_REF], the authors compared multithreading and vectorization on different processors. They showed the benefit of vectorization over multithreading for energy consumption. Special instruction set impact was also studied, like load and store AVX instruction set [START_REF] Jakobs | On the energy consumption of load/store AVX instructions[END_REF]. In [START_REF] Jakobs | Examining energy efficiency of vectorization techniques using a gaussian elimination[END_REF], the authors study the energy consumption of different implementations of Gaussian elimination on different architectures. These studies present a comparison of an execution with and without SIMD instructions, while we present a comparison of the instructions.

Studies like [START_REF] Gepner | Early performance evaluation of AVX for HPC[END_REF] and [START_REF] Inoue | How SIMD width affects energy efficiency: A case study on sorting[END_REF] can be considered complementary to our work. In [START_REF] Gepner | Early performance evaluation of AVX for HPC[END_REF], the authors studied the impact of SSE and AVX on applications performance. This work was done on older architectures. In [START_REF] Inoue | How SIMD width affects energy efficiency: A case study on sorting[END_REF], the authors compared the execution time, power and energy consumption of an AVX and SSE implementation of a sorting algorithm. They also showed the impact of varying the memory bandwidth on the performance and energy consumption of SIMD instructions. This work is complementary to our work since it studies the impact of SIMD instructions on another algorithm and studies other parameters such as memory bandwidth, while we studied the DRAM power consumption. Moreover, we focused on different applications with different profiles. [START_REF] Mazouz | An incremental methodology for energy measurement and modeling[END_REF] presents an energy model using codelets. They also study the energy consumption of scalar and vectorized instructions (using SSE and AVX2). They showed the effect of different profiles on the energy consumption of SSE and AVX2 instructions. They also showed the effect of data located on L2 and L3 cache on the energy consumption of the codelets. Our work is complementary since we study the power and energy consumption of applications for both socket and memory consumption, and study the impact of TDP on the performance and power consumption. In [START_REF] Cebrián | Performance and energy impact of parallelization and vectorization techniques in modern microprocessors[END_REF], the authors compared different vectorized instructions (SSE and AVX) while varying the number of threads using different Intel and ARM platforms. They also evaluated when turboboost improves energy consumption. This work is the closest to our work since the authors compared vectorized instructions and studied the impact of turboboost. However, the machines they used (Intel Ivy and Sandy Bridge) were such that using SIMD instructions had almost no impact on power. Moreover, we presented a study on the impact of SIMD instructions on DRAM power consumption.

VI. CONCLUSION

In this work, we studied the impact of vectorization and thermal design power on processor and DRAM power consumption. For that purpose, we used 3 different architectures and 5 applications with different behavior. Our conclusions showed that because of thermal design power, performance and power become less and less independent. As a consequence, when trying to understand an application performance, studying its power consumption and frequency can help understanding its behavior. Moreover, our study showed that although using SIMD instructions with larger buffer size improves performance and energy consumption, it has a negative impact on both DRAM and processor power consumption. However, AVX512 seem to have a different behavior where its power consumption is lower than the other instructions despite providing better performance.

Since power consumption is becoming a major problem, using power capping techniques may provide a good leverage to reduce power consumption. As a consequence, in the future, we plan to study the behavior of the processors when a power cap is applied. This will be especially interesting for application like AFF3CT and SToRM which do not reach TDP. We also plan to study how using hyperthreading in addition to vectorized instructions can impact application behavior.

Fig. 2 :

 2 Fig. 2: Frequency impact on execution time on chifflet

Fig. 3 :Fig. 4 :Fig. 5 :

 345 Fig. 3: Comparison of SIMD instructions with and without turboboost on nova

Fig. 6 :

 6 Fig. 6: Comparison of SIMD instructions power consumption when running at AVX512 frequency (1.9GHz) on skylake for HPL and SVD

TABLE I

 I

: Execution time of HPL on chifflet with AVX2 while enabling and disabling turboboost (average on 5 runs).

Table II :

 II CPUs, 8 cores per CPU and 64GB of memory. It provides SSE, AVX and AVX2 instruction sets. Note that for this processor, turboboost frequency does not depend on the SIMD instruction being used.• chifflet: Chifflet is a cluster located in Grid'5000 Lille site and is equipped with 8 Intel Xeon E5-2680 v4. We used chifflet-1 for our experiments. It is equipped with 2 CPU, 14 cores per CPU and 768 GB of memory. It provides SSE, AVX and AVX2 instruction sets. For this processor, the frequency varies according to the SIMD instruction used (as shown in TableII). • skylake: skylake is a server located in our lab at Telecom SudParis. It has 4 Intel(R) Xeon(R) Gold 6130, each equipped with 16 cores. Each NUMA node has 63 GB of memory. Thus, skylake has 64 cores and 252GB of memory in total. It provides SSE, AVX, AVX2 and AVX512 instruction sets. Like chifflet, SIMD instructions have an impact on the turboboost frequency. Moreover, AVX512 has also an impact on the base frequency as shown in TableII. All platforms run under Intel Pstate with powersave governor. Table

• nova: The Nova cluster, located in Grid'5000 Lyon site, is equipped with 23 Intel(R) Xeon(R) CPU E5-2620 v4. All experiments were run on nova-11. It is equipped with 2

TABLE II :

 II Target platforms characteristics extracted from processors documentation

Table IIIa

 IIIa

		N	NB	(PxQ)		N	NB
	nova	58912	224	(4x4)	nova	20000 200
	chifflet	100996 224	(4x7)	chifflet	28000 200
	skylake	91840	224	(8x8)	skylake	48000 200
		(a) HPL				(b) SVD
	. We used HPL
	version 2.2 compiled with OpenMPI 3.3.1 and with MKL
	library version l mkl 2018.3.22. MKL allows choosing the
	right SIMD instructions by setting the environment vari-
	able MKL ENABLE INSTRUCTIONS to SSE4 2, AVX,
	AVX2 or AVX512 (for skylake only). Note that HPL is a
	CPU-intensive application.

TABLE III :

 III HPL and SVD setup for nova, chifflet and skylake

TABLE IV :

 IV Average observed frequency, in GHz, on chifflet and skylake over all cores.

TABLE V :

 V Applications average L3 cache misses over all socket on skylake. The results are presented as a ratio of AVX512 cache misses over SSE or AVX2 cache misses.

TABLE VI :

 VI Mean memory throughput over all socket, in MByte/s, for each SIMD instruction on nova, chifflet and skylake. For SVD Bulge, only socket 0 memory throughput is presented.

LIKWID commit:

233ab943543480cd46058b34616c174198ba0459

ACKNOWLEDGMENT

Experiments presented in this paper were carried out using the Grid'5000 experimental testbed, being developed under the INRIA ALADDIN development action with support from CNRS, RENATER and several Universities as well as other funding bodies (see https://www.grid5000.fr).

We would like to thank Mathieu Faverge for his help on PLASMA. We would also like to thank Laurent Noé for his help on SToRM. Finally, we would like to thank Adrien Cassagne for his help on AFF3CT.

VII. APPENDIX

In this appendix, we provide, for each application, how it was compiled to use the desired SIMD instruction. We also provide the command line that we used for our experiments. If input files are required, we describe how we obtained them (this is especially the case for SToRM). Section III-A provides some settings (sur as matrix size for HPL).

A. HPL

HPL is compiled against MKL which allows choosing the right SIMD instruction through the environment variable MKL_ENABLE_INSTRUCTIONS. It can be set to SSE4 2, AVX, AVX2 and AVX512 (for skylake only).

HPL uses a configuration file (the values that we changed are described in Table IIIa and is simply launched with mpirun: mpirun -n $nb_cores xhpl

B. SVD

SVD is also compiled against MKL. Thus setting the desired SIMD instruction is done the same way as HPL.

The command line that we used is the following :

./time_dgesvd_tile --threads=$cores --n_range=$N:$N:1 --nb=$NB --nowarmup --nodyn --nocheck where $cores is the number of cores that we want to use, and $N and $NB are described in Table IIIb.

C. SToRM data generation and execution command line

SToRM can be compiled using the provided Makefile with the desired SIMD compilation flags (msse4.2, mavx2 and mavx512bw). We chose to run the storm-nucleotide program (and not the stormcolor). The command line that we used to run the program is:

./storm-nucleotide-sse42-x-gcc -g data/Homo_sapiens.GRCh38.d na_rm.chromosome.22.fa -r data/test_1000000.fq -N $nb_cores -i 15 -z 180 -t 200 -o /dev/null where -g takes a genome file. We downloaded the Homo sapiens file from http://bioinfo.lille.inria.fr/yass/data/ Homo sapiens/GRCh38.dna rm.chromosome.22.fa -r is the reads file to map against the genome. In order to generate the test 1000000.fq file, we used the ncbi toolkit (https://www.ncbi.nlm.nih.gov/sra/docs/toolkitsoft/). It comes with three separate tools (NGC, NCBI-VDB and SRA Toolkit). The tool provides a command fastq_dump. It generates the n first reads (-X option) from an input file (SRR7764388) and writes it to output (-Z). SRR7764388 is a file available at the ncbi database (https://www.ncbi.nlm.nih. gov/sra/SRR7764388).

Note that we used 1000000 on nova and chifflet and 10000000 on skylake.

D. AFF3CT command line

For AFF3CT, we used the compilation guideline provided in the documentation. Just like SToRM, compiling with the desired SIMD instructions is done through compilation flag. Launching the execution is done with ./bin/aff3ct --sim-type BFER --sim-cde-type TURBO -m 1.0 -M 1.0 -K 6144 --dec-type TURBO --dec-implem FAST --dec-sub-simd INTER -i 6 --itl-type LTE --dec-sf-type LTE_VEC --dec-sub-max MAX


Thomas Nussle

Pascal Thibaudeau

Stam Nicolis

Dynamic magnetostriction for antiferromagnets

In this paper we propose a Hamiltonian framework for the dynamics of magnetic moments, in interaction with an elastic medium, that can take into account the dynamics in phase space of the variables that describe the magnetic moments in a consistent way. While such a description involves describing the magnetic moments as bilinears of anticommuting variables that are their own conjugates, we show how it is possible to avoid having to deal directly with the anticommuting variables themselves, only using them to deduce non-trivial constraints on the magnetoelastic couplings. We construct the appropriate Poisson bracket and a geometric integration scheme, that is symplectic in the extended phase space and that allows us to study the switching properties of the magnetization, that are relevant for applications, for the case of a toy model for antiferromagnetic NiO, under external stresses.

I. INTRODUCTION

The interest in precision devices based on voltagecontrolled ferromagnetism is focusing the attention of theory and applications towards developing materials with strong intrinsic ferroelectric and/or ferromagnetic response; what is of particular interest is that these can be produced through an appropriate spin-lattice coupling, that can be controlled, experimentally, in real materials 1 . What is, particularly, fascinating is that such couplings lead to novel situations, where it is possible to control the magnetic properties through electric fields and vice versa. This class of phenomena is known under the name of magneto-electric effects, and there are several ways in which such magneto-electric effects can be engineered in magnetic materials.

For instance, antiferromagnets (AFs) are the magnetically ordered materials which could be the most suitable for realizing the fastest spintronic devices 2 . Recent progress in describing spin transport and spintransfert torque (STT) effect in such AF devices open a route towards multilevel memory devices with switching speeds, that could exceed those of devices made of ferromagnetic materials and semiconductors 3 . Moreover reversible strain-induced magnetization switching in ferromagnetic materials has been reported that also allows the design of a rewritable, non-volatile, non-toggle and extremely low energy straintronic memory 4 . Recently a piezolectric, strain-controlled AF memory, insensitive to magnetic fields, has been tested, as an example of controlling magnetism by electric fields in multiferroic heterostructures 5 . However a complete theoretical description of all such phenomena is, still, very much, work in progress, since, at the scales, that are relevant, it is not possible to separate the scales of the magnetic, electric and elastic responsesall must be treated together.

It is, indeed, noteworthy that no direct local coupling between electromagnetic field components is allowed in the vacuum; the magnetic response to an electric field is necessarily mediated by atomic position arrangements of magnetic moments, in a potentially ionic environment, that produce non-local mechanical strains in response 6,7 . These illustrate the different physical origins of the spin-lattice coupling-in particular, either as emerging from a macroscopic magnetostriction [START_REF] Trémolet De Lacheisserie | Magnetostriction: Theory and Applications of Magnetoelasticity[END_REF] or, at a more microscopic scale, from the usual spin-orbit coupling in materials [START_REF] Bruno | [END_REF][START_REF] Bruno | Ferienkurse Des Forschungszentrums Jülich[END_REF] , however, considerably enhanced from the effects of magneto-elastic effects.

The idea to couple the dynamics of mechanical and magnetic degrees of freedom in matter is, of course, not new, and many strategies to produce a general conceptual framework have been explored to date, covering a quite broad spectrum of approaches and applications. Illustrations of the interplay between magnetism and mechanical elasticity in matter are so diverse that a large variety of techniques and ideas have been developed, with just as many different motivations [START_REF] Bardzokas | Filshtinsky[END_REF] .

It is, also, relevant to mention magnetic molecular dynamics 12 , which emerges from the fluid mechanics framework and assumes position-dependent magnetic interactions between particles, under a tight timeevolution of their own spins. From granular fields theories, apparently radically different perspectives can be gleaned, such as Edwards field theory 13 , which focuses on the long distance correlations between mechanical stresses in large systems, assuming material properties as point-particle stress distributions and interactions.

A common focus of all these approaches is, however, not at all obvious.

In addition, until recently, sufficient experimental control over the materials that could be candidates for realizing such effects in practice was not that common. It is this last fact that has changed in recent years. This, in turn, has made the call for theoretical tools more urgent. While it has become much easier to develop computational strategies, thanks to the progress in hardware and software, theoretical concepts have tended to lag behind the events.

Magneto-elastic phenomena are, typically, viewed from the continuum mechanics perspective where the magnetic properties of materials are incorporated into constitutive nonlinear laws of electro-conductive bodies [START_REF] Maugin | Continuum Mechanics of Electromagnetic Solids[END_REF] . Even if multiscale models have been investigated for magnetoelastic couplings 15 , with good experimental comparison for single crystals and polycrystalline samples 16 , localization and homogenization procedures are still the mainstays for deducing the constitutive laws for polycrystalline media. Homogenization methods for multiscale mechanics assume the existence of well-separated scales and different relations among length scales, that lead to different effective equations, which, in turn, represent the corresponding different physical effects, appropriate at each scale [START_REF] Mei | Homogenization Methods for Multiscale Mechanics[END_REF] .

Renormalization group ideas, of course, are instrumental in providing a description of how the dynamics depends on the scale; what has been lacking, to date, is a unified construction of the dynamics of elastic and spin degrees of freedom, within a common Hamiltonian approach, at any particular scale. In particular, what has, curiously, not been fully exploited is the resolution of the constraints, that the spin degrees of freedom are subject to, since they are their own canonically conjugate variables, as their equations of motion are first order. This means that they can be most usefully expressed as multilinear combinations of anticommuting variables. As we shall see, this can lead to useful insights, in practice, about the dynamics of the coupled system.

While the generalization of the canonical formalism, that incorporates the constraints implied by spin degrees of freedom has been worked out a long time ago [START_REF] Casalbuoni | [END_REF][19][20] , it has not been applied to concrete situations of physical interest, simply because these, up to now, were not of practical relevance.

In the present paper we fill this gap. We provide a complete description of the combined dynamics of the elastic and of the spin degrees of freedom in phase space and show how the anticommuting nature of the spin degrees of freedom can be-indirectly-tested. It is the refined control over these that is the particular novelty of our approach. In future work we shall show how to incorporate stochastic effects, due to noise and/or disorder.

We consider our medium as the limit of a large number of sites; to each site is assigned a spin and a local mechanical deformation variable. We deduce the corresponding equations of motion from a formalism, that treats mechanical and spin degrees of freedom in a unified way and study the dynamics of this system, comparing it to previous studies for the purely magnetic part 21 and for the magneto-elastically coupled system 22 . We stress that the novelty of the current approach lies, on the one hand, in dealing with a consistently closed system; on the other hand, in being able to extract the information that the commuting degrees of freedom, that describe the spins are, in fact, composite operators of anticommuting fields.

The paper is organized as follows, in section II we construct a classical Hamiltonian for the mechanical, and the spin degrees of freedom and discuss their couplings. In section III, we introduce a particular Poisson bracket, which allows us to obtain the equations of motion for the time-evolution of all the dynamical variables. These are, then solved, in section IV using a symplectic integration scheme. Special attention is devoted to the stability analysis of the schemes and to the order of the different symplectic operators. This analysis is presented in the appendix. Finally in section V we illustrate the formalism for the case of a simple toy model for NiO, consisting of two spins, interacting through antiferromagnetic exchange, and subject to an external SST, as well as under an external stress.

II. THE HAMILTONIAN

In this section we propose a Hamiltonian for the combined system of elastic and spin degrees of freedom. These parts have, of course, been considered before [START_REF] Trémolet De Lacheisserie | Magnetostriction: Theory and Applications of Magnetoelasticity[END_REF] , so we shall review the salient features, before introducing our approach.

The starting point for the elastic degrees of freedom is the framework of mechanics for elastic solids, within the régime of the approximation of small deformations, supplemented by the assumption of perfect mechanical micro-reversibility [START_REF] Landau | Theory of Elasticity[END_REF] . These statements imply that an elastic material is characterized by a collection of N timedependant elastic deformations, which are described by a set of symmetric Cauchy strain tensors ǫ i IJ (t) where I, J are spatial indices, ranging from 1 to 3, and and assigned to each lattice site, i ∈ {1..N }, of the material. At equilibrium, each of these variables is independent of time and we assume that such a state not only exists, but can be relevant for the time scales of interest.

The total internal energy of this system can be described in terms of the interactions between these variables, that take into account the elastic mechanical part in the most general form 24,[START_REF]Nonlocal Continuum Field Theories[END_REF] ; to lowest, non-trivial, i.e. quadratic, order, we have

H mech V 0 = 1 2 N i=1 C i IJKL ǫ i IJ (t)ǫ i KL (t), (1)
where an implicit sum on the repeated space indices is understood and where C i IJKL are the elastic constants with V 0 as the reference volume.

In this expression we have assumed that the different sites of the material do not interact through elastic deformations (they will only interact through spin exchange forces, to be spelled out below).

For a homogeneous material, these elastic constants can be expressed in terms of the two Lamé parameters [START_REF] Landau | Theory of Elasticity[END_REF] , designated, as usual, by (λ, µ), in Cartesian coordinates as follows:

C IJKL = λδ IJ δ KL + µ (δ IK δ JL + δ IL δ JK ) (2)
This mechanical system can be excited by an external stress σ ext IJ , through the coupling

H ext V 0 = - N i=1 σ ext IJ ǫ i JI . (3)
While this stress can be, in general, a time-and spacedependent quantity, it is, nonetheless, assumed to be uniform over the sites of the material. Now we introduce the magnetic degrees of freedom, noted S i I (t) and describe how their dynamics affects the mechanical degrees of freedom, through the so-called magnetoelastic coupling. The contribution of the magnetoelastic coupling to the Hamiltonian can be deduced by assuming (global) Galilean invariance on the one hand and imposing invariance under the symmetries of the point group 26 on the other. For an isotropic medium, these requirements lead to the expressions for the total energy and the total magnetization 27 . The total internal energy is, thus, a sum over the lattice U (ǫ i IJ , S i I ) 28 . Consequently, in the approximation of small deformations, the expansion of the magnetoelastic coupling energy will contain even powers of the magnetization S i . It will, however, contain all powers of the strain tensor, ǫ.

The linear part in ǫ i IJ defines the linear, or first-order, magnetoelastic coupling, which is responsible for the magnetostriction, while the quadratic terms in ǫ i IJ define the second-order magnetoelastic coupling, that are responsible for "nonlinear response" in the elastic properties of magnetic media, that have been, also, identified as describing "anomalous effects".

Higher order terms are, usually, ignored, since their coefficients are assumed to be much smaller than the lower order terms, already discussed.

Therefore, the contribution to the Hamiltonian of the magnetoelastic coupling can be written as

H ME = N i,j=1 B (1)i,j IJKL ǫ i IJ (t)S i K (t)S j L (t) + B (2)i,j IJKLMN ǫ i IJ (t)ǫ j KL (t)S i M (t)S j N (t) + . . . (4)
If one compares the last term of equation ( 4) with equation (1), it is easy to see that a pair of effective elastic constants can be introduced that depends on the value of the spins on each site, i.e.

C eff,i,j IJLK ≡ C ij IJLK + 1 V 0 B (2)i,j IJKLMN S i M (t)S j N (t). (5)
This allows a possible interpretation of the "anomalous" temperature dependence in the elastic constants in iron single crystals as resulting from the competition between spin ordering and diffusion effects 29 .

In Appendix A we recall that the most effective "classically equivalent description" of the spin degree of freedom, that can capture such effects is not through the, commuting, variables S but their, anticommuting, "Doppelgänger" ξ, related to the S i I through

S i I ≡ - ı 2 ǫ IJK ξ i J ξ i K . (6)
The ξ can be identified with Majorana fermions, which have found many applications recently in condensed matter systems 30,31 , where new methods for controlling spin degrees of freedom have been developed. It is interesting to remark that the anticommuting variables, ξ i k are not Grassmann variables, satisfying {ξ i I , ξ j J } = 0; but, rather, {ξ i I , ξ i J } = δ ij δ IJ , i.e. that the ξ I generate a Clifford algebra on each site [START_REF] Lounesto | Clifford Algebras and Spinors[END_REF] . It is in this way that the S i I , defined through eq. ( 6), satisfy the angular momentum algebra.

If one were tempted to simply replace S by ξ in eq.( 4), it is interesting to remark that, for N = 1, and because of the symmetries of B (1) , H ME = 0, which implies that the dynamics as it is cannot be encoded by this Hamiltonian.

On the other hand, this allows us to understand the constraints on the allowed terms in the true Hamiltonian. They must, necessarily, involve more than one spins. Indeed, we can construct expressions that are multi-linear combinations of the S i I on different sites, potentially, up to order N , the number of sites, since no two identical ξ i I variables occur in the same monomial. This is, therefore, a nice way of automatically organizing the multispin terms of the Hamiltonian.

In the formalism of Atomistic Spin Dynamics, since the magnetic moments are localized 33,[START_REF] Eriksson | Atomistic Spin Dynamics: Foundations and Applications[END_REF] , it is customary to consider spatial averages, around each site, defining an effective macroscopic localized spin

S eff I = S i I i (7)
which implies that multi-linear expressions, S I S J no longer vanish identically, as was imposed by the anticommuting nature of ξ before. We can, thus, understand the relevance of this averaging procedure, in terms of the description of the spin degrees of freedom in terms of the anticommuting variables.

The microscopic mechanism behind equation ( 4) has been advocated a long time ago [START_REF] Van Vleck | [END_REF][36][37] , as a model for a two-body interaction, that is itself a pedagogical version of the quantum theory of interacting magnons and phonons 38 . In eq. ( 4) what has been left unspecified are the properties of the tensor B (1) under exchange of the indices i and j, that label the sites. What is customary is to assume that B (1) , in fact, does not depend on i and j at all-it is homogeneous across the material. By stopping the expansion at first order in the mechanical deformation in (4) and assuming homogeneity in the material, we, therefore, end up with the following expression

H ME = B (1) IJKL ǫ IJ (t)S K (t)S L (t). ( 8
)
In summary, the microscopic theory underlying Eq.( 8) describes the interaction of three particles, two of which are spin states, and one is the state of the elastic deformation. However the spin states, are, in fact, bound states of more "fundamental" entities, the anticommuting ξ.

B (1) can then be interpreted as the vertex of this interaction, that describes a spinning particle, that does not directly interact with itself-since any such interaction is inconsistent with the fact that the ξ anticommute. Such a self-interaction can, however, appear on larger scales, when spatial averages can become meaningful for describing the spin degrees of freedom. In conclusion to this section, the Hamiltonian we have constructed describes the interaction of magnetic moments through their embedding within an elastic medium. What remains to be done is to define the Poisson brackets, that can take into account the evolution in phase space of commuting, as well as anticommuting degrees of freedom. Indeed, the construction procedure for commuting canonical variables is well known, however, including anticommuting variables in a unified way has remained a rather esoteric subject-known in theory [START_REF] Casalbuoni | [END_REF][19][20] , not, however, implemented, in practice.

In the following section, we shall construct the equations of motion in the phase space of the elastic and the spinning degrees of freedom, that implements these ideas in practice.

In order to work directly with the anticommuting variables themselves, in combination with the mechanical degrees of freedom requires implementing a "graded Poisson bracket"; one way to do this is discussed in Appendix A.

III. THE POISSON BRACKETS AND THE EQUATIONS OF MOTION

To deduce from the Hamiltonian, discussed in the previous section, the equations of motion, we must define the appropriate pairs of canonically conjugate variables and consequently their Poisson bracket.

First we recognize that ε acts as a tensor and one can understand the definition of the Poisson bracket of rank-2 symmetric tensors as an application of the DeDonder-Weyl covariant hamiltonian formulation of field theory 39 . Although the context is different, the ADM procedure in general relativity also provides such a Poisson bracket 40 (with further relations, between the conjugate variables that are not relevant here).

Although no clear consensus has, in fact, emerged on the properties of Poisson bracket of rank-2 tensors 41 , if one focuses on the special case of strain tensors, that depend only on time, the following conjugate momentum can be introduced

π IJ ≡ ∂L ∂ εIJ , (9)
where L (ε IJ (t), εIJ (t)) is the unconstrained and free Lagrangian. εIJ (t) are the components of the strain-rate tensor 42 . Thus, we can build the corresponding Hamiltonian H for the time evolution with tensor variables for mechanical deformations ε and their conjugated momenta π as the corresponding Legendre transform

H(ε, π) = π IJ εIJ -L , (10)
up to a total time derivative for L . For the mechanical system only (i.e. for functions A(ε, π) and B(ε, π)), the Poisson bracket can be defined in perfect analogy to that of any particle system, that explores a given target space geometry (to which refer the indices I, J, K, L), by the usual relations

{A, B} P B = ∂A ∂ε IJ ∂B ∂π IJ - ∂A ∂π IJ ∂B ∂ε IJ . (11)
In our case, the dynamical variables are the real symmetric rank-two tensors, ε IJ and π IJ (I, J = 1, 2, 3), which are canonically conjugate in the sense that their Poisson brackets are deemed to satisfy the following properties:

{ε IJ , π KL } P B = δ IJKL , ( 12
) {ε IJ , ε KL } P B = {π IJ , π KL } P B = 0 ( 13
)
where δ IJKL is a δ "tensor", reflecting the real, symmetric, nature of the Poisson brackets 43 and defined as a product of Kronecker δs. Very schematically, we may write

δ IJKL =      δ IJ δ KL δ IJ δ LK δ JI δ LK δ JI δ KL (14)
where each choice of the RHS corresponds to a choice of indices in the Poisson brackets. This choice can be supplemented by any linear combination of these δs that enforces the symmetries of the tensors.

We now wish to include as phase space coordinates, the components of the spin vector, S. We follow reference 44 and the details are summarized in appendix A. The "generalized" Poisson bracket, for the canonical variables of our system, can be written as

{A, B} P B ≡ ∂A ∂ε IJ ∂B ∂π IJ - ∂A ∂π IJ ∂B ∂ε IJ - 1 ǫ IJK S I ∂A ∂S J ∂B ∂S K , (15)
where is introduced to restore the physical dimensions of the Poisson bracket, since S is dimensionless.

It should be stressed that this does not imply that any quantum effects are present, since the dynamics is fully classical. It is simply a bookkeeping device for a quantity that has the dimensions of angular momentumi.e. of an area in phase space-and reflects the fact that the equation of motion for S I is of first order. Quantum fluctuations will introduce the "real" .

Using this Poisson bracket, we can obtain the equations of motion for the phase space variables:

εIJ = {ε IJ , H} P B = ∂H ∂π IJ πIJ = {π IJ , H} P B = - ∂H ∂ε IJ ṠI = {S I , H} P B = 1 ε IJK S J ∂H ∂S K (16)
The consistency of this formalism can be checked by noting that these equations preserve the volume in phase space

∂ εIJ ∂ε IJ + ∂ πIJ ∂π IJ + ∂ ṠK ∂S K = 0. ( 17
)
That the dynamics preserves the volume in phase space does not, of course, imply anything about whether the system thus defined is integrable or shows Hamiltonian chaos.

The internal energy U , involving the mechanical energy for the deformed elastic medium [START_REF] Landau | Theory of Elasticity[END_REF] , the magnetic energy, defined by the Zeeman term [START_REF] Skomski | Simple Models of Magnetism[END_REF] and the magnetoelastic energy [START_REF] Trémolet De Lacheisserie | Magnetostriction: Theory and Applications of Magnetoelasticity[END_REF] , that takes into account the interaction of the magnetic moment with the medium, takes the form

U = V 0 2 C IJKL ε IJ ε KL -V 0 σ ext IJ ε IJ + B IJKL ε IJ S K S L -ω I S I , ( 18
)
where C is the fully symmetric tensor, defining the elastic response, σ ext is the external stress tensor, ω is the effective external magnetic field (expressed as a frequency) and B is the fully symmetric magnetostriction tensor.

The "kinetic" term, containing the conjugate momenta, can be written, schematically, as

H kinetic = 1 2 π IJ M -1 IJKL π KL ( 19
)
where M is a fully symmetric "mass" matrix, that describes the inertial response.

For the case of isotropic materials, it is assuming that the M tensor has the form given by Lamé, with only two characteristic constants:

M IJKL = M 0 δ IJ δ KL + M 1 (δ IK δ JL + δ IL δ JK ) . (20)
The tensors C and B are decomposed in the same way.

Consequently, the inverse of these tensors can then be deduced from

M IJKL M -1 IJMN = 1 2 (δ KM δ LN + δ KN δ LM ). Thus M -1 IJKL = -M 0 2M 1 (3M 0 + 2M 1 ) δ IJ δ KL + 1 4M 1 (δ IK δ JL + δ IL δ JK ) (21)
and the equations of motion become

               ǫIJ = M -1 IJKL π KL πIJ = -V 0 C IJKL ǫ KL + V 0 σ ext IJ -B IJKL S K S L ṠI = ε IJK ω J -2 B ABJC ǫ AB S C S K (22)
highlighting how the mechanical and magnetic subsystems are coupled. The last equation-as expected!-is a precession equation for the components of S around both the effective field ω and an additional field, that depends on the strain tensor and the spin vector.

In the following section we shall show how to solve these equations, in a way that preserves the symmetries of the extended phase space.

@Stam Stopped at 26/05/2019

IV. GEOMETRIC INTEGRATION

Solving the coupled system of equations ( 22) is the next step.

Since, in the previous section, we have shown that these equations describe a volume preserving transformation of the enlarged phase space, encompassing elastic and spin variables, it is natural to rewrite them, in terms of the action of a Liouville operator. Therefore, we shall write eqs. (22) as

ε = L ǫ ε, π = L π π, Ṡ = L S S. ( 23
)
where L is the Liouville operator. This formulation allows us to implement, manifestly, time-reversible, area preserving algorithms, for solving these equations numerically.

The general scheme is as follows: Consider an arbitrary function f of the canonically conjugate variables of our many-body system. This function, f (ε, π, S), depends on the time t implicitly; that is, through the dependence of (ε, π, S) on t. The time derivative of f is ḟ such as

ḟ = εIJ ∂f ∂ε IJ + πIJ ∂f ∂π IJ + ṠI ∂f ∂S I ≡ Lf. ( 24
)
The last line defines the (total) Liouville operator

L = εIJ ∂ ∂ε IJ + πIJ ∂ ∂π IJ + ṠI ∂ ∂S I . (25)
Equation ( 24) can be integrated formally as an initial value problem to obtain f at any time :

f (ε(t), π(t), S(t)) = e Lt f (ε(0), π(0), S(0)). (26)
It is not difficult to see that L = L ε + L π + L S . However, these single Liouville operators do not commute two-by-two, as the reader may easily check by computing

L u L v f -L v L u f =
According to the Magnus expansion 46 , however, it is always possible to express e Lt as a product of the individual operators at any given order in time, according to the so-called "splitting method" [START_REF] Hairer | Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] . This ensures that the numerical algorithm preserves phase space volumes exactly.

For instance, for a fixed timestep τ , upon expanding up to the third order in time, the following sequence of products

e Lτ = e LS τ 4 e Lπ τ 2 e LS τ 4 e Lετ e LS τ 4 e Lπ τ 2 e LS τ 4 + O(τ 3 ), (28)
can be generated; this sequence is, in fact, one of six possible, that have the property of preserving the symplectic structure of the Poisson brackets. Therefore, any one of them can be chosen. The possible combinations are presented in table I. While these schemes are free from

A τ 4 τ 2 τ 4 τ τ 4 τ 2 τ 4 1 S π S ǫ S π S 2 S ǫ S π S ǫ S 3 π ǫ π S π ǫ π 4 π S π ǫ π S π 5 ǫ π ǫ S ǫ π ǫ 6 ǫ S ǫ π ǫ S ǫ

Table I. Decomposition table of symplectic integrators

"global" errors, they are, of course, sensitive to "local" errors, due to the finite value of the timestep. It is, also, not at all obvious that all six can be implemented with comparable efficiency. It is, therefore, useful to study the numerical stability and efficiency of these different combinations, in particular, as former studies in molecular dynamics 48 and magnetic molecular dynamics 12 showed, apparently, numerical differences between them.

A sampler of such a study is presented in appendix B.

It is important to keep in mind that the one-step evolution operators for ε and π describe shifts of the corresponding tensor components, whereas the one-step evolution operator for S describes rotations. In equations e Lετ (ε(t), π(t), S(t)) = (ε(t) + τ ε(t), π(t), S(t)) (29) e Lπτ (ε(t), π(t), S(0)) = (ε(t), π(t) + τ π(t), S(t)) (30)

e LS τ (ε(t), π(t), S(t)) = (ε(t), π(t), S(t + τ ) R(τ )S(t) ) ( 31
)
where S(t + τ ) = R(τ )S(t) is given by the Rodrigues' rotation formula 49 for a spin vector around a given rotation vector ω(t) where each of its components are ωI (t) = ω I (t) -2 B JKLI ε JK (t)S L (t). These equations describe the phase space of one particle only. To describe the dynamics of a continuum, we must deduce the equations for many particles. One way to generalize eqs. (22) for the case of many particles, labeled by an index, i = 1, . . . , N , according to the conventions of the previous sections, is the following

ǫi IJ = [M i ] -1 IJKL π i KL ( 32
) πi IJ = -V 0 C i IJKL ǫ i KL + V 0 σ i ext IJ -B i IJKL S i K S i L ( 33
) Ṡi I = ε IJK ω i J - 2 B i ABJC ǫ i AB S i C S i K ( 34
)
The case of a staggered AF state is treated simply by letting N = 2. In order to simplify the mechanical part further, we can impose additional conditions pertaining to the uniformity of the external stress, mechanical constants, mass matrices and magneto-elastic constants, at each site. In what follows, we shall use the following Ansatz:

B 1 IJKL = B 2 IJKL , C 1 IJKL = C 2 IJKL , M 1 IJKL = M 2 IJKL , σ 1 ext IJ = σ 2 ext IJ . (35)
The conservative part of the precession contains a local field, which is modified to include the antiferromagnetic exchange between the sites and a single anisotropy axis n with an intensity K. One has

ω i I = 1 <ij> J ij S j I + K n J s i J n I (36)
Because of the exchange field, the Liouville operators for different spins do not commute either. A global geometric integrator, implementing the approach of Omelyan 50 , which remains accurate up to third order in the timestep expansion must, therefore, be constructed.

For any given timestep, τ , the corresponding expression for the evolution operator, reads

e LS τ = e L S 1 τ 2 e L S 2 τ e L S 1 τ 2 + O(τ 3 ). ( 37
)
For N = 2, this operator is numerically identical to the operator, obtained by permuting the site indices, 1 ↔ 2. The same reasoning is applied for the Liouville operators for the elastic variables, that enter in equations (23) and the corresponding global geometric integrators are constructed along the same lines. The system is then integrated by following one of the schemes displayed in table I.

With these tools, we can study a plethora of phenomena, that are sensitive to the coupling of magnetic, electric and elastic degrees of freedom. In the following section we shall apply this formalism for studying "switching" effects of the magnetization in antiferromagnetic media. The numerical accuracy ensured by the geometric integrators is necessary to describe picosecond switching times.

V. ANTIFERROMAGNETIC ULTRAFAST SWITCHING UNDER STRESS

In this section, we shall apply the formalism constructed above, to describe how it is possible to generate and manipulate picosecond switching of the magnetization, induced by STTs, from a short pulse of electric current, in elastic media that exhibit staggered AF order. Such a fast switching process in AFs is schematically displayed in figure 1 and can be realized, in practice, by a femtosecond laser excitation of the magnetic moments that generate a field-induced STT on the two sublattices. During the pulse, energy is transferred from conduction electrons to the sub-lattice magnetic moments via STT (this is how the electric current affects the magnetic response), and hence contributes to the exchange energy between the moments, since local moments can be canted noncollinearly. The energy due to the strong ex-change interaction between neighboring moments, which is commonly found in AFs, in particular due to magnetic anisotropy, appears as an effective inertia to the motion, leading to the appearance of a time scale, much longer than that of the excitation pulse. Afterwards, the system follows a natural path along the easy-plane to circumvent the unfavorable anisotropy barrier and finally relaxes to a new magnetic configuration. The resulting dynamics is a switching of the two sub-lattices moments to the opposite direction through the easy-plane.

This process defines the ultrafast antiferromagnetic switching effect.

During the process, the Néel vector, that probes the difference between the magnetic moments of the two sublattices, acquires a net value that can be transferred as a so-called "spin accumulation", to an adjacent nonmagnetic normal material in order to pump the produced spin current via scattering of electrons. The reverse mechanism can be realized, as well.

What was not considered before is the possibility to enhance or inhibit such a switching, depending on tensile or compression effects, that are generated by an external stress, that couples to the internal strain, produced by the intrinsic magneto-elastic interaction in such materials, as depicted in figure 1.

In order to conform to the notation used in our previous studies 22 , equation ( 34) is supplemented with a nonconservative part, labeled T , on the RHS, which includes both a transverse damping and a field-like STT torque

T i I = αǫ IJK S i J Ṡi K + G s i I s i J p J -p I s i J s i J . (38)
The resulting equations of motion, are numerically integrated, using the approach of section V. It is noteworthy that, this equation, also, describes the "backreaction" of the magnetic response on the spin transfer torque. On the other hand, we do not consider how this backreaction affects the current pulse itself, that's assumed external.

In figures 2 and 3 we report the evolution of the average magnetization m ≡ 1 2 S 1 + S 2 and the Néel vector l ≡ 1 2 S 1 -S 2 , in the presence as well as the absence of magnetoelastic coupling, for a moderate external isotropic stress of 30µ 0 M 2 s . The conditions and values of the simulation are identical to those found by Cheng et al. 21 . We start the simulations using an initial configuration, where spins are aligned along the x-axis in an antiferromagnetic configuration and apply two 10 ps electric pulses of 0.0034 rad.THz intensity, each separated by 50 ps, in the p = ẑ-direction. In addition to the exchange interaction, the spins are subjected to a global anisotropy, along the n = x-axis.

In the absence of magnetoelastic coupling, our results are identical to those by Cheng et al. 21 and to those we obtained under the same conditions in our previous work 22 .

In the presence of magnetoelastic coupling, however, as the mechanical system is undamped, the results are quite different. In figure 4 we display the evolution of the strain components over time. Because of the presence of a constant stress, a finite mass matrix and non-zero elastic constants, the mechanical system is expected to be oscillating freely, which is, indeed, observed.

These oscillations are also observed, both in the magnetization and the Néel vector in the figure 3, in the form of small wiggles, superimposed on the main component of the spin accumulation vector. A compression in one direction becomes tensile in the other one and vice-versa because of the sign of B 1 . As a consequence of such a mechanical deformation, the switching time of the magnetization is slightly faster and the amplitude of the magnetization enhanced, depending on the value of the stress and magnetoelastic constants, as observed in figure 3. On the other hand, the strain tensor itself is quasi-uniform in space on the two sub-lattices.

VI. DISCUSSION

In the previous sections, we have described how to realize a mechanism for the ultrafast switching of the magnetization, that can be induced by a STT in NiO, where both the duration and intensity can be modulated by appropriate internal strain, thanks to magnetoelastic coupling. Such a strain can be modulated, in return, by the induced magnetization process by the same magnetoelastic coupling in order to produce a modulated elastic wave. Depending on the pulse duration and transverse damping value, the terminal angle of the Néel vector evolves gradually by integer values of π as already shown 21 . In real experiments, because amplitude fluctuations of the electric pulses are inevitable, small damping values are not favorable since it may easily lead to overshoot, which is amplified by the mechanical response. When Joule heating in the normal underlying metal is taken into account, a shorter pulse with stronger current intensity should be desirable, but a modulated undamped mechanical strain is then produced. For the moment, what is not easy to predict are the mechanical damping mechanisms induced by a differential thermal stress produced by such a Joule heating, where its effect superimposes to the whole dynamics. However, at least for the case of undamped dynamics, we have developed the mathematical and computational framework, that consistently takes into account elastic and magnetic degrees of freedom and that pro-vides a starting point for designing future AFs devices for spintronics applications that combine both mechanical and magnetic coupled responses.

S I = - ı 2 ǫ IJK ξ J ξ K , (A4)
which can easily be shown to commute, as S I S J = 0. This last relation is, of course, not satisfactory. It can be easily avoided, however, by imposing that the ξ I satisfy the anticommutation relations {ξ I , ξ J } = δ ijthat they generate a Clifford algebra (up to a constant normalization) [START_REF] Lounesto | Clifford Algebras and Spinors[END_REF] . These relations imply that the S I define a spin-1/2 representation of the su(2) algebra. Higher spin representations can be defined by multilinear combinations 20 , that are relevant for describing magnetic properties of composite objects; since we can work with the S I instead of the ξ I , however, this complication will not affect us here.

This representation highlights that the spin degrees of freedom, S I , are, in fact, "composite" objects and that the "fundamental degrees of freedom" are the ξ I . Therefore, it is useful to develop the description of the dynamics, directly, in terms of the ξ I themselves. We shall recall the salient features below.

The Poisson brackets for the anticommuting variables ξ I are related to those of S I , in order that the dynamics be, indeed, equivalent, in a way that was set forth many years ago, through the construction of a corresponding graded Poisson bracket, which generalizes Poisson brackets from manifolds to super-manifolds [START_REF] Casalbuoni | [END_REF] .

In terms of any functions of anti-commuting variables ξ this bracket reads

{f (ξ), g(ξ)} PB ≡ ı f (ξ) ← - ∂ ∂ξ K -→ ∂ ∂ξ K g(ξ), (A5)
with the corresponding definition of the left and right derivative of any function of the anti-commuting variables. One can show that this bracket, also known as the "antibracket" of any two functions on a flat supermanifold, satisfies all necessary properties for a graded bracket, namely (graded) Leibniz rule, (graded) antisymmetry and (graded) Jacobi identity 53 . By taking this graded Poisson bracket for any two of these anti-commuting variables, we get {ξ I , ξ J } P B = ı δ IJ .

(A6)

This implies that any two such variables are canonically conjugate, since {ξ I , -ı ξ J } P B = δ IJ and π I ≡ -ı ξ I defines the canonical conjugate 19,20 . By using the Grassmann properties of ξ, one proves that {S I , S J } P B = 1 ǫ IJK S K (A7) which is a consistency check that the S I as defined in (A4) do realize a representation of the rotation group. In eq.(A5), if the functions of the ξ are chosen to contain only quadratic terms, then one can identify the previous bracket as a regular Poisson bracket on a Riemannian manifold for the commuting variables S {f (S), g(S)} π 2 11 (0)/4M 1 , which has to stay constant over time. We observe that this is, indeed, the case, whatever the value of the quality factor Q. We remark that, as Q → 0, the variations about the average value of the energy are suppressed, as should be expected.

Figure 6 displays some of the non-vanishing components of the strain tensor over time (here ε 22 (t) = ε 33 (t) and this third component is not reported), when using different splitting algorithms, among those displayed in Table I. We observe less than 1% of numerical relative difference between the two algorithms on the strain and its conjugate variables, and no difference (up to the machine precision) on the magnetization with a "coarse" quality factor Q, which cannot be detected by visual inspection on the figure. This difference falls to 0.1% when the quality factor is divided by 4. The same procedure can be repeated for all the splitting combinations in Table I, the conclusions previously drawn apply, also, for the magnetization, strain and its conjugate variable, depending on the frequency of appearance of the splitted operator. As expected, once a splitting algorithm is selected, the more frequently an operator is evaluated, the smaller the error, though, of course, the time required, also, grows. This opens the possibility to select optimally a proper splitting algorithm, depending on the relative intensity of εij , πij , Ṡi over time.

  0 for any function f and any combination (u, v) of the individual Liouville operator L ε , L π , L S . This means that e Lt = e Lεt+Lπt+LS t = e Lεt e Lπt e LS t .

Figure 1 .

 1 Figure 1. (Color online) Switching scheme for antiferromagnetic magneto-mechanically coupled toy model with external STT.

Figure 2 .

 2 Figure 2. (Color online) Average magnetization (upper panel) and Néel order parameter components (middle panel) for uncoupled switching with parameters in reference 22 . The lower panel displays the STT pulses. The figures agree with the reference 21 because the magnetoelastic constants are set to zero.

Figure 3 .

 3 Figure 3. (Color online) Average magnetization (upper panel) and Néel order parameter components (middle panel) for a coupled switching with parameters in reference 22 . The lower panel displays the STT pulses. Here B0 = 7.7µ0M 2 s and B1 = -23µ0M 2 s .

Figure 4 .

 4 Figure 4. (Color online) Strain components as a function of time with the magnetoelastic constants turned on with parameters in reference 22 .

Figure 5 .

 5 Figure 5. (Color online) Single-site total energy as a function of both time and variable timestep Q. Conditions of the simulation are expressed in reduced units: 2C0/µ0M 2 s = 5.1 × 10 5 , 2C1/µ0M 2 s = 3.5 × 10 5 , M1V0µ0M 2 s /2 2 = 1000, ωDC = (0, 0, 2π) rad.GHz, π11(0) = 1, s(0) = (1, 0, 0). All the other parameters not reported, included initial conditions are zero.

Figure 6 .

 6 Figure 6. (Color online)Single-site strain components as a function of time for various numerical schemes. Conditions of the simulation are identical to those for figure 5 and the results are produced for Q = 0.0025 only.

ACKNOWLEDGMENTS

All authors contributed equally to the manuscript. T.N. acknowledges financial support through a joint doctoral fellowship "CEA-Région Centre" under grant number 201600110872.

Appendix A: A model for precession through anticommuting variables

We review some issues when dealing with classical spin variables and explain some limitations which are met by using the framework of commuting variables 51,52 . A common starting point to build precessional models is choosing the appropriate representation of the su(2) algebra.

In the Heisenberg picture the Ehrenfest theorem implies that the equations of motion for the average of the spin operator Ŝ components read

where H is the Zeeman Hamiltonian operator for an external magnetic field H such that

From the commutation relations [ Ŝi , Ŝj ] = ıǫ ijk Ŝk , one can quite straightforwardly show that

which corresponds to the well-known Larmor precession of the expectation value of the spin in an external magnetic field H. Interestingly, a different-but equivalent-approach can be followed, by considering the Majorana representation, in terms of anticommuting variables, ξ I , ξ I ξ J + ξ J ξ I = 0 31 of the vector S:

Similar expressions have been found already by Casalbuoni a long time ago 19 without having attracted the attention they deserve. Following the inverse path of canonical quantization, we use this graded Poisson Bracket on any commuting quantity, which allows us to compute directly

{S, H}

thereby highlighting that the description of S in terms of ξ is an equivalent description of the dynamics.

For any time-dependent commuting functions of S(t), we can use this Poisson bracket to deduce a Liouville equation

Consequently, the equations of motion for the variables ξ read much more simple expressions as

which are known to form a non-relativistic pseudoclassical mechanics 19 . One reason why using the representation of S in terms of ξ is useful is that it is intrinsically difficult to build a Lagrangian model for the commuting spin variable S, since its canonically conjugate variable cannot be unambiguously identified. As is, by now, well known, the conjugate of the dynamical variable ξ is proportional to itself 54 . Therefore, the dynamics of spinning degrees of freedom can be described, either through a vector of commuting variables on a curved manifold, or by a vector of anti-commuting variables on a flat-though non-Riemannian-manifold. Finally, it has been found that the Majorana-fermion representation of 1/2-spin operators is also a powerful tool to straightforwardly compute spin-spin correlators 55 , which represents an advantage for computing magnetic response functions in a many-bodies formulation.

Appendix B: Numerical accuracy of the splitting algorithms

The accuracy of the numerical schemes represented in table I depends on the relative amplitude of the velocity terms ( ε, π, Ṡ), and can be monitored by checking the stability of equation ( 17) over time 46,[START_REF] Hairer | Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] .

Figure 5 displays the total energy in one domain (given by the sum of eqs.( 18) and ( 19)) upon varying the numerical precision. The splitting algorithm considered corresponds to the label A = 1 in table I. The numerical precision is controlled by the "quality factor Q", chosen at the beginning of the simulation, which produces variable timesteps τ according to the relation τ = Q/ ω . Simulation conditions produce a total energy equal to