K C Rule

G Deng

K Anand

A R Wildes

M. D J R Stewart

R A Le

Ewings

The magnetic dynamics of NiPS 3

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

NiPS 3 belongs to a family of magnetic layered van der Waals compounds. The family has been investigated relatively intensely since the late 1960s. Initial investigations focused on their low-dimensional magnetic behaviour [START_REF] Brec | Review on structural and chemical properties of transition metal phosphorus trisulfides MPS3[END_REF]. They then attracted considerable attention for the ability to intercalate other species into the host [START_REF] Grasso | Low-dimensional materials: The MPX3 family, physical features and potential future applications[END_REF]. More recently, the discovery of graphene has led to renewed interest in the family as they can be delaminated to monolayer thickness while retaining their magnetism [START_REF] Park | Opportunities and challenges of 2D magnetic van der Waals materials: magnetic graphene?[END_REF].

The family members consist of a magnetic transition metal with a 2+ ionization state, in particular Mn, Fe, Co and Ni, and the sulfur can be substituted for selenium. The crystal structures for all the thiophosphate compounds share the same C 2/m space group [START_REF] Ouvrard | Structural determination of some MnPS3, FePS3, CoPS3, NiPS3, CdPS3 layered phases[END_REF], with the transition metal atoms forming a planar honeycomb lattice in the ab planes. The metal atoms are enclosed in octahedra formed by the sulfur atoms and there is a phosphorus doublet at the centre of the honeycomb hexagons, * wildes@ill.fr with the ensemble forming layers that stack along the c axis with weak van der Waals bonding. All the compounds are antiferromagnetic, although they have different magnetic structures and anisotropies [START_REF] Brec | Review on structural and chemical properties of transition metal phosphorus trisulfides MPS3[END_REF]. NiPS 3 magnetically orders below T N ∼ 160 K with the magnetic structure shown in figure 1 [START_REF] Wildes | Magnetic structure of the quasi-two-dimensional antiferromagnet NiPS3[END_REF] along with the corresponding Brillouin zone. The magnetic moments are collinear, lying mostly along the a axis with a small component normal to the ab planes. The structure consists of ferromagnetic zig-zag chains parallel to the a axis that are antiferromagnetically coupled along the b axis, resulting in the magnetic propagation vector k M = [010]. Its paramagnetic susceptibility is isotropic [START_REF] Wildes | Magnetic structure of the quasi-two-dimensional antiferromagnet NiPS3[END_REF], suggesting that the magnetic Hamiltonian is Heisenberg-like. The Mermin-Wagner theorem states that a Heisenberg magnet cannot form long-ranged order in two dimensional materials at finite temperatures, however NiPS 3 has the highest Néel temperature of all the members in the thiophosphate family suggesting that interplanar magnetic interactions may be more important than its sister compounds. This is also suggested by recent Raman spectroscopy measurements of NiPS 3 showing that the compound retains its long-ranged magnetic structure when thinned to two layers, but that the order is lost when ACCEPTED MANUSCRIPT / CLEAN COPY FIG. 1. The magnetic structure of NiPS3 [START_REF] Wildes | Magnetic structure of the quasi-two-dimensional antiferromagnet NiPS3[END_REF], viewed along the (a) c * axis, and (b) b axis respectively. The figure was created using the VESTA software package [START_REF] Momma | VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data[END_REF]. (c) the Brillouin zone for the magnetic structure of NiPS3. Selected highsymmetry points are indicated, corresponding to the magnetic space group PC21/m (#11.57) as defined in the Bilbao Crystallographic Library [START_REF] Aroyo | Crystallography online: Bilbao Crystallographic Server[END_REF][START_REF] Aroyo | Bilbao crystallographic server: I. databases and crystallographic computing programs[END_REF][START_REF] Aroyo | Bilbao Crystallographic Server. II. Representations of crystallographic point groups and space groups[END_REF].

thinned to a monolayer [START_REF] Kim | Suppression of magnetic ordering in XXZ-type antiferromagnetic monolayer NiPS3[END_REF].

Within the flurry of recent activity concerning the family, there has been particular interest in NiPS 3 . It has the smallest spin of the family with S = 1, hence may host quantum effects, exotic magneto-lattice and magneto-electron behaviour [START_REF] Belvin | Exciton-driven antiferromagnetic metal in a correlated van der Waals insulator[END_REF][START_REF] Hwangbo | Highly anisotropic excitons and multiple phonon bound states in a van der waals antiferromagnetic insulator[END_REF][START_REF] Kang | Coherent many-body exciton in van der waals antiferromagnet NiPS3[END_REF]. Detailed knowledge of the spin dynamics is required to understand the magnetic properties of the compound. A variety of highresolution spectroscopy methods have been used to probe the dynamics of NiPS 3 , including recent efforts with light [START_REF] Kim | Suppression of magnetic ordering in XXZ-type antiferromagnetic monolayer NiPS3[END_REF][START_REF] Belvin | Exciton-driven antiferromagnetic metal in a correlated van der Waals insulator[END_REF], optical pump-probe methods [START_REF] Afanasiev | Controlling the anisotropy of a van der Waals antiferromagnet with light[END_REF] and electron spin resonance [START_REF] Mehlawat | Low-energy excitations and magnetic anisotropy of the layered van der Waals antiferromagnet Ni2P2S6[END_REF]. These methods give information on the dynamics at the Brillouin zone centre. Neutron spectroscopy is the method of choice to probe the dynamics throughout the Brillouin zone, and it has been used to measure the spin waves from a sample of powdered NiPS 3 [START_REF] Lançon | Magnetic exchange parameters and anisotropy of the quasi-two-dimensional antiferromagnet NiPS3[END_REF]. The study concluded that the spin dynamics could be modelled using a Heisenberg Hamiltonian with a uniaxial single-ion anisotropy. The nearest-neighbour exchange was found to be ferromagnetic, and the antiferromagnetic structure was stabilised by third in-plane nearest neighbour exchange which was shown to be by far the dominant term in the Hamiltonian. The singleion anisotropy gave rise to an energy gap at the Brillouin zone centre, and would further stabilise the long-ranged magnetic order.

There were a number of ambiguities and assumptions present in the powder data and their analysis. The modelling assumed that the compound was purely twodimensional, with no exchange between the ab planes. The analysis found that the low energy (∆E ≤ 10 meV) spin waves had almost equal energies at the Brillouin zone centre at Γ and the corners at C. This is ambiguous in the analysis of powder data as certain Γ and C points have approximately the same momentum magnitude, Q. Crystals of the compounds are known to twin with 120 • rotations of the ab planes about the c * axis [START_REF] Matsumoto | Nonreciprocal magnons due to symmetric anisotropic exchange interaction in honeycomb antiferromagnets[END_REF], resulting in the scattering from magnetic Γ points being potentially superimposed onto the C points. Finally, the magnitude of a spin wave gap at the Γ point was not clearly determined with the best estimate arriving at ∆E ∼ 7 meV. This is somewhat at odds with more recent measurements, some of which show that there may be two gaps in the low-energy modes at the Γ point [START_REF] Afanasiev | Controlling the anisotropy of a van der Waals antiferromagnet with light[END_REF][START_REF] Mehlawat | Low-energy excitations and magnetic anisotropy of the layered van der Waals antiferromagnet Ni2P2S6[END_REF]. Two gaps would be consistent with an easy-plane anisotropy [START_REF] Kim | Suppression of magnetic ordering in XXZ-type antiferromagnetic monolayer NiPS3[END_REF], contrasting with the analysis of the neutron powder data that used a uniaxial anisotropy [START_REF] Lançon | Magnetic exchange parameters and anisotropy of the quasi-two-dimensional antiferromagnet NiPS3[END_REF].

This article describes a neutron spectroscopy study of the spin dynamics of single crystals of NiPS 3 to resolve the ambiguities in the analysis of the powder data. The data were measured using neutron three-axis and direct-geometry time-of-flight spectrometers, and they were analysed by linear spin wave theory using the SpinW software package [START_REF] Toth | Linear spin wave theory for single-Q incommensurate magnetic structures[END_REF]. The discussion reconciles the neutron data with measurements of the spin wave energies at the Γ point using other experimental techniques and outlines limitations in using a simple Heisenberg Hamiltonian with straight-forward single-ion anisotropies for the data analysis.

II. EXPERIMENT

A. Sample preparation

Crystals of NiPS 3 were prepared using a vapour transport method that has previously been described in detail [START_REF] Wildes | Magnetic structure of the quasi-two-dimensional antiferromagnet NiPS3[END_REF]. A total of 8 of the largest crystals with a total mass of ≈ 0.5 g were chosen for the neutron experiments, being platelets with typical dimensions 5 × 5 × 0.2 mm 3 .

This family of compounds is known to be prone to characteristic twinning involving the rotation of the a and b axes by 120 • about the c * axis [START_REF] Murayama | Crystallographic features related to a van der Waals coupling in the layered chalcogenide FePS3[END_REF], which is normal to the platelets. The twinning results in diffraction patterns that appear to have 6-fold symmetry. The degree of twinning in each of the crystals was determined by comparing the integrated intensities of neutron diffraction peaks that do not have something approaching 6-fold symmetry, such as {022} or, in the magnetically ordered state, {010}, for the three 120 • rotations.

The magnetic properties of NiPS 3 are affected if crystals are glued directly to supports in preparation for a measurement [START_REF] Wildes | Magnetic structure of the quasi-two-dimensional antiferromagnet NiPS3[END_REF], possibly due to strong magnetostric-ACCEPTED MANUSCRIPT / CLEAN COPY tion or intercalation. Consequently, all the crystals were wrapped in aluminium foil prior to mounting, and glue was never applied directly to a crystal.

B. Three-axis spectrometry

Neutron three-axis spectrometry was used to study the low energy dispersion from the best crystal. The crystal had dimensions ∼ 10×10×0.3 mm 3 and a domain ratio of 0.58:0.26:0.16. The foil-wrapped sample was glued to an appropriate support using a small amount of GE varnish. Two supports were used: one to access the 0kl scattering plane, and one to access hk0, defined with respect to the dominant domain. The sample was simply removed from its aluminium envelope between experiments, then wrapped in a new envelope before being realigned and mounted on the new support.

Two three-axis spectrometers were used.

The IN8 instrument at the Institut Laue-Langevin, France, was used to study the scattering in the 0kl plane [START_REF] Wildes | The magnon dispersion in single crystals of NiPS3[END_REF]. The instrument was configured with a pyrolytic graphite (PG) 002 monochromator and analyser. The wavenumber for the analyser was fixed to k f = 2.662 Å-1 for the experiment. Higher order contamination was suppressed using a PG filter before the analyser. Both monochromator and analyser were vertically focussed and horizontally flat, and 40 ′ collimators were used before and after each of them. A liquid helium cryostat was used to cool the sample to 1.8 K.

The Taipan instrument at the Australian Nuclear Science and Technology Organisation, Australia, was used to study the scattering in the hk0 plane. The instrument was configured in the same manner as IN8, with the exceptions that 40 ′ collimators were only used before and after the sample and that temperature control down to 5 K was achieved using a closed-cycle refrigerator.

The resolution on both instruments was calibrated by mapping the elastic scattering from a nuclear Bragg peak, and by measuring the energy width of the incoherent scattering at an equivalent position in Q. The 200 peak from single crystals of MnPS 3 and NiPS 3 were used to calibrate IN8 [START_REF] Wildes | The nature of the anisotropy[END_REF] and Taipan respectively. The measured energy full-width half-maximum of the elastic line was 0.79(3) meV on IN8 and 0.82(3) meV on Taipan. The scattering was subsequently compared to calculations using the Popovici method [START_REF] Popovici | Resolution of slow-neutron spectrometers. 4. Triple-axis spectrometer resolution function, spatial effects included[END_REF] in the RESCAL5 library [23] for MATLAB R . The comparison was satisfactory using the nominal instrument parameters, and a sample mosaic of 25 ′ .

C. Time-of-flight spectrometry

Time-of-flight spectrometry experiments were conducted using the MERLIN direct geometry spectrometer at the ISIS spallation neutron source, UK [START_REF] Bewley | MERLIN, a new high count rate spectrometer at ISIS[END_REF].

These measurements [START_REF] Wildes | The spin wave scattering from the quasi-two dimensional antiferromagnet[END_REF] gave an overview of the dynamic structure factor over the entire Brillouin zone. All eight crystals were combined and co-aligned for the experiment. The crystals were fixed to aluminium plates using a small amount of GE vanish and some Fomblin R oil. The c * axis was normal to the plates and the samples were co-aligned such that the 0kl plane was horizontal, defined by the majority domain for the individual crystals. The plates were then mounted on a common "toaster rack" aluminium block. Temperature control to 5 K was achieved using a cryorefrigerator.

The scattering was measured for a series of 1 • rotations about the normal to the 0kl plane. A total of 84 rotations were measured. The scattering could then be combined into a 4-dimensional volume and cut or sliced along desired reciprocal space directions using the HO-RACE software [START_REF] Ewings | HORACE: Software for the analysis of data from single crystal spectroscopy experiments at time-of-flight neutron instruments[END_REF]. MERLIN can profit from rep-rate multiplication, whereby the scattering from multiple incident neutron energies can be measured in a single neutron pulse [START_REF] Russina | First implementation of Repetition Rate Multiplication in neutron spectroscopy[END_REF][START_REF] Nakamura | First demonstration of novel method for inelastic neutron scattering measurement utilizing multiple incident energies[END_REF]. Data with E i = 16.9, 30.4 and 70 meV were simultaneously collected. The resolution for the three incident energies had a full-width half maximum at ∆E = 0 of 0.5, 1.1 and 3.4 meV respectively.

III. RESULTS

A. Neutron three-axis spectroscopy, and the low energy spectra

Three-axis scans show details that were not apparent in the analysis of the data from powdered samples [START_REF] Lançon | Magnetic exchange parameters and anisotropy of the quasi-two-dimensional antiferromagnet NiPS3[END_REF]. Figure 2 shows energy scans at Q = 120, corresponding to Γ in figure 1c, and at Q = 3 2 1 2 0 and 3 2 1 2 0 corresponding to C points in Brillouin zone corners. The data at the Γ point shows a clear increase at ∆E ∼ 6 meV, while a similar increase is seen at ∆E ∼ 9 meV for the data at the C points. These features were expected and correspond to gapped modes observed in the powder data. The gap at the Γ point is consistent with terahertz spectroscopy [START_REF] Belvin | Exciton-driven antiferromagnetic metal in a correlated van der Waals insulator[END_REF] and optical pump-probe [START_REF] Afanasiev | Controlling the anisotropy of a van der Waals antiferromagnet with light[END_REF] measurements that showed features at 5.3 meV and 3.8 meV respectively. The intensity does not decrease at larger energy transfer as the line shape is highly asymmetric due to the convolution of the instrument resolution with the very steep spin wave dispersion.

Close inspection at smaller energy transfers shows that the count rates for ∆E ≤ 5 meV are systematically larger at Γ than at either of the C points, hinting at extra spectral weight. This is due to a second, lower energy spin wave mode. Figure 3 shows constant energy cuts at ∆E = 2 meV, which is approximately the lowest energy transfer where data can be confidently taken to be clean from contamination due to the elastic scattering. Scans along h20 and 1k0 are shown in figure 3 (a) and (c) respectively and a weak, but clear and consistent, peak is observed at the 120 Brillouin zone centre in each. Similar scans through the C points did not consistently have The weak mode is steeply dispersive along a * and b * , which may be expected, however it is also surprisingly dispersive along c * considering that NiPS 3 was thought to be quasi-two dimensional. Figure 3 (e) and (f) show scans along 0k1 and 01l. Figure 3 (c) and 3 (e) follow the same trajectory through the Brillouin zone and may be compared directly, both showing weak peaks at the respective Brillouin zone centres. A peak is also observed at 011 in figure 3 (f). A true two-dimensional layered compound would have no magnetic exchange between the planes. Any spin waves normal to the planes would be dispersionless, and scans in this direction would show no peaks. That there is a peak in figure 3 (f) shows that there is a magnetic exchange between the layered planes in NiPS 3 which appears to be considerably larger than that for the sister compounds MnPS 3 [START_REF] Wildes | Spin waves and the critical behaviour of the magnetization in MnPS3[END_REF] and FePS 3 [START_REF] Lançon | Magnetic structure and magnon dynamics of the quasi-two-dimensional antiferromagnet FePS3[END_REF].

Similar scans at larger energy transfers give an indication of the dispersion of the weak mode. Q scans along h 1 2 0, h20, 1k0 and 01l are shown in figures 4 (a)-(d) respectively. Figure 4 (a) shows that no clear peaks are observed in the Brillouin zone corner until an energy transfer of ∆E ∼ 8 meV, however peaks are present in all the scans at the 120 and 011 Brillouin zone centres. The peak intensities are very weak until ∼ 6 meV, where they increase dramatically.

Dispersive modes would give two peaks in the constant energy cuts shown in figure 4. Unfortunately, the instrumental resolution is too coarse to separate the two peaks, however the fitted widths of a single peak function give an indication of the dispersion. Gaussian fits to the peaks are shown in the figure, and the Q widths in h, k and l widths are plotted as a function of energy transfer in figure 5 (a)-(c). Data determined from fits around the 120 and 011 positions are shown as black and red points respectively. Estimates for the instrumental resolution are also shown as black and red dashed lines, correlated to the data points of the same colour. The estimates give the expected widths if the mode was a delta function in Q, corresponding to a spin wave that was infinitely steep. The data at small energy transfers are close to the resolution estimates, but increase approximately linearly with energy, as expected from linearly dispersive excitations.

The measurements had insufficient resolution to determine whether the lowest energy mode at the Γ point was gapped. A dedicated measurement with a cold neutron spectrometer would be required for this task. However, it is highly likely to be gapped as the ordered moments are ultimately collinear with a well-defined crystallographic axis, suggesting a non-zero magnetic anisotropy. Measurements using optical pump-probe [START_REF] Afanasiev | Controlling the anisotropy of a van der Waals antiferromagnet with light[END_REF] and electron spin resonance [START_REF] Mehlawat | Low-energy excitations and magnetic anisotropy of the layered van der Waals antiferromagnet Ni2P2S6[END_REF] techniques have found evidence for a low-energy gap, estimated to be 1. [START_REF] Bewley | MERLIN, a new high count rate spectrometer at ISIS[END_REF]

B. Time-of-flight spectroscopy

The MERLIN measurements performed with E i = 16.9 meV and 30.4 meV confirm the presence of the low energy mode observed in the three axis data. Examples of the scattering for a selection of slices through the fourdimensional (Q, ∆E) data are shown in figure 6. The scattering was very weak, and statistics were improved by slicing and combining the data along symmetrically identical directions. In particular, data at the same magnitude of l were combined.

Figures 6(a), (d) and (g) show the nominally elastic (∆E = 0) scattering, measured with E i = 16.9 meV, in the (hk1), (h1l) and (0kl) planes respectively. Nuclear Bragg peaks in the C 2 m space group have systematic absences for h + k = 2n, where n is an integer. Clear magnetic Bragg peaks, representing Γ points in figure 1(c), are visible at the expected [010] propagation vector positions, such as 011 and 121. Weaker peaks are also visible at the 1 2 1 2 1 positions in figures 6(a). These are due to 011 magnetic Bragg peaks from the 120 • rotational domains. They are not exactly in the scattering plane as c * is not orthogonal to hk1 and a 120 • rotation about this axis will map 011 from one of the domains onto a lattice position with l = 1 6 or 5 6 in the figure. However, the peaks are extended along the c * axis due to their quasi-two-dimensional nature, as is apparent in figures 6(d) and (g). Furthermore, the data in figure 6(a The energy integration limits were chosen to be as large as possible while avoiding contamination from either elastic scattering or contributions from the stronger modes at higher energies. Very weak intensity is observable at the Γ points 011 in all the figures, and arguably at 012 in figures 6(e) and (h). The peaks are very sharp in Q in all three figures, reflecting the dispersion of the spin waves in three dimensions. No extra intensity is observed at any of the Brillouin zone corners, or at the 1 are only present at the Brillouin zone centres and that the instrumental resolution is sufficient that the slices are free of similar contributions from crystal twins.

Higher energy slices, integrating over 9 ≤ ∆E ≤ 11 meV, are shown in figures 6(c), (f) and (i). Figure 6(c) shows data at this energy that are consistent with those measured using three-axis spectroscopy, with clear intensity spots appearing at the Γ point at Q = 011 and at the C points at Q = ± 1 2 5 2 1. The extent of the dispersion along c * is apparent in figures 6(f) and (i), where the sharp spots at lower energy have become streaks along (01l) due to the appearance of the more intense and larger-gapped mode. The streaks show that this energy is close to the maximum in the dispersion along this direction. The dispersion is apparent in figure 6(j), which shows a (Q, ∆E) slice along (01l) measured with E i = 30.4 meV neutrons. The higher- intensity mode is the most clearly visible, with minima at l = 1 and 2 and an energy of ∼ 13 meV at the B point at l = 3 2 . Individual modes are not observed as the slice was integrated over ±0.1 rlu in both h and k, and the dispersion in these directions is so steep that the spin waves in these directions are also included in the integration.

(Q, ∆E) slices through the data taken with E i = 70 meV show the dispersion along high-symmetry directions at higher energies. Figures 7 show representative slices, with (a), (d) and (g) showing slices parallel to b * for selected values of h and integer values for l and (j) showing a slice along (04l). Qualitative inspection showed that slices in equivalent reduced lattice units were identical, and the data were subsequently combined for better statistics.

Inspection of the data shows the presence of a third ACCEPTED MANUSCRIPT / CLEAN COPY mode at the Γ point, in addition to the low-energy modes shown in the three-axis measurements and the data in figure 6. It appears most clearly as a minimum at ∆E ∼ 34 meV in a steeply dispersive mode at k = 4 in figure 7(a) and at k = 5 in figure 7(g). A slice along (04l) reveals that this mode is dispersionless within the resolution of the measurement, as shown in figure 7(l). On the other hand, while the spectral weight varies as a function of the reduced lattice units, only one clear mode is observed in the slices along the Y-C Brillouin zone boundary, as shown in figure 7(d).

The data indicate that there are at least three spin waves in the compound and that some branches are close to, or potentially exactly, degenerate. The observation is consistent with linear spin wave theory which predicts four spin waves, one for each magnetic atom in the unit cell, some of which may be degenerate.

IV. ANALYSIS WITH LINEAR SPIN WAVE THEORY

Within the resolution of the measurements, the spin dynamics of NiPS 3 consists of well-defined spin waves and can be analysed using linear spin-wave theory with an appropriate Hamiltonian. NiPS 3 has four moments in the magnetic unit cell, and in principle this results in the Hamiltonian being expressed as an 8×8 matrix. Previous analysis of neutron inelastic scattering from a powdered sample [START_REF] Lançon | Magnetic exchange parameters and anisotropy of the quasi-two-dimensional antiferromagnet NiPS3[END_REF] used a Heisenberg Hamiltonian with a singleion anisotropy,

Ĥ = 1 2 ij J ij S i • S j + D i (S x i ) 2 , ( 1
)
where J ij is the exchange interaction between neighbours i and j, and D is the anisotropy strength. The previous analysis assumed a uniaxial anisotropy with the x axis collinear to the aligned moment directions. This simplifies the Hamiltonian, allowing it to be expressed as a 4×4 matrix which can be analytically diagonalised to give expressions for S (Q, ∆E). Many of the spin wave branches become degenerate. Notably, there are only two doublydegenerate modes at the Γ point, one at low and one at high energy. Also, only one four-fold degenerate mode is present along the Y-C Brillouin zone boundary. Subsequent experiments using Raman spectroscopy [START_REF] Kim | Suppression of magnetic ordering in XXZ-type antiferromagnetic monolayer NiPS3[END_REF] suggested that the compound is better described using an XXZ Hamiltonian:

Ĥ = 1 2 ij J ij S x i S x j + S y i S y j + αS z i S z j + D x i (S x i ) 2 + D z i (S z i ) 2 , (2)
with x again parallel to the aligned moment direction and z normal to both this direction and b. The parameter α permits anisotropy in the magnetic exchange parameters.

This Hamiltonian must be expressed as an 8 × 8 matrix which is not easily diagonalised analytically, but can be solved numerically using the SpinW software package [START_REF] Toth | Linear spin wave theory for single-Q incommensurate magnetic structures[END_REF] written for MATLAB R . The program is also compatible with HORACE [START_REF] Ewings | HORACE: Software for the analysis of data from single crystal spectroscopy experiments at time-of-flight neutron instruments[END_REF] and the combined software can be used to calculate the expected neutron scattering, including a convolution with the instrumental resolution. The more general Hamiltonian in equation 2 lifts the degeneracy of the spin wave branches in a manner consistent with the observation of two low-energy modes at the Γ point, and hence equation 2 was used to analyse the data. Manipulation of the parameters shows that the anisotropy must be dominated by an easy-plane singleion anisotropy largely confined to the ab planes, i.e. by a relatively large and positive D z . The value for D x must be small, and α must be very close to one. Departures from these conditions led to large and dramatic discrepancies with the data, in particular with the appearance of multiple modes along the Y-C Brillouin zone boundaries. A value of α = 1 was subsequently adopted to reduce the number of free parameters for all subsequent analysis.

SpinW was not explicitly used to fit the data as they had relatively poor statistics and a complicated background that made fitting two-dimensional data sets problematic. Instead, data were extracted at given highsymmetry points in the Brillouin zone where peaks due to spin waves were least ambiguous. The peaks were fitted to determine the spin wave energies, and the energies could then be fitted using analytical expressions.

The spin wave energies determined from the data are listed in table I. A total of 7 data points were extracted. The low-energy excitations at the Γ and C points were extracted from the three-axis data. The fits included a flat background and a convolution with the instrumental resolution, and examples are shown in figure 2. The spin wave energy at the B point was extracted from the MERLIN data collected at E i = 30.4 meV, shown in figure 6(j), and energies for ∆E > 15 meV were extracted from the data at E i = 70 meV, shown in figure 7. All the MERLIN data were fitted with Gaussians. Table I also lists their corresponding reciprocal and reduced-reciprocal lattice positions. An additional point was added to account for a gap in the lowest energy mode at Γ position. The energy was estimated to be at 1.1 ± 0.1 meV from the average of the measurements using optical pump-probe [START_REF] Afanasiev | Controlling the anisotropy of a van der Waals antiferromagnet with light[END_REF] and electron spin resonance [START_REF] Mehlawat | Low-energy excitations and magnetic anisotropy of the layered van der Waals antiferromagnet Ni2P2S6[END_REF] techniques.

SpinW is capable of calculating an analytical expression for the spin waves from small unit cells using the MATLAB R symbolic toolbox, but it was not capable of returning a general expression from the 8 × 8 matrix required for NiPS 3 . However, it was able to produce expressions for the eigenvalues at the high-symmetry points. The expressions were calculated in a manner consistent with previous analysis of the transition metal-PS 3 compounds [START_REF] Kim | Suppression of magnetic ordering in XXZ-type antiferromagnetic monolayer NiPS3[END_REF][START_REF] Wildes | Spin waves and the critical behaviour of the magnetization in MnPS3[END_REF][START_REF] Lançon | Magnetic structure and magnon dynamics of the quasi-two-dimensional antiferromagnet FePS3[END_REF], with magnetic exchange up to the third nearest-neighbour in the ab planes and to the nearest neighbours between the planes. The spin wave ACCEPTED MANUSCRIPT / CLEAN COPY The data point marked with an asterisk is an average taken from optical pump-probe [START_REF] Afanasiev | Controlling the anisotropy of a van der Waals antiferromagnet with light[END_REF] and electron spin resonance [START_REF] Mehlawat | Low-energy excitations and magnetic anisotropy of the layered van der Waals antiferromagnet Ni2P2S6[END_REF] measurements. The corresponding eigenvalue attribution is defined with respect to the order in table III in the appendix, and the fitted energy for that eigenvalue is also listed. All energies are in meV. expressions therefore contained only six free parameters: J 1..3 for the intraplanar exchange between first to third neighbours, J ′ for the interplanar exchange, D x for a uniaxial anisotropy along the collinear moment axes and D z for an easy-plane anisotropy, with the spin taken to be S = 1. The expressions are listed in the appendix. The spin wave energies were then assigned to one of the up-to-four modes calculated at the appropriate wavevector position with the corresponding analytical expression, and the ensemble was then fitted with a Levenberg-Marquardt regression algorithm. The fitted energies are listed in table I, and the output parameters are listed in table II. The assigned eigenvalues used for the fitting, defined in table III in the appendix, are also listed in table I. The choice was somewhat arbitrary, especially at higher energy transfers, as many modes are almost doubly-or, at the Y-point, quadruply-degenerate. Choosing another mode made little difference to the final fit parameters.

Q k-
Having relatively few data points with respect to the number of free parameters and the quasi-degeneracy of many of the modes gave unreasonable uncertainties in the FIG. 8. The spin wave dispersions calculated with the fitted exchange parameters for NiPS3 using equation 2 with the parameters listed for the current analysis in table II.

fitted parameters when they were all permitted to be free. Their stabilities were therefore tested by initially fitting with all the parameters free to find the global solution, and then fixing all but one of them in turn for subsequent fits. The parameters did not change in these fits, and the resulting uncertainties are listed in table II.

The fitted parameters are shown in table II, alongside the previously-published parameters determined from neutron spectroscopy on a powdered sample [START_REF] Lançon | Magnetic exchange parameters and anisotropy of the quasi-two-dimensional antiferromagnet NiPS3[END_REF] and from modelling the two-magnon spectra from Raman spectroscopy [START_REF] Kim | Suppression of magnetic ordering in XXZ-type antiferromagnetic monolayer NiPS3[END_REF]. All values are consistent with the Hamiltonian convention used by SpinW.

The exchange parameters determined by neutron spectroscopy on powdered and single crystal samples are consistent showing that J 1 is ferromagnetic, 2 is much smaller and is antiferromagnetic, and that the magnetic dynamics are dominated by a very large antiferromagnetic third nearest-neighbour interaction, J 3 . The magnitude of J 1 is smaller in the current analysis than that determined from the powder data, however this is compensated by the inclusion of a relatively large ferromagnetic exchange J ′ . The magnitude of the anisotropy is approximately the same, but the identification of the new low-energy mode at the Γ point establishes that it is dominantly co-planar with a very small collinear contribution along the aligned moment direction. The spin wave energies throughout the Brillouin zone have been calculated using the new set of parameters and is shown in figure 8.

By contrast, the exchange parameters determined from Raman spectroscopy are all approximately the same magnitude and are antiferromagnetic, and the anisotropies are also considerably larger in magnitude. The spin wave dispersion calculated using these parameters does not match the neutron data, and hence these parameters do not appear to be representative of NiPS 3 .

The energies determined from using linear spin wave theory with the parameters in table II are listed in table I, and they match the data acceptably well. The parameters were subsequently used with the combined SpinW/HORACE software to calculate the neutron scattering both with and without a convolution with the instrumental resolution. The calculations are shown in fig-

ACCEPTED MANUSCRIPT / CLEAN COPY ures 6 and 7, and are directly comparable to the experimental data on the corresponding row. The agreement between measured data and the calculations are mostly satisfactory, confirming that the fitted parameters accurately describe the spin wave dynamics. However, there are some discrepancies, hinting at limitations in the analysis method to describe the scattering.

The clearest discrepancies are visible when comparing the measured data in figure 6(j) with the corresponding calculations in 6(k) and (l). The calculations give the lower energy mode significantly more spectral weight than the higher energy mode at l = 1 and 2, whereas the neutron data in figures 2 and 6 show the opposite. According to linear spin wave theory, the lowest energy mode is due to S y S y fluctuations while the higher energy mode is due to S z S z fluctuations. This is expected as the relatively large positive value of D z favours the moments to lie in the xy planes, hence the in-plane fluctuations will cost less energy than those out-of-plane. The S y S y fluctuations give no intensity at Q = 010 as y b * , and thus they have no component perpendicular to the momentum transfer Q. However, they become readily visible, and more intense than the S z S z fluctuations, with increasing l. That the measured intensities do not match the calculations suggests that the modes are not so simply polarized.

The three-axis data at Q = 120, shown in figure 2, hint at another discrepancy concerning the spectral weight at the Γ point. The data were fitted with a function to determine the energy gap associated the larger of the two low-energy modes, and the fit is shown in the figure .  The fit included a convolution with the estimated resolution function, hence the line shape is asymmetric due to the steepness of the dispersion, which was set from inspection of the MERLIN data. The fit matches the increase of the intensity at ∆E ∼ 6 meV, but decreases more rapidly with increasing energy than the data. The extra intensity in the data at ∆E 10 meV may be due to spurious scattering or contamination from other twinned domains, although an overall view of the experimental results suggest that the latter contribution is minor. Alternatively, the extra intensity may be the result of more exotic effects, such as quantum fluctuations or hybridization.

Other discrepancies include the fitted value for the fitted dispersion along 00l being slightly steeper than observed, as may be seen by comparing figures 6(j) and (k), and indications of spectral weight at ∆E ∼ 50 meV at the Q = 1 2 40 and 1 2 50 positions, hinted at in 7(d). The former may be within the uncertainty of the fit, however, this is not the case for the latter. The calculation in 7(e) shows that some extra intensity at higher energies will result from the instrumental resolution, but not enough to account for the observed data. The Hamiltonian in equation 2 will not give two modes that are separated by ∆E ∼ 10 meV at the Y point without severely compromising the energies at other points in the Brillouin zone. These discrepancies hint at the need to allow for effects beyond the theory used here, or for extra terms in the Hamiltonian.

V. DISCUSSION

there are limitations in the ability of the analysis to describe the measured neutron intensity, linear spin wave theory using equation 2 and the parameters listed in table II allow for significant insight into the spin dynamics of NiPS 3 .

The strength of the interplanar exchange, J ′ is particularly noteworthy. It is the strongest of all the transition metal-PS 3 determined to date, amounting to ∼10% the value of J 1 , hence NiPS 3 is more three-dimensional than its sister compounds. This helps explain the magnitude of the Néel temperature which, at T N ∼ 155 K, is considerably larger than MnPS 3 (78 K), FePS 3 (123 K) and CoPS 3 (120 K), especially when combined with the very strong antiferromagnetic J 3 which, again, is much stronger than any of the exchanges determined in the sister compounds. It also explains the observation from Raman spectroscopy that NiPS 3 maintains a long-ranged magnetic structure when thinned down to two monolayers, but does not order when thinned to a monolayer [START_REF] Kim | Suppression of magnetic ordering in XXZ-type antiferromagnetic monolayer NiPS3[END_REF]. Without a relatively strong interplanar exchange, the magnetic anisotropy is insufficient to stabilise longranged order.

The relatively weak anisotropy is predominantly XYlike, which was also a conclusion from Raman spectroscopy measurements [START_REF] Kim | Suppression of magnetic ordering in XXZ-type antiferromagnetic monolayer NiPS3[END_REF], however a small uniaxial anisotropy must be present to force the moments to align almost parallel to the a axis. The analysis allows for this anisotropy through the inclusion of D x . The fitted magnitude may be compared with observations of a spinflop transition in NiPS 3 when a field is applied in the ab planes, given as B SF ≈ 6 Tesla and B SF ≈ 14 Tesla based on magnetometry [START_REF] Basnet | Highly sensitive spin-flop transition in antiferromagnetic van der Waals material MPS3 (M = Ni and Mn)[END_REF] and 31 P NMR [START_REF] Bougamha | 31 P NMR investigation of quasi-two-dimensional magnetic correlations in T2P2S6 (T=Mn, Ni)[END_REF] respectively. Mean field theory gives the spin flop transition at [START_REF] Blundell | Magnetism in Condensed Matter[END_REF]:

B SF = D x M 2 (2AM 2 -D x ), ( 3
)
where M is the ordered moment, which is is 1.05 µ B from neutron diffraction [START_REF] Wildes | Magnetic structure of the quasi-two-dimensional antiferromagnet NiPS3[END_REF], and

A = - 4zJ ef f N (gµ B ) 2 ( 4
)
where z is the number of nearest neighbours and J ef f is an effective nearest-neighbour exchange. The fitted parameters in table II give D x = -0.01 meV and a weighted mean of J ef f = 2.825 meV. Using these with z = 3 for the honeycomb lattice gives a spin flop field of B SF = 7.3 Tesla, which is close to the values determined from experiment and giving confidence that D x is a believable value.

ACCEPTED MANUSCRIPT / CLEAN COPY

As a final point, the quantitative nature of neutron scattering for studying lattice and magnetic dynamics must be reinforced. Neutron spectroscopy allows access to dynamics over the whole Brillouin zone. The first Born approximation holds very well for neutron spectroscopy, thus the cross-section is proportional to the dynamic structure factor S (Q, ∆E) and the measured intensity distribution may be directly compared with calculations from first principles. The parameters presented in table II represent the best description of the magnetic dynamics to date, however the discrepancies between the measured and calculated intensities for the low-energy modes show that elements are lacking in the relatively straightforward analysis and further work is needed. Any future discussion, though, must be able to produce quantitatively the spin wave energies and neutron intensities apparent in the data presented here.

VI. CONCLUSIONS

Neutron spectroscopy has been used to measure the one-magnon spin dynamics of NiPS 3 . The have been analysed with linear spin wave theory using a Heisenberg Hamiltonian with two single-ion anisotropies. The analysis successfully reproduces the measured spin wave energies, and the resulting exchange parameters and anisotropies represent the best estimates to date for reproducing the dynamical structure factor. The analysis does not quite reproduce the measured neutron intensities, however, particularly for the spin fluctuations at the lowest energies, hinting at other more exotic dynamics in the compound. showing the analytic expressions for the eigenvalues at specific points in the Brillouin zone, as given by SpinW. The parameters D x and D z are taken to be always negative, for uniaxial anisotropy, and positive, for planar anisotropy, respectively. The parameters U and V for the eigenvalues at the Z point are defined in equations A1 and A2 Γ 0 0 0 1: 2S (D x (D x -D z -J1 -4J2 -3J3)) 2: 2S ((D x -D z -J3 + J ′ ) (D x + J1 -2J3 + J ′ )) 2: 2S ((D x -D z + J ′ ) (D x + J1 -3J3 + J ′ ))

1 2 Z 0 1 2 0 1: S 2U -2 (V )
1 2 B 00 1 2 1: 2S ((D x -D z + 4J ′ ) (D x -J1 -4J2 -3J3 + 4J ′ )) 1 2 2: 2S ((D x + 4J ′ ) (D z -D x -J1 -4J2 -3J3 + 4J ′ )) 1 2 3: 2S ((D x -D z + 2J1 -4J2 + 2J ′ ) (D x + J1 -3J3 + 2J ′ )) 1 2 4: 2S ((D x + 2J1 -4J2 + 2J ′ ) (D x -D z + J1 -3J3 + 2J ′ )) 1 2

FIG. 2 .

 2 FIG. 2. Neutron three-axis spectroscopy data measured using Taipan at Q = 120, Q = 3 2

4 FIG. 3 . 2 1 2

 4322 FIG. 3. Neutron three-axis scans performed at an energy transfer of ∆E = 2 meV along different trajectories in reciprocal space. (a) and (c) respectively show scans parallel to the a * and b * axes through the Brillouin zone centre at Q = 120. (b) and (d) respectively show scans parallel to the a * and b * axes through the Brillouin zone corner at Q = 3 2

1 2

 1 FIG. 4. Neutron three-axes scans performed at constant energy transfers at different reciprocal lattice points, together with fits of a Gaussian function where appropriate. Data at equivalent reciprocal lattice points and different energies have been separated vertically for clarity. (a) Scans parallel to the a * axis through the Brillouin zone corner at Q = 3 2

FIG. 5 .

 5 FIG. 5. Fitted Gaussian full-widths-at-half-maximum for neutron three-axes scans performed at constant energy transfers. The scan trajectories are parallel to the (a) a * , (b) b * and (c) c * axes respectively. The black points correspond to measurements from Taipan centred at Q = 120, and red points to measurements from IN8 centred at Q = 011. The black and red dashed lines show the calculated widths for the corresponding black and red points if the dispersion was infinitely steep, corresponding to a delta function in Q .

2 1 2

 22 FIG. 6. Neutron spectroscopy measurements from MERLIN, focusing on the low energy modes and the dispersion along c * . Experimental data at positive and negative l have been combined. (a) -(c) Constant energy slices in hk |1|, with -0.1 ≤ |l| ≤ 0.1 rlu. (d) -(f) Constant energy slices in h1 |l|, with -0.1 ≤ k ≤ 0.1 rlu. (g) -(i) Constant energy slices in 0k |l|, with -0.1 ≤ h ≤ 0.1 rlu. Data in (a), (d) and (g) are centred at ∆E = 0 ± 1 meV; in (b), (e) and (h) at ∆E = 3.5±1.5 meV; and (c), (f) and (i) at ∆E = 10±1 meV. (j) Q, ∆E slice along 01l, measured with Ei = 30.4 meV. (k) Combined SpinW/HORACE calculation of S ⊥ (Q, ∆E) for the same slice shown in (j), accounting for the instrument resolution. The data have had a small constant added (5 counts) to resemble background. (l) SpinW calculation for the same slices as shown in (j) and (k), without the instrument resolution.

FIG. 7 .

 7 FIG. 7. Neutron spectroscopy Q, ∆E slices from MERLIN. All data were collected with Ei = 70 meV. (a) Combined slices from 0k0 and 0k |1|, with corresponding SpinW/HORACE calculations including (b) and without (c) convolution with the instrument resolution. (d) Combined slices from 1 2 k0 and 1 2 k |1|, with corresponding SpinW/HORACE calculations including (e) and without (f) convolution with the instrument resolution. (g) Combined slices from |1| k0 and |1| k |1|, with corresponding SpinW/HORACE calculations including (h) and without (i) convolution with the instrument resolution. (j) Slice from 04 |l|, with corresponding SpinW/HORACE calculations including (b) and without (c) convolution with the instrument resolution with corresponding SpinW/HORACE calculations including (k) and without (l) convolution with the instrument resolution. The SpinW/HORACE plots in (b), (e), (h) and (k) have had a small constant added (5 counts) to resemble background.

1 2 2 : 2 3: 2 4:

 222 2S ((D x -D z ) (D x -J1 -4J2 -3J3)) 1 2S ((D x -D z + 2J1 -4J2 + 2J ′ ) (D x + J1 -3J3 + 2J ′ )) 1 2S ((D x + 2J1 -4J2 + 2J ′ ) (D x -D z + J1 -3J3 + 2J ′ )) ((D x -J3 + J ′ ) (D x -D z + J1 -2J3 + J ′ )) 1 2

  ((D x + J ′ ) (D x -D z + J1 -3J3 + J ′ )) 1 2

TABLE I .

 I Table showing the spin wave energies at specific reciprocal lattice points, extracted from the experimental data.

TABLE II .

 II Tableshowingthe fitted exchange parameters from the current data, along with the previously published results from Lançon et al.[START_REF] Lançon | Magnetic exchange parameters and anisotropy of the quasi-two-dimensional antiferromagnet NiPS3[END_REF] and Kim et al.[START_REF] Kim | Suppression of magnetic ordering in XXZ-type antiferromagnetic monolayer NiPS3[END_REF]. Negative J denotes a ferromagnetic exchange. D x and D z are the magnitudes of the single ion anisotropies along the a and c * axes respectively. A negative D denotes a uniaxial anisotropy, and a positive value denotes a planar anisotropy.

			vector Attributed Energies from	Fit to
			eigenvalue	experiment		eigenvalue
	0 0 0 *	Γ	1	1.1(1)		1.3
	1 2 0	Γ	2	6.5(2)		5.8
	0 4 0	Γ	3	35.1(3)		34.7
	1 2 30 01 3 2 1 7 2 0 0 9 2 0	Y B Z Z	1 1 1 2	41.1(3) 13.8(2) 39.4(5) 49.0(2)		41.0 15.6 37.5 49.4
	3 2	1 2 0	C	1	8.5(1)		7.7
			Reference [16]	Reference [6]	Current work
	J1		-3.8		3.18		-2.6 (2)
	J2		0.2		4.82		0.2 (1)
	J3		13.8		9.08		13.5 (3)
	α		1		0.66		1
	J ′		-		-		-0.32 (3)
	D x	-0.3		-0.89		-0.010 (9)
	D z	0		2.85		0.21 (8)

TABLE III .

 III Table

ACCEPTED MANUSCRIPT / CLEAN COPY

ACKNOWLEDGMENTS

We thank the Institut Laue-Langevin, experiment number 4-01-1608 [20], and Australian Nuclear Science and Technology Organisation (ANSTO), experiment number DB 13377, for the allocation of neutron beam time. Experiments at the ISIS Neutron and Muon Source were supported by a beamtime allocation RB2010329 [25] from the Science and Technology Facilities Council. We also thank Dr. A. Piovano and Dr. A. Ivanov, F. Charpenay and V. Gaignon for assistance with the IN8 experiment, and J. Shalala, the sample environment group, and Dr. A. Stampfl at ANSTO for assistance with the Taipan experiment. ARW and KA thank Dr. L. Testa and Prof. H. Rønnow for assistance with the preliminary identification of suitable single crystals using x-ray Laue diffraction.

Appendix A: Eigenvalues for high-symmetry points

ACCEPTED MANUSCRIPT / CLEAN COPY


Simon Cadiou
email: sicadiou@gmail.com

François Robin

Romain Lecigne

Pascal Guggenbuhl

Raphaël Guillin

Raphael Guillin

Keywords: sacroiliac, bone marrow edema, anatomical variation, MRI, CT

    Bone marrow edema of the sacroiliac joint and anatomical variations: Look back the posterior bony excess. Comment on: "Spine abnormalities associated with bone edema on sacroiliac joints MRI in patients with non-inflammatory chronic back pain" by Marnas et al.

J o u r n a l P r e -p r o o f

We read with much interest the original article about bone marrow edema (BME) on magnetic resonance imaging (MRI) on sacroiliac joints (SIJs) in patients with non-inflammatory chronic back pain [START_REF] Marnas | Spine abnormalities associated with bone edema on sacroiliac joints MRI in patients with non-inflammatory chronic back pain[END_REF]. Marnas et al. found that static disorders and degenerative lesions of spine were associated to BME SIJs in patients with chronic back pain. Usually, in case of BME located in the upper and anterior part of the SIJs, the origin is mainly degenerative. If, BME is located in the inferior and central to posterior part of the SIJs, origin is mostly inflammatory [START_REF] O'shea | Inflammatory and degenerative sacroiliac joint disease in a primary back pain cohort[END_REF][START_REF] Baillet | Bone Marrow Edema in the Sacroiliac Joint-Degenerative Sacroiliac Joint Disease Might Be More Likely Than Spondyloarthritis: Comment on the Article by Turina et al[END_REF]. Marnas et al. showed in an original study that BME of the SIJs in patients with mechanical back pain could be seen in the inferior (56.3%) and posterior parts (52.1%) too.

Anatomical variations of the SIJs are now well-identified causes of inflammatory sacroiliitis diagnosis false positive on MRI [START_REF] El Rafei | Sacroiliac joints: anatomical variations on MR images[END_REF]. However they can be truly symptomatic and associated with pain [START_REF] Ziegeler | Joint anatomy in axial spondyloarthritis: strong associations between sacroiliac joint form variation and symptomatic disease[END_REF]. As recently underlined by Braun et al. [START_REF] Braun | What's new on the sacroiliac joint ? Comment on: "Joint anatomy in axial spondyloarthritis: strong associations between sacroiliac joint form variation and symptomatic disease[END_REF], anatomical variation of SIJs in healthy population is a very frequent condition (about one third of patients). In case of specific anatomical conditions other that degenerative spine, mechanical stresses and anatomical variations may also lead to abnormal mobility of SIJs. The mobility of SIJs is mainly represented by a movement of nutation and count-nutation in the sagittal plane, like a "swing" movement centered by the interosseous sacroiliac ligament, with an estimated range of motion of 2.27 ± 1.35° [START_REF] Lindsey | Evaluation of a minimally invasive procedure for sacroiliac joint fusionan in vitro biomechanical analysis of initial and cycled properties[END_REF] and even more in case of ligament injuries. For example, pregnancy can cause SIJs pain due to asymmetric laxity [START_REF] Damen | Pelvic pain during pregnancy is associated with asymmetric laxity of the sacroiliac joints[END_REF]. Hyperlaxity or excessive stiffness with repeated micro-trauma can also cause degenerative and mechanical damages of SIJs [START_REF] Vleeming | The sacroiliac joint: an overview of its anatomy, function and potential clinical implications[END_REF].

Here, we described a case of a 27-y.o. female with chronic inflammatory low back pain suggestive of axial spondyloarthritis (SpA). A first MRI showed BME of SIJs with normal spine MRI. Given the lack of effectiveness of NSAIDs, the patient was referred to our center to start biological therapy. Clinical examination suggested a sacroiliac origin of the pain with the lack of HLA B27 antigen, a normal protein C reactive level (0.3 mg/dL) and no extra-articular feature of SpA. A new MRI of the SIJs (Figure 1) showed a slight BME on the left sacroiliac joint with dysmorphic and mechanical changes in the inferior part of SIJ on CT. A contralateral BME on the iliac and sacral bone was also noticed, without structural changes. We proposed that these bilateral BME was due to a mechanical stress on the right SIJ resulting from a mobility of SIJs favored by the left anatomical variation in the posterior part. In Figure 2, we can see a dysmorphic sacroiliac joint of the inferior part of the SIJs. It can be identified as a "circumflex accent" in semicoronal view on CT. The nutation, in a swinging movement, may lead to repeated microtrauma of sacrum on ilium, resulting in sclerosis and BME.

Many authors have described anatomical variations of the SIJ [START_REF] El Rafei | Sacroiliac joints: anatomical variations on MR images[END_REF][START_REF] Umay | Frequency of anatomical variation of the sacroiliac joint in asymptomatic young adults and its relationship with sacroiliac joint degeneration[END_REF][START_REF] Kiil | Anatomical variation at the sacroiliac joints in young adults: estimated prevalence by CT and concomitant diagnostics by MRI[END_REF] but it can be difficult to apply this observation in clinical practice, especially on MRI, which visualize poorly the cortical bone and do not allow a semi-axial and semi-coronal reconstruction. As in our case and previous studies [START_REF] Baillet | Bone Marrow Edema in the Sacroiliac Joint-Degenerative Sacroiliac Joint Disease Might Be More Likely Than Spondyloarthritis: Comment on the Article by Turina et al[END_REF], a BME may be associated with SIJs pain. Patients with low back pain and posterior anatomical variations of the SIJs with BME can be grouped as a "posterior bony excess": accessory sacroiliac joint, iliosacral complex and dysmorphic change, which are conditions leading to an excess of contact between ilium and sacrum and then cause BME in the articular part [START_REF] El Rafei | Sacroiliac joints: anatomical variations on MR images[END_REF][START_REF] Kiil | Anatomical variation at the sacroiliac joints in young adults: estimated prevalence by CT and concomitant diagnostics by MRI[END_REF], like a "swinging" movements. These posterior bony excess variations could be then a specific cause of chronic mechanical low back pain, perhaps accessible to infiltration or local treatment as in Bertolotti's syndrome or facet syndrome.

Our reports illustrate the usefulness of identifying anatomical variation of SIJs in case of BME, notably with a double assessment by CT and MRI. We suggest that the BME of SIJs in patients with noninflammatory back pain described by Marnas et al. may not be due only to static disorders and spinal pathologies, but also sometimes to specific anatomical conditions of the SIJs in the posterior part.

Figure 1 .

 1 Figure 1. Sacroiliac joint on CT showing semi-axial view (A and B) and semi-coronal view (C) in CT, and axial view (D and E) and coronal view (F) in MRI. The white arrow show the posterior dysmorphic change of the left ilium in its inferior, almost as an accessory SIJ. The star show the BME of the right SIJ.

Figure 2 .

 2 Figure 2. A and B : the coronal and axial view (respectively) on MRI, showing bilateral BME of the inferior and posterior part of the SIJs (arrows). C and D: semi-coronal (dashed line) and semi-axial (full line) on CT, showing in the same area than image A and B, bilateral sclerosis, bone cysts and gas (circumflex accent), suggesting mechanical lesions.


T P Perrin

C Y M Morio

T Besson

H A Kerhervé

G Y Millet

Dr Jérémy Rossi
email: jeremy.rossi@univ-st-etienne.fr

Comparison of skin and shoe marker placement on metatarsophalangeal joint kinematics and kinetics during running

Keywords: motion capture, footwear, markerset, metatarsophalangeal joint

This study investigated the effects of marker placement (skin-vs shoe-mounted) on metatarsophalangeal joint (MTP) kinematics and kinetics during running. Fifteen trained men ran on a 15-m track at 10 and 13 km/h with three (low, standard and high stiffness) shoe longitudinal bending stiffnesses (LBS). Reflective markers were fixed on the shoe upper, and on the skin using holes cut in the shoe. Three-dimensional marker positions and ground reaction forces were recorded at 200 and 2,000 Hz, respectively. Kinematic and kinetic parameters were analyzed using one-dimensional metrics (statistical parametric mapping). MTP joint was less dorsiflexed at midstance ([57% to 100%] of braking phase and [0% to 48%] of pushing phase), and the MTP joint plantarflexion moment was higher ([22% to 55%] of pushing phase) with the shoe markerset in comparison with the skin markerset. The effect of LBS on MTP angle was found to be significant for a larger percentage of each stride using the shoe markerset compared to the skin markerset. However, the effect of LBS on plantarflexion moment was significant with the shoe markerset only. The effect of running speed on MTP angle was significant for a larger percentage of each stride with the skin markerset. This study demonstrates that the placement of markers influences the measurement of MTP kinematics and kinetics and that these effects are mediated by other variables such as LBS or running speed.

It is concluded that the shoe markerset does not fully reflect the movement of the MTP joint.

Introduction

Motion capture systems and force plates are standard tools used to measure joint kinematics and kinetics with applications in sports medicine and sports performance. Biomechanical studies commonly focus on ankle, knee and hip kinematics to describe locomotion, considering the foot as a unique segment. The rearfoot and the metatarsophalangeal joint (MTP) are usually neglected as they contribute little to the running propulsion. Although some studies have highlighted the elastic energy-saving characteristics of the foot [START_REF] Ker | The spring in the arch of the human foot[END_REF][START_REF] Stearne | The Foot's Arch and the Energetics of Human Locomotion[END_REF], only 0.5 to 2% of the whole-body's positive work is produced by the MTP joint [START_REF] Cigoja | Does increased midsole bending stiffness of sport shoes redistribute lower limb joint work during running[END_REF][START_REF] Hoogkamer | The biomechanics of competitive male runners in three marathon racing shoes: A randomized crossover study[END_REF]. The actual role of foot motion has been reassessed with the advent of carbon plate inserts in racing shoes due to the perturbation induced by the plate on the MTP joint [START_REF] Hoogkamer | The biomechanics of competitive male runners in three marathon racing shoes: A randomized crossover study[END_REF][START_REF] Roy | Shoe midsole longitudinal bending stiffness and running economy, joint energy, and EMG[END_REF][START_REF] Stefanyshyn | Mechanical energy contribution of the metatarsophalangeal joint to running and sprinting[END_REF].

Two types of marker placements can be used for non-invasive estimation of in-shoe joint kinematics -as MTP joint and rearfoot kinematics -using motion capture. The simplest is to fix markers directly on the shoe sole or upper, as in seminal studies analyzing MTP joint motion [START_REF] Stefanyshyn | Mechanical energy contribution of the metatarsophalangeal joint to running and sprinting[END_REF]. A more precise method involves fixing the markers directly on the skin, using holes or 'windows' made in the shoe upper (for a review, see [START_REF] Arnold | Quantifying foot kinematics inside athletic footwear: A review[END_REF]. Previous research suggests that these two placements do not yield the same results regarding rearfoot kinematics. Some studies have shown that shoe markersets significantly underestimated i) the calcaneus range of motion (ROM) across all planes of motion [START_REF] Alcantara | Calcaneus range of motion underestimated by markers on running shoe heel[END_REF][START_REF] Trudeau | The calcaneus adducts more than the shoe's heel during running[END_REF], and ii) the external rotation velocity and peak ankle eversion velocity [START_REF] Sinclair | Differences in Tibiocalcaneal Kinematics Measured with Skin-and Shoe-Mounted Markers[END_REF], concluding that the rearfoot moves inside the shoe. Prior investigations also showed that shoe markersets over-estimate foot pronation [START_REF] Stacoff | The movement of the heel within a running shoe[END_REF], eversion, and tibio-calcaneal rotation across all planes [START_REF] Reinschmidt | Tibiocalcaneal motion during running, measured with external and bone markers[END_REF], suggesting that shoe markersets do not only measure foot kinematics but also the deformation of the shoe [START_REF] Reinschmidt | Tibiocalcaneal motion during running, measured with external and bone markers[END_REF]. Overall, this suggests that foot motion may not be accurately captured by measuring shoe motion. However, to our knowledge, no study has compared the effect of marker placement (foot vs shoe) on the measurement of MTP joint movement during running.

The effect of varying longitudinal bending stiffness (LBS) on MTP joint has been investigated using either skin markersets [START_REF] Cigoja | Does increased midsole bending stiffness of sport shoes redistribute lower limb joint work during running[END_REF](Cigoja et al., , 2020[START_REF] Cigoja | Can changes in midsole bending stiffness of shoes affect the onset of joint work redistribution during a prolonged run[END_REF][START_REF] Day | Does running speed affect the response of joint level mechanics in non-rearfoot strike runners to footwear of varying longitudinal bending stiffness?[END_REF][START_REF] Flores | The stiff plate location into the shoe influences the running biomechanics[END_REF][START_REF] Oh | The bending stiffness of shoes is beneficial to running energetics if it does not disturb the natural MTP joint flexion[END_REF][START_REF] Willwacher | Does specific footwear facilitate energy storage and return at the metatarsophalangeal joint in running[END_REF][START_REF] Willwacher | The gearing function of running shoe longitudinal bending stiffness[END_REF], or shoe markersets [START_REF] Beck | Adding carbon fiber to shoe soles does not improve running economy: A muscle-level explanation[END_REF][START_REF] Farina | Creating footwear for performance running[END_REF][START_REF] Healey | Longitudinal bending stiffness does not affect running economy in Nike Vaporfly Shoes[END_REF][START_REF] Hoogkamer | The biomechanics of competitive male runners in three marathon racing shoes: A randomized crossover study[END_REF][START_REF] Madden | Forefoot bending stiffness, running economy and kinematics during overground running[END_REF][START_REF] Roy | Shoe midsole longitudinal bending stiffness and running economy, joint energy, and EMG[END_REF][START_REF] Stefanyshyn | Mechanical energy contribution of the metatarsophalangeal joint to running and sprinting[END_REF]. Overall, higher LBS reduces dorsiflexion of the MTP joint [START_REF] Day | Does running speed affect the response of joint level mechanics in non-rearfoot strike runners to footwear of varying longitudinal bending stiffness?[END_REF][START_REF] Healey | Longitudinal bending stiffness does not affect running economy in Nike Vaporfly Shoes[END_REF][START_REF] Hoogkamer | The biomechanics of competitive male runners in three marathon racing shoes: A randomized crossover study[END_REF][START_REF] Madden | Forefoot bending stiffness, running economy and kinematics during overground running[END_REF][START_REF] Oh | The bending stiffness of shoes is beneficial to running energetics if it does not disturb the natural MTP joint flexion[END_REF][START_REF] Willwacher | Does specific footwear facilitate energy storage and return at the metatarsophalangeal joint in running[END_REF], reduces MTP joint dorsiflexion angular velocity [START_REF] Day | Does running speed affect the response of joint level mechanics in non-rearfoot strike runners to footwear of varying longitudinal bending stiffness?[END_REF][START_REF] Healey | Longitudinal bending stiffness does not affect running economy in Nike Vaporfly Shoes[END_REF][START_REF] Hoogkamer | The biomechanics of competitive male runners in three marathon racing shoes: A randomized crossover study[END_REF][START_REF] Willwacher | Does specific footwear facilitate energy storage and return at the metatarsophalangeal joint in running[END_REF], decreases MTP joint negative work ( ) (Cigoja et al., 𝑊 - 2019; Day & Hahn, 2021; Farina et al., 2019; Healey & Hoogkamer, 2022; Hoogkamer et al., 2018), and increases positive work ( ) [START_REF] Cigoja | Does increased midsole bending stiffness of sport shoes redistribute lower limb joint work during running[END_REF](Cigoja et al., , 2020;;[START_REF] Willwacher | Does specific footwear facilitate energy storage and return at the metatarsophalangeal joint in running[END_REF].

𝑊 +

Although similar conclusions might seem to be drawn from the two methods regarding changes in LBS, a direct comparison of the two methods still remains to be done. Additionally, MTP joint biomechanics are also affected by running speed: MTP joint is more involved when we go faster with an increase of MTP joint range of motion, maximal plantarflexion moment, and moment at time of peak dorsiflexion [START_REF] Day | Optimal footwear longitudinal bending stiffness to improve running economy is speed dependent[END_REF]). Yet, it is not known whether the effect of running speed on MTP joint kinematics and kinetics is consistent across the two marker placements (Day & Hahn, 2019;[START_REF] Day | Does running speed affect the response of joint level mechanics in non-rearfoot strike runners to footwear of varying longitudinal bending stiffness?[END_REF][START_REF] Flores | The stiff plate location into the shoe influences the running biomechanics[END_REF].

The present study aimed to compare MTP joint kinematics and kinetics during running using skin vs shoe markersets at various LBSs and speeds, to determine whether these two methods could be used interchangeably. We hypothesized that the shoe markerset would not represent MTP joint kinematics and kinetics as accurately as the skin markerset. We also hypothesized that the direction of changes elicited by both LBS and running speed would be similar across markersets, i.e. there would be no effect of marker placement on the interpretation of the influence of LBS and running speed on MTP joint biomechanics.

Methods

Participants

Fifteen recreational male runners were recruited (age 25 ± 5 years, mass 74.3 ± 6.9 kg, height 1.77 ± 0.06 m). They were all rearfoot strikers and free from lower limb injury in the past six months. They provided their written informed consent, in accordance with the Declaration of Helsinki prior to the experiment.

Experimental design

Participants initially ran on a treadmill for 5 min at their preferred running speed for warm up and shoe familiarization purposes [START_REF] Divert | Stiffness Adaptations in Shod Running[END_REF] for each of the three shoe conditions (HIGH: high stiffness; STA: standard stiffness; LOW: low stiffness). Participants then performed 5 trials on a 15-m runway at two different speeds (10 and 13 km/h) for each shoe condition. Running speed was measured using photocells (Whitty Tab, Microgate, Mahopac, NY) placed at the beginning and at the end of the measurement area at ⅔ of the runway. A trial was considered valid if actual speed was within ± 5% of the required speed. Ground reaction forces (GRF) were measured at 2,000 Hz using six force platforms (Kistler, Winterhour, Switzerland) embedded in the runway. Kinematics were recorded at 200 Hz using a motion capture system with 12 optoelectronic cameras (Motion analysis corporation, Raptor 4, Santa Rosa, California) tracking the position of 27 markers placed on the right lower limb of the participants using the markerset described in Besson et al. study (2019) for all markers proximal to the ankle. To measure MTP joint motion, 6 markers were placed on the first and fifth metatarsals, on the hallux and on the posterior, medial and lateral faces of the calcaneus, using holes cut in the shoe (skin markerset, abbreviated FTM for "foot markerset"), in accordance with [START_REF] Bishop | A method to investigate the effect of shoe-hole size on surface marker movement when describing in-shoe joint kinematics using a multi-segment foot model[END_REF] and [START_REF] Shultz | Determining the maximum diameter for holes in the shoe without compromising shoe integrity when using a multi-segment foot model[END_REF]. To measure simultaneously FTM and Shoe-mounted markers (shoe markerset: SHM), SHM were fixed on the shoe upper just below each hole (Figure 1), and were therefore not placed perfectly on the anatomic landmark, but very precisely and accurately placed to respect the alignment of the MTP and avoid shift between foot and shoe MTP joint axis.

INSERT FIGURE 1 NEAR HERE

Shoe conditions

Three visually similar shoe prototypes differing exclusively in LBS were used in the study. The standard shoe had an intermediate level of LBS (STA: 8.4 N/mm), a model with 9 mm grooves cut on the forepart of the midsole provided a lower LBS (LOW: 5.1 N/mm) and a model with a carbon plate embedded between the midsole and the outsole provided a higher LBS (HIGH: 16.2 N/mm). A bespoke three-point flexion test was used to measure shoe LBS without the contribution from the shoe upper, using a frame supporting the shoe on 70% of the length of the sole, and a slow compression dynamometer bending the sole at 1 mm/s until a 5 mm displacement was achieved. For the comparison with other studies, one should be aware that with the addition of the shoe upper, the LBS can be considered slightly greater than the reported values (Fraser et al. 2014). Force was measured at the stabilization point, and then released at 1 mm/s. Measurements were repeated 3 times over 5 consecutive cycles for each shoe and averaged to obtain the LBS of each shoe.

Data Analysis

Visual3D software (C-Motion, Germantown, Maryland) was used to compute MTP joint angle, angular velocity and moment on the sagittal plane. Kinematic marker trajectories and GRF data were low-pass filtered at 30 Hz using a zero-lag 4 th order Butterworth.

Stance phase was detected using a 30 N threshold on the vertical component of the GRF. It has been shown that LBS variation induces different distribution of the braking and pushing phases [START_REF] Willwacher | Does specific footwear facilitate energy storage and return at the metatarsophalangeal joint in running[END_REF]. As we planned to conduct curve analyses using statistical parametric mapping (SPM), we divided the trajectory between braking and pushing phase according to the negative and positive part of the anteroposterior component of the GRF, respectively. This allowed us to avoid temporal shifts between the different phases within the analyzed curves.

Each phase was time normalized from 0-100%. MTP joint was modeled as an oblique axis from the first to the fifth metatarsal [START_REF] Smith | Measurement procedures affect the interpretation of metatarsophalangeal joint function during accelerated sprinting[END_REF]. The joint center was non-fixed, perpendicular to the center of pressure (COP) (Day & Hahn, 2019a). MTP joint movement was defined as the movement of the hallux around the MTP joint axis in the rearfoot coordinate system. MTP joint moment, power, and work were studied only during the period of the stance phase when the COP was ahead of the MTP joint axis on the transverse plane, that we called MTP joint pushing time (Day & Hahn, 2019). Negative ( and positive MTP joint work ( ) were calculated 𝑊 -) 𝑊 +

as the integral of the negative and positive parts of the MTP joint power curve, respectively.

The MTP joint movement variability was assessed by calculating, for each condition, the standard deviation of the 5-trials for minimal and maximal dorsiflexion angle (°) during stance phase.

Statistics

Zero-dimensional analyses were performed on Rstudio (R-Studio Inc., Boston, Massachusetts), for minimal, maximal, and mean values of MTP joint angle, angular velocity, angular moment, positive and negative work, and MTP joint movement variability using the mean value of the five trials for each participant and condition. Normality was checked using Shapiro-Wilk tests.

Three-way repeated-measures ANOVA (markerset speed LBS), followed when × × applicable by Bonferroni post-hoc tests, were performed to compare the markerset effect and interactions. Effect sizes were evaluated using partial eta-squared ( ) and interpreted according 𝜂 2 𝑝 to Cohen's scale (small effect: 0.01 < < 0.06; moderate effect: 0.06 < < 0.14; large effect: [START_REF] Cohen | Statistical Power Analysis for the Behavioral Sciences[END_REF].

𝜂 2 𝑝 𝜂 2 𝑝 > 0.14,
𝜂 2

𝑝

Curve analyses with 1-D statistical parametric mapping (SPM) of the entire kinematics and kinetics time series were performed using Matlab (The Matworks, Natick, Massachusetts) and the open-source spm1d code (www.spm1d.org, [START_REF] Pataky | Zero-vs. One-dimensional, parametric vs. Non-parametric, and confidence interval vs. Hypothesis testing procedures in one-dimensional biomechanical trajectory analysis[END_REF]. Data normality was checked with the d'Agostino test. SPM three-way repeated-measures ANOVA (markerset × speed LBS), followed when applicable by Bonferroni post-hoc tests, were performed to × compare the markerset effect and interactions. The permutation method was applied to SPM ANOVA if normality assumptions were not met. SPM two-way repeated-measures ANOVA (speed LBS) were computed for each markerset, followed when applicable by Bonferroni × post-hoc tests. Alpha level was set at 0.05.

Results

One participant was excluded from the analysis as shoe and skin marker trajectories could not be accurately distinguished during the reconstruction process. Therefore, data from 14 subjects were analyzed.

Markerset influence on MTP joint kinematics and kinetics

Contact, braking, and pushing durations and MTP joint pushing duration are presented in Table 1. There was a significant markerset effect on MTP joint pushing duration (F=86.4; p<0.001; =0.87). MTP joint pushing duration was 4.5 ± 1.0% higher using skin markerset (FTM) than 𝜂 2 𝑝 using shoe markerset (SHM). There was also a significant effect of speed on MTP joint pushing duration (F=89.6; p<0.001; =0.87), with higher speed leading to shorter MTP pushing 𝜂 2 𝑝 duration. There was no significant interaction between LBS, speed and markerset on MTP joint pushing duration.

INSERT TABLE 1 NEAR HERE

There was a main effect of markerset on MTP joint angle during the second half of the braking phase [57 to 100%] and the first half of the pushing phase [0 to 48%] (Figure 2A). MTP joint angle was higher during peak plantar flexion of the stance (around midstance) phase when calculated from SHM compared to FTM (+3.5 ± 4.2°, p<0.001, =0.60). None of the markerset 𝜂 2 𝑝 speed, markerset LBS and markerset speed LBS interactions were significant for

× × × ×
MTP joint angle (respectively F=4.2, p=0.14; F=1.4, p=0.51; F=2.2 p=0.36). There was a main effect of markerset on MTP joint plantarflexion moment (Figure 2B), which was higher with SHM compared to FTM between 22 and 55% of the moment production phase, with a 7 ± 11%

INSERT FIGURE 2 NEAR HERE

Markerset influence on the effect of LBS

The main effect of LBS on the MTP joint angle was significant with SHM (p<0.05 [0 to 83%] of braking phase and [43 to 100%] of pushing phase) and with FTM (respectively [9 to 23%] and [88 to 100%] of braking and pushing phases). A higher dorsiflexion angle was observed with increased LBS (Figure 3A). The LBS effect on peak MTP joint dorsiflexion angle was significant with SHM (F=34.92, p<0.001, =0.73), but not with FTM (F=3.12, p=0.061,

𝜂 2 𝑝 𝜂 2 𝑝 =0.19
). There was a significant main effect of LBS on MTP joint moment only with SHM and at the end of pushing phase ([65 to 99%]). Higher LBSs yielded higher plantarflexion moments (Figure 3B).

INSERT FIGURE 3 NEAR HERE

Markerset influence on the effect of running speed

The effect of speed on MTP joint angle was significant with both markersets in the braking phase (FTM: [39 to 100%], SHM: [0 to 21%] and [40 to 78%]). During the pushing phase, the effect of speed on MTP joint angle was only significant with FTM ([51 to 85%], Figure 4A).

The effect of speed on MTP joint moment was significant over a large part of the stance phase, independently of the markerset ([2 to 92%] for FTM; [2 to 80%] for SHM, Figure 4B).

INSERT FIGURE 4 NEAR HERE

Discussion

This study is the first to investigate MTP joint markers position during running. The primary aim was to compare MTP joint kinematics and kinetics between shoe-and skin-mounted markers during running. We hypothesized that the shoe markerset would not represent the movement of the MTP joint as accurately as the skin-mounted markerset. The main results of this study are that, as compared to skin markerset, shoe markerset overestimated plantarflexion angle at midstance, and peak dorsiflexion moment. The second aim of this study was to investigate the influence of marker placement on the effects of LBS and running speed. We hypothesized that the interpretation of LBS and speed effects would be similar between markersets. Our results do not support this hypothesis as the effects of LBS and speed were dependent on the markerset. For instance, the effect of LBS on MTP joint moment was only significant with the shoe-mounted markerset.

Markerset influence on MTP joint kinematics and kinetics

This study provides evidence that SHM are not representative of MTP joint motion measured using FTM, peak plantarflexion angle and moment were overestimated, compared to FTM.

Such large effects of markerset on MTP joint kinematics and kinetics can be explained by the movement of the foot inside the shoe, and by the additional deformation of the shoe measured with SHM. For instance, we found that MTP joint measurements using SHM were less dorsiflexed than with FTM at midstance (during peak plantarflexion). This difference could originate i) from the slight depression of the metatarsals at midstance pressing in the compliant shoe insole and thereby limiting MTP joint plantarflexion, and ii) in shoe upper overflexion at midstance [START_REF] Miura | Deformation and energy dissipation behaviors of shoe upper material under practical using condition[END_REF]. Combined, these results indicate that peak plantarflexion may be overestimated with SHM. After midstance, the difference between markersets ceases to be significant as the COP continues to move forward of the MTP joint towards the toes.

Overestimation of the MTP joint peak plantarflexion moment by SHM can be explained by differences in MTP joint angle measurements between markersets. The MTP joint lever arm is inversely related to MTP joint dorsiflexion angle. A reduced MTP joint dorsiflexion with SHM implies that a higher MTP joint plantarflexion moment is required to overcome the increased lever arm. Differences in joint moment, and angular velocity (Appendix 1) between markersets leaded to a markerset effect on the measurement of negative and positive work (Appendix 2).

With SHM, the higher dorsiflexion velocity (see Appendix 1. Figure 1) associated with higher plantarflexion moment during the pushing phase led to overestimating compared to FTM.

𝑊 -

At the end of the stance phase, when MTP joint is plantarflexing (approximately 90 to 100% of stance phase), SHM was underestimated because plantarflexion velocity was lower than 𝑊 + with FTM and the plantarflexion moment was similar to FTM.

The shorter MTP joint pushing duration measured using SHM may be due to a backward shift of the foot inside the shoe during running compared to the static standing position, creating a misalignment between MTP joint axes measured with skin or shoe markersets.

Finally, MTP joint angle variability was noticeably lower with SHM than FTM (see Figure 2), especially at midstance (when dorsiflexion angle reaches its minimum). Specifically, the foot appears to move significantly more from trial to trial than the shoe. This could be explained by the foot movement inside the shoe which is not recorded with SHM. In this sense, reducing the movement variability by using SHM could artificially increase the statistical power compared to the same analysis performed using FTM, suggesting that MTP joint kinematics and kinetics should ideally be measured on the skin and not on the shoe, to be more in line with real foot motion.

Markerset influence on the effect of running speed and LBS

Increases in running speed led to higher MTP joint dorsiflexion, higher angular velocities and joint moments in the current study, as it has previously been shown in the literature [START_REF] Day | Optimal footwear longitudinal bending stiffness to improve running economy is speed dependent[END_REF], 2021;[START_REF] Flores | The stiff plate location into the shoe influences the running biomechanics[END_REF][START_REF] Fu | Effect of the Construction of Carbon Fiber Plate Insert to Midsole on Running Performance[END_REF], showing that the involvement of the MTP joint is higher when running speed increase [START_REF] Day | Optimal footwear longitudinal bending stiffness to improve running economy is speed dependent[END_REF]. Our findings indicate that the effect of speed on MTP joint moment is similar using either FTM or SHM, but the magnitude of the effect when studying MTP joint angle is greater with FTM, especially at midstance and during the pushing phase. This difference can be due to the slight depression of the metatarsals in the compliant midsole at midstance which may be greater at higher speeds due to higher ground reaction forces [START_REF] Fu | Effect of the Construction of Carbon Fiber Plate Insert to Midsole on Running Performance[END_REF][START_REF] Nigg | The influence of running velocity and midsole hardness on external impact forces in heel-toe running[END_REF]. Therefore, the increased dorsiflexion at midstance at 13 compared to 10 km/h was only detected with FTM.

Both markersets showed that increased LBS reduced MTP joint dorsiflexion, increased contact time, pushing time, in line with the literature [START_REF] Hoogkamer | The biomechanics of competitive male runners in three marathon racing shoes: A randomized crossover study[END_REF][START_REF] Madden | Forefoot bending stiffness, running economy and kinematics during overground running[END_REF][START_REF] Oh | The bending stiffness of shoes is beneficial to running energetics if it does not disturb the natural MTP joint flexion[END_REF][START_REF] Rodrigo-Carranza | The effects of footwear midsole longitudinal bending stiffness on running economy and ground contact biomechanics: A systematic review and meta-analysis[END_REF][START_REF] Willwacher | Does specific footwear facilitate energy storage and return at the metatarsophalangeal joint in running[END_REF]. However, our findings indicate that the effects of LBS on MTP joint kinematics and kinetics depend on the markerset placement. Indeed, the LBS effect i) was significant for the MTP joint angle for a longer relative duration of the stance phase when using SHM vs FTM, ii) was significant for the MTP joint moment with SHM only. It has previously been shown that a higher LBS led to lower MTP joint peak dorsiflexion angle during propulsion (because the shoe is harder to bend) with both SHM [START_REF] Hoogkamer | The biomechanics of competitive male runners in three marathon racing shoes: A randomized crossover study[END_REF][START_REF] Madden | Forefoot bending stiffness, running economy and kinematics during overground running[END_REF] and FTM [START_REF] Oh | The bending stiffness of shoes is beneficial to running energetics if it does not disturb the natural MTP joint flexion[END_REF][START_REF] Willwacher | Does specific footwear facilitate energy storage and return at the metatarsophalangeal joint in running[END_REF]. Our results showed an effect of LBS for a longer relative duration with SHM compared to FTM, with a significant LBS effect on peak dorsiflexion angle detected only with SHM. Such differences between markersets could be due to: i) the larger angle variability when measured with FTM, which could affect the detection of a statistical effect, ii) the higher shoe deformation compared to foot deformation with increased LBS. However, [START_REF] Oh | The bending stiffness of shoes is beneficial to running energetics if it does not disturb the natural MTP joint flexion[END_REF] and [START_REF] Willwacher | Does specific footwear facilitate energy storage and return at the metatarsophalangeal joint in running[END_REF] detected LBS effect on peak dorsiflexion with FTM, but the range of LBS used on these studies is larger than in the present study, with carbon plate thickness two times that of the highest shoe LBS in this study. Overall, differences between markersets are important enough to lead to different conclusions when comparing LBS. Even if the effect of LBS on contact and pushing duration, and dorsiflexion velocity were not different between markersets, SHM was more sensitive to changes in MTP joint dorsiflexion and joint moment as a function of LBS than FTM. Caution is therefore advised when interpreting LBS effects using shoe markersets as the actual motion of the foot may be occulted.

One limitation of the present study is that SHM could not be placed exactly on the same anatomical landmarks as the skin markerset and were placed just below. One cannot ruled out the fact that if may induced some discrepancies with the perfect anatomical position [START_REF] Mcdonald | Unholey shoes: Experimental considerations when estimating ankle joint complex power during walking and running[END_REF]. However, it was postulated that the movement of the foot inside the shoe was more determinant of the markerset differences. Moreover, this study design was preferred over performing separate trials, as it allowed comparing shoe and skin markersets during the same running strides and avoided repeating running trials that may add variability [START_REF] Mann | Association of previous injury and speed with running style and stride-to-stride fluctuations[END_REF].

Our three-way repeated-measures ANOVA was underpowered to detect Markerset ×LBS × speed interactions due to a small sample size. Even if it is not the center of our analysis, interactions part of the results should be taken with caution. For clarity, the effect of using threeway instead of a two-way repeated measures ANOVA is discussed in Appendix D. Our results are similar between the two analyses (three vs. two-factor RMANOVAs). Therefore, our interpretation and discussion are valid.

Conclusion

Our results indicate that markers fixed on the external surface of the shoe vs directly on the skin are not interchangeable for the measure of MTP joint kinematics and kinetics. Shoe markersets do not reflect the movement of the foot inside the shoe as it cannot detect the slight sliding of the MTP joint into the midsole at midstance, and the effects of LBS and speed on MTP joint kinematics were found to be dependent on the markerset used. Although the skin markerset increased variability, it allowed measuring actual foot motion independent of shoe deformation. Therefore, we recommend using skin markersets when possible, especially when the aim of the study is to analyze LBS. Furthermore, when appraising the contribution of previous literature relating the effects of LBS on running biomechanics, attention must be given to the methodologies that were used.

Conflict of interest statement

No potential conflict of interest was reported by the authors.

FIGURES CAPTIONSFigure 1 .

 1 FIGURES CAPTIONS

Figure 3 .

 3 Figure 3. MTP joint angle (panel A) and joint angular moment (panel B) with shoe markerset (SHM) or skin markerset (FTM) for each LBS. Black horizontal bars denote periods of significant effects of LBS. Dotted black vertical line separates braking and pushing phases. HIGH: high stiffness; STA: standard stiffness; LOW: low stiffness.

Figure 4 .

 4 Figure 4. MTP joint angle (panel A) and joint angular moment (panel B) with shoe markerset (SHM), or skin markerset (FTM) for each speed. Black horizontal bars denote periods of significant effects of speed. Dotted black vertical line separates braking and pushing phases (10: 10 km/h; 13: 13 km/h).

Table 1 . Mean values ± standard deviations for the duration of each part of the contact duration measured for all speed and LBS conditions.

 1 Shoe markerset; FTM: skin markerset; LOW: low LBS; STA: standard LBS; HIGH: high LBS; P: p-value; : Partial eta-squared; n.s: non-significant; * indicates a significant main 𝜂 2 𝑝 effect of the Anova, *: p<0.05, **: p<0.01, *** p<0.001.

				LBS		LBS effect	Speed effect
		Speed (km.h -1 )	LOW	STA	HIGH	P	𝜂 2 𝑝	P	𝜂 2 𝑝
	Contact duration (ms)	10 13	268 ± 17 233 ± 14	268 ± 17 234 ± 13	272 ± 14 236 ± 14	*	0.24	***	0.97

Acknowledgements

The authors would like to thank Decathlon SportsLab and Kalenji for providing the shoes. The authors would like to thank Dr Crane Rogers for the English editing of the manuscript.