I Auffret

L Joumard

R C Vila

L D Sousa

I L Buda-Prejbeanu

Prejbeanu

M Mansueto

A Chavent

S Auffret

I Joumard

L Vila

R C Sousa

L D Buda-Prejbeanu

I L Prejbeanu

B Dieny

Spintronic memristor for neuromorphic circuits based on the angular variation of tunnel magnetoresistance

Keywords: Magnetic tunnel junction, tunnel margnetoresistance, memristor, exchange bias, neuromorphic, dry friction

Spintronic memristors for neuromorphic circuits based on the angular variation of tunnel magnetoresistance

Introduction

The exponential increase of the amount of data produced worldwide in the last decade pushes our technologies toward their limits. The total associated power consumption steadily increases every year. In this context, bio-inspired technologies seem to provide a promising solution to face these trends. In particular, neuromorphic computing can be used for fast and energy-efficient applications through the implementation of networks of artificial neurons and synapses. Some concepts of such spintronic components were proposed in the last few years, exploiting different physical mechanisms to achieve memristor functionality. Memristors are non-volatile memory devices mimicking artificial synapses, whose resistance can take several intermediate values between minimum and maximum, the resistance varying monotonously as a function of the polarity of the current pulses flowing through the device. Memristor concepts also exist based on resistive oxides [START_REF] Kim | Comprehensive Physical Model of Dynamic Resistive Switching in an Oxide Memristor[END_REF] or ferroelectric materials [START_REF] Chanthbouala | A ferroelectric memristor[END_REF] . However, based on the knowledge acquired during MRAM development, the expected advantages of spintronic memristors are a better write endurance and a lower cell to cell variability than for their RRAM or ferroelectric-based counterparts. Their drawback is the lower ROFF/RON ratio which should be at least partially compensated by their lower cell to cell variability. A first proposed concept of spintronic memristor is based on the controlled domain wall (DW) nucleation and propagation in the free layer of a magnetic tunnel junction (MTJ) [START_REF] Lequeux | A magnetic synapse: multilevel spin-torque memristor with perpendicular anisotropy[END_REF][START_REF] Kim | Multilevel anomalous Hall resistance in a single Hall cross for the applications of neuromorphic device[END_REF][START_REF] Azam | Voltage control of domain walls in magnetic nanowires for energy-efficient neuromorphic devices[END_REF] . The pinning and depinning of the DW at different positions in MTJ track define the intermediate resistance states. The main drawback is the reliability of the intermediate states. Moreover, the number of pinning sites, directly related to the number of intermediate resistance states, is linked to the size of the device, making it hardly scalable. Other propositions rely on the use of several binary MTJs connected in series or in parallel configurations 6, [START_REF] Ostwal | Novel Compound Synapse Using Probabilistic Spin-Orbit-Torque Switching for MTJ-Based Deep Neural Networks[END_REF] . The probabilistic switching of each individual MTJ due to Spin Transfer Torque (STT) under pulses of current of properly chosen duration and amplitude allows to realize the memristor function. In this case, the number of resistance levels strictly depends on the number of interconnected MTJs. This limits again the scalability despite the good reliability of the intermediate states. A third proposal exploits the spin orbit torque (SOT) in a SOT-MRAM cell modified to include an antiferromagnetic material constituting the heavy metal line, exchange biased with a perpendicularly magnetized storage layer 8, [START_REF] Kurenkov | Artificial Neuron and Synapse Realized in an Antiferromagnet/Ferromagnet Heterostructure Using Dynamics of Spin-Orbit Torque Switching[END_REF] . A gradual switching occurs in µm-sized dots under pulses of current flowing through the heavy metal line. Besides, spike-timingdependent plasticity (STDP) is observed in this device. By reducing the size of the device, the number of intermediate states reduces up to a binary behavior at few hundreds nm of diameter. Commonly to all those devices, the size is a limiting parameter because directly linked to the number of intermediate states.

In the present study, we realized a nano-sized memristive device whose functionalities do not depend directly on the lateral size. This, together with the other advantages of the MRAM technology, allows the creation of a scalable, reliable and easily integrable device. The Please do not adjust margins Please do not adjust margins idea is to exploit the transport properties of an in-plane magnetized tunnel junction. It is known that the conductance of such device depends on the relative orientation of the magnetizations of the two magnetic electrodes (analyzer and free layer) [START_REF] Slonczewski | Conductance and exchange coupling of two ferromagnets separated by a tunneling barrier[END_REF][START_REF] Jaffrès | Angular dependence of the tunnel magnetoresistance in transition-metal-based junctions[END_REF] . It varies as the cosine of the angle between the two vectors. The stabilization of the magnetization of the free layer along different in-plane directions yields the possibility to obtain intermediate resistance values between the lowest and the highest (respectively for parallel and antiparallel alignment of the two magnetizations). Concerning the writing mechanism, in order to discretely rotate the magnetization in the plane, the spin transfer torque (STT) coming from an out-of-plane polarizing layer can be exploited through the application of properly designed pulses [START_REF] Mansueto | Realizing an Isotropically Coercive Magnetic Layer for Memristive Applications by Analogy to Dry Friction[END_REF][START_REF] Mansueto | Isotropically coercive free layer integration in a magnetic tunnel junction for neuromorphic applications[END_REF] . In a previous study, we have shown how to realize a free layer able to stabilize the magnetization along different in-plane directions [START_REF] Mansueto | Realizing an Isotropically Coercive Magnetic Layer for Memristive Applications by Analogy to Dry Friction[END_REF][START_REF] Mansueto | Isotropically coercive free layer integration in a magnetic tunnel junction for neuromorphic applications[END_REF] . Such isotropic properties can be obtained in a trilayer structure composed of two ferromagnets sandwiching a thin antiferromagnetic material. The thickness of the antiferromagnetic layer is adjusted to provide enhanced coercivity but no exchange bias at the temperature of operation (typically around 2nm for IrMn at 300K) [START_REF] Mansueto | Realizing an Isotropically Coercive Magnetic Layer for Memristive Applications by Analogy to Dry Friction[END_REF] . In such full sheet film sample, the coercive field and remanence are independent on the in-plane direction in which the magnetic field is applied. Those properties can be phenomenologically modeled with the Landau-Lifshitz-Gilbert (LLG) equation including a dry friction-like term. The addition of this dissipative term led to a good agreement with the experimental results. Moreover, the integration of such layer in a nano-sized pillar preserves its isotropic properties [START_REF] Mansueto | Isotropically coercive free layer integration in a magnetic tunnel junction for neuromorphic applications[END_REF] as confirmed by rotating field experiment and measurement of the tunnel magnetoresistance (TMR) variation versus in-plane angle. To achieve a memristor behavior which requires a monotonous variation of the resistance versus current pulse polarity, it is important to limit the angular excursion of the free layer magnetization between 0° and 180° corresponding respectively to the parallel and antiparallel magnetic configuration. In our earlier study [START_REF] Mansueto | Realizing an Isotropically Coercive Magnetic Layer for Memristive Applications by Analogy to Dry Friction[END_REF] , we showed that this could be achieved by the application of a static transverse field of suitable amplitude acting on the free layer magnetization. This field is applied in-plane perpendicular to the reference layer magnetization. However, applying such field is not very convenient in practical devices which motivates for instance the search for field free implementation in the context of SOT-MRAM [START_REF] Garello | Manufacturable 300mm platform solution for Field-Free Switching SOT-MRAM[END_REF] . To circumvent the difficulties associated with application of a static field, we show in this study that, by properly adjusting the two STT relative contributions (from the perpendicular polarizer and the inplane reference layer), it is possible to gradually rotate the magnetization of the free layer clockwise or counterclockwise depending on the pulse polarity, within a 0°-180° excursion, thus achieving memristor functionality. In the framework of our model, we distinguish different working regions for our device as a function of the two STT contributions. Guided by the simulations, we realize a nano-sized device exhibiting the expected memristive properties.

Model

The proposed concept is actually an extension of a previously studied memory device investigated for the realization of sub-ns binary switching thanks to the phenomenon of precessional switching induced by a perpendicular polarizer on an in-plane magnetized free layer 15-20 . This device was sometime called Orthogonal Spin Transfer Torque MRAM (OST-MRAM) [START_REF] Liu | Ultrafast switching in magnetic tunnel junction based orthogonal spin transfer devices[END_REF] . In the present study, the original device is modified by introducing an isotropically coercive free layer instead of a bistable free layer with uniaxial anisotropy. In fact, the idea is to exploit the STT coming from the perpendicular polarizer to pull the magnetization of the free layer out-of-plane and make it rotate around the demagnetizing field. By properly adjusting the pulse amplitude and duration, the magnetization can be switched by 180°. Besides, it was shown that by playing with the relative contributions of the two STT (from perpendicular polarizer and from planar reference layer), one can achieve 180° switching independently of the pulse duration which makes the concept much easier to implement in practical devices 17, [START_REF] Lacoste | Modulating spin transfer torque switching dynamics with two orthogonal spin-polarizers by varying the cell aspect ratio[END_REF] . We will show that the same is true for the present memristor device.

Our model based on a modified LLG equation with the addition of a dissipative dry friction term [START_REF] Mansueto | Realizing an Isotropically Coercive Magnetic Layer for Memristive Applications by Analogy to Dry Friction[END_REF] is used here to investigate the free layer magnetization (M) dynamics in such device:

𝜕𝑴 𝜕𝑡 = -𝛾 (𝑴 × µ 0 𝑯 𝒆𝒇𝒇 ) + 𝛼 𝐺 𝑀 𝑆 ( 𝑴 × 𝜕𝑴 𝜕𝑡 ) + 𝛽𝑀 𝑆 𝑴× 𝜕𝑴 𝜕𝑡 |𝑴× 𝜕𝑴 𝜕𝑡 | (1)
where 𝛾 is the gyromagnetic ratio, 𝛼 𝐺 the Gilbert damping, 𝑀 𝑆 the saturation magnetization and 𝛽 the dry friction parameter. In this equation, the first term describes the magnetization precession around the effective field, the second term represents the Gilbert damping dissipation and the third term the dry friction term associated with the exchange like coupling with the thin antiferromagnetic layer [START_REF] Mansueto | Realizing an Isotropically Coercive Magnetic Layer for Memristive Applications by Analogy to Dry Friction[END_REF][START_REF] Mansueto | Isotropically coercive free layer integration in a magnetic tunnel junction for neuromorphic applications[END_REF] . From the magnetization dynamics, the resistance variation can be straightforwardly derived from the expression of the conductance of a magnetic tunnel junction [START_REF] Slonczewski | Conductance and exchange coupling of two ferromagnets separated by a tunneling barrier[END_REF] :

𝐺 = 1 𝑅 = 𝐺 90°+ ∆𝐺 2 𝑴. 𝑴 𝑹 𝑀 𝑠 .𝑀 𝑅𝑠 (2)
Where 𝐺 90° is the junction conductance in the configuration when the free layer magnetization is orthogonal to that of the reference layer (also called analyzer), ∆𝐺 is the absolute magnetoconductance defined as the variation of conductance between the parallel magnetic configuration (free layer magnetization parallel to that of the reference layer, high conductance state) and antiparallel configuration (low conductance state). ∆𝐺 is related to the usual (1+𝑇𝑀𝑅) . In (2), 𝑴 𝑹 is the reference layer vector magnetization and 𝑀 𝑅𝑠 its amplitude.

We first analyze the contribution of torque coming from the in-plane analyzer (i.e. reference layer) by including in the equation the torque 𝑻 𝒂𝒏 = 𝛾𝑎 𝑎𝑛// 𝑉𝑴 × (𝑴 × 𝑳 ) where 𝐿 = (1 0 0), 𝑉 is the applied voltage and 𝒂 𝒂𝒏// is the Slonczewski prefactor. This last can be expressed as

𝑎 𝑎𝑛// = ℏ 𝟐𝒆 𝜼 𝒂𝒏 𝒍 𝑴 𝒔 𝑹×𝑨
, where 𝑙 is the thickness of the free layer, 𝑅 × 𝐴 the resistance x area product and 𝜂 𝑎𝑛 the spin torque efficiency. All the following results are obtained for a layer of cylindrical shape, 100nm in diameter, and at T = 0K. Looking for in-Please do not adjust margins Please do not adjust margins plane (Mz = 0) stationary solutions ( 𝜕𝑴 𝜕𝑡 = 𝟎), we find that the threshold current density can be expressed as

𝐽 𝑎𝑛 = 𝛽 𝑀 𝑆 𝑙 2 𝑒 𝛾 ℏ 𝜂 𝑎𝑛 (3)
In a similar way, considering the STT from the perpendicular polarizer only: 𝑻 𝒑𝒆𝒓𝒑 = 𝛾𝑎 𝑝𝑒𝑟𝑝// 𝑉𝑴 × (𝑴 × 𝑷 ) with P= (0 0 1), it is possible to compute the critical current required to trigger the outof-plane precession (OPP) :

𝐽 𝑝𝑒𝑟𝑝 = 𝛽 𝑀 𝑆 𝑙 2 𝑒 𝛾 ℏ 𝜂 𝑝𝑒𝑟𝑝 (4)
If now both torques are inserted in the equation, the critical current density for which in-plane static solutions become unstable writes:

𝐽 𝑑𝑜𝑢𝑏𝑙𝑒 = 𝛽 𝑀 𝑆 𝑙 2 𝑒 𝛾 ℏ 1 √𝜂𝑝𝑒𝑟𝑝 2 +𝜂 𝑎𝑛 2 (5)
In Fig. 1, the three critical lines are drawn with respect to the spin torque efficiency of the perpendicular polarizer 𝜂 𝑝𝑒𝑟𝑝 for a fixed value of 𝜂 𝑎𝑛 = 0.3 (obtained from the measured TMR value of the device with the formula 𝜂 𝑎𝑛 = √𝑇𝑀𝑅(𝑇𝑀𝑅+2) 2(𝑇𝑀𝑅+1) [START_REF] Sun | Spin-torque switching efficiency in CoFeB-MgO based tunnel junctions[END_REF] ).

Considering a generic initial configuration of the magnetization of the free layer, the dotted line 𝐽 𝑎𝑛 represents the limit of stability of such state. Below this threshold, there is the stability region where the magnetization is at equilibrium whatever its in-plane direction, while, above it, the magnetization becomes unstable. The dashed line 𝐽 𝑝𝑒𝑟𝑝 is the critical line for the excitation of out-of-plane precession of the magnetization (OPP). Above this line only OPP can occur. The green dashed dotted line 𝐽 𝑑𝑜𝑢𝑏𝑙𝑒 represents the correction to the stability line 𝐽 𝑎𝑛 due to the other STT contribution originating from the perpendicular polarizer. As for the precessional switching mechanism, a magnetic stack comprising two similar tunnel junctions would be ideal to fully exploit the OPP region. In this case, close to the top part of Fig. 1, the OPP would start as soon as the magnetization is no more in the stability region because of the dominant effect of the STT contribution from the perpendicular polarizer. However, as we will see further, this would not allow to keep the angular excursion of the free layer magnetization between 0° and 180° and therefore to achieve the required monotonous variation of resistance for a given pulse polarity.

In the unstablity region, where the two STT contributions are competing, it is not possible to exactly define analytically the magnetization dynamics. For values of 𝜂 𝑝𝑒𝑟𝑝 close to zero, the STT torque from the analyzer is dominating, leading to a binary switching of the free layer. In this case, the addition of the dry friction term blocks the stochastic behavior at the beginning of the switching mechanism, making the magnetization switching almost impossible at low current densities (blue line in Fig. 2a). By increasing the torque originating from the perpendicular polarizer, the sum of the two contributions leads to a binary switching but in a time significantly shorter than in the previous case (orange line in Fig. 2a). This is due to the initial torque impulse induced by the perpendicular polarizer which pulls the free layer magnetization away from the analyzer magnetization direction. In this region, the torque 𝑇 𝑎𝑛 favors a switching towards parallel or antiparallel configuration through a precessional trajectory around the demagnetizing field under the influence of 𝑇 𝑝𝑒𝑟𝑝 . For higher values of 𝜂 𝑝𝑒𝑟𝑝 , the magnetization starts to precess because of the dominant influence of the perpendicular STT with a trajectory distorted by the STT term from the in-plane analyzer (green line in Fig. 2a). For the final purpose of realizing a memristive behavior, the two regions of OPP and unstability could be exploited. The first case has been already simulated in ref [START_REF] Mansueto | Realizing an Isotropically Coercive Magnetic Layer for Memristive Applications by Analogy to Dry Friction[END_REF] leading to a memristive behavior with pulses in the ps-range. Despite this, the monotonicity of the resistance variation with the polarity of the voltage can be achieved only with the application of an in-plane transverse field which limits the angular excursion of the free layer magnetization only in a half plane (between 0° and 180°). Moreover, because of the very short pulse width, due to the relatively high frequency of oscillations, the write margin in terms of pulse duration and amplitude is very small (as for the precessional switching case). For these reasons, working in the unstability region and adjusting the relative contributions of the two torques to limit the angular excursion of the free layer magnetization in a half plane seems to be a better option for this application. In particular, the fast binary switching, as the orange line in Fig. 2a, represents an ideal switching mechanism for our case. Indeed, it is then possible to gradually rotate the magnetization step by step under pulse of currents (See Fig. 2b) while limiting the magnetization excursion between the P and AP states. This ensures the monotonicity of the resistance variation for a given polarity of the voltage even without the application of an in-plane field. Because of those advantages, we chose to fabricate a device designed to operate in the unstability region of Fig. 1.

Magnetic stack

To exploit the effects of the combined STT in the region described above with a torque efficiency from the perpendicular polarizer being about one sixth of the torque efficiency from the in-plane analyzer (corresponding to the orange curve in Fig. 2a), the isotropically coercive free layer is sandwiched between a Cu metallic spacer separating it from the perpendicular polarizer and a tunnel barrier separating it from the in-plane analyzer (i.e. reference layer) as depicted in Fig. 3. Indeed, it is known that the spin-transfer torque efficiency in metallic spin-valves is significantly lower than in magnetic tunnel junctions [START_REF] Liu | Ultrafast switching in magnetic tunnel junction based orthogonal spin transfer devices[END_REF] . This allows the device to operate in the unstablity region of the diagram in Fig. 1 with 𝜂 𝑝𝑒𝑟𝑝 significantly lower than 𝜂 𝑎𝑛 and, at the same time, have a reading mechanism based on the TMR of the top junction.

In the magnetic stack, the polarizer and the analyzer are implemented as synthetic antiferromagnets with an RKKY antiferromagnetic coupling through a thin Ru spacer 22-24 . In the polarizer, the perpendicular magnetization is induced by the (Co/Pt) multilayers of composition [Co(0.5)/Pt(0.25)]6 and [Co(0.5)/Pt(0.25)]3 (in parenthesis, the thicknesses are in nm) [START_REF] Zeper | Perpendicular magnetic anisotropy and magneto-optical Kerr effect of vapor-deposited Co/Pt-layered structures[END_REF] . The reference is a Co(1) layer followed by the spacer Cu(3) with a lamination of Cu(0.2) to enhance the spin-polarization thanks to a shortening of the spin-diffusion length [START_REF] Houssameddine | Spin-torque oscillator using a perpendicular polarizer and a planar free layer[END_REF] . The isotropic free layer Py(1)/IrMn(2.1)/Py (1) is inserted between a layer of Co(0.6) and a layer of FeCoB(1.6) to optimize the STT coming from both sides and the TMR of the top junction. This last is based on an MgO(1.2) tunnel barrier with an RxA around 6Ωµm [START_REF] Chanthbouala | A ferroelectric memristor[END_REF] . As in conventional STT-MRAM, the MgO barrier is surrounded by two FeCoB electrodes themselves separated from the rest of the stack by two thin Ta layers. These Ta laminations aim at absorbing the B from the FeCoB electrodes upon annealing of the structure and insuring a structural transition between the FeCoB/MgO/FeCoB body centered cubic (bcc) part of the stack and the bottom and top parts which have a face centered cubic (fcc) structure [START_REF] Apalkov | Magnetoresistive Random AccessMemory[END_REF] . The top SAF is based on a Co(1) and Co(2.5)/IrMn (7) structure with a FeCoB(1.2) in contact with the MgO. The overall stack may look complex. However, its structure as well as the materials used in the different layers comprised in the stack are very similar to those used in STT-MRAM cells as well as in magnetoresistive heads implemented in hard disk drives [START_REF] Fullerton | Spintronics, Magnetoresistive Heads, and the Emergence of the Digital World[END_REF] . Therefore, this type of stack is manufacturable without major difficulties using state-of-the art fabrication tools Fig. 3 Schematic of the memristive device with the detailed composition.

The devices are fabricated in arrays of pillars with diameters between 80nm and 200nm. After nanofabrication, the isotropic properties of the free layer have been confirmed [START_REF] Mansueto | Isotropically coercive free layer integration in a magnetic tunnel junction for neuromorphic applications[END_REF] .

Single pulse experiments

The first experiments aimed to find the pulse characteristics (amplitude and duration) for which the STT yields a discrete rotation of the free layer magnetization and consequently discrete variation of resistance value. To do this, the device was first set in the parallel state through the application of a strong in-plane field. Then, a single pulse was applied followed by the read out of the final resistance state. The device was reset in the P state before each pulse. This procedure was repeated for different pulse amplitudes and durations similarly to a switching probability experiment on a binary device. In this case, instead of the switching probability, we expect to see the actual resistance variation dependence on the pulse parameters. The resulting experimental points (each averaged over 10 measurements) are shown in Fig. 4a for a junction with a nominal diameter of 100nm and a TMR of 30%. The magnetization angle is directly extracted from the measurement through the cosine variation of the conductance [START_REF] Jaffrès | Angular dependence of the tunnel magnetoresistance in transition-metal-based junctions[END_REF] and shown in Fig. 4b. As expected, a monotonicity of the resistance variation is obtained with respect to both the pulse amplitude and duration. For each amplitude of the voltage, a threshold pulse length separates a region in which the magnetization is not moving from a region in which discrete rotations are obtained. Those threshold values are decreasing with the increase of the voltage amplitude, as expected. The experimental points are fitted with a line that represents an average behavior assuming that the magnetization angle varies linearly with time (constant angular velocity). We note that the oscillation frequencies extracted turn out to be significantly lower than the GHz range. This is likely the result of the important increase of the Gilbert damping already observed in FM/AF bilayers [START_REF] Spenato | Study of dynamic properties and magnetic anisotropies of NiFe/MnPt in the critical thickness range[END_REF][START_REF] Mccord | Magnetization Dynamics of Coupled Ferromagnetic-Antiferromagnetic Thin Films[END_REF] as well as FM/AF/FM trilayers [START_REF] Xu | Influence of Thickness on Magnetic Properties and Microwave Characteristics of NiFe/IrMn/NiFe Trilayers[END_REF] (6) from ref [START_REF] Mansueto | Realizing an Isotropically Coercive Magnetic Layer for Memristive Applications by Analogy to Dry Friction[END_REF] , where 𝜃 𝑀 is the magnetization polar angle. The linear behavior of this last equation with respect to the applied voltage can be used to fit the experimental points as in Fig. 4c. From this fit, the extracted threshold voltage for the triggering of the OPP is around 0.6V. From this value, that corresponds to the critical line 𝐽 𝑝𝑒𝑟𝑝 in eq. 4, one can estimate a value of 𝜂 𝑝𝑒𝑟𝑝 ≈ 0.15 ± 0.02 that leads, through eq. 6, to a value of 𝛼 𝐺 ≈ 0.07 ± 0.01. The experimental points can be compared to the analytical diagram in Fig. 1. The comparison indicates that the device operates indeed in the unstability region as expected. Note that, in order to achieve the limited angular excursion of the magnetization (as in the orange curve in Fig. 2a), a ratio 𝜂 𝑎𝑛 /𝜂 𝑝𝑒𝑟𝑝 ≈ 2 seems to be sufficient experimentally.

Despite the fact that the values of 𝜂 𝑝𝑒𝑟𝑝 and 𝛼 𝐺 perfectly match the expectations, they can be affected by an error coming from the fact that the equations used to extract them are valid in the OPP region and in the long pulses approximation while the experimental points are taken in the unstability region where no analytical expressions can be computed. Also, the importance of the STT coming from the perpendicular polarizer is confirmed by the fact that, in devices as in Fig. 3 but without perpendicular polarizer, it is observed that no switching occurs for pulse parameters as in Fig. 4a, exactly as predicted by the blue line in Fig. 2a.

Multiple pulses experiments

In the following experiment, a train of pulses was sent to the device to observe the evolution of its resistance. The device was set in the P state only at the beginning and not between two consecutive pulses as in the previous experiment. The pulses were designed in a way to not switch completely the magnetization but to rotate it so as to reach several intermediate resistance values. For this reason, the voltage amplitude was set to 0.9V with a varying pulse duration for the results shown in the following. An important parameter chosen for this experiment was the delay between two consecutive pulses that was set to 1s to get an idea of the thermal stability of the intermediate resistance states (discussed in the next section). The resistance variation to a series of 50 pulses of positive polarity is shown in Fig. 5a for different pulse durations. The magnetization, that is not much affected by the 1ns pulses, gradually rotates towards the AP state proportionally to the pulse duration. The resistance is then saturating to the AP state confirming that the effect of the double STT acts as in the unstable region of Fig. 1 and not as in OPP (for which the resistance is supposed to oscillate with the same polarity of the voltage). In both insets of Fig. 5a, it can be observed that the intermediate resistance states turn out to be thermally stable along the 1s between the two pulses confirming the possibility to stabilize the magnetization along different in-plane directions at least during a 1s period. One observes that although the overall trend of the resistance variation under a train of pulses of same polarity is monotonous, the individual resistance jumps are not always varying in a monotonous way (as in the two insets of Fig. 5a). A random telegraph type of noise in the resistance variation seems to be superposed to the general trend. For this reason, the number of pulses needed to completely switch from the P to the AP state is larger than what would be predicted from Fig. 4a (for which an average over 10 measurements reduces the random effect). This effect, that will be discussed in the next paragraph, is not dependent on the voltage polarity as shown in Fig. 5b. Even if the STT effect is visible in the monotonicity of the average resistance with the polarity of the voltage, the presence of random resistance jumps associated with voltage pulses persists even in the P or AP states where the STT from the in-plane analyzer is supposed to be zero (but not the one from the perpendicular polarizer). We emphasize that what we note as P and AP states here are not fully saturated P and AP states but rather remanent state of the free layer magnetization along the P and AP directions. We know already from earlier sheet film studies [START_REF] Mansueto | Realizing an Isotropically Coercive Magnetic Layer for Memristive Applications by Analogy to Dry Friction[END_REF] that the hysteresis loops of the free layer are not square. The remanence is about 80% of the saturation magnetization (Fig. 4b in ref 12). At the microscopic scale, due to the randomness induced by the interfacial coupling between the ferromagnetic and antiferromagnetic layer, the magnetization of the free layer is viewed as being distributed in a cone around an average direction. This explains why the resistance jumps around the nominally AP state can reach higher resistance values than those represented by the dotted lines in Fig. 5a andb which correspond to the resistance value in the Please do not adjust margins Please do not adjust margins AP remanent state and not AP saturated state. The same reasoning applies around the P state for resistance jumps reaching resistance values below that in the remanent P state.

Note also that the P-AP and AP-P transitions are asymmetric because of the residual dipolar field from the uncompensated out-of-plane SAF and/or the in-plane one that favors precession only towards one direction.

Fig. 5 a) Resistance variation due to a train of pulses of 0.9V with varying pulse durations and 1s delay between two consecutive pulses. In the zooms, the blue curve shows an almost monotonous behavior while the green curve exhibits an average monotonous trend of the resistance variation superposed with random jumps. b) Resistance variation under a train of pulses of both polarities. These measurements were performed on two different MTJs of nominal diameter 100nm exhibiting a TMR of ~30%.

Temperature effect

The presence of a thin IrMn layer in the free layer, with a blocking temperature of the order of 100°C, makes this device sensitive to temperature. For this reason, the effect of temperature was investigated in terms of stability of the intermediate resistance states and Joule heating induced by the voltage pulses.

A first glimpse of the thermal stability of the resistance levels reached through STT is shown in the insets of Fig. 5a. An example of resistance distribution associated with such curve is shown in Fig. 6a. The distribution of intermediate states depends on the pulse amplitude and duration since the angular step that the storage layer magnetization makes at each pulse depends on these parameters. The distribution shown in Fig. 6a was derived from the resistance variation shown in Fig. 5a for pulse 3ns long. An almost continuous range of resistances is achieved between the P and AP states. To analyse the stability of single resistance levels, only few levels were selected as in Fig. 6b. The selection criteria in order to make them distinguishable considers that the distance between two peaks has to be greater than the largest FWHM of the two distributions. The distribution of such resistances show a FWHM below 3Ω, meaning a resistance variation below 0.1% in a time range of 1s. To determine the thermal stability of the resistance states at RT, the magnetization was saturated, through the application of a strong inplane field, in different directions and the resistance was measured for 3 minutes. As can be seen from Fig. 6c, the FWHM is around 17Ω with a resistance variation below 2%. The results suggest that, for a device with TMR 100% and same 𝑅 × 𝐴, about 50 intermediate states could be used for applications in which the read/write access would occur at typical frequency of a few Hz. Moreover, this confirms that the resistance jumps, observed in the previous experiment and that occur during the pulses, are mainly due to Joule heating and correspond to thermally activated switching of small domains in the free layer. To further investigate this, we tried to separate the effects of STT and Joule heating on the resistance variation. For each polarity of the voltage, to minimize the role of the STT compared to that of Joule heating, we set the system close to the P or AP state. This strongly reduces the STT from the in-plane analyzer but not that of the perpendicular polarizer. A statistical analysis of the resistance jumps due to a train of pulses in the saturated AP state is shown in Fig. 7a for different pulse lengths in a device of 150nm of diameter. Interestingly, the distributions of jumps turn out to be peaked Please do not adjust margins Please do not adjust margins around an average value corresponding to a small fraction of the normalized TMR. This means that, when the magnetization is in the AP state, the temperature effect is to switch the magnetization of a relatively small part of the free layer with a certain characteristic size. Moreover, this area is observed to increase with the pulse length which can be ascribed to a more significant heating of the device with longer pulses (Fig. 7a). From the knowledge of the amplitude of the resistance jumps (plotted as horizontal scale in Fig. 7a), one can derive the size of the switching domains, the distribution of which is plotted in Fig. 7b. The Joule heating activates the switching of small portions of the free layer as it would do in a granular layer. To further investigate this, the same analysis was performed on devices with different diameters with a fixed pulse length of 3ns. As shown in Fig. 7c, the resistance jumps appear in a similar peaked distributions whatever the size of the device (the noise limit is taken at 2% according to the previous measurement). Furthermore, the amplitude of the resistance jumps decreases as the diameter of the device increases (Fig. 7c, shift towards smaller normalized TMR jumps as the diameter of the device increases). However, if we now recalculate the size of the switching domains taking into account the area of the various devices, we find that the three distributions present the same average value for the diameter of the switching domains of the order of 45nm as shown in Fig. 7d. This tends to confirm the idea that the free layer as a sort of granular behavior with grain size of the order of 45nm. This phenomenon can be explained with the presence of the IrMn layer in the free layer. Its granular structure and its spin glass-like magnetic state can create, through exchange interaction at the interface with the two ferromagnetic layers, a granular type of magnetic behavior [START_REF] Mansueto | Realizing an Isotropically Coercive Magnetic Layer for Memristive Applications by Analogy to Dry Friction[END_REF] . These grains can switch even around the P and AP states when the voltage pulses are applied by thermal activation. The temperature reached in the free layer during a voltage pulse can be estimated in a MTJ within a unidimensional model through the expression

𝑇(𝑡) = 𝑇 0 + 𝑃 𝑑 𝐴 𝐾 [1 -𝑒𝑥𝑝 (- 𝑡 𝜏 )] (7)
where 𝑇 0 is the initial standby temperature (RT), 𝑃 𝑑 the heating power density during the pulse (𝑉² 𝑅𝐴 ⁄ ), 𝐴 the device area, K the effective heat conductance between the tunnel barrier and the bottom and top electrodes considered to act as thermal bath at 𝑇 0 , 𝜏 the characteristic heating time which depends on the pillar specific heat and thermal conductance and t the pulse length [START_REF] Prejbeanu | Thermally assisted MRAM[END_REF] . From earlier studies on thermally assisted MRAM which used relatively similar material stacks [START_REF] Spenato | Study of dynamic properties and magnetic anisotropies of NiFe/MnPt in the critical thickness range[END_REF] , the following statements can be made : For 𝑡 < 20ns, the temperature increases rapidly with the pulse length t and almost independently on the device size. For 𝑡 > 20ns, the temperature increase tends to slow down due to heat diffusion towards the top and bottom metallic electrodes and gradually the temperature reaches an asymptotic regime which corresponds to a balance between the heat production around the tunnel barrier and the heat flow towards the electrodes. In our experiments using pulses between 1ns and 4ns long, we are still in a regime where the temperature rise in the MTJ is roughly proportional to the pulse duration. Considering that in thermally assisted MRAM, a temperature rise of about 150°C was achieved in about 5ns with a heating power of 100mW/m² and taking into account the dependence on the MTJ diameter, we can roughly estimate that in the present junction, a 2ns pulse of 0.9V can yield an increase of temperature of the order of 40°C while a 4ns long pulse yields an increase of temperature of about 80°C. Since the coercivity of the free layer measured in quasistatic mode by VSM tends to vanish at 100°C, the Joule heating induced by the voltage pulses can indeed trigger the switching of small domains in the free layer.

Conclusions

In conclusion, we proposed a device based on the angular variation of TMR in an in-plane MTJ. Thanks to the property of isotropic coercivity of a free layer based on a FM/AF/FM structure, we are able to discretely rotate the magnetization by STT and stabilize it at different in-plane angles. We found that in single pulse experiments, we could adjust the relative amplitude of the STT originating from the in-plane analyzer and perpendicular polarizer to limit the angular excursion of the free layer magnetization in a half-plane which is mandatory to achieve a monotonous variation of resistance for each pulse polarity as required for memristor functionality. In the multiple pulses experiments, we demonstrated an average monotonous behavior of the resistance variation with the voltage polarity. An almost continuous variation of resistance could be achieved between a minimum resistance state corresponding to parallel remanent magnetic configuration and a maximum resistance corresponding to antiparallel remanent magnetic configuration. By analyzing few of these levels, resistance variations below 0.1% in 1s timescale and below 2% in 3 minutes were found. For both voltage polarities, random resistance jumps are observed on top of the general Please do not adjust margins Please do not adjust margins monotonous trend in resistance variation. These random jumps are associated with thermally activated switching of small domains in the free layer magnetization due to Joule heating occurring during the current pulses. These domains are about 45nm in dimension whatever the diameter of the pillar and seems to be associated with a granular structure of the free layer.

The device here described has intrinsically a good down-size scalability since its working principle (based on angular variation of the TMR) is not directly linked to its size. In particular, a reduction in the diameter of such device would first reduce the Joule heating [START_REF] Prejbeanu | Thermally assisted MRAM[END_REF] thanks to additional heat dissipation though the sidewalls of the pillar and, at the same time, bring the device to a more macrospin behavior thus reducing significantly the random contributions occurring during application of the voltage pulses. Such memristor could be implemented in crossbar architectures to perform the cumulate and multiply function, key in neuromorphic computing [START_REF] Grollier | Neuromorphic spintronics[END_REF] .

Fig. 1

 1 Fig.1 Analytical diagram of the working regions for a device with 𝜂 𝑎𝑛 = 0.3. The experimental points in black, red and blue correspond to an applied voltage respectively of 0.8V, 1V and 1.2V. The parameters used for all the simulations presented in this work are: 𝑀 𝑠 = 10 6 A/m, 𝛼 𝐺 = 0.04, 𝑅 × 𝐴 = 10 Ωµm², 𝑙 = 4 nm, 𝛽 = 5 x 10 8 rad/s, Ku = 0 J/m 3 , N = (0.05 0.05 0.899).

Fig. 2 a

 2 Fig.2 a) Magnetization dynamics under STT in the unstability region for different values of 𝜂 𝑝𝑒𝑟𝑝 and for current density J = 2.5 x 10 11 A/m². b Resistance variation (top panel) induced by a train of pulses

Ta0. 2 /

 2 FeCoB (1.6) FeCoB (1.2)/Ta0.2 Please do not adjust margins Please do not adjust margins

Fig. 4 a

 4 Fig.4 a) Resistance variation with the pulse duration for different pulse amplitudes. b) In-plane angle of the magnetization varying with the pulse duration as in a). The black, red and blue fitting lines correspond respectively to angular velocity of the magnetization. c) Linear fitting of the angular velocity variation with the pulse amplitude.

Fig. 6 a

 6 Fig.6 a) Resistance distribution in a multiple pulse experiment as in Fig. 5a. b) Resistance distribution for few distinguishable levels. c) Resistance distribution for intermediate resistance values induced by field saturation at different in-plane angles for a 3min timescale. The measurements shown in this figure were performed on a 100nm MTJ with a TMR around 26%.

Fig. 7 a

 7 Fig.7 a) Resistance jumps distributions for different pulse lengths in a device of diameter 150nm. b) Resistance jumps distributions as in a) expressed versus diameter of the switching domains. c) Resistance jumps distributions for different device diameters and pulse length of 3ns. d) Resistance jumps distributions as in c) normalized with respect to the corresponding diameters.

  and expected from the expression

	𝜔 𝑀 = -	𝛾𝑎 𝑝𝑒𝑟𝑝// 𝑉 𝛼 𝐺	-	𝛽 𝛼 𝐺 𝑠𝑖𝑛(𝜃 𝑀 )	𝜔 𝑀 |𝜔 𝑀 |

This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 5

This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 7

Acknowledgements

This work was funded under the ERC Adv grant MAGICAL N° 669204.


Laurent Proville

Dynamical structure factor of a nonlinear Klein Gordon lattice

Keywords: numbers: 63.20.Ry, 61.10.Dp, 61.12.Bt, 61.14.Dc

The quantum modes of a nonlinear Klein Gordon lattice have been computed numerically [L. Proville, Phys. Rev. B 71, 104306 (2005)]. The on-site nonlinearity has been found to lead to phonon bound states. In the present paper, we compute numerically the dynamical structure factor so as to simulate the coherent scattering cross-section at low temperature. The inelastic contribution is studied as a function of the on-site anharmonicity. Interestingly, our numerical method is not limited to the weak anharmonicity and permits to study thoroughly the spectra of nonlinear phonons.

I. INTRODUCTION

The metal hydrides are typical compounds whose the technological importance 1 aroused many spectroscopic inquiries. The spectroscopy of hydrogen modes allows to work out the occupied interstitial gaps, e.g. octahedral or tetrahedral and thus the phase structure of the alloy. Further, the determination of the anharmonicity of hydrogen sites 2,3 permits to evaluate the on-site potential landscape 3 and so the hydrogen ability to migrate or to dimerise. The inelastic neutron scattering (INS) revealed the on-site anharmonicity in the model hydrides as NbH, TaH and PdH (for a review see Ref. 4). In the compounds TiH and ZrH, that anharmonicity was found [START_REF] Kolesnikov | [END_REF]6 to lead to the optical phonon bound states, i.e., some well-known excitations in molecular crystals 7,[START_REF] Agranovich | Spectroscopy and Excitation Dynamics of Condensed Molecular Systems[END_REF][START_REF] Califano | Lattice dynamics of molecular crystals[END_REF][START_REF] Gellini | [END_REF]11 (also quoted as biphonon and bivibron) where the anharmonicity of internal covalent bonds overpasses the interaction between molecules. As in the metal hydrides, the proton dynamics has been studied by INS in molecular crystals as polyglycine 12 and 4-methyl pyridine 13 . In the latter, the bound states of the methyl group rotational modes proved to last several days 13,14 (see also Refs. 15,16). The trapping of energy, whether it is in proton vibrations or in the intrinsic modes of molecules is a consequence of the emergence of breather, a theoretical paradigm that has been matter of intensive research in different contexts (for instance see Refs. 17,18,19,20,21,22,23, 24,25,26,27,28,29). The distinctive property of breather is to gather the spatial localization and time periodicity. For the case of a translational-invariant quantum lattice, the phonon bound states can be thought as the siblings of breathers 19,20,30,31 , as phonon bound states and breathers both stem from the lattice anharmonicity. The formers appeared, though, earlier in literature (see Refs. 25,32 and references therein). The link between breather and phonon bound states has been studied in different works [START_REF] Scott | Nonlinear science[END_REF]24,25,31,34 . Consequently, the experiments that enhanced the existence of phonon bound states (for instance see Refs. 5,7,10) are equally some concrete evidences of the emergence of quantum breathers. In addition, several experiments 35,36,37 were recently dedicated to breather and its important role in the energy storage and transport. As an overview, one may retain that the studies on breather aim at improving the standard harmonic treatment of lattice dynamics which flaws are well-known in solids [START_REF] Ashcroft | Solid State Physics[END_REF] .

The scattering of neutrons [START_REF] Lovesey | Theory of neutron scattering from condensed matter[END_REF] , X-rays [START_REF] Burkel | Inelastic scattering of X-rays with very high energy resolution[END_REF] and electrons [START_REF] Crescenzi | Electron scattering and related spectroscopies[END_REF] provided numerous insights into condensed matter and molecular physics. In particular, the inelastic scattering allows to probe the dynamics of crystals and molecules. The dynamical structure factor, S(q, ω) which is proportional to the scattering cross-section can be calculated by solving the Schrödinger equation for the lattice modes (within Born approximation), so this quantity is essential for bridging theory to experiment. However the lattice eigen-modes may be figured out exactly only in the harmonic lattice model (see for instance Ref. 39,42) and thus the anharmonic contributions to the energy are neglected, although they stem from the atomic interaction potential. In order to examin the effect of anharmonicity on S(q, ω), some approaches have been attempted in different quantum lattice models, e.g., the Hubbard model for bosons [START_REF] Agranovich | [END_REF]44 and the Klein Gordon model with a weak on-site anharmonicity 30 . The latter is often quoted as nonlinear Klein Gordon lattice (KG) and, noteworthy, it may account for the intrinsic anharmonicity of molecular crystals 7,[START_REF] Califano | Lattice dynamics of molecular crystals[END_REF][START_REF] Gellini | [END_REF]11 , molecule bonds 33,45 or metals where interstitial gaps are filled with light particles, as metal hydrides [START_REF] Eckert | [END_REF]3,[START_REF] Kolesnikov | [END_REF]6 . It is the reason why we proposed a numerical method 31 that is tractable for different type of nonlinearity, to compute the nonlinear quantum modes. In the present paper, our previous developments are used to evaluate the coherent Dynamical Structure Factor (DSF) of the KG lattice, at low temperature. The contribution to the DSF of the nonlinear quantum modes, known as either the phonon bound states [START_REF] Kolesnikov | [END_REF]7,[START_REF] Agranovich | Spectroscopy and Excitation Dynamics of Condensed Molecular Systems[END_REF][START_REF] Califano | Lattice dynamics of molecular crystals[END_REF][START_REF] Gellini | [END_REF] or else the quantum breathers 20,23,30 , is the central purpose of our work. As the DSF standard derivation is obtained in the harmonic approximation [START_REF] Lovesey | Theory of neutron scattering from condensed matter[END_REF]42 , treating the anharmonicity as a perturbation, we propose a different scheme where this approximation is not required. For instance, the Bloch identity, which is a key for the conventional calculation, is not invoked in our theory. We simply introduce a Taylor expansion of the DSF with respect to the atomic displacements, computing numerically the coefficients of the series.

Our method is tested by comparison to the standard analytical calculation on the purely harmonic lattice. Our approach allows us to deal with different strength of the anharmonicity. When the on-site nonlinearity dominates the intersite coupling, the nonlinear quantum modes lead to some anharmonic resonances well separated from the multi-phonon continua, in the inelastic spectrum. These anharmonic resonances have same order as for fundamental phonons whereas the unbound multi-phonons yield a DSF at least one order below. On the other hand, their dispersion is found to decrease dramatically with the value of the energy transfer, which differentiates them from the broad phonon resonance. By way of contrast, when the intersite coupling is larger than the on-site nonlinearity, the anharmonicity of the spectrum is weak compared to the phonon dispersion. In such a quasi-harmonic lattice, the dynamical response of the biphonon (a two phonon bound state) may yet exhibit a significant magnitude, still much larger than the two-phonon, even though the biphonon resonance occurs at an energy transfer that approaches the unbound phonons band. Actually, the biphonon signature is found to dominate the two-phonon DSF provided that the lattice nonlinearity is not strictly zero. Alongside our numerics, a perturbation theory is developed and proves accurate for strongly nonlinear lattices. We also show how the model parameters may be adjusted so as to simulate the inelastic scattering of a material and thus to work out the energy landscape of the inner particles as well as their interactions. Finally, the present study will serve as a basis for future simulations of the KG incoherent scattering cross-section.

After a brief introduction to the nonlinear KG lattice, the DSF is derived in Sec. II. in Sec. III. Results and comments are detailed in the same section while the Sec. IV summarizes them and draws some perspectives.

II. COHERENT DYNAMICAL STRUCTURE FACTOR OF THE NONLINEAR KG LATTICE

We assume that some light particles form a regular network, whether it is inside a molecule or a crystalline solid. Were the particle dynamics independent, the Hamiltonian would have read:

H 0 = l [ p 2 l 2m + V (x l )], (1)
where x l and p l are displacement and conjugate momentum (i.e., [x l , p l ] = i) of a particle at site l. For sake of simplicity, we choose to work on a one-dimensional lattice with a single direction for the motion of atoms, denoted by the unit vector u, which is named as polarization in the followings. The single dimension proved relevant in several concrete studies [START_REF] Kolesnikov | [END_REF]6,13,14 . The local potential V (x l ) may then be expanded as a fourth order series:

V (x l ) = a 2 x 2 l + a 3 x 3 l + a 4 x 4
l . Some higher order terms could be added with no difficulties for our numerics, so as to simulate eventually a specific shape of V . In truth, the interaction between nearest neighbors, that may be direct or mediated by the heavy atoms that compose the skeleton of a material, involves a displacement coupling that is modeled by a quadratic term. So the Hamiltonian of the particle ensemble now reads:

H = H 0 -c l,j=<l> (x l -x j ) 2 , ( 2
)
where c is the coupling parameter and the symbol (< l >) describes the first neighbors of the site labeled by l. Introducing the dimensionless operators P l and X l , the Hamiltonian can be reformulated as follows:

H = Ω l [ P 2 l 2 + A 2 X 2 l + A 3 X 3 l + A 4 X 4 l + C 2 X l j=<l> X j ], (3)
where the fundamental frequency Ω, as well as the dimensionless coefficients A 3 , A 4 and C have been defined in Ref. 31. The total number of sites is denoted by N, the lattice parameter by a 0 and the orientation of the chain is defined by a unit vector v. We introduce ψ, the eigenstates of H, and particularly the groundstate ψ GS , the phonons ψ ph (k) and the biphonons ψ bi (k). These states are computed numerically 31 as well as their eigenvalues E GS , E ph (k) and E bi (k), respectively. The wave momentum, denoted by k, verifies k = 2πk N a 0 v where k is an integer. The groundstate wave momentum k 0 takes its value in the reciprocal lattice, k 0 = 2πk 0 a 0 v where k 0 is an integer. Apart when it will be noted, k 0 is fixed to zero. The lattice is assumed to be maintained at very low temperature, i.e., kT < E ph , so the scattering induces some excitations solely from ψ GS . The transition probability is proportional to the DSF 42 :

S(q, ω ψ(k) ) = 1 N | j < ψ GS |e i[q.r j ] |ψ(k) > | 2 (4)
where ω ψ(k) = (E ψ(k) -E ψ GS )/ , q is the scattering vector and r j represents the atomic position at site j. That position can be expressed as r j = ja 0 v+LX j u, where L = /(mΩ)

is the length scale of vibrations, fixed to a reasonable value of 3% of a 0 . The polarization vector u is the principal axis of atomic displacement around the equilibrium position (ja 0 v).

We now spell out our method for computing the DSF. As, strictly speaking the Bloch identity does not hold for a nonlinear lattice 42 , we propose a derivation which differs from the conventional one. Since L is much smaller than a 0 , the exponential function in Eq. ( 4) can be expanded as a Taylor series with respect to the weakest of its arguments, the term proportional to L:

S(q, ω ψ(k) ) = 1 N | j e i[q.v]ja 0 × p (i[q.u]L) p p! < ψ GS |X p j |ψ(k) > | 2 . ( 5
)
This development holds provided that the order of the expansion is large enough. The series convergence has been tested by increasing the order up to reach the desired precision. The agreement to the standard analytical calculation on the purely harmonic lattice (see Sec.

III) is a necessary condition of validity that has been fulfilled too. In the left hand side of Eq. ( 5), the bracket < ψ GS |X p j |ψ(k) > can be replaced by < ψ GS |X p 0 |ψ(k) > multiplied by a phase factor e -i([k.v]a 0 j) , because of the translational invariance. The sum over the subscript j in Eq. (5) gives zero for all q, aside from the wave vectors that match the momentum conservation [(qk).v] = 0. This condition, as well as the conservation of the energy are supposed to be satisfied for a given state ψ(k). The calculation of S(q, ω ψ ) can then be reduced to the determination of the bracket D(ψ(k), p) =< ψ GS |X p 0 |ψ(k) >, since one may check that:

S(q, ω ψ(k) ) = Nδ [(q-k).v] | p (i[q.u]L) p p! D(ψ(k), p)| 2 . ( 6
)
We study the variation of the DSF along a single direction which the angles with respect to u and v are α u and α v , respectively. Then S(q, ω ψ ) depends only on the magnitude of the momentum transfer |q| along that direction and the conservation law fixes

|q|cos(α v ) = [k.v].
This condition can be achieved whether |α v | = π/2 which may be reasonably assumed for a low-angle scattering. In a same manner, when |α u | = π/2 the inelastic part of the structure factor is identically null because [q.u] = 0 (excepted for the elastic DSF which then equals N). Consequently, we assume that |α v | = |α u | = π/2 which keeps the physics of our problem, avoiding the unimportant cases. We denote by q the scalar product [q.u] for a vector q that matches the conservation law, so that we have also q = [k.v]. Then, the DSF reads:

S(q, ω ψ(q) ) = N| p (iqL) p p! D(ψ(q), p)| 2 . ( 7
)
Since the optical modes are characterized by a weak value of C in Eq. ( 3), the first step of our treatment is concerned with a perturbation theory with respect to the intersite coupling.

At C = 0, namely the anti-continuum (e.g. , uncoupled, atomic or molecular) limit 19 , the phonon states and the phonon bound states can all be written as some Bloch waves 31 :

B α (k) = 1 √ N j e -i[k.v]a 0 j × φ α,j Π l =j φ 0,l (8)
where φ α,j is a on-site wave function that depends only on X j . Actually, φ α,j is the αth eigenstate of the on-site Hamiltonian:

h l = P 2 l 2 + A 2 X 2 l + A 3 X 3 l + A 4 X 4 l . (9)
To a first order in C, the lattice groundstate ψ GS is simply given by the product B 0 = Π l φ 0,l while Eq. ( 8) gives a phonon for α = 1 and a biphonon for α = 2. Let us denote by V α,n , the projection of the state φ α onto the nth on-site harmonic oscillator eigen-state |n > and

develop the bracket G(α, 0, p) =< φ 0 |X p |φ α >: G(α, 0, p) = n≥0,m≥0 V 0,m V α,n < m|X p |n > ( 10
)
where the subscript has been dropped for ease. We point out that for C = 0, we have

D(B α (k), p) = G(α, 0, p)/ √ N and D(B 0 , p) = G(0, 0, p).
One more thing has to be done before reaching our goal, is to compute < m|X p |n >. To that purpose, the Bose-Einstein operators [i.e. a + = √ 2(X -iP ) and a = √ 2(X + iP )] are used to expand X p . We wrote a fortran program which realizes the expansion, respecting the commutation rule, [a, a + ] = 1.

For example, the output for p = 10 is:

X 10 = 1 32 [a +10 + a 10 + 10(a + a 9 + a +9 a) + 45(a +2 a 8 + a +8 a 2 + a +8 + a 8 )
+ 120(a +3 a 7 + a +7 a 3 ) + 360(a + a 7 + a +7 a)

+ 210(a +4 a 6 + a +6 a 4 ) + 1260(a +2 a 6 + a +6 a 2 ) + 630(a +6 + a 6 ) + 252a +5 a 5 + 2520(a +5 a 3 + a +3 a 5 )

+ 3780(a +5 a + a + a 5 ) + 3150(a +4 + a 4 + a +4 a 4 ) + 4725(2a +2 a 4 + 2a +4 a 2 + 4a +2 a 2 + a +2 + a 2 + 2a + a) + 12600(a + a 3 + a +3 a + a +3 a 3 ) + 945]. (11)
The writing of X p would take more than one full page for values of p larger than 30. Our program allows us to compute the bracket < m|X p |n > up to the power p = 60, for different integers m and n. The association of this program with the numerical diagonalization of h l (Eq. ( 9)) which fixes the coefficients V α,n in Eq. ( 10) 31 , permits to compute the coefficients G(α, 0, p) that can be tabulated for different model parameters. Finally, to a first order in C, one obtains:

S(q, ω Bα(q) ) = | p (iqL) p p! G(α, 0, p)| 2 (12)
for the inelastic scattering where the momentum conservation imposes q = 2πk/(Na 0 ), whereas the elastic response is given by:

S(q, 0) = N| p (iqL) p p! G(0, 0, p)| 2 (13)
where q = 2πk 0 /a 0 [k and k 0 range over integers]. Formally, these results are given for a one-dimensional lattice but they can be extended to higher dimensional lattice by summing over the coordinates and polarizations. To a first order in C, the frequency ω Bα(q) can also be evaluated by:

ω Bα(q) = Ω(γ α -γ 0 -2C × G(α, 0, 1) 2 cos(qa 0 )) ( 14
)
where the coefficient γ α is the eigen-energy of h l (Eq. ( 9)) associated to Φ α . For the perfect harmonic lattice, the first order in C in the standard calculation [START_REF] Lovesey | Theory of neutron scattering from condensed matter[END_REF]42 of S(q, ω) is equivalent to Eqs. ( 12) and ( 13).

When C is large, the Hamiltonian in Eq. ( 3) must be diagonalized numerically after expanding in a suitable basis. We worked with a Bloch wave basis 31 given by:

B [Π j α j ] (k) = 1 A [Π j α j ] j e -i[k.v]a 0 j × Π l φ α l ,l+j (15)
where A [Π j α j ] ensures the normalization. Each eigenstate, characterized by a wave vector k, can be written as a linear combination of those Bloch waves:

ψ(k) = Π j α j W ψ,Π j α j B [Π j α j ] (k), (16)
where the subscript Π j α j identifies a single Bloch wave in Eq. ( 15). The numerical diagonalization has been carried out for different lattice sizes with no noticeable discrepancy in the eigenspectrum in increasing N. In Fig. 1, the same energy cutoff on the Bloch wave basis has been fixed for both N = 23 and N = 33. Apart on their wave vector, the eigenvalues are found to be independent of N. In Fig. 1 (a), the optical phonon branch is plotted whereas, in Fig. 1 (b), the energy region of the first overtone is plotted for same parameters as in (a). One notes clearly the biphonon branch and the two-phonon band, thoroughly described elsewhere [START_REF] Agranovich | Spectroscopy and Excitation Dynamics of Condensed Molecular Systems[END_REF][START_REF] Scott | Nonlinear science[END_REF][START_REF] Eilbeck | Localization and Energy Transfer in Nonlinear Systems[END_REF]30,31 . Interestingly, for small enough nonlinear parameters, the anharmonicity of the lattice is negligible compared to the phonon branch width and the biphonon branch disappears (see Fig. 1 (c)). The bracket D(ψ, p) can now be written as follows:

D(ψ(k), p) = (Π j α j ,Π j β j ) W * ψ GS ,Π j β j W ψ,Π j α j < B [Π j β j ] (0)|X p 0 |B [Π j α j ] (k) > ( 17
)
where the last term in the right hand side can be detailed further:

< B [Π j β j ] (0)|X p 0 |B [Π j α j ] (k) > = 1 A [Π j α j ] A [Π j β j ] n 1 ,n 2 e -i[k.v].a 0 n 1 ×Π l =j δ(α l+n 1 , β l+n 2 )G(α l+n 2 , α l+n 1 , p). (18)
This equation concludes our computation task on the KG lattice dynamical response to the scattering of a external beam of light or particles. Our numerical treatment, including the diagonalization of H and the computations of Eqs. 7, 17 and 18 takes few hours on a conventional desktop PC. The convergence of the series in Eq. ( 7) has been tested by comparing our results for different development orders, e.g. 30, 40, 50 and 60. As expected, the smaller is the scattering vector, the better is the precision. For a momentum transfer that ranges over 30 lattice Brillouin zones, the difference between the DSF computed with series orders 50 and 60 is less than 1 %.

III. RESULTS AND DISCUSSIONS

The Fig. 2 shows our typical result that is a 3D plot of the inelastic S(q, ω) for a nonlinear lattice which model parameters are those of Figs. 1 (a) and (b). Three noteworthy resonances emerge with same order of magnitude. These three resonances correspond to the eigen-energies of the nonlinear phonon states, i.e., the phonons, biphonons and triphonons.

The binding energy of the biphonon, evaluated at the center of the Brillouin zone (BZ) is around a few percent of the fundamental phonon excitation, as measured in different metal hydrides [START_REF] Kolesnikov | [END_REF]6 . Several model parameters have been tested and lead to qualitatively similar results. For instance, besides the different values of the energy transfer, the biphonon and triphonon resonances occur in a lattice with a quadratic-quartic on-site potential, i.e., A 3 = 0 in Eq. 3. In the biphonon and triphonon resonances, one may recognize some peaks that will be examined further. In order to analyze our results more quantitatively it is, though, easier to work with a 2D plot that represents the S(q, ω) profile versus the projection q = [q.u] of the scattering vector q. In order to verify the accuracy of our computations, presented in Sec. II, we compare our numerical results to the exact analytical ones that are achieved in a purely harmonic lattice. In that case, the profile of the Debye-Waller factor, i.e., S(q, 0)/N is plotted in Fig. 3 and shows a convincing agreement since our numerical results scarcely differ from the analytical ones. It ensures us that our series expansion of S(q, ω) has converged. Only the non-zero values of the Debye-Waller have been retained in the plot, i.e., q = 2πk 0 /a 0 . A similar accord is obtained for the inelastic part of S(q, ω) (see further in the same section). The calculation made in a standard harmonic approximation 42

gives a Gaussian dependency of the Debye-Waller factor S(q, 0)/N = exp(-2W (q)) where

2W (q) = (qL) 2 | < ψ GS |X 2 j |ψ GS > |.
In this expression, the bracket is simply given by 1/(2N) K 1/ω ph (K) where ω ph (K) = 1 + 2Ccos(K.a 0 ) and K ranges in the first BZ.

When the nonlinear parameters are no longer negligible in the Eq. ( 3), we note, in Fig. 3 that the Debye-Waller of a nonlinear lattice differs from the harmonic one. The Gaussian form is yet roughly conserved so there is no qualitative changes involved by the on-site anharmonicity. Our perturbation theory proves sufficient under the condition that the intermolecular coupling is not too large. In agreement with our results, the earlier study of B.V. Thompson 46 concluded that a cubic-anharmonic term involved a negligible correction upon the Debye-Waller factor. The calculation of Thompson is, in fact, a second order perturbation in the cubic term, which is thus assumed to be weak enough for the perturbation theory to be valid. Not surprizingly, it leads to a correction that is also weak. However, as we shall see, the most significant contribution of the nonlinearity to the DSF spectra proves to be inelastic and is due to the phonon bound states.

In a same manner as for the elastic scattering, our numerical computation of the inelastic DSF is compared to the exact calculation 42 in the case of a harmonic lattice. The contribution of phonon states is then given by: S(q, ω ph (q)) = q 2 2ω ph (q) exp (-2W (q)), (19) with same notations as previously used in the same section. The obtained agreement, demonstrated in Fig. 4, confirms the validity of our theoretical approach and consequently the convergence of our series development. The DSF profile appears as being continuous because the lattice size is large enough to blur the discreteness of the Fourier space at the scale of Fig. 4, which ranges over 60 BZ. Apart from the S(q, ω) ripple, the perturbation theory captures quite well the main variations of the envelop of S(q, ω ph ). This ripple stems from the dispersion of ω ph and so, its amplitude increases with C to become some peaks for the acoustic phonons in a harmonic lattice 48 . When the sign of C is inverted in our model (C becomes negative), the local minima of S(q, ω ph ) are shifted at the edges of the lattice Brillouin zones instead of being in the middle when C > 0. In Fig. 4, the DSF resonance involved by the two-phonon states is also plotted. It is usually termed multiple-scattering and is one order of magnitude smaller than the one-phonon process. If the anharmonicity of V can not be neglected, the profile of the one-phonon resonance in S(q, ω) shows no qualitative changes compared to the harmonic lattice (see Figs. 4 and 5 (a)). However, the response of the biphonon, which profile is denoted by S(q, ω bi ), appears clearly at higher energy transfer than one-phonon (see Figs. 2 and5 (a)). The resonance associated to the biphonon is one order of magnitude larger than for the unbound phonon states and has same order as S(q, ω ph ). A gap is opened between the biphonon branch and the two-phonon band (Fig. 1 (b)) so the biphonon resonance occurs for an energy transfer which does not match the harmonics of the fundamental phonons. As mentioned above, that energy gap is also called the binding energy of the biphonon 8 . It is also worth noting that the maximum of the S(q, ω bi ) envelop is reached for a scattering vector that differs from the envelop maximum of S(q, ω ph ). This maximum occurs, indeed, at larger q for the biphonon. The S(q, ω bi ) envelop can be evaluated within the perturbation theory (see Figs. 5 (a)) which provides a rather satisfactory approximation, although the S(q, ω bi ) ripple has a larger amplitude than for phonons. To interpret our results, it is needed to dwell on the ripple of the biphonon contribution to the DSF. As for the phonon, this ripple is related to the dispersion of the biphonon branch. In the case where the binding energy of biphonon is not very large, say no more than a few percent of the phonon energy, the different bottoms of the S(q, ω bi ) ripple occur at center of the successive BZ while the tops are reached at the BZ edges.

According to several tests, it is a systematic feature which, in contrast to phonon, does not depend on the sign of C. Moreover, as shown in the insets of Figs. 5 (a) and (b), the minima of S(q, ω bi ) and the maxima of the two-phonon resonance correspond one to one. That may be explained as follows. In our perturbation theory the biphonon states are approximated by the Bloch waves that bear a single on-site excitation, α = 2 in Eq. ( 8).

In case C = 0, since the biphonon binding energy is not very large, these Bloch waves are hybridized with some other Bloch waves, given by Eq. ( 15), and particularly those bearing two on-site excitations α = 1, at different sites. The degeneracy-lifting of the latest states yields the unbound two-phonon band. The ripple of the DSF can not be analyzed through our first order perturbation theory (see Figs. ) because here the width of the two-phonon band is maximum. The smaller the energy gap is the larger the hybridization is, so as a result, one finds a larger contribution of the two on-site excitation Bloch waves into the biphonon eigenstate at the center of each BZ. The Bloch waves with multiple excitations at distinct sites yield a zero response to scattering, S(q, ω) = 0. One deduces that the biphonon response S(q, ω bi ) is minimum when the contribution of the two on-site excitation Bloch waves is maximum, i.e., at the center of each BZ. On the other hand, the unbound two phonons resonance comes to a maximum at the center of each BZ because of the contribution of the Bloch waves with a single on-site excitation α = 2 to the two-phonon eigenstates. Conversely, at the edge of the lattice Brillouin zones, the Bloch wave hybridization is minimum (because the energy gap is maximum) so that the biphonon response is maximum and may even reach the value computed within the perturbation theory (Fig. 5 (b)). Another interesting point that is made clear within the previous discussion is why the S(q, ω bi ) ripple is not shifted when the sign of C is changed, unlike S(q, ω ph ). Indeed, reversing the C sign leaves unchanged the two-phonon band shape, i.e., the band width is still maximum at the BZ center. So it is for the biphonon branch which the dispersion is mainly due to the contribution of the two on-site excitation Bloch waves 31 (the dispersion computed within Eq. ( 14) for α = 2 is much smaller than the effective biphonon band width). So the gap between the biphonon and two-phonon bands is equally unchanged under the inversion of the coupling sign, i.e., the energy gap is still maximum at the BZ edges and minimum at the center. According to our analysis of the S(q, ω bi ) ripple, it makes clear why the maxima and minima of the biphonon DSF are independent of the sign of C. Our argumentation holds provided that the biphonon binding energy remains inferior to a few percent of the fundamental excitations.

When the anharmonicity of V is weak compared to the phonon band width, the biphonon gap closes and a pseudogap forms 47 at the edge of the lattice Brillouin zones. Then the biphonon branch hardly appears nearby the two-phonon band (see Fig. 1 (c)). The lattice is said quasi-harmonic since the spectrum anharmonicity vanishes. There is, indeed, no anharmonic resonances, neither in the density of state nor, accordingly, in the DSF. That may happen for higher order phonon bound states as the hardening of the quartic on-site term in V (Eq. ( 3)) could compensate the softening of the cubic potential in a certain range of energy. Although that seemingly perfect harmonicity, the S(q, ω) ripple may yet betray the nonlinearity of V . In truth, at the BZ center, the DSF associated to the biphonon, referred to as S(q, ω bi ), has same order as the DSF of two-phonon states, which magnitude is one order below the phonon resonance (see Fig. 5(b)). In contrast, the profile S(q, ω bi ) reaches its maxima where the pseudogap opens (see Fig. 1 (c)), i.e., at the BZ boundaries and the ripple alongside each BZ can be seen as a series of peaks. In Fig. 5 (b), the largest peak remains much larger than the two-phonon response and it has same order as the maximum of the phonon DSF profile, S(q, ω ph ). The maxima of the biphonon DSF decrease as the nonlinear parameters tends to zero. However, before to reach that point, we have seen that the anharmonicity of the energy spectrum vanishes as in Fig. 1 (c). Consequently, we propose to dub the S(q, ω bi ) peaks by nonlinear Bragg peaks, since they gather two features that are opposite to the usual Bragg scattering: first, they are inelastic peaks and second, they satisfy the Bragg reflection condition but shifted of a half reciprocal lattice vector (i.e., π/a 0 for the one-dimensional chain). According to Ref. 47, the pseudogap may occur in lattices with a higher dimension so that similar results as those presented here can be expected in different lattices. In spite of the relative simplicity of our model compared to the case encountered in practice, it indicates that the biphonon may still emerge and contribute significantly to the scattering, even though the spectrum seems perfectly harmonic. In Figs. without opening a substantial gap which would permit to identify the anharmonicity in the energy density of state. In fact, the larger is the phonon dispersion the larger is the range of nonlinear parameters that leads to a quasi-harmonicity. In Fig. 6 (a) the profile of S(q, ω)

is plotted for different values of the energy transfer ω, e.g., ω ph for phonon, ω bi for biphonon and ω tri for triphonon. The triphonon resonance is found to reach a significant magnitude that is comparable to phonon and biphonon. Such a high order phonon bound state has been identified in different materials, e.g., TiH 5 and ZrH 6 and HCl 10 . It is worth noting the triphonon spectrum ripple, which originates from the hybridization between the bound and unbound phonon states, as for the biphonon spectrum. This triphonon ripple is shifted of π/a 0 compared to the biphonon DSF profile.

In what follows, we spell out how our theory could be linked to some concrete cases.

To simulate a scattering experiment, we must account for the resolution of the technics. If one assumes that the energy resolution is not sufficient to differentiate the inelastic lattice resonances (as it may be the case, for instance, in the inelastic X-rays scattering), we have to integrate our simulation with respect to the energy transfer ω. The elastic contribution is skipped from that integral so as to distinguish the inelastic part of the DSF. We computed the corresponding integral and reported the result as a dashed line in Fig. 6 (a). The amplitude and the position of the maximum of the DSF integral are related to the nonlinear parameters. In a strictly harmonic lattice, that integral scarcely differs from the phonon response, whereas in Fig. 6 (a) we note a clear difference between the dashed line and the profile S(q, ω ph ). The envelop of the DSF integral has a form that might be fitted by a superposition of two Gaussian functions, centered at different momenta. It is a standard treatment in the interpretation of the experimental S(q, ω) profile (see for instance Fig. 8 in Ref. 13). If the nonlinear Bragg peaks are large (i.e., when the anharmonicity is weak compared to the phonon dispersion but that is yet not negligible) they may emerge in the integral as shown for large q in Fig. 6 (a). In a strongly nonlinear lattice, a large gap separates the biphonon branch from the two-phonon energy band so the nonlinear Bragg peaks shrinks to a ripple. Then it becomes rather difficult to depict the contribution of the nonlinear states in the DSF integral, aside from the shift of the envelop to larger q.

On the other hand, if the energy resolution is sufficient to resolve the energy dispersion of the phonon branch but the accuracy over the scattering vector is not to separate the BZ (as it may be the case, for instance, in the neutron scattering of a powder specimen, since u and v are randomly distributed), then our simulation must account for the contribution of the distinct BZ. As in a one-dimensional lattice the only symmetry is the inversion, we have to sum the factor 2S(q + k 0 , ω(q)) over the reciprocal vector k 0 , where q and ω(q) are both assumed to match the momentum and the energy conservation, respectively. Since S(q + k 0 , ω(q)) decreases exponentially for a large enough value of (q + k 0 ) (see Fig. 6 (a)), the sum is ensured to be finite. The result of our treatment is shown in Fig. 6 (b). We note the three resonances due to the nonlinear states, i.e., the biphonon, the triphonon and the quadriphonon that are separated from the broad resonances of either the phonon or the unbound phonons. Here, we must enhance that whether the inelastic scattering is coherent or incoherent, the lattice resonances occur at same energy transfer. Thus, from that point of view it is relevant to compare our simulation to the experimental spectrum of incoherent INS in metal hydrides. Actually, the model parameters in Fig. 6 (b) has been adjusted so as to obtain the same energy resonances for the phonon, biphonon and triphonon as in the spectrum shown in Fig. 2 of Ref. 49. The energy unit is Ω ≈ 107meV , the dimensionless coupling is C = 0.06 and the nonlinear parameters are A 3 = 0.17 and A 4 = 0.0231. In that manner, the on-site nonlinearity and thus the energy landscape of the light particles, here the deuterium, may be worked out to a better accuracy than in a purely quadratic model Although our work is concerned with a one-dimensional lattice, the results enhanced in the present study should hold in higher dimension where the DSF depends on the lattice symmetries and polarizations. The extension to higher dimension could be achieved in our perturbation theory with no particular difficulty. This approach proves adequate under the condition that the intersite coupling is very weak, as expected in the metal hydrides as NbH and TaH, from the study of J. Eckert et al. [START_REF] Eckert | [END_REF] . Our numerics, more accurate, could also be extended to higher dimension but that would require either to use a powerful computer or else a dramatic restriction on the site number, which could yet be relevant for some molecules as benzen 33 .

IV. CONCLUSION

The dynamical structure factor (DSF) of the nonlinear Klein Gordon chain (KG) has been calculated for different strength of the on-site anharmonicity. This has been possible thanks to our numerics that permit to diagonalized accurately the non-quadratic Hamiltonian 31 .

The DSF has been expanded as a Taylor series of the atomic displacements, which avoids the use of the conventional harmonic approximations. We found that the on-site nonlinearity leads to phonon bound states and consequently to some anharmonic resonances in the DSF spectrum, which somehow confirms other works on different quantum lattices 30,[START_REF] Agranovich | [END_REF]44 . The treatment of the intersite coupling in a perturbation theory proved satisfactory to tackle a lattice which nonlinearity is stronger than the dispersion. In contrast, when the interaction between first neighbors dominates the on-site nonlinearity, which remains however non-zero, the amplitude of the dynamical structure factor makes visible the biphonon, although the energy transfer is almost harmonic. In such a case, the variation of the biphonon DSF with respect to the transfer momentum q exhibits some peaks that are much larger than the multiple-scattering. This would hint that, provided that the DSF could be resolved accurately in q, the nonlinear behavior of certain materials could be worked out, even though their spectra show no apparent energy anharmonicity.

In certain metal hydrides whose hydrogen anharmonicity is significant, we proposed a scheme to work out from the INS spectra, the potential landscape of the interstitial gaps.

To approach some concrete cases, the theory requires, thought, further developments to consider, for instance, the three dimensionality of a realistic sample or the incoherent scattering.

In a first step, this could be achieved within our perturbation theory. In addition to the possibility of simulating the incoherent INS in certain hydrides and molecular crystals, as those quoted in the present paper, it might be worth studying the nonlinear surface modes that could be investigated practically by different technics, e.g., electron scattering or infrared adsorption. The low dimensionality of the surface could give an advantage to achieve a simulation that would be physically accurate, particularly upon the lattice geometry and various polarizations. As an example, we note the case of the CO molecules adsorbed on a Ru(111) surface that has been studied both theoretically and experimentally 53,54,55 . As a low dimensional system, a vicinal surface, which might exhibit a regular one-dimensional nano-structure over several micrometers, would be the ideal substrate to explore the coherent scattering of a appropriate nonlinear surface modes. The behavior of molecules adsorbed on a vicinal surface presents some specificity (see for instance Refs.56,57 and references therein)

and thus the feasibility of such a study remains under question. Though, the aim would be to carry out a momentum resolved experiment, as those attempted on the PtCl ethylene diamine chlorate 20,58 , to check whether the predictions established within the KG model are confirmed in practice, especially about the biphonon contribution to the DSF. Another possibility to test our theory, would be on the O 2 solid 20,50 where the oxygen cross-section is mainly coherent. Although, in the β-phase of the O 2 crystal, the inter-molecular coupling overpasses the anharmonicity of the O 2 stretch 50 , our study showed that in such a case the inelastic scattering due to biphonon could yet be larger in magnitude than the two-phonon response. Further, the phonon bound states could emerge at higher energy transfer as found for the triphonon and quadriphonon in metal hydrides 49 . At least from a theory viewpoint, the β-phase of the O 2 solid might be worth reexamining with recent experimental devices.

Finally, according to the authors of Ref.3, the potential landscape of hydrogen in PdH, that has been worked out from the INS spectra 3 agrees well with a first-principles calculation made independently 59 . On the basis of that successful comparison, one may expect that a standard first-principles theory could be used to calibrate the KG parameters so as to simulate ab-initio the dynamical structure factor. The X axis, which unit is π/a 0 bears the scattering vector q over a range of 60 Brillouin zones.

Fig. 4: A plot of the coherent inelastic S(q, ω) profile versus q = [q.u], for phonon in a harmonic lattice: C = 0.05, A 3 = A 4 = 0. The thin dashed line has been obtained by a perturbation theory. The thin solid line corresponds to our numerical treatment and the thick dashed line to a plot of the exact formula 42 (the two latest curves are hardly distinguished in the graph). In the inset, the area in the dot-dashed rectangle is magnified.

The X axis unit is π/a 0 . Fig. 5: A plot of the coherent inelastic S(q, ω) profile versus q = [q.u], for phonon, perturbation theory and the solid lines to our numerical treatment. The X axis unit is π/a 0 and the scattering vector q ranges over 60 Brillouin zones. The inset shows a magnification of the two-phonon contributions.

Fig. 6: In (a), the same as in Fig. 5(a) but for different model parameters: C = 0.06, A 3 = 0.17 and A 4 = 0.0231. In addition, have been plotted the S(q, ω) profile for the triphonon and the integral of S(q, ω) over the energy transfer ω (dashed line). The elastic contribution has been skipped. In (b), a plot of the profile of the S(q, ω) integral over the reciprocal lattice (see the text), versus the energy transfer ω, for same parameters as in (a). The inset in (b) shows the corresponding potential landscape V (X j ) and its quantum levels. Our energy unit Ω has been fixed to 107 meV.

  Fig. 1 (b) and Ref. 31) because here the width of the two-phonon band is maximum. The

1

  (c) and 5 (b), we have chosen intentionally a set of parameters that yields a eigen-spectrum which is almost harmonic. A pseudogap hardly opens between the biphonon branch and the two-phonon band. The parameter A 3 could be doubled [as in Fig.2(c) inRef.47]

(

  see the inset in Fig.6 (b)). The same simulation could have been achieved in other alloys as TiH and ZrH[START_REF] Kolesnikov | [END_REF]6 . In their earlier studies, A.I. Kolesnikov and co-workers[START_REF] Kolesnikov | [END_REF]6,49 recognized the phonon bound states resonances in the INS spectra of TiH-D and ZrH-D. To analyse their experimental spectra, the authors compared their data to a simulation made with a model proposed initially by V.M. Agranovich 8 , i.e., a one-dimensional boson Hubbard lattice. In such a model it is, though, difficult to figure out the potential landscape of the interstitial gaps. The single dimension of the model was shown to be reasonable because of the strong anisotropy of the H-H interaction along the metal c-axis. Since we are studying the coherent DSF, a accurate simulation of the spectra obtained by Kolesnikov et al. for the incoherent INS is out of purpose. We shall attempt such a exercise in a future work, devoted to the computation of the incoherent DSF. The amplitudes of various resonances as well as their widths might then be inferred. In Fig. 6 (b), our simulation has been realized with a typical energy resolution function, i.e., a triangle function centered at ω with a width of 0.02ω, which corresponds to the crystal analyser TFXA, ISIS Rutherford Appleton Laboratory, described in Ref. 6. The dynamical trajectory of the TFXA machine has not been simulated in Fig.6 (b). It is noteworthy that the amplitude of the nonlinear resonances has same order as for phonons. Similar calculations with different model parameters even showed that, for a stronger anharmonicity the nonlinear resonances may be even larger in magnitude. The reason for such a behavior is that the density of state enhances as the eigenstates band width decreases. This effect is sharp around the narrow branches as those of phonon bound states 51 (see Fig.1 (b)) which explains why the biphonon resonance may dominate the phonon one. The increase in the density of state may be recognized too in the single phonon response because the phonon branch is flat at the edge and at the center of the lattice BZ (see Fig.1 (a)). Thus one sees two side-bands in the phonon resonance, at the upper and lower boundary in energy. These side-bands may be identified as phonon wings. According to our reading of the above cited works, the experimental INS spectra of the metal hydrides do not exhibit that property but, once again, one must left that point aside, until we study the incoherent DSF. On the other hand, it is worth noting that the form of the lattice response around the first and the second overtone energy regions, in Fig.6(b), is similar either to the infrared adsorption spectra 7 in crystals as CO 2 , N 2 O and OCS or to the Raman spectrum of H 2 solid 11 . Indeed the emergence of a sharp peak, associated to a biphonon (or to a bivibron 11 ) occurs near a small hill-like resonance that is due to the unbound phonon states. In our simulation, a gradual variation of the model parameters so as to decrease the anharmonicity shows the same behavior as the pressure-induced bound-unbound transition, at 25 GPa in H 2 solid 11 , or at 34 GPa in D 2 solid52 . Around that transition, the biphonon (or bivibron) peak broadens and decreases in magnitude. In our model, the pressure variation can be simulated by a change of the coupling parameter C due to the fact that neighbouring molecules are moved closer together because of the external pressure.

Fig. 1 :

 1 Fig. 1: Energy spectrum of a chain composed of N unit cells, N = 33 (full circles) and N = 23 (empty circles): (a) the optical phonon branch for model parameters C = 0.05, A 3 = 0.12 and A 4 = 0.01; (b) the two-phonon energy region for the same parameters as (a); (c) the same as in (b) but for C = 0.05, A 3 = 0.05 and A 4 = 0.01. The Y axis unit is Ω. The wave vector is reported on the X axis, whose unit is (π/a 0 ) and ranges over the lattice first Brillouin zone.

Fig. 2 :

 2 Fig.2: (Color online) A 3D plot of the inelastic S(q, ω) as a function of the dimensionless energy transfer ω and the scalar product q = [q.u] of the transfer momentum q and the polarization u. The parameters are the same as in Figs.1 (a) and (b). The figure is easier to depict in color.

Fig. 3 :

 3 Fig.3:A plot of the Debye-Waller factor [S(q, ω = 0)/N] versus the scalar product q = [q.u], for different model parameters. Each symbol corresponds to a Bragg peak. Our numerical results are reported for C = 0.05, A 3 = A 4 = 0 (diamonds) and C = 0.05, A 3 = 0, A 4 = 0.2 (triangles). The exact formula derived in a harmonic model (dashed line) confirms the former while the latter has been computed within a perturbation theory (circles) too.

  biphonon and two-phonon states. The model parameters are: (a) C = 0.05, A 3 = 0.12 and A 4 = 0.01; (b) C = 0.05, A 3 = 0.05 and A 4 = 0.01. The dashed lines correspond to our

FIG. 1 :

 1 FIG. 1: (2005) L. Proville

  FIG. 2: (2005) L. Proville

Acknowledgments

I gratefully acknowledge Robert S. MacKay who invited me to attend a neutron scattering experiment 58 at Rutherford Appleton Laboratory, which motivated the start of the present study. I also thanks M. Guttmann who pointed the metal hydrides as a possible concrete example for the theory.