Matteo Pancaldi
email: matteo.pancaldi@elettra.eu

Francesco Guzzi
email: francesco.guzzi@elettra.eu

Charles S Bevis

Michele Manfredda

Jonathan Barolak

Iuliia Bykova

Dario De Angelis

Giovanni De Ninno

Mauro Fanciulli

Luka Novinec

Stefano Bonetti

Emanuele Pedersoli

Arun Ravindran

Benedikt Rösner

Christian David

Thierry Ruchon

Alberto Simoncig

Marco Zangrando

Daniel E Adams

Paolo Vavassori

Maurizio Sacchi

George Kourousias

Giulia F Mancini

Flavio Capotondi

High-resolution ptychographic imaging at a seeded free-electron laser source using OAM beams

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Maxwell's equations show that electromagnetic radiation carries energy, polarization and angular momentum. While the first two properties have extensively been used in several experimental approaches, only recently the latter quantity, i.e. the specific spatial distribution of the electromagnetic field phase related to the orbital angular momentum (OAM) of light, has emerged as a valuable experimental tool. The first observation of OAM carried by the electromagnetic field was reported by Allen et al. [1] in 1992. In their work, the authors showed that light described by an exp(𝑖ℓ𝜙) spiraling phase structure intrinsically possesses an OAM equivalent to ℓħ per photon, where the parameter ℓ defines the so-called "topological charge order" of the beam (ℓ ∈ ℤ). After this groundbreaking work, OAM light quickly unlocked new areas of research due to its peculiar applications in both fundamental and applied physics [2][3][4]. Nowadays, OAM beams are harnessed in a great variety of applications covering different fields, such as the study of dichroic and chiroptical materials [5][6][7][8][START_REF] Van Veenendaal | Prediction of Strong Dichroism Induced by X Rays Carrying Orbital Momentum[END_REF][START_REF] Fanciulli | Electromagnetic theory of helicoidal dichroism in reflection from magnetic structures[END_REF][START_REF] Fanciulli | Observation of Magnetic Helicoidal Dichroism with Extreme Ultraviolet Light Vortices[END_REF][START_REF] Rouxel | Hard X-ray helical dichroism of disordered molecular media[END_REF], manipulation and trapping of particles [START_REF] Garcés-Chávez | Observation of the Transfer of the Local Angular Momentum Density of a Multiringed Light Beam to an Optically Trapped Particle[END_REF][START_REF] Adachi | Orbital motion of spherical microparticles trapped in diffraction patterns of circularly polarized light[END_REF], optical communication systems [START_REF] Barreiro | Beating the channel capacity limit for linear photonic superdense coding[END_REF], quantum technologies [START_REF] Molina-Terriza | Twisted photons[END_REF][START_REF] Erhard | Twisted photons: new quantum perspectives in high dimensions[END_REF], and imaging [START_REF] Ritsch-Marte | Orbital angular momentum light in microscopy[END_REF].

Regarding microscopy applications, OAM beams have been proposed as a tool for enhancing the edges detection in phase contrast microscopy [START_REF] Fürhapter | Spiral phase contrast imaging in microscopy[END_REF][START_REF] Sakdinawat | Soft-x-ray microscopy using spiral zone plates[END_REF][START_REF] Sharma | Fourier spectrum analysis of spiral zone plates[END_REF], and for increasing the spatial resolution [START_REF] Tamburini | Overcoming the Rayleigh Criterion Limit with Optical Vortices[END_REF][START_REF] Li | Beating the Rayleigh limit: Orbital-angular-momentum-based super-resolution diffraction tomography[END_REF][START_REF] Wei | Sub-100nm resolution PSIM by utilizing modified optical vortices with fractional topological charges for precise phase shifting[END_REF][START_REF] Tham | Beating Rayleigh's Curse by Imaging Using Phase Information[END_REF]. It has been proven that, at optical wavelengths, OAM beams enhance the image quality in phase-sensitive techniques, providing a uniform contrast at the interface between different elements, dislocations or morphological height variations [START_REF] Ritsch-Marte | Orbital angular momentum light in microscopy[END_REF][START_REF] Li | Complex wavefront sensing based on coherent diffraction imaging using vortex modulation[END_REF][START_REF] Aspden | Heralded phase-contrast imaging using an orbital angular momentum phase-filter[END_REF]. Basically, in such experiments a spiral phase plate in a 4-f optical system acts as a bidimensional Hilbert filter in Fourier space, generating a contrast enhancement at the edges of the target structure. This is of particular importance in low-transmission contrast samples like biological cells, where topological height information can be retrieved through the analysis of the phase contrast [START_REF] Ritsch-Marte | Orbital angular momentum light in microscopy[END_REF][START_REF] Jesacher | Shadow Effects in Spiral Phase Contrast Microscopy[END_REF]. At shorter wavelengths, diffraction imaging techniques based on beam coherence are routinely employed to obtain high-resolution diffraction-limited images of nanostructured samples [START_REF] Miao | Coherent X-Ray Diffraction Imaging[END_REF][START_REF] Prosekov | Methods of Coherent X-Ray Diffraction Imaging[END_REF][START_REF] Karl | Full-field imaging of thermal and acoustic dynamics in an individual nanostructure using tabletop high harmonic beams[END_REF]. Among them, ptychography represents a complete in-situ characterization tool, since it allows retrieving both the transmission and phase contrast image of the sample, as well as wavefront-sensing information of the illumination function [START_REF] Thibault | Probe retrieval in ptychographic coherent diffractive imaging[END_REF][START_REF] Maiden | An improved ptychographical phase retrieval algorithm for diffractive imaging[END_REF][START_REF] Pfeiffer | X-ray ptychography[END_REF][START_REF] Guizar-Sicairos | Ptychography: A solution to the phase problem[END_REF]. In combination with OAM beams, ptychography has been exploited in the characterization of the vortex beam structure, and in the determination of the topological charge order [START_REF] Vila-Comamala | Characterization of x-ray phase vortices by ptychographic coherent diffractive imaging[END_REF][START_REF] Dzhigaev | Ptychographical imaging of the phase vortices in the x-ray beam formed by nanofocusing lenses[END_REF][START_REF] Loetgering | Generation and characterization of focused helical x-ray beams[END_REF][START_REF] Baluktsian | Mode purity and structural analysis of x-ray vortices generated by spiral zone plates[END_REF]. More recently, beneficial effects given by the structured illumination of OAM beams for the ptychographic imaging of extended samples have been demonstrated by Eschen et al. [START_REF] Eschen | Material-specific high-resolution table-top extreme ultraviolet microscopy[END_REF] and Wang et al. [START_REF] Wang | High-fidelity ptychographic imaging of highly periodic structures enabled by vortex high harmonic beams[END_REF]. In such cases, the enriched spatial frequency content of the OAM illumination function optimizes the distribution of scattered photons at the detector plane, relaxing the dynamic range requirements and improving the signal-to-noise ratio (SNR). Moreover, due to the reduced scattering symmetry and the increased diversity of collected diffraction patterns, the convergence and the robustness of ptychographic reconstruction algorithms are improved with respect to the case of conventional Gaussian illumination.

In this work, we report on high quality ptychography image reconstructions achieved by exploiting zone-plate-generated OAM beams [START_REF] Rösner | High resolution beam profiling of X-ray free electron laser radiation by polymer imprint development[END_REF][START_REF] Rebernik Ribič | Extreme-Ultraviolet Vortices from a Free-Electron Laser[END_REF] at an extreme ultraviolet (EUV) seeded free-electron laser (FEL) source [START_REF] Allaria | Highly coherent and stable pulses from the FERMI seeded free-electron laser in the extreme ultraviolet[END_REF]. With respect to previous studies [START_REF] Sala | Pulse-to-pulse wavefront sensing at free-electron lasers using ptychography[END_REF][START_REF] Kharitonov | Flexible ptychography platform to expand the potential of imaging at free electron lasers[END_REF][START_REF] Kharitonov | Single-shot ptychography at a soft X-ray free-electron laser[END_REF], we experimentally observe that the improved pointing stability introduced by the zone plate optics, combined with the high intensity and wavelength stability of a seeded FEL source, allow us to obtain high quality images of a nanostructured Siemens star test object with single-shot per position illumination. Moreover, an increase in the achievable spatial resolution of the reconstructed ptychography images is observed at high topological charge orders, due to the structured illumination of OAM beams. This result extends the ability of coherent diffraction imaging techniques to perform time-resolved high-resolution imaging of large area samples. Finally, we observe that the retrieved illumination function can be used for determining the main residual aberrations induced by the focusing optical system. Besides the importance for imaging experiments, our analysis is also relevant for dichroic experiments involving OAM beams [START_REF] Fanciulli | Electromagnetic theory of helicoidal dichroism in reflection from magnetic structures[END_REF][START_REF] Rouxel | Hard X-ray helical dichroism of disordered molecular media[END_REF], since modifications in the expected intensity pattern for opposite ℓ values can introduce artifacts in the evaluation of the dichroic signal.

Ptychography and experimental setup

The ptychographic technique

Ptychography merges the main advantages of lensless imaging [START_REF] Chapman | Coherent lensless X-ray imaging[END_REF], since it combines the use of phase retrieval algorithms with a scanning probe approach typical of x-ray microscopy [START_REF] Gianoncelli | Soft X-ray Microscopy Techniques for Medical and Biological Imaging at TwinMic-Elettra[END_REF], in order to reconstruct diffraction-limited morphological information of extended samples. An essential requirement of ptychography is the spatial overlap (typically > 70%) among the areas of the sample generating single diffraction patterns. This redundancy in collected data information allows for lifting the isolated object constraint of traditional coherent diffraction imaging [START_REF] Pfeiffer | X-ray ptychography[END_REF]. Being a computational imaging technique, ptychography requires an iterative reconstruction procedure, which has to be carefully checked and tuned in order to avoid artifacts. The net result is a high-resolution complex-valued image (transmission amplitude and phase) of the sample and a refinement of the probing beam illumination function [START_REF] Guizar-Sicairos | Ptychography: A solution to the phase problem[END_REF]. As a consequence, ptychography finds applications not only as an imaging technique, but also as a wavefront-sensing metrological diagnostic tool [START_REF] Vila-Comamala | Characterization of high-resolution diffractive X-ray optics by ptychographic coherent diffractive imaging[END_REF][START_REF] Schropp | Full spatial characterization of a nanofocused x-ray free-electron laser beam by ptychographic imaging[END_REF][START_REF] Moxham | Aberration characterization of x-ray optics using multi-modal ptychography and a partially coherent source[END_REF]. Figure 1(a) shows a sketch of a typical experimental setup: the scattered beam is measured in transmission geometry, recording the diffraction patterns at the detector plane jointly with the sample coordinates. After scanning the desired set of sample positions, an iterative procedure solves the inverse problem of finding the object and the average illumination that better approximate the set of measured scattering patterns. The accuracy and convergence of a ptychography reconstruction rely on the precise knowledge of the input parameters which describe the wave propagation in the reconstruction algorithm (i.e., pulse's central wavelength and bandwidth, beam size, pixel size, propagation distances and sample positions).

Our ptychography algorithm, described in Ref. [START_REF] Guzzi | A Parameter Refinement Method for Ptychography Based on Deep Learning Concepts[END_REF], automatically factorizes the underlying object and the common illumination shed on it [START_REF] Thibault | Probe retrieval in ptychographic coherent diffractive imaging[END_REF][START_REF] Maiden | An improved ptychographical phase retrieval algorithm for diffractive imaging[END_REF], which can also include optical aberrations, like astigmatism or coma. Moreover, the trade-off between photon flux and spatial coherence of the illumination can be relaxed, considering partially coherent beam conditions, e.g. when multiple propagating modes incoherently add at the detector plane, which would result in a blurred diffraction pattern [START_REF] Spence | Coherence and sampling requirements for diffractive imaging[END_REF]. In this case, to solve the inverse problem, probe decomposition techniques (such as the one presented in Ref. [START_REF] Thibault | Reconstructing state mixtures from diffraction measurements[END_REF]) are employed to computationally disentangle the weighted contribution of each individual mode to the recorded diffraction pattern dataset. For single-shot per position experiments, the decomposition of the illumination function in multiple modes intrinsically allows for estimating a lower bound for the degree of transversal coherence of the beam [START_REF] Vartanyants | Origins of decoherence in coherent X-ray diffraction experiments[END_REF][START_REF] Enders | Ptychography with broadbandwidth radiation[END_REF].

In order to improve the object self-consistency during the iterative process [START_REF] Cop | A multiresolution approach to the 3D reconstruction of a 50S ribosome from an EMtilt series solving the alignment problem without gold particles (electron microscopy)[END_REF], we impose a refinement of the sample positions associated with each diffraction pattern in the scan [START_REF] Guzzi | A Parameter Refinement Method for Ptychography Based on Deep Learning Concepts[END_REF][START_REF] Guzzi | Automatic Differentiation for Inverse Problems in X-ray Imaging and Microscopy[END_REF]. This additional optimization procedure allows for the evaluation of spurious vibrations on the beam pointing in terms of position offsets, providing a quantitative measure of the setup stability, and it is of paramount importance for the final image quality. Similar iterative sampleposition refinement approaches have also been employed in Refs. [START_REF] Sala | Pulse-to-pulse wavefront sensing at free-electron lasers using ptychography[END_REF][START_REF] Kharitonov | Flexible ptychography platform to expand the potential of imaging at free electron lasers[END_REF] for enforcing a common illumination, even in the extreme case of strong pointing instability due to a vibrating beam.

Description of the experimental setup

The imaging capabilities of ptychography with OAM beams have been experimentally investigated at the DiProI beamline [START_REF] Capotondi | Invited Article: Coherent imaging using seeded free-electron laser pulses with variable polarization: First results and research opportunities[END_REF] of the FERMI FEL source [START_REF] Allaria | Highly coherent and stable pulses from the FERMI seeded free-electron laser in the extreme ultraviolet[END_REF]. The end-station is equipped with an active focusing system based on Kirkpatrick-Baez (K-B) mirrors [START_REF] Raimondi | Kirkpatrick-Baez active optics system at FERMI: system performance analysis[END_REF], which were set in flat-mirrors configuration for compensating the natural divergence of the FEL source and producing a nearly collimated beam. The focusing of the EUV radiation and the generation of OAM beams were performed by means of nanofabricated spiral zone plates (SZPs) [START_REF] Rösner | High resolution beam profiling of X-ray free electron laser radiation by polymer imprint development[END_REF][START_REF] Rebernik Ribič | Extreme-Ultraviolet Vortices from a Free-Electron Laser[END_REF] directly placed inside the main vacuum chamber of the end-station. The FEL light was linearly polarized in the horizontal plane with a wavelength 𝜆 of 18.9 nm (65.6 eV). Taking into consideration the beamline transmission and the 500 nm-thick Zr solid-state filter used for attenuating the total source intensity, we estimate a pre-focusing fluence of about 0.1 µJ/cm 2 irradiating the SZPs, well below their damage threshold and guaranteeing long-term stability.

A sketch of the experimental setup is shown in Fig. 1(a). Motorized piezo stages were used for precisely switching between seven SZPs for the generation of ℓ = 0, ±1, ±2, ±3 OAM beams [62]. All SZPs were fabricated with nominally identical physical parameters [START_REF]SZP parameters: diameter D of 1.92 mm, outermost zone width dr of 1.62 μm[END_REF], in order to provide a focal length 𝑓 of 164.7 mm for 𝜆 = 18.9 nm radiation, with an estimated first order focusing efficiency of about 12% [START_REF] Rebernik Ribič | Extreme-Ultraviolet Vortices from a Free-Electron Laser[END_REF]. According to the average input fluence and to the theoretical full width at half maximum spot size (1.67 μm), we calculate a fluence of about 15 mJ/cm 2 at the focal plane for the Gaussian illumination case (ℓ = 0). An order sorting aperture and a central beam stop were used for spatially filtering the higher diffraction orders of the SZPs and for blocking the unfocused direct beam, respectively.

Ptychography scans were performed on a 100 nm-thick hydrogen silsesquioxane (HSQ) Siemens star structure [START_REF] Horstmeyer | Standardizing the resolution claims for coherent microscopy[END_REF] patterned on top of a 200 nm-thick silicon membrane. The outer diameter of the Siemens star is approx. 500 μm, while the smallest features close to the center have a size of 100 nm, as shown in Fig. 1(b). The sample was mounted on a five-axis piezo goniometer, which allows for controlled movements with 100 nm of resolution. For each ℓ value, the center of the Siemens star was scanned in a 7×7 grid, where the grid points were chosen according to a Fermat spiral pattern with a step size of 1 μm [START_REF] Huang | Optimization of overlap uniformness for ptychography[END_REF]. At each point of the scan, a single FEL shot was used for illuminating the sample, and the diffracted intensity was recorded by an in-vacuum Princeton MTE2048 CCD camera (2048 × 2048 pixels, with a 13.5 µm pixel size), placed at a distance 𝑑 of 140 mm from the sample plane and operated in 4×4 binning mode.

Results

Ptychography reconstructions at a seeded FEL

The ptychographic reconstructions were carried out using the M-rPIE algorithm [START_REF] Guzzi | A modular software framework for the design and implementation of ptychography algorithms[END_REF] implemented in the SciComPty software package [START_REF] Guzzi | A Parameter Refinement Method for Ptychography Based on Deep Learning Concepts[END_REF], by including the simultaneous propagation of three independent illumination modes [START_REF] Thibault | Reconstructing state mixtures from diffraction measurements[END_REF]. For each ℓ value, three ptychographic scans were recorded, and representative reconstructions are shown in Fig. 2. Each panel is composed by amplitude and phase of the main illumination mode, and by the corresponding object amplitude. The reconstructed pixel size can be precisely calibrated in order to retrieve the beam and object sizes at the sample plane, as reported in Supplement 1, Sec. S1. Since beam coherence is a fundamental property of FEL sources, Fig. 3(a) shows the modal weights associated to the three illumination modes for each reconstruction. The modal weight 𝛽 𝑖 is obtained by normalizing the integrated intensity of the i-th mode with respect to the total intensity. Most of the weight (82% in average) is associated with the first decomposition mode, while each one of the higher order modes contributes to less than 12% to the total illumination function, hence allowing us to consider only the reconstruction of the main illumination modes in Fig. 2. As discussed in Ref. [START_REF] Enders | Ptychography with broadbandwidth radiation[END_REF], the decomposition of the illumination function in orthogonal modes provides insight into the factors affecting the speckle visibility in the diffraction patterns, which can be related to, e.g., partial transverse and longitudinal coherence of the radiation, vibration of the beam around an average position, detector point spread function. It is worth pointing out that, in single-shot per position experiments, the speckle visibility degradation is a representative measurement of the beam partial coherence since, in contrast to multiple-shot exposure, the blurring of a single-shot diffraction pattern cannot be affected by beam pointing fluctuations. As a consequence, it is possible to evaluate the degree of transversal coherence 𝜁 for the beam through the modal weights, according to

𝜁 = ∑ 𝛽 𝑖 2 𝑖 (∑ 𝛽 𝑖 𝑖 ) 2 ⁄ [67]
. Since the modal decomposition can favor the formation of unphysical modes that spread part of the diffraction-uncorrelated intensity outside the object area [START_REF] Enders | Ptychography with broadbandwidth radiation[END_REF] (see Supplement 1, Sec. S2), we can only provide a lower boundary for 𝜁. Considering the average < 𝛽 𝑖 > values in Fig. 3(a), we obtain that 𝜁 = (68 ± 2)%. Such lower boundary for the degree of transversal coherence is in good agreement with previous estimates obtained on a seeded FEL source by means of the Young's double slit interference experiment [START_REF] Allaria | Highly coherent and stable pulses from the FERMI seeded free-electron laser in the extreme ultraviolet[END_REF][START_REF] Capotondi | Invited Article: Coherent imaging using seeded free-electron laser pulses with variable polarization: First results and research opportunities[END_REF]. In a recent work, Kharitonov et al. [START_REF] Kharitonov | Flexible ptychography platform to expand the potential of imaging at free electron lasers[END_REF] performed a similar analysis using a self-amplified spontaneous emission (SASE) FEL source [START_REF] Pellegrini | The history of X-ray free-electron lasers[END_REF]. They reported, in a comparable spectral range, a transversal coherence of about 54%, and the need for nine independent illumination modes to obtain reliable and stable reconstructions of the object. Our mode decomposition analysis shows favorable conditions in terms of transversal coherence quality for seeded FEL radiation in comparison to a SASE FEL source. Moreover, the seeding scheme also provides a higher degree of longitudinal coherence with respect to SASE, with the generated radiation behaving statistically as a laser [START_REF] Yu | Seeded X-ray free-electron laser generating radiation with laser statistical properties[END_REF]. As a consequence, the narrow spectral bandwidth of seeded FEL radiation allows for a well-defined focal length for the chromatic SZP optical element, so increasing the number of photons available inside the coherent scattering volume and improving the contrast visibility of the speckle pattern. These properties make the seeded FEL an ideal source for performing time-resolved high-resolution ptychography experiments.

Besides beam coherence, ptychography reconstructions are very sensitive to positioning errors. Beam positioning errors can be a severe challenge for stable and reliable reconstruction algorithms, especially for ptychographic imaging experiments at FEL sources due to intrinsically larger pointing instabilities with respect to synchrotron radiation facilities [START_REF] Sala | Pulse-to-pulse wavefront sensing at free-electron lasers using ptychography[END_REF][START_REF] Kharitonov | Flexible ptychography platform to expand the potential of imaging at free electron lasers[END_REF]. Moreover, in single-shot per position experiments the preservation of a precise positioning during the sample scan is crucial, since fluctuations are not mitigated by averaging over many shots at the expense of a worse image resolution. An estimate of our system stability can be obtained by analyzing the output of the position-refinement step of our algorithm [START_REF] Guzzi | A Parameter Refinement Method for Ptychography Based on Deep Learning Concepts[END_REF]. In this procedure, the position of each diffraction pattern of the scan is optimized from its starting point to iteratively extract the offset along the scanning axis for improving the convergence of the ptychographic reconstruction. Figure 3(b) shows the distribution of the final refinement offsets in the X-Y sample plane for all sample positions considered during the experiment (approx. 1000 points). The full widths at half maximum of both the X-and Y-projection are about 300 nm. This analysis clearly shows that, due to the large numerical aperture, the use of diffractive optical elements to focus a collimated beam strongly mitigate the detrimental effects of pointing instabilities at the imaging plane by more than one order of magnitude with respect to previously performed experiments at FEL sources, where conventional reflective K-B optics were used [START_REF] Sala | Pulse-to-pulse wavefront sensing at free-electron lasers using ptychography[END_REF][START_REF] Kharitonov | Flexible ptychography platform to expand the potential of imaging at free electron lasers[END_REF].

Analysis of the reconstructed illumination function

In the previous section we have shown that most of the intensity of the illumination function is concentrated in the first mode of the modal decomposition. For each ℓ value, Fig. 2 reports the corresponding principal mode illumination, showing that the ptychographic reconstructions successfully identify the correct phase winding direction and topological charge order of the selected SZP. Indeed, by looking at the central part of the figure panels reporting the illumination profiles, a 2𝜋 phase winds |ℓ| times around the propagation axis, and the direction is opposite for opposite signs of the topological charge order. However, a careful inspection of all the retrieved spots shows two principal distortions with respect to the ideal "donut-like" shape of a Laguerre-Gaussian mode [1,[START_REF] Rösner | High resolution beam profiling of X-ray free electron laser radiation by polymer imprint development[END_REF]: (i) an azimuthal modulation of the beam intensity and (ii) a separation of the expected vortex singularity in multiple poles for topological charges |ℓ| > 1. With the help of simulations (see Supplement 1, Sec. S3), we were able to identify the 45-degree astigmatism induced by the beamline K-B optics as the major source of the abovementioned probe beam deformations. The effect of this optical aberration is to split the central phase singularity of an ideal Laguerre-Gaussian beam into |ℓ| single-order singularities aligned along a line [START_REF] Baluktsian | Mode purity and structural analysis of x-ray vortices generated by spiral zone plates[END_REF][START_REF] Ya | Astigmatic telescopic transformation of a high-order optical vortex[END_REF][START_REF] Wang | Generalised Hermite-Gaussian beams and mode transformations[END_REF]. Despite the presence of astigmatism, the total OAM topological charge order is conserved, as shown by the proper winding of the overall reconstructed phase and as discussed in Ref. [START_REF] Abramochkin | Generalized Gaussian beams[END_REF].

The presence of 45-degree astigmatism was confirmed by numerically propagating the reconstructed ℓ = 3 illumination function along the Z-axis by means of the free-space propagator for a monochromatic wave [START_REF] Paganin | Coherent X-Ray Optics[END_REF]. Figure 4(a),(b) shows sections of this free-space propagation corresponding to the XZ-plane and to the YZ-plane, respectively. In these projected views, the beam waist is at the same position, approximately at 𝑍 = -200 μm, showing the perfect matching between tangential and sagittal focusing planes. However, along planes rotated by 45° with respect to the X-Y coordinate system, a relative difference of 400 μm in the beam waist position is observed (see Fig. 4(c),(d)), confirming the presence of a residual 45-degree astigmatism not fully compensated by the K-B optics.

It is worth noting that, since distortions induced by astigmatism have a different impact on the illumination function for opposite ℓ values, experiments based on OAM chiral dichroism [START_REF] Fanciulli | Electromagnetic theory of helicoidal dichroism in reflection from magnetic structures[END_REF][START_REF] Rouxel | Hard X-ray helical dichroism of disordered molecular media[END_REF] must handle with great care the quantitative comparison between measurements obtained for opposite ℓ values. Indeed, spurious chiral signals can be generated (e.g.) by spatial intensity variations due to optical aberrations. Ptychography, with its wavefront-sensing capabilities, could represent the ideal metrological technique for in-situ beam characterization, in order to minimize optical aberrations before performing the target chiral experiments. Fig. 4. Free-space propagation of the ℓ = 3 main illumination mode. (a) Section of the freespace propagation corresponding to the XZ-plane (top view). The white line highlights the position of the beam waist. (b) Section of the free-space propagation corresponding to the YZplane (side view). The white line highlights the position of the beam waist, which is the same as for (a). (c) Section of the free-space propagation corresponding to the X'Z-plane. (d) Section of the free-space propagation corresponding to the Y'Z-plane. The sections in panels (c) and (d) correspond to planes rotated by 45° with respect to the X-Y coordinate system. The fact that in these views the beam waist is in different positions, as marked by the white dashed lines, confirms the presence of 45-degree astigmatism.

Analysis of the reconstructed object image

The intrinsic structuring of an OAM beam and the additional contribution to the illumination diversity provided by astigmatism are instrumental in enhancing the quality of the reconstructed object images. Figure 2 already shows that the fine details in the central part of the Siemens star can readily be identified, allowing to estimate a resolution comparable with the pixel size (96 nm), an order of magnitude better with respect to a SASE FEL source [START_REF] Kharitonov | Flexible ptychography platform to expand the potential of imaging at free electron lasers[END_REF][START_REF] Kharitonov | Single-shot ptychography at a soft X-ray free-electron laser[END_REF] and comparable with more stable sources, like table-top lasers for high-harmonic generation [START_REF] Eschen | Material-specific high-resolution table-top extreme ultraviolet microscopy[END_REF][START_REF] Wang | High-fidelity ptychographic imaging of highly periodic structures enabled by vortex high harmonic beams[END_REF][START_REF] Tadesse | Wavelength-scale ptychographic coherent diffractive imaging using a high-order harmonic source[END_REF] and synchrotron radiation facilities [START_REF] Guizar-Sicairos | Role of the illumination spatial-frequency spectrum for ptychography[END_REF][START_REF] Maiden | Soft X-ray spectromicroscopy using ptychography with randomly phased illumination[END_REF]. To obtain a quantitative evaluation of image resolution, Fig. 5(a) shows the experimental average Fourier ring correlation (FRC) curves calculated averaging the FRC results obtained by permuting pairs of different reconstructed images collected at a given topological charge order |ℓ|. As described in Ref. [START_REF] Van Heel | Fourier shell correlation threshold criteria[END_REF], the image resolution can be extracted by defining a threshold for the FRC values, which allows getting the maximum spatial frequency to achieve a predefined SNR. By considering a SNR of 0.41 (corresponding to the common 0.5-bit threshold), we obtained a resolution of (118 ± 7) nm, (128 ± 8) nm, (96 ± 5) nm, and (85 ± 5) nm for |ℓ| = 0, 1, 2, 3, respectively. Under similar experimental conditions, a clear improvement of about 30% in the image resolution is observed for the SZP with the highest topological charge order (|ℓ| = 3) with respect to conventional Gaussian beam illumination (ℓ = 0). Similar beneficial effects in imaging resolution going from Gaussian illumination to OAM beams have recently been reported by Tan et al. [START_REF] Tan | Enhancing critical resolution of a ghost imaging system by using a vortex beam[END_REF] at optical wavelengths and by Eschen et al. [START_REF] Eschen | Material-specific high-resolution table-top extreme ultraviolet microscopy[END_REF] in the EUV regime. We ascribe this resolution improvement to the beam structuring and to its consequent enhancement of the probe spatial-frequency spectrum. On one hand, this effect increases the SNR due to a more efficient distribution of scattered photons at the detector plane; on the other, it increases the diversity of the recorded diffraction patterns during the scan, speeding up the reconstruction convergence.

Besides the resolution improvement due to the higher frequency content of structured beams [START_REF] Guizar-Sicairos | Role of the illumination spatial-frequency spectrum for ptychography[END_REF][START_REF] Tan | Enhancing critical resolution of a ghost imaging system by using a vortex beam[END_REF], the effect of an OAM beam illumination is also visible at low spatial frequencies (below 0.1 in our case, see Fig. 5(a)), where the FRC values grow along with the value of |ℓ|. According to Odstrčil et al. [START_REF] Odstrčil | Towards optimized illumination for high-resolution ptychography[END_REF], the observed enhancement is due to the improvement in the reconstruction accuracy of the object phase, which is determined by the minor sensitivity to unwrapped phase errors with respect to constant phase probes when the size of the investigated object is larger than the illumination spot size (i.e. at low reconstructed frequencies). In other words, the OAM phase gradient at length scales comparable to the beam spot size allows to unequivocally determine the object phase among successive scanning positions during the reconstruction process, so improving the image quality at low spatial frequencies. This enhancement of image quality at spatial frequencies comparable with the beam spot size agrees with the results reported by Li and Li [START_REF] Li | Beating the Rayleigh limit: Orbital-angular-momentum-based super-resolution diffraction tomography[END_REF]. Indeed, their work suggests that OAM light can exceed the Rayleigh limit and achieve super-resolution when the beam is tightly focused to a spot size comparable to the radiation's wavelength. Despite the highest resolution being associated with the highest ℓ value, the |ℓ| = 1 resolution falls out of the general trend, suggesting the existence of a trade-off between the increase of the high frequency spectral content induced by non-ideality of the illumination function and the minimum resolvable separation between two adjacent structures [START_REF] Tan | Enhancing critical resolution of a ghost imaging system by using a vortex beam[END_REF]. We ascribe this evidence to the fact that, even if the resolution should increase with |ℓ| due to the higher frequency content of the illumination function, the variations in the illumination function profile should match the peculiar object features to be efficient in terms of SNR. In other words, when the characteristic spatial frequency of the probing beam does not match with the minimum sample feature, the reduced speckles contrast of the far field diffraction pattern produces a deterioration of image quality. To corroborate this point, Fig. 5(b) shows the amplitude of the Fourier transform of the product between the reconstructed illumination functions in Fig. 2 (for ℓ = 0, 1, 2, 3) and the average reconstructed object. From this quantity, which in logarithmic scale is proportional to Fraunhofer diffraction [START_REF] Paganin | Coherent X-Ray Optics[END_REF], it is indeed possible to compare the contrast of the far field speckle pattern for different ℓ values. It is apparent that, at high spatial frequencies (e.g. above 0.75, see the black dashed circle in Fig. 5(b)), the details of the fine features in the speckle are comparable for ℓ = 0 and 1, while the structuring is more defined for ℓ = 2 and 3, so confirming the trend emerging from the FRC curves.

Conclusions

In this work, we report on a series of single-shot per position ptychography experiments at a seeded FEL source, which were designed to evaluate the effect of the OAM topological charge order ℓ on the reconstructed image quality. Our analysis demonstrates an improvement of more than one order of magnitude in spatial resolution and beam stability with respect to previously reported experiments at FEL sources [START_REF] Sala | Pulse-to-pulse wavefront sensing at free-electron lasers using ptychography[END_REF][START_REF] Kharitonov | Flexible ptychography platform to expand the potential of imaging at free electron lasers[END_REF]. We observe a clear trend relating the achieved image resolution with ℓ: the higher the |ℓ| value, the better the image resolution and quality (also considering the improvement of image contrast at low spatial frequencies). The higher topological charge orders (|ℓ| ≥ 2) show an improvement of about 30% in image resolution with respect to conventional Gaussian beam illumination, paving the way for achieving superresolution with OAM beams in the EUV regime. Moreover, the beam characterization, performed through the ptychographic reconstruction of the illumination functions, confirms that OAM beams are very sensitive to optical aberrations, which distort the beam transversal profile with respect to ideal Laguerre-Gaussian intensity patterns. For this reason, the experimentally achieved resolution shows a minimum for |ℓ| = 1, and we were able to ascribe this apparent anomaly to the reduced speckles contrast at high spatial frequencies, which is then improved by increasing ℓ. Our findings are relevant for the design of fully single-shot FELbased ptychographic imaging experiments exploiting OAM light, since the use of high ℓ values can maximize the beneficial effects of the non-ideality of the illumination function determined by residual beam aberrations.

Fig. 1 .

 1 Fig. 1. (a) Schematics of the experimental setup for ptychography. The OAM beam is generated and focused on the sample plane (at a distance 𝑓) using a spiral zone plate (SZP). The order sorting aperture (OSA) allows only light from the first diffraction order to reach the sample, while the higher orders are blocked. Upstream to the SZP optics, a central beam stop (CBS) blocks the direct beam that does not carry OAM. The sample is then scanned in the XY-plane in order to record a set of diffraction patterns with a CCD camera (placed at a distance 𝑑), to be used for the ptychography reconstruction of both the object and the illumination functions. (b) Scanning electron microscopy image of the Siemens star used as a test transmitting sample, and magnification of the central part, where the smallest features have a size of 100 nm.

Fig. 2 .

 2 Fig. 2. (a)-(g) Ptychographic reconstructions of the Siemens star for ℓ = 0, ±1, ±2, ±3. For each panel, the main illumination amplitude, the main illumination phase and the object amplitude are shown. Both the white markers in (a) correspond to a length of 5 μm, and they differ by a factor of two.

Fig. 3 .

 3 Fig. 3. (a) Modal weights 𝛽 𝑖 associated to the three illumination modes for each reconstruction. The dashed lines highlight the average 𝛽 𝑖 values, while the shaded areas represent the standard deviations. (b) Offsets in the scan positions as determined by the position-refinement algorithm along the X-and Y-axis for all the reconstructions. The full widths at half maximum of both the projections in the lateral plots are 300 nm.

Fig. 5 .

 5 Fig. 5. (a) Average Fourier ring correlation (FRC) curves calculated on the experimental ptychography reconstructions. The upper horizontal axis shows the corresponding half-period resolution calibrated with the reconstructed pixel size. The shaded areas represent the standard deviation calculated over all the datasets illuminated with the same |ℓ|. The black dashed line indicates the 0.5-bit threshold, which is a common criterion for defining the image resolution and corresponds to a SNR of 0.41. (b) Amplitude of the Fourier transform of the product between the main mode (𝑝 1 ) of the reconstructed illumination functions in Fig. 2 and the average reconstructed object 𝑜. For comparison, each quadrant is associated to a different ℓ value and it is normalized to its average amplitude in the central part of the pattern. The black dashed circle has a radius of 0.75.

Acknowledgments. We thank Maya

Kiskinova for the critical reading of the manuscript. Part of this work has been supported under the European Research Council (ERC) grant agreement No. 866026 (S-BaXIT). G.F.M. acknowledges funding by the European Union's Horizon 2020 research and innovation program, through the grant agreement No. 851154 (ULTRAIMAGE), from Fondazione Cariplo (NanoFast 2020.2544) and from MUR FARE PiXiE (R207A8MNNJ). C.S.B. was supported by Marie Curie Postdoctoral Fellowship DECIPHER, grant [101067016 HORIZON MSCA-2021-PF-01]. P.V. acknowledges support from the Spanish Ministry of Science and Innovation under the Maria de Maeztu Units of Excellence Programme CEX2020-001038-M and the project PID2021-123943NB-I00 (MICINN/FEDER). T.R., M.S. and M.F. acknowledge support from the French "Agence Nationale de la Recherche" under contract ANR-21-CE30-0037 (HELIMAG).

Disclosures.

The authors declare no conflicts of interest.

Data availability. Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Supplemental document.

See Supplement 1 for supporting content.

High-resolution ptychographic imaging at a seeded free-electron laser source using OAM beams: supplemental document S1. Retrieved beam and object sizes from ptychography reconstructions

The ptychography reconstructions allow obtaining a quantitative measurement of the beam and object sizes at the sample plane, since the reconstructed pixel size 𝛿𝑢 can be precisely evaluated:

))

where 𝜆 is the wavelength (18.9 nm), 𝑑 is the sample-to-CCD distance (140 mm), and 𝑙 is the size of the CCD sensor (27.648 mm). Hence, the corresponding reconstructed pixel size is 𝛿𝑢 ≈ 95.7 nm.

For fully characterizing the probe beam, it is interesting to have an estimate of the illuminating spot size. Figure S1 shows the reconstructed illumination amplitudes for ℓ = ±1, ±2, ±3. The radius of the white dashed-line circles corresponds to the radius of maximum amplitude 𝜌 ℓ,𝑚𝑎𝑥 generated at focus by a spiral zone plate (SZP) and calculated according to [1]

where 𝛾 |ℓ|-1,1 is the first root of the (|ℓ| -1)th-order Bessel function of the first kind, 𝜆 is the wavelength, 𝑓 is the focal length (164.7 mm), and 𝐷 is the diameter of the SZP (1.92 mm). It turns out the lobes of maximum amplitude in the ptychography reconstructions fall close to the calculated ideal radius of maximum amplitude: 1.24 μm for ℓ = ±1, 1.98 μm, for ℓ = ±2, and 2.65 μm for ℓ = ±3.

S2. Illumination modes from ptychography reconstructions

To highlight the relevance of the three 𝑝 𝑖 illumination modes for each ptychography reconstruction, Fig. S2(a) shows the modes' amplitude for ℓ = 0, 1, 2, 3, while Fig. S2(b) shows the numerical propagation of those modes to the detector plane via Fraunhofer integral [3]. Evidently, the main modes (𝑝 1 ) contain both the orbital angular momentum illumination and the signature of the optical elements upstream of the SZP (i.e., the central beam stop and the 2 mm input aperture of the DiProI experimental chamber). On the other hand, the intensities of the 𝑝 2 and 𝑝 3 modes are confined to the outer part of the reconstructed region, as also observed by Enders et al. [2]. This result can be understood by looking at the propagation of the 𝑝 2 mode, where the footprint of the used square-shaped order sorting aperture (downstream of the SZP) is visible. By transferring amplitude to the outer mode areas, the reconstruction algorithm can properly handle the presence of diffraction-uncorrelated intensity in the recorded diffraction patterns, without having to correlate it with features in the object transmission function [2].

S3. Effect of astigmatism on simulated datasets

To have a qualitative understanding of the impact of astigmatism in our experiment, we simulated the focal spot generated by a ℓ = -1 SZP with different illuminating wavefronts [4][5][6][7][8]. In order to mimic the experimental conditions, a SZP of focal length 𝑓 = 164.7 mm at 𝜆 = 18.9 nm was considered. To speed up the computation time, the SZP area was limited to a total diameter 𝐷 = 0.8 mm (corresponding to an outermost zone width 𝑑𝑟 = 4 μm). The central part of the illuminating beam was masked by a central beam stop (300 μm in diameter), and the emerging electric field from the SZP was numerically propagated via Fresnel integral [3] to the sample plane.

The SZP was illuminated with five different fields, as summarized in the table of Fig. S3: (a) a spherical wave with curvature radius 𝑟 𝑡 = 𝑟 𝑠 = 17 m (comparable with the Hartmann wavefront sensor measurements performed during the experiment on the direct beam), (b) an elliptical wavefront with curvature radii 𝑟 𝑡 = 17 m and 𝑟 𝑠 = 12.75 m, (c) an elliptical wavefront with curvature radii 𝑟 𝑡 = 17 m and 𝑟 𝑠 = 8.5 m, and the same elliptical wavefront rotated by (d) 𝜗 = 22.5° and by (e) 𝜗 = 45° around the propagation direction. To be noticed that the usage of two different curvature radii 𝑟 𝑡 , 𝑟 𝑠 accounts for an astigmatic wavefront, and 𝜗 is the astigmatism angle with respect to the optical system (i.e. the detector, which defines the x-and y-axis). The curvature radii are not referred to the x-and y-axis when 𝜗 ≠ 0°, hence the "tangential" (𝑡) and "sagittal" (𝑠) labelling. The amplitude and phase of these configurations at a distance 𝑓 from the SZP are shown in Fig. S3(a)-(e). In case of a spherical wave, Fig. S3(a), the amplitude displays the ideal "donut-like" shape, and the phase wraps continuously around 2𝜋. In case of an astigmatic wavefront, two bright spots appear in the intensity distribution, and, as it is apparent by comparing Fig. S3(b) and (c), their relative intensity depends on the 𝑟 𝑡 𝑟 𝑠 ⁄ ratio of the illuminating wavefront. Figure S3(b),(c) also shows that, for 𝜗 = 0°, the axis passing through the spots is rotated by 45° with respect to the (𝑥, 𝑦) reference system. On the other hand, a non-zero astigmatism angle further affects such a rotation angle, as shown in Fig. S3(d),(e). In particular, for 𝜗 = 45° (45-degree astigmatism), the two spots lay along the xaxis. This is consistent with the behavior of the reconstructed illumination amplitudes shown in Fig. 2 of the main article.

We conclude with a remark on the role of the topological charge order ℓ. Figure S4 displays the simulated intensity distributions generated by SZPs as a function of ℓ for 𝑟 𝑡 = 17 m, 𝑟 𝑠 = 8.5 m, and 𝜗 = 45°. Both a variation of the intensity distribution of the focused spots (with the appearance of |ℓ| singularities) and an increase in the distance of the peak intensity from the pattern center are visible. In particular, the latter quantity approximately grows as the ideal maximum-amplitude radius 𝜌 ℓ,𝑚𝑎𝑥 defined in Eq. (S2) even in the presence of astigmatism, as also confirmed by the experimental measurements in Fig. S1.


Nasser Rezzoug
email: nasser.rezzoug@univ-poitiers.fr

Antun Skuric

Vincent Padois

David Daney

Simulation Study of the Upper-limb Wrench Feasible Set with Glenohumeral Joint Constraints

Keywords: wrench feasible set, musculoskeletal model, glenohumeral joint, dislocation, stability

The aim of this work is to improve musculoskeletal-based models of the upper-limb Wrench Feasible Set i.e. the set of achievable maximal wrenches at the hand for applications in collaborative robotics and computer aided ergonomics. In particular, a recent method performing wrench capacity evaluation called the Iterative Convex Hull Method is upgraded in order to integrate non dislocation and compression limitation constraints at the glenohumeral joint not taken into account in the available models. Their effects on the amplitude of the force capacities at the hand, glenohumeral joint reaction forces and upper-limb muscles coordination in comparison to the original iterative convex hull method are investigated in silico. The results highlight the glenohumeral potential dislocation for the majority of elements of the wrench feasible set with the original Iterative Convex Hull method and the fact that the modifications satisfy correctly stability constraints at the glenohumeral joint. Also, the induced muscles coordination pattern favors the action of stabilizing muscles, in particular the rotator-cuff muscles, and lowers that of known potential destabilizing ones according to the literature.

26

Introduction

One fundamental assumption of human centered robot control is the ability to measure or estimate the physical capacity of humans to continuously adapts the robot assistance level based on the realtime need of its human counterpart [START_REF] Carmichael | Estimating physical assistance need using a musculoskeletal model[END_REF][START_REF] Figueredo | Human comfortability: Integrating ergonomics and muscular-informed metrics for manipulability analysis during human-robot collaboration[END_REF][START_REF] Pehlivan | Minimal assist-as-needed controller for upper limb robotic rehabilitation[END_REF]. Moreover, within the framework of computer-aided design of workstations, the knowledge of human force capacities enables to evaluate to what extent a task is in adequacy with the capacities of the operators [START_REF] Figueredo | Human comfortability: Integrating ergonomics and muscular-informed metrics for manipulability analysis during human-robot collaboration[END_REF][START_REF] Perdeaux | Investigating the role of the shoulder musculature during maximum unilateral isometric exertions[END_REF], to define ergonomic criteria of discomfort [START_REF] Maurice | Human-oriented design of collaborative robots[END_REF] and to implement models of muscular fatigue [START_REF] Ma | A new muscle fatigue and recovery model and its ergonomics application in human simulation[END_REF][START_REF] Savin | Movement variability and digital human models: Development of a demonstrator taking the effects of muscular fatigue into account[END_REF][START_REF] Savin | Evidence of movement variability patterns during a repetitive pointing task until exhaustion[END_REF]. Different methodologies have been proposed to assess human force capacities.

Data based approaches rely on the availability of experimental data [START_REF] Khalaf | A normative database of isokinetic upperextremity joint strengths: Towards the evaluation of dynamic human performance[END_REF][START_REF] Guenzkofer | Elbow torque ellipses: investigation of the mutual influences of rotation, flexion, and extension torques[END_REF][START_REF] Lannersten | Isometric strength in flexors, abductors, and external rotators of the shoulder[END_REF][START_REF] Hall | A comparison between measured female linear arm strengths and estimates from the 3d static strength prediction program (3dsspp)[END_REF][START_REF] Kotte | Normative values of isometric elbow strength in healthy adults: a systematic review[END_REF][START_REF] Delfa | The 'arm force field' method to predict manual arm strength based on only hand location and force direction[END_REF] which are numerous but limited to particular postures. Whitin the context of model based approaches, several algorithms using musculoskeletal models, predominantly of the upper-limb, are available [START_REF] Skuric | On-line feasible wrench polytope evaluation based on human musculoskeletal models: An iterative convex hull method[END_REF][START_REF] Valero-Cuevas | A Mathematical Approach to the Mechanical Capabilities of Limbs and Fingers[END_REF][START_REF] Carmichael | Estimating physical assistance need using a musculoskeletal model[END_REF][START_REF] Ingram | Improving anterior deltoid activity in a musculoskeletal shoulder model -an analysis of the torque-feasible space at the sternoclavicular joint[END_REF].

From posture, muscular forces and muscles moment arms, the set of achievable wrench at the hand called the Wrench Feasible Set (WFS) can be computed under the form of a convex polytope. Its main interest is that it allows an exhaustive evaluation of force capacities whatever the direction of force application or the posture. The base assumption of all these algorithms is that muscular forces are the limiting factor. However, other factors can limit force capacities such as joint stability (Blache et al., 2017a,b;[START_REF] Ingram | Improving anterior deltoid activity in a musculoskeletal shoulder model -an analysis of the torque-feasible space at the sternoclavicular joint[END_REF], safety mechanisms implemented by the central nervous system for preventing joint injuries [START_REF] Assila | Glenohumeral joint and muscles functions during a lifting task[END_REF][START_REF] Maurice | Human-oriented design of collaborative robots[END_REF] or balance constraints [START_REF] Perdeaux | Investigating the role of the shoulder musculature during maximum unilateral isometric exertions[END_REF]. When considering the upper-limb and the glenohumeral (GH) joint, the stability issue is of paramount importance [START_REF] Lippitt | Glenohumeral stability from concavity-compression: A quantitative analysis[END_REF][START_REF] Dickerson | Experimental evaluation of a computational shoulder musculoskeletal model[END_REF][START_REF] Labriola | Stability and instability of the glenohumeral joint: The role of shoulder muscles[END_REF]Blache et al., 2017a;[START_REF] Sinha | Use of botulinum a toxin in irreducible shoulder dislocation caused by spasm of pectoralis major[END_REF]. Indeed, this joint has the greatest range of motion of all the joints of the human body. This is due to the relative laxity of the shoulder ligaments, but especially to the shallow conformation of the humeral head on the glenoid fossa, which surface is relatively flat. The side effect is that the GH joint is prone to dislocation when the shear component of the GH Joint Reaction Force (JRF) exceeds a certain percentage of its compressive component [START_REF] Lippitt | Glenohumeral stability from concavity-compression: A quantitative analysis[END_REF][START_REF] Dickerson | Experimental evaluation of a computational shoulder musculoskeletal model[END_REF]. To prevent this phenomenon, several muscles including the rotator cuff muscles act to stabilize the GH joint by orienting the JRF toward the middle of the glenoid fossa while pressing the humeral head against it [START_REF] Lippitt | Glenohumeral stability from concavity-compression: A quantitative analysis[END_REF]. Therefore, the induced muscle coordination pattern may have a significant effect on strength capacities and, if not accounted for, may lead to a potential overestimation of the WFS. However, to our knowledge, the only works implementing GH non dislocation constraints considered the shoulder torque feasible set only [START_REF] Ingram | Improving anterior deltoid activity in a musculoskeletal shoulder model -an analysis of the torque-feasible space at the sternoclavicular joint[END_REF] or were used to asses shoulder muscles forces with inverse dynamics and static optimization [START_REF] Dickerson | Experimental evaluation of a computational shoulder musculoskeletal model[END_REF]Blache et al., 2017b,a) or with EMG assisted optimization [START_REF] Assila | Emg-assisted algorithm to account for shoulder muscles co-contraction in overhead manual handling[END_REF][START_REF] Assila | Glenohumeral joint and muscles functions during a lifting task[END_REF].

To fill this gap and to improve models of WFS computed from musculoskeletal models, the recent iterative convex hull method (ICHM) [START_REF] Skuric | On-line feasible wrench polytope evaluation based on human musculoskeletal models: An iterative convex hull method[END_REF] is upgraded in order integrate two important limiting factors of force capacities i.e. GH joint non dislocation and limitation of GH JRF compression component. The core of the ICHM consists in finding the vertices of the WFS by solving a linear programming problem. To ensure GH joint non dislocation, additional constraints take the form of inequalities involving the admissible ratio between shear and compressive components of the GH JRF [START_REF] Lippitt | Glenohumeral stability from concavity-compression: A quantitative analysis[END_REF][START_REF] Dickerson | Experimental evaluation of a computational shoulder musculoskeletal model[END_REF]Blache et al., 2017b;[START_REF] Ingram | Improving anterior deltoid activity in a musculoskeletal shoulder model -an analysis of the torque-feasible space at the sternoclavicular joint[END_REF]. Also additional constraints bound the GH JRF compression. In this context, the objective of this work is to test in silico the effects on the WFS of the ICHM implementation incorporating the GH constraints.

In particular, the following points are verified:

• P 1 : the decrease of the range of maximal achievable forces,

• P 2 : the orientation of the GH JRF toward the middle of the glenoid fossa,

• P 3 : the decrease of the amplitude of the GH JRF compression,

• P 4 : the modifications of muscle force coordination pattern with an expected increase in the contribution of the rotator cuff (Blache et al., 2017a) and biceps muscles [START_REF] Labriola | Stability and instability of the glenohumeral joint: The role of shoulder muscles[END_REF] and a decrease in that of the deltoids and pectoralis major muscles [START_REF] Labriola | Stability and instability of the glenohumeral joint: The role of shoulder muscles[END_REF].

They are tested by considering a validated musculoskeletal model of the upper-limb [START_REF] Saul | Benchmarking of dynamic simulation predictions in two software platforms using an upper limb musculoskeletal model[END_REF] in three different postures.

Methods

Wrench Feasible Set definition

Let us consider a serial kinematic chain C with n rotational degrees of freedom (dofs) representing the upper-limb. Its configuration is specified by the joint angle vector q ∈ R n and the joint velocity vector q ∈ R n .τ ∈ R n refers to the joint torque vector. The task space is of dimension m = 3, the end effector position is defined by x ∈ R m and the end-effector wrench by F ∈ R m . C is actuated by d muscles (d > n) producing muscular tension force t ∈ R d . The joint torque τ ∈ R n , generated by the muscle force t ∈ R d at joint configuration q ∈ R n , can be calculated as follows:

τ = -L T (q)t (1)
L(q) T ∈ R n×d is the transpose of the moment arm matrix, which is defined as the muscle length Jacobian relating the space of joint and muscle length velocities

l = L(q) q, L ij = ∂l i ∂q j (2)
where l ∈ R d is the muscle lengthening velocity and q ∈ R n the joint angular velocity. The negative sign in equation ( 1) makes the force applied in the length shortening direction of the muscle positive.

The muscles have a limited force capacity and each component t i = 1..d of t is to lie within an interval t i , t i . Also, 0 ≤ t i < t i because a muscle can only pull and not push. The tension t i developed by muscle i is defined according to a Hill type model [START_REF] Saul | Benchmarking of dynamic simulation predictions in two software platforms using an upper limb musculoskeletal model[END_REF].

A general achievable set of muscle-tendon forces t ∈ t, t forms a d-dimensional hyper-rectangle with side lengths equal to the ranges of each of the d muscle forces. Using the moment arm matrix L(q)

and equation ( 1) this hyper-rectangle can be projected into the n-dimensional space of joint torques τ , forming the convex polytope of achievable joint torques P τ

P τ = τ ∈ R n | τ = -L T (q)t, t ∈ t, t (3)
Furthermore, the dual relationship of the m-dimensional Cartesian wrenches F and the generalised joint torques τ is given through the Jacobian transpose matrix J T (q)F = τ and defines the achievable Cartesian feasible wrench:

P F = F ∈ R m | J T (q)F = τ , τ ∈ P τ (4)
Combining (3) and (4) the implicit definition of P F is

P F = F ∈ R m | J T (q)F = -L T (q)t, t ∈ t, t (5)

Outline of the ICHM algorithm

The ICHM is a recent algorithm for the feasibility set analysis of a generic class of linear algebra problems. It is suited for the determination of the feasible Cartesian WFS P F associated to a musculoskeletal model of the human upper limb (5) [START_REF] Skuric | On-line feasible wrench polytope evaluation based on human musculoskeletal models: An iterative convex hull method[END_REF]. It is based on two processes. The first one consists in finding the vertices of P F by solving a linear programming (LP) problem for a chosen set of directions c ∈ R m in the Cartesian task space. The second one consists in computing the convex hull of the current vertices in order to define a new directions set c orthogonal to the current polytope faces. This process is iterated until the distance between the new vertices and the current polytope faces are below a predefined threshold. In the case of the WFS P F determination, the LP problem solved during the vertex search for appropriately chosen c directions is the following: r) forms an orthonormal basis of Ker(J(q)) = Im(J T (q)) ⊥ . It is obtained from a singular value decomposition of

max t -c T J T + (q)L T (q)t s.t. V T 2 L T (q)t = 0 t ≤ t ≤ t (6) V 2 ∈ R n×(n-
J (q) = U ΣV T where V = [V 1 , V 2 ] T .
The null space basis of J(q) is spanned by the columns of V corresponding to singular values equal to zero and stored in V 2 . r is the rank of J(q) and J T + (q) is the Moore-Penrose pseudo-inverse of the J (q) T .

The LP problem (6) finds a vertex or face of P F furthest to a face of the current polytope which normal is c by finding the maximum of the projection of a wrench F = -J T + (q)L T (q)t on the c direction. The equality constraints V T 2 L T (q)t = 0 ensure that L T (q)t ∈ Im(J T ) and therefore that F = -J T + (q)L T (q)t is an exact solution for the wrench F in the equation J T (q)F = -L T (q)t. In other words, the equality constraints guarantee that the wrench F can be generated by the joint torque τ = -L T (q)t such that F = J T + (q)τ . The LP formulation is very flexible and allows to add new constraints relative to the GH non dislocation involving the JRF.

Adding GH joint non dislocation constraints to the vertex search algorithm of the ICHM

Formulation in terms of GH Joint Reaction Force

Due to the particular conformation of the humerus and the glenoid fossa of the scapula, GH joint dislocation can occur if the JRF exerted by the humerus head on the the glenoid fossa is outside a particular irregular conic region (Figure 1). Its boundaries can be defined by considering the admissible ratio µ j between the shear and compression components of the JRF which are parallel and orthogonal to the plane of the glenoid fossa [START_REF] Dickerson | Experimental evaluation of a computational shoulder musculoskeletal model[END_REF][START_REF] Lippitt | Glenohumeral stability from concavity-compression: A quantitative analysis[END_REF], respectively. The interior of the polyhedral conic region is defined by a serie of affine constraints represented by half planes (Figure 2). For a direction j defined by an angle θ j within the glenoid fossa plane (y GH , z GH ), the constraint is expressed as following:

µ j cosθ j sinθ j GF F H/S ≤ 0 (7)
where GF F H/S is the Joint Reaction Force exerted by the humerus head (H) on the scapula glenoid fossa (S) expressed in the frame of the glenoid fossa (superscript GF). The first component of GF F H/S = [F n , F t1 , F t2 ] is compression and the other two are shear forces. F n along x GF is always negative due to the choice of the glenoid fossa frame (Figure 1).

For a set of directions, a linearized version of the cone boundaries is obtained by stacking the constraints (7) in the matrix C DIS :

C DIS GF F H/S ≤ 0 (8)
--please insert Figure 1 around here ----please insert Figure 2 around here --

Incorporating muscular and external forces

Let consider the static equilibrium of the humerus expressed in the global frame (denoted by the superscript 0). The weight of the segments is neglected.

0 F S/H + 0 F F/H + Σ 0 F M/H = 0 (9)
0 F S/H and 0 F F/H represents the forces exerted on the humerus by the scapula and forearm, respectively, Σ 0 F M/H represents the action of forces (expressed in the global frame) of muscles acting on the humerus and crossing the GH joint. An external force 0 F EXT exerted at the hand is transmitted to the humerus according to the third Newton's law. Then, 0 F F/H = 0 F EXT . 0 F H/S the Joint Reaction

Force exerted by the humerus on the scapula in the global frame can be expressed as following:

0 F H/S = Σ 0 F M/H + 0 F EXT (10)
The contribution of muscular forces to the GH JRF, Σ 0 F M/H , can be replaced by 0 C JRF (q)t with 0 C JRF ∈ R m×d . Each column of 0 C JRF (q), called the force direction matrix [START_REF] Ingram | Improving anterior deltoid activity in a musculoskeletal shoulder model -an analysis of the torque-feasible space at the sternoclavicular joint[END_REF], stores the contribution n i of a 1N force of muscle i to the GH JRF:

0 C JRF (q) = n 1 ... n i ... n d i = 1..d (11)
Combining the action of muscles and that of 0 F EXT , the GH JRF expression in the global frame is :

0 F H/S = 0 C JRF (q)t + 0 F EXT (12)
It is considered that 0 F EXT acting at the hand is counterbalanced by one element F ∈ P F . Since the equality constraints in the LP (6) ensure that L T (q)t ∈ Im(J T ); for a certain t, 0 F EXT can be expressed as:

0 F EXT = J T + (q)L T (q)t (13)
Then ( 12) becomes:

0 F H/S = ( 0 C JRF (q) + J T + (q)L T (q))t (14)
Expressed in the frame of the glenoid fossa ( 14) becomes:

GF F H/S = GF R 0 0 F H/S GF F H/S = GF R 0 ( 0 C JRF (q) + J T + (q)L T (q))t (15)
where GF R 0 is the rotation matrix from the glenoid fossa frame to the global frame.

Constraints and LP formulation

The GH non dislocation constraints are obtained by substituting GF F H/S in ( 8) by its expression in (15):

C DIS GF R 0 ( 0 C JRF (q) + J T + (q)L T (q))t ≤ 0 Ct ≤ 0 (16)
In addition, the amplitude of the first component of GF F H/S corresponding to compression is kept under a maximum value by considering the following constraint:

Dt -f max ≤ 0 ( 17
)
where D is the first row of GF R 0 ( 0 C JRF (q) + J T + (q)L T (q)) and f max the maximum allowed compressive force.

Finally, the formulation of the LP problem incorporating GH non dislocation constraints and limitation of GH JRF compression is:

max t -c T J T + (q)L T (q)t s.t. V T 2 L T (q)t = 0 Ct ≤ 0 Dt -f max ≤ 0 t ≤ t ≤ t (18)

Implementation details

The seven degrees of freedom and 50 muscles Opensim upper-limb model proposed in [START_REF] Saul | Benchmarking of dynamic simulation predictions in two software platforms using an upper limb musculoskeletal model[END_REF] was used. The Jacobian matrix, muscles moment arm matrix, maximal and minimal muscle forces were obtained with the Matlab Opensim API.

Using the calcReactionOnParentExpressedInGround method of the Opensim joint object, each column n i of C JRF (q) was obtained by dividing 0 F H/Si = t i n i , the individual contribution of muscle i = 1..d to the GH JRF, by the value of the muscle tension t i :

n i = 0 F H/Si t i (19)
The contribution n i to the GH JRF was computed for each muscle by setting its activation a to a random value in the interval ]0, 1], to ensure that t i > 0, and by setting the force of the other muscles to 0 (both passive and active components). In order to quantify the static GH JRF, the rotational inertia was set to huge values to prevent joint rotational acceleration and acceleration of the humerus center of gravity allowing to obtain individual contribution of a single muscle to the GH JRF in a static condition.

Since (8) was expressed in the frame of the glenoid fossa, the rotation matrix from this frame to the global frame was needed. Three markers were attached to the scapula of the chosen Opensim model (Figure 3) and the glenoid fossa frame was defined according to the method proposed in (Blache et al., 2017b):

--please insert Figure 3 around here --

x GF = (gf 2 -gf 3 ) × (gf 1 -gf 3 ) ∥(gf 2 -gf 3 ) × (gf 1 -gf 3 )∥ y GF = gf1+gf2 2 -gf 3 ∥ gf1+gf2 2 -gf 3 ∥ z GF = x GF × y GF 0 R GF = x GF y GF z GF GF R 0 = 0 R GF T (20)
The x GF axis corresponded to the direction of the compressive component of the GH JRF which was always negative because the humeral head was pressing against the glenoid fossa. The y GF and z GF axes corresponded to shear forces because it was assumed that the glenoid fossa was almost flat.

The µ i ratios were defined such that the constraint boundary was an ellipse with a ratio of 0.5 in the cranial-caudal directions and 0.3 in the anterior-posterior directions [START_REF] Dickerson | Experimental evaluation of a computational shoulder musculoskeletal model[END_REF][START_REF] Lippitt | Glenohumeral stability from concavity-compression: A quantitative analysis[END_REF]. f max was set to 1500N according to [START_REF] Assila | Emg-assisted algorithm to account for shoulder muscles co-contraction in overhead manual handling[END_REF]Blache et al., 2017a).

Results

The WFS was computed using the ICHM algorithm with and without the GH constraints. Three postures with different levels of arm abduction and elbow flexion were considered (Table 1 and Figure 4). They were chosen such that a range of GH abduction was covered and also because the third posture corresponded to the so-called "apprehension test" for shoulder luxation testing [START_REF] Mcmahon | Muscles may contribute to shoulder dislocation and stability[END_REF][START_REF] Mcmahon | A novel cadaveric model for anterior-inferior shoulder dislocation using forcible apprehension positioning[END_REF]. The forearm pronation-supination, wrist flexion-extension and radial-ulnar deviation were set to 0°.

--please insert Table 1 around here ----please insert Figure 4 around here ----please insert Figure 5 around here ----please insert Figure 6 around here --

The results were divided in three parts. Firstly, the WFS with and without GH non dislocation constraints were compared to verify P 1 . In order to assess the global difference between WFS with and without GH constraints, the following procedure was used. The two WFS were superimposed. A set of rays originating from the center of the WFS with GH constraints was constructed. They were defined by a set of direction v(r, θ, φ), with r = 1, θ ∈ [0, 360°] and φ ∈ [0, 180°], which spherical coordinates vary with an increment of 1°. Each ray intersects [START_REF] Tuszynski | Triangle/ray intersection[END_REF] with both WFS and the "length" (in N) of the segment between the two points of intersection was computed for each ray. The graphical representation consisted of a sphere projected in 2D whose surface was colored according to the "length" of the segments corresponding to the v direction (Figure 6). In addition, the difference between the minimal and maximal forces for each polytope and the segments "length"

for some particular directions in the global reference frame (front, back, up, down, right and left) was also assessed (Table 2).

In the second part, the GH JRF corresponding to each vertex of the WFS preventing GH dislocation was obtained from ( 14) and the constraints satisfaction was checked by computing the ratio between shear and compressive components (P 2 ). The limit of compression forces was also checked (P 3 ).

In the third part, the coordination of the shoulder muscles was investigated to verify P 4 . A particular attention was paid to the forces of the deltoid and pectoralis major muscles on one side and rotator cuff and biceps muscles on the other side, which were known to have a destabilizing and stabilizing action on the GH joint, respectively [START_REF] Labriola | Stability and instability of the glenohumeral joint: The role of shoulder muscles[END_REF][START_REF] Sinha | Use of botulinum a toxin in irreducible shoulder dislocation caused by spasm of pectoralis major[END_REF].

Effect of GH joint contraints on the WFS

The GH constraints reduced the WFS amplitude. The maximal forces were reduced by 66, 45 and 71 N for Posture 1, 2, and 3 respectively. Also, the maximal reduction were noticed when both non dislocation and compression constraints were ON. For posture 1; the main reductions were for front, up and left (55, 28 and 20 N respectively). For posture 2 it was front, left and down (124, 99 and 45 N, respectively) and for posture 3 front and down (66 and 24N). The global difference, of the sphere projected in 2D is given in Figure 5.

--please insert Table 2 around here --

GH Constraints satisfaction

The ratio of the GH JRF shear forces divided by the compressive component superimposed with the boundary of the GH non dislocation constraints for the three chosen postures are presented in Figure 7. Without the non dislocation constraints, the developed muscular forces conducted to a dislocation of the GH joint for 61%, 55%, and 76% of the vertices of the corresponding WFS for postures 1, 2, and 3, respectively (ratios depicted by red circle lying outside of the GH boundaries). When GH constraints were ON, these percentages dropped to 5%, 2%, and 3%. For posture 1, the GH JRF vectors are directed upward (Figure 7.A), anteriorly for postures 2 and in the anterior-inferior direction for posture 3 (Figures 7.B and 7.C). When the GH constraints are active, the ratios between shear and compressive lye inside the constraints boundaries validating P 2 . The JRF depicted in Figure 8 confirms the orientation of the dislocating JRF and show the effect of the constraints limiting the amplitude of compressive forces (blue thick arrows in Figure 8).

--please insert Figure 7 around here ----please insert Figure 8 around here ----please insert Figure 9 around here ----please insert Figure 10 around here --

Muscle force coordination

The effect of GH constraints on shoulder muscles force coordination was investigated. The pectoralis major, anterior, middle and posterior parts of deltoid were considered for their destabilizing effect and the rotator cuff muscles (subscapularis, intraspinatus, supraspinatus and teres minor) and biceps (for high GH abduction values) for their stabilizing contribution [START_REF] Labriola | Stability and instability of the glenohumeral joint: The role of shoulder muscles[END_REF][START_REF] Sinha | Use of botulinum a toxin in irreducible shoulder dislocation caused by spasm of pectoralis major[END_REF].

Firstly, the maximum compression constraint was relaxed in order to compare the sole effect of the non dislocation constraint i.e. the direction of the GH JRF within the glenoid fossa. For the WFS without GH constraints, the mean and standard deviation of muscular forces was computed separately for the vertices conducting to GH dislocation (red circles in Figure 7) and those not conducting to GH dislocation (green stars in Figure 7). The mean and standard deviation of muscular forces was assessed for all the vertices of the WFS with GH non dislocation constraints. The results are depicted in Figure 9 for the three postures.

Discussion

Force amplitudes

The consideration of GH contraints for avoiding dislocation and too excessive joint load had an effect on the WFS at the hand. A decrease of the maximal forces was noticed particularly in the anterior and mediolateral left directions and depended on the posture.

The non dislocation constraints had a moderate effect on the maximal forces at the hand but resulted in a modification of muscle coordination with a modulation of tension that ensured GH non dislocation. This illustrate the muscular redundancy implemented in the musculoskeletal model allowing a redistribution of muscular forces conducting stabilize the GH joint without affecting significantly the maximal forces. On the other hand, the maximum compression constraint contributed more substantially to the lowering of forces at the hand by a global decrease of all muscle tensions (Figure 10).

Muscle coordination

Apart from generating joint movements the muscles action is also directed toward the stabilization of the joints. This is especially true for the GH joint. Our results are coherent with those of the literature with a decrease of the force amplitude of the deltoids and pectoralis major in the three postures (Figure 9) with a particular implication of pectoralis muscle exemplified by the effect of its spasms in GH dislocation [START_REF] Sinha | Use of botulinum a toxin in irreducible shoulder dislocation caused by spasm of pectoralis major[END_REF]. In our results, it is particularly noticeable for this muscle with a decrease of the mean tension of 216 N, 575N and 746N corresponding to 27%, 54% and 77% when comparing the WFS with GH non dislocation constraints to the WFS vertices conducting to dislocation for posture 1, 2, and 3 respectively. Moreover, in posture 3, corresponding to the apprehension test, it can be seen that without the proposed constraints, the JRF is directed toward the anterior-inferior direction. This observation is consistent with the literature indicating that GH dislocation occurs in this direction when performing forceful apprehension test on cadaveric upper-limbs [START_REF] Mcmahon | Muscles may contribute to shoulder dislocation and stability[END_REF][START_REF] Mcmahon | A novel cadaveric model for anterior-inferior shoulder dislocation using forcible apprehension positioning[END_REF].

Limitations

The present model has several limitations. Firstly, only muscular contribution to joint non dislocation was considered but not the static restraints such as capsulolabrum, articular surface geometry and intracapsular pressure. However, it is hypothesized that the non dislocation ratios consider their contribution in some way [START_REF] Lippitt | Glenohumeral stability from concavity-compression: A quantitative analysis[END_REF][START_REF] Dickerson | Experimental evaluation of a computational shoulder musculoskeletal model[END_REF]. Secondly, the WFS is computed in an isometric condition neglecting the segments weight and movement dynamics. To relax these hypotheses the equation relating muscle tension t and F should be replaced by the following one:

J T (q)F = -L(q)t -τ b (21)
where τ b is the bias joint torque produced by the gravity and the model dynamics. Results on muscles tension are global and are averaged for all the WFS vertices and therefore their standard deviation is high. It would be valuable to differentiate the results by considering the Cartesian directions.

However, the global presentation was chosen for better synthesizing the results. Another limitation is that the results heavily depend on the quality of the musculoskeletal model in terms of muscle force and moment arm matrix. The tested model was generic without any scaling and was used to provide a proof of concept for the proposed method. Finally, an experimental validation is mandatory and will be conducted on the basis of the experiment carried out in [START_REF] Rezzoug | Upper-limb isometric force feasible set: Evaluation of joint torque-based models[END_REF].

Conclusion

The aim of this work was to propose an extension of the ICHM to assess the upper-limb WFS taking into account the limiting factor of GH joint stability. The adaptations of the ICHM were described and implemented. The results confirmed the lowering of the WFS due to the association of GH constraints on the JRF. The non dislocation constraints had an effect on the coordination of muscular forces while compression constraints limited muscle tensions. The change in muscle tension were investigated and in global agreement with the literature, the activity of some muscular groups was lowered or increased. Taking into account biomechanical constraints such as those implemented on the GH joint allow to improve models of musculoskeletal capacities for better simulating human behavior in various application such as collaborative robotics and ergonomics.

Figure 2 :

 2 Figure 2: Glenohumeral dislocation constraints boundaries defined by a set of half-planes.

Figure 3 :

 3 Figure 3: Glenoid fossa reference frame.

Figure 4 :

 4 Figure 4: Considered upper-limb postures.

Figure 5 :

 5 Figure 5: WFS and skeletal model representation for the three postures from two different points of view. Light (cyan): WFS without GH constraints, Dark (magenta): WFS with GH constraints.

Figure 6 :

 6 Figure 6: Difference (N) between polytopes according to Cartesian directions represented by their longitude and latitude (rad) for the three postures A.: posture 1, B. posture 2, and C. posture 3.

Figure 7 :Figure 8 :

 78 Figure 7: GH constraints verification: for one WFS vertex, the abscissa and ordinate of each 2D point correspond to the ratio of the shear component of the GH JRF around the Z and Y axis of the glenoid fossa frame, respectively, divided by the compressive component around the X axis. The red ellipse in each panel represent the GH non dislocation constraints boundaries. The red circles correspond to vertices of the WFS computed without GH constraints conducting to GH dislocation, the green stars to vertices of the WFS without GH constraints not conducting to GH dislocation, and finally the blue crosses correspond to vertices of a WFS with GH constraints on. Panels A, B, and C present the results for postures 1, 2, and 3 respectively

Figure 9 :

 9 Figure 9: Mean and standard deviation of muscles tension according to posture and GH non dislocation constraints status. The GH maximum compression constraints are off. A, B, and C stand for posture 1, 2, and 3, respectively. PCEM: pectoralis major, DELT: sum of anterior, middle and posterior deltoid, ROTCUFF: sum of rotator cuff muscles: supraspinatus, infraspinatus, teres minor, and subscapularis, TMAJ: teres major, LAT: latissimus dorsi, CORB: coracobrachialis, TRIlong : long

Figure 10 :

 10 Figure 10: Muscles tension according to posture and GH compression constraints: blue (dark)constraint is off, orange(light) constraint is on. GH non dislocation constraints status is ON for both cases. A, B, and C stands for posture 1, 2, and 3, respectively. PCEM: pectoralis major, DELT: mean of anterior, middle and posterior deltoid, ROTCUFF: mean of rotator cuff muscles: supraspinatus, infraspinatus, teres minor, and subscapularis, TMAJ: teres major, LAT: latissimus dorsi, CORB: coracobrachialis, TRIlong : long head of the triceps brachii, BIC: mean of biceps short and long heads.

Figure 3 :

 3 Figure 3: Glenoid fossa reference frame

Figure 6 :Figure 7 :

 67 Figure 6: Difference (N) between polytopes according to Cartesian directions represented by their longitude and latitude (rad) for the three postures A.: posture 1, B. posture 2, and C. posture 3

Figure 8 :

 8 Figure 8: Each arrow depicts the GH JRF corresponding to one vertex of a WFS. The thin magenta arrows correspond to a WFS with GH constraints off and the thick blue arrows to a WFS with GH constraints on. Panels A, B, and C present the results for posture 1, 2, and 3 respectively. The left, middle and right figures present the results in the frontal, sagittal, and coronal planes, respectively.

Figure 9 :

 9 Figure 9: Mean and standard deviation of muscles tension according to posture and GH non dislocation constraints status. The GH maximum compression constraints are off. A, B, and C stand for posture 1, 2, and 3, respectively. PCEM: pectoralis major, DELT: sum of anterior, middle and posterior deltoid, ROTCUFF: sum of rotator cuff muscles: supraspinatus, infraspinatus, teres minor, and subscapularis, TMAJ: teres major, LAT: latissimus dorsi, CORB: coracobrachialis, TRIlong : long head of the triceps brachii, BIC: mean of biceps short and long heads.

Figure 10 :

 10 Figure10: Muscles tension according to posture and GH compression constraints: blue (dark) constraint is off, orange (light) constraint is on. GH non dislocation constraints status is ON for both cases. A, B, and C stands for posture 1, 2, and 3, respectively. PCEM: pectoralis major, DELT: mean of anterior, middle and posterior deltoid, ROTCUFF: mean of rotator cuff muscles: supraspinatus, infraspinatus, teres minor, and subscapularis, TMAJ: teres major, LAT: latissimus dorsi, CORB: coracobrachialis, TRIlong : long head of the triceps brachii, BIC: mean of biceps short and long heads.

Table 1 :

 1 Glenohumeral and elbow joint angles in degrees for the selected postures. EL: Elevation angle, SE: shoulder elevation, SR: Shoulder rotation, and EF: Elbow flexion.

	Posture	EL (°) SE (°) SR (°) EF (°)
	Posture 1	10	20	0	70
	Posture 2	0	70	-30	90
	Posture 3	0	90	-90	90

Table 2 :

 2 Force difference (N) between the WFS without and with GH constraints according to posture, constraints and directions. disloc : dislocation constraints, comp : compression constraints, min and max: directions of the minimal and maximal force of the WFS without GH constraints, respectively, front, back, up, down, right and left correspond to directions in the global reference frame.

	Posture disloc comp min max front back up down right left
	1	ON	OFF	0	0	0	0	0	0	0	4
	1	ON	ON	9	66	55	3	28	0	10	20
	2	ON	OFF	1	0	22	2	0	13	1	46
	2	ON	ON	2	45	124	8	1	45	2	99
	3	ON	OFF	0	7	48	0	-2	0	3	-1
	3	ON	ON	1	41	66	2	7	24	7	2

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work has been partially funded by the BPI France Lichie project.