Julien Deparday

Patrick Bot

Benoit Augier

Marc Rabaud

Dario Motta

Modal Analysis of Pressures on a Full-Scale Spinnaker

ρ Density of air (1.25 kg/m 3 ) S Sail Area of the asymmetrical Spinnaker (65 m 2 )

While sailing offwind, the trimmer typically adjusts the downwind sail "on the verge of luffing", letting occasionally the luff of the sail flapping. Due to the unsteadiness of the spinnaker itself, maintaining the luff on the verge of luffing needs continual adjustments. The propulsive force generated by the offwind sail depends on this trimming and is highly fluctuating. During a flapping sequence, the aerodynamic load can fluctuate by 50% of the average load. On a J/80 class yacht, we simultaneously measured timeresolved pressures on the spinnaker, aerodynamic loads, boat and wind data. Significant spatio-temporal patterns are detected in the pressure distribution. In this paper we present averages and main fluctuations of pressure distributions and of load coefficients for different apparent wind angles as well as a refined analysis of pressure fluctuations, using the Proper Orthogonal Decomposition (POD) method. POD shows that pressure fluctuations due to luffing of the spinnaker can be well represented by only one proper mode related to a unique spatial pressure pattern and a dynamic behavior evolving with the Apparent Wind Angles. The time evolution of this proper mode is highly correlated with load fluctuations. Moreover, POD can be employed to filter the measured pressures more efficiently than basic filters. The reconstruction using the first few modes allows to restrict to the most energetic part of the signal and remove insignificant variations and noises. This might be helpful for comparison with other measurements and numerical simulations.

1 julien.deparday@ecole-navale.fr

INTRODUCTION

In research and development in sail aerodynamics, full-scale testing, wind tunnel testing and numerical simulation have always been complementary. Numerical simulation allows to efficiently investigate different designs without the cost of creating sails [START_REF] Chapin | A Multidisciplinary Computational Framework for Sailing Yacht Rig Design & Optimization through Viscous FSI[END_REF][START_REF] Durand | FSI investigation on stability of downwind sails with an automatic dynamic trimming[END_REF][START_REF] Ranzenbach | Mainsail Planform Optimization for IRC 52 Using Fluid Structure Interaction[END_REF][START_REF] Viola | A CFD-based wing sail optimisation method coupled to a VPP[END_REF]. Nowadays, advanced computational resources have enhanced numerical simulation and have allowed to couple fluid and structural solvers to create Fluid-Structure Interaction simulations [START_REF] Chapin | A Multidisciplinary Computational Framework for Sailing Yacht Rig Design & Optimization through Viscous FSI[END_REF][START_REF] Durand | FSI investigation on stability of downwind sails with an automatic dynamic trimming[END_REF][START_REF] Ranzenbach | Mainsail Planform Optimization for IRC 52 Using Fluid Structure Interaction[END_REF][START_REF] Augier | Numerical study of a flexible sail plan submitted to pitching: Hysteresis phenomenon and effect of rig adjustments[END_REF][START_REF] Lombardi | A strongly coupled fluid-structure interaction model for wind-sail simulation[END_REF][START_REF] Renzsch | Fluid Structure Interaction Simulation of Spinnakers-Towards Simulation Driven Sail Design[END_REF][START_REF] Trimarchi | Wrinkle development analysis in thin sail-like structures using MITC shell finite elements[END_REF]. However wind tunnel testing and full-scale testing are required for comparison and validation [START_REF] Hansen | Comparison of Wind Tunnel and Full-Scale Aerodynamic Sail Force measurements[END_REF][START_REF] Renzsch | An experimental validation case for fluid-structure-interaction simulations of downwind sails[END_REF][START_REF] Viola | Sail pressures from fullscale, wind-tunnel and numerical investigations[END_REF]. Wind tunnel testing has the advantage to be in a controlled environment where a balance can be used to measure the forces created by the sails on the boat frame (Campbell, 2014a[START_REF] Flay | A twisted flow wind tunnel for testing yacht sails[END_REF][START_REF] Graf | Investigation of the Flying Shape of Spinnakers in a Twisted Flow Wind Tunnel, 19th Cheasapeake Sailing yacht Symposium[END_REF][START_REF] Zasso | Twisted Flow Wind Tunnel Design for Yacht Aerodynamic Studies[END_REF]. Those results can easily be used to create a Velocity Prediction Program (Campbell, 2014b, Le Pelley and[START_REF] Le Pelley | Effective Wind Tunnel Tesing of Yacht Sails Using a Real-Time Velocity Prediction Program[END_REF]. Nevertheless with wind tunnel testing, some rules of similitude are violated as the Reynolds number, or the ratio of fabric weight to wind pressure or the ratio of membrane stress to wind pressure. Full-Scale testing does not have those issues, and permits to determine yacht performance in real sailing conditions. Those experiments need complex set-up in a harsh environment but actual aerodynamic loads can be assessed in a variety of ways. Sail boat dynamometers [START_REF] Herman | A sail force dynamometer: design, implementation and data handling[END_REF][START_REF] Hochkirch | Fullscale hydrodynamic force measurement on the Berlin sailing dynamometer[END_REF][START_REF] Masuyama | The work achieved with the sail dynamometer boat "Fujin", and the role of full scale tests as the bridge between model tests and CFD[END_REF] measured forces from upwind sails transmitted to the boat frame. [START_REF] Fossati | A novel full scale laboratory for yacht engineering research[END_REF] created a sail boat dynamometer with the possibility of measuring aerodynamic forces of downwind sails. [START_REF] Augier | Experimental validation of unsteady models for fluid structure interaction: Application to yacht sails and rigs[END_REF] carried out experiments where loads on the rigging lines and sails were measured. They contributed to a better comprehension on interaction between the wind, the rigging and the sails. Le Pelley et al. (2015) measured the forces and the directions on the three corners of spinnakers. Le Pelley et al. (2012), [START_REF] Lozej | Pressure distribution on sail surfaces in real sailing conditions[END_REF], [START_REF] Motta | Experimental investigation of asymmetric spinnaker aerodynamics using pressure and sail shape measurements[END_REF], [START_REF] Viola | On-water pressure measurements on a modern asymmetric spinnaker[END_REF] measured pressures on sails for upwind and downwind sails. However downwind sails are more complex to study than upwind sails mainly due to their non-developable 3D shape with highly cambered sections and massively detached flow around a thin and very flexible membrane. Due to the dynamic behavior of this unsteady fluid-structure interaction, the pressures on the sail vary quickly. Even in stable conditions, offwind sails have an inherent unsteadiness. One key feature of spinnaker unsteadiness comes from the flapping at the leading edge, also called luffing. We have previously investigated pressure evolution during luffing [START_REF] Deparday | measurements of pressures, sail shape and forces on a full-scale spinnaker, 23rd HISWA Symposium on Yacht Design and Yacht Construction[END_REF][START_REF] Motta | An investigation of the dynamic behaviour of asymmetric spinnakers at full-scale[END_REF]. In [START_REF] Deparday | measurements of pressures, sail shape and forces on a full-scale spinnaker, 23rd HISWA Symposium on Yacht Design and Yacht Construction[END_REF], we showed an example where flapping of spinnaker creates pressure peaks at the leading edge increasing the aerodynamic force dynamically by 50%. Due to the non-stationarity of the environment while sailing, spatio-temporal pressure data are complex to analyze and therefore to simulate. However significant and different spatio-temporal patterns can be spotted [START_REF] Motta | An investigation of the dynamic behaviour of asymmetric spinnakers at full-scale[END_REF] and might be produced by different physical causes (Fluid-Structure Interaction, wind variations, boat motions, etc.). In this paper we present an approach to decompose complex pressure evolutions into simpler modes. It would then allow easier analysis and comparison with simulations. This paper presents results of full-scale experiments of an instrumented J/80 class yacht in offwind conditions were loads, pressures on the spinnaker, boat and wind data were measured. After presenting the experimental apparatus, average and fluctuations of pressures and loads are presented. The next section is the use of the Proper Orthogonal Decomposition (POD) method on pressures to create a simpler model of the complex variations of pressure distribution in time. We show then that the method also helps to highlight the correlation between the main evolution of pressures and the variations of loads.

EXPERIMENTAL SETUP

An instrumented J/80 class sailing yacht, an 8 meter onedesign cruiser racer was used during those experiments. A tri-radial asymmetrical spinnaker with a surface of about 65 m 2 with a 12 meter long rounded luff was hoisted as well as a mainsail of 17 m 2 . Boat and wind data, loads on the standing rigging and on the sails were recorded. Moreover pressure taps, developed by the Yacht Research Unit from the University of Auckland were stuck on the spinnaker to acquire the dynamic pressure distribution. They were synchronized with the other data thanks to an acquisition software, RTMaps developed by Intempora which received every signal at their own rate and timestamped them "on the flow". A resampling was applied during the post processing to obtain synchronous data for easier analysis. Figure 1 shows the arrangement of all the sensors set onboard. This setup for downwind navigation is a further development of the experimental system described in [START_REF] Augier | Experimental validation of unsteady models for fluid structure interaction: Application to yacht sails and rigs[END_REF] which was used for measurements in upwind navigation.

Loads

The standing rigging (shrouds, forestay and backstay) is fitted with custom-made turnbuckles and shackles equipped with strain gages. The running rigging (the corners of the mainsail and of the spinnaker -head, tack and clew-) is equipped with instrumented shackles too. For the standing rigging and the mainsail the sensors are connected to a load acquisition system Spider8 from HBM. Voltages are received from all strain gages and amplified. They are then converted in digital data at a rate of 25 Hz. Thereafter they are transferred to the real-time acquisition software, RTMaps. Due to the high displacements of the spinnaker -in the order of magnitude of 1 to 5 meters for the spinnaker used in those experiments-, the instrumented shackles on the three corners of the spinnaker communicate wirelessly to the acquisition system. The clew sensor is connected via a wire running along the foot of the sail to a small box located near the tack point of the spinnaker (see Figure 2). This box contains two strain gage amplifiers, one for the tack sensor and one for the clew. A microcontroller receives and transmits data at a sampling frequency of 25 Hz to the receiver inside the boat via a wireless and low consumption ZigBee network. Another box is located at the head position for the head instrumented shackle. The delay between the emission and reception of data is insignificant compared to dynamics in sailing. The errors of measurement are less than 2% of the measurement range (10 000 N for the shrouds, forestay and for the mainsail sheet, 5000 N for the backstay and other instrumented shackles on the mainsail and spinnaker).

Pressures on spinnaker

On the spinnaker, 44 low range differential pressure sensors (Honeywell XSCL04DC) are located on the surface along 4 horizontal stripes: at 1/4, 1/2, 3/4 and 7/8th height of the spinnaker (see Figures 1 and2). 12 transducers are used on each of the first 3 stripes and 8 for the top one. There is a higher concentration of pressure taps near the leading edge to be able to record potential leading edge suction peaks. Those sensors measure a difference of pressure between the suction side and the pressure side using the piezoresistive effect. There is no need for a measure of a reference pressure, a complicated task in full-scale experiments. Those sensors are stuck on one side -at the pressure side when sailing on portside tack-and are positioned facing 2mm-diameter holes on the sail to measure the pressure jump across the sail without significant air leak. Punctured light sail cloth patches are applied on the pressure taps to profile the sensors. This custom built pressure system was designed by the Yacht Research Unit at the University of Auckland. These pressure transducers are connected by wires to the receiver inside the boat and thus are synchronized with the other data. The pressure sensors have a sampling frequency of approximatively 10 Hz with a maximum range of ±1 kPa and a resolution of 0.5 Pa. The pressure acquisition system is more described in [START_REF] Motta | Experimental investigation of asymmetric spinnaker aerodynamics using pressure and sail shape measurements[END_REF].

Procedure

Sea trials were performed in the bay of Brest, France, offshore Ecole Navale. During those experiments the weather conditions were stable:

• average true wind speed: 6 m/s (12 kn) • gust: 8 m/s (16 kn) • wind direction: 270 • (westerly wind). Stable. Flat water.

Even in conditions considered as "stable" (with no gust, no wind shift, on flat water and fixed trimming), offwind sails have an inherent unsteadiness, like luffing (flapping at the leading edge). To keep "stable" conditions, a standard procedure must be followed. During those experiments, the controlled inputs were the apparent wind angle and the trim of the spinnaker. Trimmer and helmsman were kept the same for the whole test. The trimmer adjusted the spinnaker at the optimum trim (i.e. on the verge of luffing at the leading edge). The helmsman kept the apparent wind angle as constant as possible.

During the post-processing routine, periods of 5 seconds minimum were labelled "stable" when the standard deviation of the apparent wind angle (AWA) was below 4 • and the standard deviation of the apparent wind speed (AWS) was below 10% of the average. Those periods were extended in time as long as those criteria were met. A large range of AWA (between 55 • and 140 • ) is swept by the "stable" periods found, with a certain redundancy for most of the AWA. Each stable period is processed individually.

AVERAGES AND FLUCTUATIONS Pressures

The average pressure distribution and loads are compared according to the apparent wind angle. A large range of apparent wind angles (AWA) has been met in a rather constant true wind speed (TWS) between 5.8 m/s and 7.1 m/s thus between 11.2 kn to 13.8 kn. The apparent wind angle is measured at the mast head. One should be aware this measure is affected by the twist of the wind, the upwash effect from the sails and the heel.

To display the pressure distribution on the whole sail from discrete measurement points, a linear Radial Basis Function interpolation has been used. On following figures where the pressure distribution is displayed on the spinnaker, blue crosses show where the pressure measurement sensors were located on the sail. Thus pressures at those blue crosses are actual measured values, when the pressure distribution is interpolated between the stripes and pressure taps. With no information on the sail boundaries, values at the top (above 7/8th) and at the bottom (below 1/4th) are extrapolated. The shapes used to display the pressure distributions come from other experiments with the same spinnaker where photogrammetric measurements were carried out to acquire the flying shapes. Time-Averaged pressure distributions for similar apparent wind angles have good repeatability. Moreover the pressure distribution evolves clearly with the AWA. Figure 3a presents 3 characteristic pressure spatial distributions at 66 • , 118 • and 140 • . It shows the coefficient of the difference of pressure as commonly defined in aerodynamics:

∆C P = P leeward -P windward 1 2 ρ(AWS) 2 .
At tight angles as AWA 66 • , a bulb of high suction is found at the leading edge in the top half of the spinnaker (∆C P ≈ -3) which produces high aerodynamic force. ∆C P on the rest of the sail is around -2 increasing to -1/-0.5 at the trailing edge.

At AWA around 110 • -120 • , the area where the peak of suction occurs is smaller around half of the spinnaker height and the absolute value lower. On the rest of the spinnaker, the pressure coefficient on the spinnaker is rather constant around -1.2 and increasing to -0.5 at the trailing edge. At AWA 140 • , the decrease of suction is even more visible on the whole sail with almost no suction peak at the leading edge and with a reduction of |∆C P | along the flow up to a positive ∆C P at the trailing edge even on the actual measured points. Positive pressure coefficient means a collapse of the sail at the trailing edge and thus an unstable flying shape. It is consistent with what the authors have noticed during experiments: at large AWA, the spinnaker starts collapsing first at the leech and not at the luff. While the AWA is increased, not only is a clear decrease of absolute differential pressure coefficient (from -3 down to 0 about), but also the AWS decreases (from 7 m/s to 3.5 m/s about). So the absolute values of ∆P decrease even more dramatically: At tight AWA, around 65 • , the order of magnitude of differential pressure is -40 Pa, and only -4 Pa at large AWA -around 140 • -.

Figure 3b shows the standard deviation for the corresponding AWA. Standard deviation on the whole sail is interpolated from the standard deviations calculated on the pressure taps only. Higher standard deviations mean bigger variations of pressure during a "stable" period. Despite a clear difference for the pressure distribution on the whole spinnaker depending on the AWA, pressure variations during the "stable" periods have similar spatial patterns. Strong variations are found at the leading edge, around 1, on the whole height for 66 • and 118 • while the rest of the spinnaker has a standard deviation of about 0.2. However, while the order of magnitude of standard deviation of ∆C P is similar for every AWA, the relative variation of pressures compared with the average pressure coefficient varies. Variations are more significant for large AWA (around 120 • -140 • ) than for tight AWA. For tight AWA, the standard deviation is around 30 Pa thus 75% of the average pressure. For large AWA, the standard deviation is around 8 Pa thus 2 times bigger than the average pressure.

Loads

Figure 4 displays the load coefficients on the three corners of the spinnaker according to the apparent wind angle:

C F = Load
1 2 ρ S(AWS) 2 with S the sail area of the spinnaker. In Figure 4, only "stable" periods of 10 seconds minimum are taken. Even though the periods chosen are "stable", loads can vary significantly. Therefore each period for a specific average apparent wind angle is displayed as a box plot. The central red mark is the median, and the edges of the box are the lower and upper quartile. The lower quartile (Q1) splits off the lowest 25% of data from the highest 75%. The upper quartile (Q3) splits off the highest 25% of loads from the lowest 75%. The box represents the interquartile range (IQR = Q3 -Q1). It contains 50% of the loads recorded during one "stable" period. The whiskers show the maximum and minimum loads recorded.

Figure 4 shows also the general trend of the load coefficients on the three corners according to the apparent wind angle. Head and tack have similar evolution with a decrease especially between 110 • and 140 • respectively from 0.8 to 0.5 and from 0.7 to 0.3. Whereas clew load coefficient is approximatively constant around 0.4. As explained previously, when the AWA is increased, the AWS decreases. While at clew point, the absolute loads mostly decrease only due to the decrease of the AWS, at tack and head points the absolute loads decrease even more significantly. To confirm this, Figure 5 displays the evolution of load coefficient using the True Wind Speed (TWS) for the non-dimensional coefficient:

C F = Load
T W S = (AWS cos(AWA) -BS) 2 + (AWS sin(AWA)cos(heel))
2 TWS is rather constant for every AWA (around 12 kn). A higher decrease is noticed when the AWA is increased for the head and tack loads than for the clew load which varies only a little.

To analyze variations of load coefficients from an aerodynamic point of view the AWS is used for the non dimensional coefficient C F . Q3 can be seen as an arbitrary separation between the small variations of loads around the median (inside the IQR) and the peaks of loads (in the top quarter). In Figure 4, the IQR has the same relative range for every AWA, between -10% and +10% of the median for the most loaded corners (head and tack) and between -5% and +5% for the clew. Since the IQR is relatively constant and small, one can conclude that for every AWA, without taking into account peaks of loads, the averaged load coefficient is rather stable and varies slightly.

For the most loaded corners (head and tack) the upper whisker (maximum load) is about 20% higher than the median for tight angles (AWA < 100 • ) and 30% for large angles (AWA > 120 • ). For the clew point the upper whisker is always about 15% higher than the median for all AWA. Relative variations of loads seem fairly constant for every AWA, and slightly bigger for large AWA at the head and tack points. It is interesting to note that most variations of loads are present at the head and tack points, the closest points of the leading edge where the highest variations of pressure occur. While at clew point, the relative variation of loads is smaller.

Those variations of loads and pressures are unsteady even during "stable" periods. However specific patterns might be spotted and might be linked to different causes. The yacht motion and its influence on the apparent wind (pitching and rolling of the boat), gusts (pure aerodynamic cause), vortex shedding, or a change of the spinnaker shape as luffing (unsteady fluid-structure interaction) could make the spinnaker forces vary. Therefore, we would like to extract patterns in order to decompose complex pressure evolutions into simpler modes. Those pressure modes could help to describe a temporal global behavior in a better way than analyzing each pressure sensor signal, and could be correlated with other recorded data. We decided to use the Proper Orthogonal Decomposition method to characterize the spatial pattern of pressure variations.

DECOMPOSITION INTO MODES Proper Orthogonal Decomposition method

The Proper Orthogonal Decomposition (POD), is based on the Karhunen-Loeve expansion and also called Principal Component Analysis, PCA. It was first introduced in the context of Fluid Mechanics by [START_REF] Lumley | The structure of inhomogeneous turbulence[END_REF]. The input data (in our case ∆C P (x, t) ) can be expanded into orthogonal basis functions φ i (x) with time coefficient a n (t):

U (x, t) = n a n (t) • φ n (x).
As proper modes are derived from the data itself (data driven decomposition), there is no need of a-priori knowledge or education scheme. Moreover, each basis function has its own amount of fluctuation energy different from each other. These functions are statistically optimal in the least mean-square sense. As a result, fluctuation energy drops down quickly which means a low number of modes is needed in the expansion to reproduce the main variations of the field. POD is a powerful tool for generating lower dimensional models of dynamical systems.

Most of the time, POD is used on the fluctuations of the input data only. After subtracting the average component (seen as the zeroth mode) from the data, a matrix U is created as a set of N observations (commonly called snapshots) of M records. Each column contains all fluctuating input data (M values) from a specific snapshot and each row contains all snapshots (N snapshots) from a specific measurement point.

U =      u 11 u 12 • • • u 1N u 21 u 22 • • • u 2N . . . . . . . . . . . . u M 1 u M 2 • • • u M N     
Then the auto covariance matrix C (MxM) is calculated as:

C = U * U T
because M >> N . We have M= 44 measurement points and N ≈ 20000. However in fluid mechanics, it is common to have N >> M when using PIV or CFD results for example. For those cases the so-called "Snapshot POD" introduced first by [START_REF] Sirovich | Turbulence and the dynamics of coherent structures part i: coherent structures[END_REF] is used. For our experiments, the "Direct POD" has been applied. The corresponding eigenvalue problem of the auto covariance matrix is solved:

C * Φ = λ * Φ
The eigenvectors Φ(i) are the POD modes. POD modes are sorted in descending order according to their corresponding eigenvalue λ(i) which represent their energy. The POD mode with the highest corresponding eigenvalue is mode 1. The expansion coefficient (or mode time coefficient) is calculated as follows:

a = U T * Φ

POD results

Following results presented here are for a "stable" period with an average AWA of 69 • , but is representative to what we have observed for different periods at different AWA. This point will be discussed further in the article. Figure 6 shows the energy distribution for each POD mode. The first mode contains almost 45% of the fluctuation energy. Mode 2 and 3 represent only 15% each. And most of the time other modes have less than 5% of the fluctuation energy.

It is clear that the first mode is dominant compared to the others.

The pressure distribution evolution can be simplified by taking only the first terms of the expansion. The reconstruction using the first modes allows to only keeping the most energetic part of the signal and removing insignificant variations and noises. The precision of the reconstruction has been calculated according to the number of modes. With mode 0 (the average) and mode 1, 85% of the signal is already reconstructed. With 3 modes, the error of reconstruction of the pressure signals is 10%. About 10 modes are required to achieve a reconstruction with less than 5% of difference.

Figure 7 presents the first 3 modes. The scale value is arbitrary. To represent a fluctuation of Cp, they must be multiplied by the corresponding mode time coefficient depending on the time -which can be positive or negative -(presented in Figure 9). Mode 1 has a bulb of pressure on the top half of the spinnaker at the leading edge and a smaller bulb of opposite sign on the bottom half height of the spinnaker. Mode 2 is similar with the standard deviation pattern presented before. Mode 3 and further modes display smaller coherent patterns and may change with the period used. Table 1 shows the energy distribution for different "stable" periods at different AWA. The ratio of energy of each mode number is rather constant for every AWA. Moreover each mode number has a similar pattern of pressure distribution, even though it happens that mode 2 and mode 3 are inverted in a few cases. To compare modes, the maximum value of time coef- ficient is taken to be multiplied by the spatial mode: max(a 1 (t)) • φ 1 (x). Maximum value of time coefficient is used since we want to analyze and compare dominant variations. Moreover, by mathematical definition: ∀n, a n (t) = 0 and thus is irrelevant to be used for comparison. Maximum values of mode 1 for different AWA are presented in Figure 8 at 4 different stripes where the pressures are measured (1/4, 1/2, 3/4, 7/8th height of the spinnaker). Differential pressure coefficients have comparable shapes for every AWA, except that |∆C P | on the bottom half is slightly smaller for deeper AWA. The bulb of suction at the leading edge at 7/8 and 3/4 of the spinnaker height is always present and a smaller bulb of positive ∆C P at 1/2 and 1/4 height is also spotted. Even if the POD method is a data driven decomposition (i.e. modes are derived from the data itself), there is a good repeatability of POD modes for "stable" periods when the spinnaker has a fixed trim. Moreover mode 1, which plays an important role in the fluctuation of pressures, could be defined as a unique mode whatever the AWA. POD modes evolve in time. When the time coefficient of a corresponding mode is at an extremum, the corresponding mode is then preponderant. Analyzing time coefficients would then help to link pressure variations with other recorded data. Figure 9 shows the evolution of the time coefficient for the first 3 modes. Amplitudes of mode 1 are bigger than the other modes as expected due to the larger energy it possesses. Mode 1 and mode 2 are slightly correlated; mode 2 is shifted of a quarter of a pseudo-period. It means when mode 1 is maximal, mode 2 is null. A typical pseudo period for mode 1 stands out for this AWA 69 • . Furthermore, for different "stable" periods, at different AWA -not displayed here-, similar variations of the temporal coefficient of mode 1 are detected. However the dynamics change with the AWA. The pseudo-period is measured, and the corresponding frequency is displayed in Table 2 with the corresponding average AWS of the "stable" period. The reduced frequency is calculated as follows:

f r = f s • √ S AW S
with S the sail area, thus √ S = 8.3 m. When the AWA is increased, the typical pseudo frequency of the time coefficient of mode 1 is reduced by a factor of 2.5. The AWS decreases at a similar rate. A ratio of 2.1 is found between AWA 57 • and 140 • . Therefore the reduced frequency is nearly constant with a small decrease. Figure 10 displays the pseudo-frequency of the time coefficient of mode 1 as a function of the AWS. There is a linear dependence of the pseudo-frequencies with the AWS. It demonstrates that mode 1 is mostly driven by aerodynamic phenomena as expected, and not by mechanical resonance of the rigging or of the membrane of the sail. In conclusion, the spatial pattern of mode 1 does not change with the AWA, but only the temporal dynamics with an increase of the pseudo-period when the AWA increases. Table 3 presents the cross-correlation of all data measured for the specific period presented in [START_REF] Deparday | measurements of pressures, sail shape and forces on a full-scale spinnaker, 23rd HISWA Symposium on Yacht Design and Yacht Construction[END_REF] where loads and flapping were strongly correlated. The normalized cross-correlation is calculated with the time coefficients of the first three modes. Cross correlation between two signals X and Y is defined as follow:

C xy (τ ) = E[(X(t 2 ) -µ X(t2) ) -(Y (t 1 ) -µ Y (t1) )]
where E[•] is the expected value operator, τ = t 2 -t 1 is the shift applied between two signals. µ X and µ U Y are the mean functions. The cross correlation matrix is calculated to determine the correlations of every signal with each other. The values are between 0 when not correlated at all -in white in the table-and 1 when signals have the same dynamics -in red in the table. Colors in Table 3 highlight the correlations between experimental data. The diagonal represents the auto-correlation of every signal. A very strong correlation is present between the loads except with the forestay and the shroud D1 leeward as they are not loaded and slack. Loads are also slightly correlated with the yaw. The shift between the two signals is about 1 second in advance for the loads. A peak of aerodynamic loads might modify the aerodynamic center of effort and thus change the equilibrium of the sailing yacht and make the course vary. The time coefficient of mode 1 is also well correlated with the loads (around 0.9). There is no delay between mode 1 and the spinnaker aerodynamic loads. Spinnaker aerodynamic loads are 0.1 s in advance of the standing rigging loads. In this case the pressure evolution is instantaneously transmitted to the corners of the spinnaker which then transmit this increase of loads to the shrouds and backstay. Mode 1 is a good parameter to define peaks of loads due to flapping of the spinnaker. Here mode 2 is reasonably correlated with mode 1 with a coefficient of 0.72. The shift between the two modes is 0.4 s. As explained previously with Figure 9, mode 2 is shifted of a quarter of a pseudo-period. Thus the pseudo-frequency for this specific period should be around T = 4 • 0.4 s = 1.6 s. This pseudo-period corresponds to what was presented in [START_REF] Deparday | measurements of pressures, sail shape and forces on a full-scale spinnaker, 23rd HISWA Symposium on Yacht Design and Yacht Construction[END_REF]. In conclusion, mode 1 describes well the flapping of the luff. It can be represented as a unique spatial pattern for every AWA and thus is a good indicator of the dynamics of the spinnaker. Its temporal pseudo-period slows down as the AWS decreases when the AWA is increased. Head 0.59 0.96 0.94 0.95 0.69 1.00 1.00 0.98 0.90 0.63 0.47 0.58 0.57 0.71 0.63 0.53 0.57 Tack 0.59 0.95 0.94 0.95 0.70 1.00 1.00 0.98 0.88 0.61 0.45 0.54 0.55 0.71 0.62 0.51 0.58 Clew 0.59 0.93 0.94 0.95 0.76 0.98 0.98 1.00 0.85 0.65 0.44 0.63 0.58 0.70 0.67 0.61 0.60 mode 1 0.65 0.91 0.85 0.88 0.66 0.90 0.88 0.85 1.00 0.72 0.43 0.52 0.63 0.74 0.54 0.47 0.38 mode 2 0.62 0.66 0.66 0.68 0.53 0.63 0.61 0.65 0.72 1.00 0.46 0.64 0.84 0.60 0.68 0.48 0.29 mode 3 0.59 0.45 0.41 0.41 0.52 0.47 0.45 0.44 0.43 0.46 1.00 0.66 0.41 0.45 0.56 0.75 0.51 roll 0.60 0.58 0.48 0.54 0.49 0.58 0.54 0.63 0.52 0.64 0.66 1.00 0.83 0.83 0.77 0.88 0.47 pitch 0.69 0.61 0.61 0.61 0.52 0.57 0.55 0.58 0.63 0.84 0.41 0.83 1.00 0.76 0.85 0.68 0.28 yaw 0.51 0.75 0.68 0.71 0.65 0.71 0.71 0.70 0.74 0.60 0.45 0.83 0.76 1.00 0.66 0.56 0.46 AWA 0.73 0.58 0.63 0.63 0.60 0.63 0.62 0.67 0.54 0.68 0.56 0.77 0.85 0.66 1.00 0.58 0.52 AWS 0.59 0.49 0.50 0.54 0.54 0.53 0.51 0.61 0.47 0.48 0.75 0.88 0.68 0.56 0.58 1.00 0.50 BS 0.40 0.57 0.44 0.46 0.44 0.57 0.58 0.60 0.38 0.29 0.51 0.47 0.28 0.46 0.52 0.50 1.00 CONCLUSIONS Full-Scale experiments were carried out on an instrumented J/80 sailing yacht where loads, pressure distribution on the spinnaker boat and wind data were measured at different downwind angles. We show pressure coefficients and load coefficients decrease when the AWA is increased, whereas variations of load and pressure coefficients are mainly constant. Therefore variations relative to the average loads or pressures are bigger for larger AWA. Moreover we found that most of the pressure and load variations are mainly at the luff for every AWA. A POD analysis has been used on pressure signals in order to identify the most energetic patterns. The first mode of the POD method is uniquely identified for every AWA, and is well associated with the flapping of the luff producing correlated variations of loads on the 3 corners of the spinnaker. It has a unique spatial mode of pressure fluctuations for all AWA. Only the temporal behavior differs with the AWA. An identified typical pseudo-period of flapping has been calculated thanks to this first mode. It shows a linear decrease of the pseudo-frequency with the AWS. Flapping of the luff might be the result of an aerodynamic phenomenon. Moreover POD enables to characterize a global unsteady behavior instead of analyzing all local pressure time series This paper presented a way to characterize the pressure evolution due to flapping. It permits to show to sail designers where the highest variations of pressure occur when flapping. Moreover thanks to the POD method, the first modes allow to reconstruct a signal with the main variations (i.e. with the most energetic part) and remove noises. Therefore comparison with other measurements or numerical simulations is simplified.

Figure 1 :

 1 Figure 1: General arrangement of the experimental set-up on the J/80. 16 load sensors (green discs), 44 pressure taps (red circles), and wind and boat sensors (blue squares).

Figure 2 :

 2 Figure 2: Photograph of the tack of the spinnaker used.

Figure 3 :

 3 Figure 3: Pressure distributions, time average pressure coefficient ∆C P (a) and the fluctuations (b) for 3 typical AWA (66 • , 118 • , 140 • ). Blue crosses show the positions of the pressure taps.

1 2 ρ

 2 S(TWS) 2 for Figure5only. with the TWS formula from[START_REF] Fossati | Aero-hydrodynamics and the performance of sailing yachts[END_REF]:

Figure 4 :

 4 Figure 4: Boxplot for load coefficients C F at the three corners of the spinnaker for different AWA. Central red mark is the median, the box represent the interquartile range between the lower and upper quartiles. It contains 50% of the loads. The upper and lower whiskers indicate the minimum and maximum values.

Figure 5 :

 5 Figure 5: Boxplot for load coefficients C F with the TWS as reference at the three corners of the spinnaker for different AWA. Central red mark is the median, the box represent the interquartile range between the lower and upper quartiles. It contains 50% of the loads. The upper and lower whiskers indicate the minimum and maximum values.

Figure 6 :

 6 Figure 6: Energy distribution of the fluctuations for the first 10 modes for AWA 69 • .

Figure 7 :

 7 Figure 7: First 3 spatial orthogonal modes from POD method for AWA 69 • . Mode 1 is the most energetic mode.

Figure 8 :Figure 9 :

 89 Figure 8: Mode 1 for maximum time coefficient max(a 1 (t)• φ 1 (x) for different AWA.

  Figure 10: Pseudo-frequency of the time coefficient of the first POD mode (red dots) for different AWA as a function of the AWS. A linear interpolations fits the experimental data.

Table 1 :

 1 Energy distribution of the fluctuations for the first 10 modes for different AWA.

	Modes 66 deg 98 deg 120 deg 140 deg
	1	43%	46%	45%	45%
	2	21%	16%	19%	15%
	3	8%	9%	10%	11%
	4	7%	6%	6%	7%
	5	4%	5%	3%	5%
	6	3%	4%	3%	3%
	7	2%	3%	2%	2%
	8	2%	2%	2%	2%
	9	2%	2%	2%	2%
	10	1%	1%	1%	1%

Table 2 :

 2 Pseudo-frequency and reduced frequency for different AWA.

	AWA( • ) f s (Hz) AWS(m/s) f r (-)
	57	1.5	7.0 1.78
	66	1.3	6.4 1.69
	69	1.3	6.2 1.73
	98	1	5.4 1.53
	118	0.7	3.7 1.57
	140	0.6	3.3 1.52

Table 3 :

 3 Cross-correlation between different signals recorded during experiments, and first 3 modes of the POD. Correlation values (and their corresponding colors) vary between 0 (in white) meaning no correlation and 1 (in red), same dynamics. .62 0.88 1.00 0.99 0.79 0.94 0.94 0.94 0.85 0.66 0.41 0.48 0.61 0.68 0.63 0.50 0.44 D1 windw'd 0.63 0.92 0.99 1.00 0.79 0.95 0.95 0.95 0.88 0.68 0.41 0.54 0.61 0.71 0.63 0.54 0.46 D1 leew'd 0.78 0.62 0.79 0.79 1.00 0.69 0.70 0.76 0.66 0.53 0.52 0.49 0.52 0.65 0.60 0.54 0.44

	forestay	backstay	V1 windw'd	D1 windw'd	D1 leew'd	Head	Tack	Clew	mode 1	mode 2	3 mode	roll	pitch	yaw	AWA	AWS	BS
	forestay 1.00 0.58 0.62 0.63 0.78 0.59 0.59 0.59 0.65 0.62 0.59 0.60 0.69 0.51 0.73 0.59 0.40
	backstay 0.58 1.00 0.88 0.92 0.62 0.96 0.95 0.93 0.91 0.66 0.45 0.58 0.61 0.75 0.58 0.49 0.57
	V1 windw'd 0

ACKNOWLEDGEMENTS

This project has received funding from the European Union's Seventh Programme for research, technological development and demonstration under grant agreement No PIRSES-GA-2012-318924, and from the Royal Society of New Zealand for the UK-France-NZ collaboration project SAILING FLUIDS.


MD, PhD Nicolas Bonnevialle
email: nicolasbonnevialle@yahoo.fr

MD, PhD Philipe Clavert

MD a , Marine Arboucalot

MD David Bahlau

MD, PhD e , Thomas Bauer

MD, PhD f , SOFCOT Matthieu Ehlinger

Contribution of arthroscopy in the treatment of anterior glenoid rim fractures: a comparison with open surgery

published or not. The documents may come

Contribution of arthroscopy in the treatment of anterior glenoid rim fractures: a comparison with open surgery

Level of evidence: Level III; Retrospective Cohort Design; Treatment Study Keywords: Glenoid; fracture; shoulder; Ideberg; arthroscopy; shoulder instability Scapula fractures with extension to the glenoid are rare. In the Ideberg classification of intra-articular fractures, type IA corresponds to an anterior glenoid rim fracture, which occurs with glenohumeral dislocation in two-thirds of cases. [START_REF] Ideberg | Epidemiology of scapular fractures. Incidence and classification of 338 fractures[END_REF] When the articular fragment makes up more than 20% of the glenoid surface, reduction and stabilization may be proposed to prevent glenohumeral instability and osteoarthritis (OA) from developing over time. [START_REF] Goss | Fractures of the glenoid cavity[END_REF][START_REF] Tauber | Arthroscopic screw fixation of large anterior glenoid fractures[END_REF] Arthroscopy techniques offer less invasive and, theoretically, more precise solutions than standard open surgery through the anterior deltopectoral approach for treating these intra-articular fractures. [START_REF] Bauer | Arthroscopic treatment of glenoid fractures[END_REF][START_REF] Jiang | Do reduction and healing of the bony fragment really matter in arthroscopic bony Bankart reconstruction?: a prospective study with clinical and computed tomography evaluations[END_REF][START_REF] Kavanagh | Open reduction and internal fixation of displaced intra-articular fractures of the glenoid fossa[END_REF][START_REF] Mayo | Displaced fractures of the glenoid fossa. Results of open reduction and internal fixation[END_REF][START_REF] Nakagawa | Bone fragment union and remodeling after arthroscopic bony Bankart repair for traumatic anterior shoulder instability with a glenoid defect: influence on postoperative recurrence of instability[END_REF][START_REF] Plath | Arthroscopic suture anchor fixation of bony Bankart lesions: clinical outcome, magnetic resonance imaging results, and return to sports[END_REF][START_REF] Schandelmaier | Fractures of the glenoid treated by operation. A 5-to 23-year follow-up of 22 cases[END_REF][START_REF] Scheibel | Arthroscopic reduction and fixation of large solitary and multifragmented anterior glenoid rim fractures[END_REF][START_REF] Scheibel | Open reconstruction of anterior glenoid rim fractures[END_REF][START_REF] Sugaya | Arthroscopic osseous Bankart repair for chronic recurrent traumatic anterior glenohumeral instability[END_REF][START_REF] Tauber | Arthroscopic screw fixation of large anterior glenoid fractures[END_REF][START_REF] Yamamoto | Osseous defects seen in patients with anterior shoulder instability[END_REF] Arthroscopy also provides an opportunity to assess and treat any associated rotator cuff tendon injuries or greater tuberosity fracture. [START_REF] Rouleau | Surgical treatment of displaced greater tuberosity fractures of the humerus[END_REF][START_REF] Tauber | Arthroscopic screw fixation of large anterior glenoid fractures[END_REF] There is currently no evidence that one surgical technique is superior for treating Ideberg IA fractures. This study assessed the clinical and radiologic outcomes of Ideberg type IA glenoid fractures treated using conventional open surgery compared with those treated with arthroscopic surgery. We hypothesized that arthroscopic surgery would reduce the complication rate and result in better clinical outcomes than open surgical treatment. During the study period, 74 glenoid fractures were treated surgically, and 56 were Ideberg IA fractures. One patient was lost to follow-up before the 12-month follow-up visit. Thus, 10 patients in group O and 45 in group A were available for analysis.

Materials and methods

This

Surgical technique

The technical choice of open or arthroscopic surgery depended on surgeon preference.

Patients in group O were operated on in the beach chair position under general anesthesia with an interscalene nerve block. A deltopectoral approach was used. Subscapular tenotomy was performed to access the glenohumeral joint. After exploring the joint and locating any osteochondral lesions, the surgeon reduced the fragment and held it in place temporarily with Kirschner wires. The fragment was secured with one or two 3.5-mm diameter cannulated screws (Fig. 1). If a proximal humeral fracture was also present, additional fixation with screws or a plate was done through the same approach.

Patients in group A were operated on in the beach chair position (n = 33) or in lateral decubitus (n = 13). One or 2 instrument portals (anterosuperior, anteroinferior) were established along with an optical portal in the posterior soft point. After draining the hemarthrosis, the joint was explored to identify any osteochondral lesions. The fracture was reduced and stabilized with anchors according to the Sugaya [START_REF] Sugaya | Arthroscopic osseous Bankart repair for chronic recurrent traumatic anterior glenohumeral instability[END_REF] technique (n = 36) or with a combination of anchors and screw fixation (n = 10) introduced via an accessory lower anterior portal through the subscapularis (Fig. 2).

The patient's operated-on arm was immobilized postoperatively with a sling in internal rotation for 4 to 6 weeks, no matter which surgical technique was used. Pendulum exercises were initiated within a few days of surgery, and active movements were allowed starting at 6 weeks, under the control of a physiotherapist.

Clinical evaluation

Patients were reviewed for clinical and radiologic follow-up assessments and to provide consent to participate in this study. Active range of motion was measured with a goniometer with the patient seated and performing anterior elevation, external rotation with the elbow at the side, and internal rotation (maximal vertebral level reached by the patient's thumb) movements. The outcomes were assessed objectively with the Constant Score and subjectively with an analog scale to determine satisfaction (excellent, good, average, or poor outcome). [START_REF] Constant | A clinical method of functional assessment of the shoulder[END_REF] Complications and reoperations were noted.

Radiographic evaluation

Preoperatively, plain anteroposterior x-ray images and a computed tomography scan were used to classify the fracture based on the Ideberg classification. [START_REF] Ideberg | Epidemiology of scapular fractures. Incidence and classification of 338 fractures[END_REF] The largest width of the fracture was measured on the computed tomography scan. Fracture healing was assessed at 6 months postoperatively on anteroposterior and lateral views. Signs of osteoarthritis were assessed at the latest follow-up according to Samilson and Prieto. [START_REF] Samilson | Dislocation arthropathy of the shoulder[END_REF]

Statistical methods

Measurements are expressed as the mean, range, and standard deviation. The D'Agostino-Pearson test was used to analyze the data distribution. Paired values were compared using a paired t test, and unpaired results were compared using the Mann-Whitney test. The χ 2 or Fisher exact test was used to compare categoric data. The significance level was set at P < .05.

Results

Study cohort

Patients were a mean age of 46 years (range, 20-77 years), and there were 13 women and 42 men. The fracture was caused by a sports-related injury in 38% (n = 17). Glenohumeral dislocation was associated in 87% (n = 48). Ten patients also had a proximal humeral fracture: undisplaced greater tuberosity (n = 7), displaced greater tuberosity (n = 2), and 3-part fracture (n = 1). Characteristics of the patients in groups O and A are summarized in Table I. Group O had more patients with (8-18) .29 Time to surgery, d 2.5 ± 3 (0-10) 5 ± 5 (1-22) .01

Data are presented as mean ± standard deviation (range), as number of patients, or as number (%).

an associated humeral fracture, and they were operated on more quickly.

Surgical data

The

Clinical outcomes

At a mean follow-up of 30 months (range, 12-115 months), the mean Constant score was 82 ± 16 points (range, 20-100 points) in the entire cohort. Active anterior elevation was 162°± 30°(range, 45°-180°), external rotation was 51°± 19°(range, 10°-80°), and internal rotation was 8 ± 2 points (range, 2-10 points). The results for the 2 groups are compared in Tables II andIII. Patients in group A had better strength, anterior elevation, and external rotation than those in group O. There were no significant differences between the 2 groups in the absolute and adjusted Constant Score.

The presence of osteochondral glenoid lesions detected intraoperatively had a negative effect on the overall Constant Score at the final follow-up (76 points vs. 88 points; P = .002). When humeral fractures were excluded from the analysis to make the groups more even, no significant difference was found in the average absolute and adjusted Constant Score between group O (82 patients [120%], P = .78) and group A (84 patients [27%], P = .36). The subjective outcome was good or excellent in 90% of group O patients and in 91% of group A patients (P = .19).

Radiologic outcomes

At the last follow-up visit, all fractures had healed, regardless of the type of surgery performed. The incidence of OA is reported in Table IV. There were no significant differences between groups.

Complications

There were 5 complications (9%) in the cohort: 3 in group O (30%) and 2 in group A (4%; P = .03). In group O, there were 2 infections and 1 case of joint stiffness requiring 3 reoperations: 2 lavage/débridement procedures with removal of hardware and 1 arthrolysis procedure with a poor result at 12 months. In group A, 2 cases of joint stiffness required 2 arthrolysis procedures with removal of hardware. These 2 patients recovered full range of motion uneventfully.

Discussion

Our initial hypothesis was partially validated: arthroscopic and conventional open surgical treatment of Ideberg IA glenoid fractures results in similar clinical outcomes at the early to midterm follow-up. However, the postoperative complication rate was significantly higher with open surgical treatment.

No matter which type of surgical technique was used, the outcome was negatively affected by the presence of posttraumatic osteochondral lesions. Kavanagh et al [START_REF] Kavanagh | Open reduction and internal fixation of displaced intra-articular fractures of the glenoid fossa[END_REF] reported the outcomes of a case series of 9 patients with a scapular fracture with glenoid extension treated surgically with an open procedure. One case of heterotopic ossification and 1 case of glenohumeral OA were identified after a mean of 4 years of follow-up. The largest study was reported by Mayo et al. [START_REF] Mayo | Displaced fractures of the glenoid fossa. Results of open reduction and internal fixation[END_REF] Of the 27 treated patients, 89% had anatomic reduction and 82% had good or excellent clinical outcomes. Only 1 infection was reported. Infections occurred in 9% of patients in the Schandelmaier et al [START_REF] Schandelmaier | Fractures of the glenoid treated by operation. A 5-to 23-year follow-up of 22 cases[END_REF] study, in which 22 glenoid fractures were fixed through a posterior or anterior approach. This was also the most frequent complication (2 of 10) in our series of patients treated by open surgery, whereas no infections occurred in the arthroscopy group. This difference has been reported with other types of injuries and appears to be an advantage common to every arthroscopic and minimally invasive treatment. [START_REF] Baker | Arthroscopic versus open rotator cuff repair: which has a better complication and 30-day readmission profile?[END_REF][START_REF] Day | Comparison of short-term complications after rotator cuff repair: open versus arthroscopic[END_REF] During the open anterior surgical approach, subscapular tenotomy allowed us to verify the articular reduction of the glenoid fracture. As with surgery for chronic instability, subscapularis tenotomy can-by itself-cause the postoperative stiffness and strength deficit that we observed. [START_REF] Bonnevialle | Selective capsular repair for the treatment of anterior-inferior shoulder instability: review of seventy-nine shoulders with seven years' average follow-up[END_REF][START_REF] Samilson | Dislocation arthropathy of the shoulder[END_REF][START_REF] Sato | Architectural and biochemical adaptations in skeletal muscle and bone following rotator cuff injury in a rat model[END_REF] Simply splitting the subscapularis muscle may have reduced the consequences of open surgery, as recommended for a Latarjet bone block procedure. [START_REF] Paladini | Long-term subscapularis strength assessment after Bristow-Latarjet procedure: isometric study[END_REF] The arthroscopic stabilization technique described by Sugaya et al [START_REF] Sugaya | Arthroscopic osseous Bankart repair for chronic recurrent traumatic anterior glenohumeral instability[END_REF] was used most often in our cohort. In patients with chronic instability, Sugaya et al [START_REF] Sugaya | Arthroscopic osseous Bankart repair for chronic recurrent traumatic anterior glenohumeral instability[END_REF] reported the outcomes after tendon-to-bone suture repair with anchors of the anteroinferior bone fragment that ranged in size from 2% to 21% of the glenoid's diameter. We used a mixed technique to leverage the advantages associated with anchors to reduce the glenoid fragment and stiffer fixation through screwinduced compression, to achieve optimal primary stability.

Tauber et al [START_REF] Tauber | Arthroscopic screw fixation of large anterior glenoid fractures[END_REF] described the outcomes of 10 patients who underwent arthroscopic screw fixation at more than 2 years of follow-up. As in our study, they reported no serious intraoperative complications. Nevertheless, the screws had to be removed in 1 patient because of impingement with the humeral head during internal rotation. [START_REF] Scheibel | Open reconstruction of anterior glenoid rim fractures[END_REF] This problem is not specific to the surgical technique but rather to screw fixation and its orientation. Absorbable screws can be used to avoid this complication; however, their mechanical properties are controversial. [START_REF] Balestro | Graft osteolysis and recurrent instability after the Latarjet procedure performed with bioabsorbable screw fixation[END_REF][START_REF] Scheibel | Arthroscopic reduction and fixation of large solitary and multifragmented anterior glenoid rim fractures[END_REF] Scheibel et al [START_REF] Scheibel | Arthroscopic reduction and fixation of large solitary and multifragmented anterior glenoid rim fractures[END_REF] found OA in 28% (6 of 21) of their patients at mean follow-up of 33 months. The presence of OA was correlated to patient age at the time of the surgery, not to the quality of the articular reduction. This correlation had previously been established in studies of chronic anterior shoulder instability. [START_REF] Bonnevialle | Selective capsular repair for the treatment of anterior-inferior shoulder instability: review of seventy-nine shoulders with seven years' average follow-up[END_REF][START_REF] Buscayret | Glenohumeral arthrosis in anterior instability before and after surgical treatment: incidence and contributing factors[END_REF][START_REF] Mizuno | Long-term results of the Latarjet procedure for anterior instability of the shoulder[END_REF] In our study, the type of surgical technique did not affect OA development, despite the reduction being superior when an open procedure was performed (90% vs. 67%).

Our study has certain limitations related to its retrospective nature. First, the groups were not fully comparable in the number of patients and associated injuries. Glenoid fractures associated with a proximal humeral fracture were preferentially treated with open surgery. Nevertheless, the subgroup analysis found no significant effect of this parameter on the final Constant Score.

Second, reduction and healing was only assessed on plain x-ray images, which could have underestimated step-off or nonunion.

Our study's greatest strengths are (1) comparison of 2 types of surgical techniques and (2) a larger sample size than in other published studies.

Conclusion

Surgical treatment of Ideberg IA glenoid fractures has similar early to midterm outcomes, whether performed with a conventional open technique or with arthroscopic surgery. However, range of motion seems to be better and the complication and reoperation rates are lower after the arthroscopic technique. The presence of osteochondral lesions negatively affects the outcome.

  was a retrospective study of patients treated between December 1, 2006, and January 1, 2016. Patients were included who had (1) a recent Ideberg IA 9 glenoid fracture (<3 weeks), (2) had been treated with standard open surgery (group O) or arthroscopic surgery (group A), and (3) had at least 12 months of follow-up consisting of clinical and radiologic assessments. Patients with Ideberg IB, II, III, IV, or V glenoid fractures were excluded.

Figure 1 (

 1 Figure 1 (A) Axial computed tomography slice shows an Ideberg IA glenoid fracture in a left shoulder. (B) Radiographs taken at 5 years after open reduction and internal fixation (anteroposterior view).

Figure 2 (

 2 Figure 2 (A) Computed tomography with 3-dimensional reconstruction and humeral head subtraction showing an Ideberg IA glenoid fracture in the left shoulder. (B) Sagittal computed tomography slice shows the arthroscopic screw fixation.

Table I

 I Comparison of epidemiologic data for group O (open) and group A (arthroscopy)

	Variable	Group O	Group A	P value
		(n = 10)	(n = 45)
	Age, yr	47 ± 12 (20-65)	41 ± 17 (22-77)	.87
	Sex			.75
	Male	8	34
	Female	2	11
	Sports-related injury	4	13	>.99
	Work compensation	1	1	.23
	Concurrent dislocation	8	40	.59
	Concurrent proximal humeral fracture	4 (40)	6 (13)	.04
	Glenoid fragment width, mm	9.5 ± 2.7 (7-17)	10.2 ± 2.5

Table II

 II Comparison of Constant score for group O (open) and group A (arthroscopy)

	Constant score	Group O (n = 10)	Group A (n = 45)	P value
	Pain (/15 points)	12.5 ± 4 (5-15)	12.6 ± 3 (5-15)	.6
	Activity (/20 points)	16.1 ± 6 (3-20)	18.3 ± 3 (6-24)	.4
	Mobility (/40 points)	32.6 ± 11 (8-40)	34.3 ± 6 (20-40)	.4
	Strength (/25 points)	13.2 ± 7 (2-25)	20.2 ± 4 (10-25)	.001
	Total (/100)	74 ± 6 (20-100)	84 ± 12 (52-98)	.07
	Adjusted, %	120 ± 38 (29-144)	133 ± 13 (91-142)	.06
	Data are presented as mean ± standard deviation (range).

Table III

 III Comparison of active range of motion for group O (open) and group A (arthroscopy)

	Range of motion	Group O	Group A (n = 45)	P value
	Anterior elevation, °144 ± 30 (180-45)	164 ± 25 (180-90)	.04
	External rotation, °34 ± 18 (10-50)	49 ± 7 (10-80)	.006
	Internal rotation, points	8.5 ± 3 (2-10)	8.2 ± 2 (2-10)	.84

Data are presented as mean ± standard deviation (range).

Table IV

 IV

	Presence of glenohumeral osteoarthritis in group O
	(open) and group A (arthroscopy)
	Variable	Group O	Group A	P value
		(n = 10)	(n = 45)
	Follow-up, mo	25 (12-72)	33 (12-115)
	Samilson OA stage 22
	1	1 (10)	4 (9)	NR
	2	0	3 (7)
	3	0	1 (2)
	Total	1 (10)	8 (18)	.52

OA, osteoarthritis; NR, not reported. Data are presented as mean (range) or number (%).

Disclaimer

All the authors, their immediate families, and any research foundations with which they are affiliated have not received any financial payments or other benefits from any commercial entity related to the subject of this article.


Debarati Chatterjee
email: dchatterjee@lpccaen.in2p3.fr

Jérôme Novak
email: jerome.novak@obspm.fr

Micaela Oertel
email: micaela.oertel@obspm.fr

Magnetic field distribution in magnetars

Keywords: numbers: 97.60.Jd, 26.60.-c, 26.60.Dd, 04.25.D-, 04.40.Dg

Using an axisymmetric numerical code, we perform an extensive study of the magnetic field configurations in non-rotating neutron stars, varying the mass, magnetic field strength and the equation of state. We find that the monopolar (spherically symmetric) part of the norm of the magnetic field can be described by a single profile, that we fit by a simple eighth-order polynomial, as a function of the star's radius. This new generic profile applies remarkably well to all magnetized neutron star configurations built on hadronic equations of state. We then apply this profile to build magnetized neutron stars in spherical symmetry, using a modified Tolman-Oppenheimer-Volkov (TOV) system of equations. This new formalism produces slightly better results in terms of massradius diagrams than previous attempts to add magnetic terms to these equations. However, we show that such approaches are less accurate than usual, non-magnetized TOV models, and that consistent models must depart from spherical symmetry. Thus, our "universal" magnetic field profile is intended to serve as a tool for nuclear physicists to obtain estimates of magnetic field inside neutron stars, as a function of radial depth, in order to deduce its influence on composition and related properties. It possesses the advantage of being based on magnetic field distributions from realistic self-consistent computations, which are solutions of Maxwell's equations.

I. INTRODUCTION

The macroscopic structure and observable astrophysical properties of neutron stars depend crucially on its internal composition and thus the properties of dense matter. The Equation of State (EoS) determines global quantities such as observed mass and radius. Transport properties such as thermal conductivity and bulk viscosity have an effect on cooling observations as well as emission of gravitational waves. As we enter an era of multimessenger astronomy, it is crucial to construct consistent microscopic and macroscopic models in order to correctly interpret astrophysical observations.

There are a large number of astrophysical observations, e.g. soft-gamma repeaters (SGR) or anomalous X-ray pulsars (AXP), that indicate the existence of ultramagnetized neutron stars or magnetars [1]. While such observations only probe the surface magnetic field, there is no way to measure directly the maximum magnetic field in the interior. Using the simple virial theorem, one may estimate the maximum interior magnetic field to be as high as 10 18 G. If such large fields exist in the interior, they may strongly affect the energy of the charged particles by confining their motion to quantized Landau levels and consequently modify the particle population, transport properties as well as the global structure [2][3][4][5][START_REF] Gomes | 14th International Workshop on Hadron Physics (Hadron Physics[END_REF][START_REF] Gomes | [END_REF][8][9][START_REF] Gomes | Compact Stars in the QCD Phase Diagram III[END_REF][START_REF] Dexheimer | [END_REF][12][13][14]. However, it is necessary to know the magnetic field amplitude at a given location in the star, i.e. a magnetic field distribution, in order to determine its effect on the internal composition and EoS.

The ideal way to tackle that problem would of course be to self-consistently solve the neutron star structure equations endowed with a magnetic field, i.e. combined Einstein, Maxwell and equilibrium equations, together with a magnetic field dependent EoS, as done by Chatterjee et al. [15]. This solution is complicated by the fact that in presence of a magnetic field, the neutron star structure strongly deviates from spherical symmetry and the spherically symmetric Tolman-Oppenheimer-Volkov (TOV) equations are no longer applicable for obtaining the macroscopic structure of a the neutron star [15][16][START_REF] Dexheimer | Compact Stars in the QCD phase diagram V (CSQCD V)[END_REF]. For small magnetic fields, perturbative solutions have been developed [18], but can no longer be applied for field strengths which might influence matter properties.

There have been several attempts to determine neutron star structure assuming an ad hoc profile of the magnetic field, without solving Maxwell's equations, within the TOV system (see e.g. [19][20][21]). To that end, many authors employ the parameterization introduced twenty years ago by Bandyopadhyay et al. [22], where the variation of the magnetic field norm B with baryon number density n B from the centre B c to the surface B s of the star is given by the form

B(n B /n 0 ) = B s + B c [1 -exp(-β(n B /n 0 ) γ )] , (1)
with two parameters β and γ, chosen to obtain the desired values of the maximum field at the centre and at the surface. This is an arbitrary profile, which possesses the same symmetries as the baryon density distribution in the star. Parameters (β, γ) are chosen such that the surface field is consistent with observations and the max-imum field prevailing at the center conforms to the virial theorem.

Lopes and Menezes [23] later introduced a variable magnetic field, which depends on the energy density rather than on the baryon number density:

B = B c ǫ M ǫ 0 γ + B s , (2)
where ǫ M is the energy-density of the matter alone, ǫ 0 is the central energy density of the maximum mass nonmagnetic neutron star and a parameter γ > 0, arguing that this formalism reduces the number of free parameters from two to one. The authors put forward as additional motivation the fact that it is the energy density and not the number density that is relevant in TOV equations for structure calculations. To account for anisotropy in the shear stress tensor, they applied the above field profile in a formalism [24], where the different elements containing the pressure are "averaged", leading to shear stress tensor of the form diag(B 2 /24π, B 2 /24π, B 2 /24π) [START_REF] Dexheimer | [END_REF]25]. Nevertheless, this approach within the TOV system is still spherically symmetric and the parameter γ not related to any experimental or observational constraint. There have also been suggestions of the magnetic field profile being a function of the baryon chemical potential [START_REF] Dexheimer | [END_REF] as:

B(µ B ) = B s + B c 1 -exp(b (µ B -938) a 938 ) , (3)
with a = 2.5, b = -4.08 × 10 -4 and µ B given in MeV.

In contrast to the profiles in Eqs. (1,2), such a formula avoids that a phase transition induces a discontinuity in the effective magnetic field. However, it was subsequently pointed out by Menezes and Alloy [26] that any of the above ad hoc formulations for magnetic field profiles are physically incorrect since they do not satisfy Maxwell's equations. In particular, it is obvious that assuming such a magnetic field profile in a spherically symmetric star implies a purely monopolar magnetic vector field distribution, which is incorrect. The inconsistency of this type of approach can be seen, too, by inspecting the most general solution of the equations of hydrostatic equilibrium in general relativity for a spherically symmetric star. In Schwarzschild coordinates, t, r, θ, φ , the line element reads:

ds 2 = -e -2Φ d t2 + 1 - 2Gm r -1 dr 2 +r 2 d θ2 + sin 2 θd φ2 , (4)
where m(r) and Φ(r) are the two relativistic gravitational potentials defining the metric (at the Newtonian limit, m represents the total mass enclosed in the sphere of radius r, and Φ/c 2 becomes the Newtonian gravitational potential). The resulting coupled system of equations for the star's structure has been derived by Bowers and Liang [27] and reads

dm dr = 4πr 2 ε dΦ dr = 1 - 2Gm rc 2 -1 Gm r2 + 4πG p r c 2 r dp r dr = -ε + p r c 2 dΦ dr + 2 r (p ⊥ -p r ) , (5)
with an energy-momentum tensor of the form T µν = diag(ε, p r , p ⊥ , p ⊥ ), where p r and p ⊥ are the radial and tangential pressure components. This is the most general energy-momentum tensor one can use assuming spherical symmetry and it goes beyond the perfect-fluid model, for which p r = p ⊥ . One may be tempted to cast a general electromagnetic energy-momentum tensor assuming a perfect conductor and isotropic matter, and for a magnetic field pointing in z-direction (see e.g. [15]) into this form. However, in the case of the electromagnetic energymomentum tensor T θθ = T φφ (look at Eqs. (23d)-(23e) of [15]), in clear contradiction with the assumption of Bowers and Liang (5) in spherical symmetry. Another problem arises from the fact that lim r→0 (T rr -T θθ ) = 0 and thus, the last term in Eq. ( 5) diverges at the origin (from the first line in this equation, one sees that the quantity m(r) ∼ r3 and therefore m/r 2 does not diverge). This discussion shows that there cannot be any correct description of the magnetic field in spherical symmetry.

Starting from two-dimensional numerical models, Dexheimer et al. [28,29] performed a fit to the shapes of the magnetic field profiles following the stellar polar direction as a function of the chemical potentials (as in [START_REF] Dexheimer | [END_REF]) by quadratic polynomials instead of exponential ones as

B(µ B ) = (a + bµ B + cµ 2 B ) B 2 c µ , (6)
where a, b, c are coefficients determined from the numerical fit. Unfortunately, no check of the validity of this fit has been shown in these works for other directions. In Ref. [30], a density dependent profile is applied within a perturbative axisymmetric approach à la Hartle and Thorne [31], but without solving Maxwell's equations. It remains, however, that the star's deformation due to the magnetic field implies that such a density (or equivalent) dependent profile depends on the direction, thus will be different looking e.g. in the polar or the equatorial direction.

In view of all these intrinsic difficulties, we will not propose here a simple scheme for solving structure equations of magnetized stars -to that end we refer to the publicly available numerical codes assuming axial symmetry [15,[START_REF] Bucciantini | XNS: Axisymmetric equilibrium configuration of neutron stars[END_REF]. Instead, since in many cases it might be sufficient to have an idea of the order of the value of the magnetic field strength to test its potential effect on matter properties, our aim is to provide a "universal" magnetic field strength profile from the surface to the interior obtained from the field distribution in a fully selfconsistent numerical calculation from one of these codes.

Further, we probe the applicability of this profile for determining the structure of magnetized neutron stars in an approximate way in spherical symmetry compared with full numerical structure calculations. As we will show, qualitatively the correct tendency can be reproduced for some NS properties, but to reproduce quantitatively correct results, the full solution has to be applied.

The paper is organized as follows. Sec. II describes our physical models, including the EoSs we use in this manuscript, together with the numerical techniques applied to solve the models. Sec. III provides the magnetic field profiles derived numerically by varying certain physical parameters, to achieve a generic profile for the monopolar part of the norm of the magnetic field. This profile is then applied in Sec. IV to a modified TOV system, to see its effect on NS masses and radii. Finally, Sec. V gives a summary of our work, together with some concluding remarks.

II. FORMALISM AND MODELS

In this section, we summarize the numerical approach for self-consistently modelling magnetized neutron stars. More details can be found in in [15,16,33].

A. Non-rotating magnetized neutron stars in general relativity

Due to the high compactness of neutron stars, we consider models within the theory of general relativity and solve coupled Einstein-Maxwell partial differential equations. We follow the scheme described in Bonazzola et al. [37], who considered the general case of rotating neutron stars, with the assumptions of stationarity, axial and equatorial symmetry, and circular spacetime, where the metric is given in the quasi-isotropic gauge, different from that used in TOV systems (4), by:

ds 2 = -N 2 dt 2 + C 2 r 2 sin 2 θ (dϕ -N ϕ dt) 2 +A 2 dr 2 + r 2 dθ 2 , (7)
where N, N ϕ , A and C are the relativistic gravitational potentials which are, as all other fields in Secs. II-III, functions of the coordinates (r, θ) only (independent from the ϕ-coordinate).

In this work, we shall restrict ourselves to the case without rotation, which in particular implies that there is no electric field in the models (perfect conductor). Nevertheless, as said in the introduction, the presence of a magnetic field induces a distortion of the stellar structure, which cannot remain spherically symmetric. Due to spacetime symmetries and circularity condition, only two magnetic field geometries can be described within this framework: a purely poloidal magnetic field (see Bocquet et al. [16]) or a purely toroidal one (see Kiuchi and Yoshida [38] and later Frieben and Rezzolla [39]).

In this work, we consider only purely poloidal magnetic fields, meaning that the only non-trivial components are B r (r, θ) and B θ (r, θ). This choice results in an asymptotically dipolar magnetic field distribution.

Matter is supposed to be composed of a perfect fluid coupled to the magnetic field. In Chatterjee et al. [15], it has been shown that the use of a magnetic field dependent EoS and inclusion of magnetization in the equations have negligible effects on neutron star structure, at least up to a polar magnetic field B pole ∼ 5 × 10 17 G (roughly corresponding to to a central magnetic field value b c ∼ 2 × 10 18 G) with a simple quark model EoS. We therefore neglect magnetic field dependency of the EoS and magnetization here, but they could be included in a straightforward way [15]. Matter is also assumed to be perfectly conducting and the magnetic field originates from free currents, moving independently from the perfect fluid. Equilibrium equations are obtained from the divergence-free condition of the energy-momentum tensor, and can be written as a first integral of motion, [15,37]. It is mostly the Lorentz force term in this equilibrium equation which distorts the stellar structure and makes it deviate from spherical symmetry. To summarize, given an equation of state (EoS) for nuclear matter (see Sec. II C hereafter), we thus solve the system of coupled Einstein-Maxwell equations, together with magnetostatic equilibrium. These models are then characterized by their gravitational mass (M G , see [37] for a definition), their EoS (see Sec. II C) and the central magnetic field, b c .

B. Numerical methods

The equations to be solved to get axisymmetric solutions form a set of six non-linear elliptic (Poisson-like) partial differential equations, coupled together with noncompact support (sources for gravitational field extend up to spatial infinity). These equations are solved using the same procedure as described in Bocquet et al. [16], employing the numerical library lorene [START_REF] Gourgoulhon | LORENE: Spectral methods differential equations solver[END_REF] based on spectral methods for the representation of fields and the resolution of partial differential equations (see Grandclément and Novak [41]).

Numerical accuracy of the axisymmetric solutions is checked through an independent test, the so-called relativistic virial theorem (Bonazzola and Gourgoulhon [42] and Gourgoulhon and Bonazzola [43]). This gives an upper bound on the relative accuracy of the obtained numerical solution, and we checked that it always remained lower than 10 -4 for the axisymmetric models presented in Sec. III.

C. Equations of state

The system of equations described above is closed by the EoS for nuclear matter relating the pressure p to the I. Saturation density , nsat, the binding energy at saturation EB, the compression modulus K, the symmetry energy, Esym and its slope, L, of symmetric nuclear matter are listed for the different EoS models employed. The neutron star maximum gravitational mass and radius at a fiducial mass of MG = 1.4M⊙ are given for cold spherical stars. Λ is the tidal deformability during inspiral of a binary neutron star merger, calculated with the chirp mass as measured for GW170817 [34,35]. Λ is only very weakly dependent on the mass ratio of the two stars, q = M1/M2 and is given here for a reference value of q = 0.8 (the same EoS is assumed for both stars). Note that the STOS EoS is as well excluded by the constraints on Λ obtained by GW170817 [35] as by nuclear physics experiments indicating much lower values for the slope of the symmetry energy [36]. We will nevertheless keep this EoS model as representative of an extreme case with large neutron star radius and symmetry energy slope.

Model nsat EB K Esym L Mmax R1.4 Λ(q = .8) (fm -3 ) (MeV) (MeV) (MeV) (MeV) (M⊙) (
baryon density n B . Our selection of EoSs for the present work has been guided by the idea to represent a large variety of different neutron star compositions and nuclear properties, derived from completely different nuclear physics formalisms. This was done in order to achieve an unbiased universal parameterization applicable to any realistic nuclear EoS. We consider one EoS model resulting from a microscopic calculation ("BL EOS") [44] 1 . It uses the Brueckner-Hartree-Fock formalism to tackle the many-body-problem and employs chiral interactions for the basic two-and three-body nuclear interactions.

In addition, we consider several phenomenological mean field models. These models are unified models in the sense that the crust EoS has been obtained with the same nuclear interaction than the one for the core guaranteeing consistency at the crust-core transition. They include two non-relativistic Skyrme parameterizations ("SLy9"and "SLy230a") [45,46] with the crust model from Gulminelli and Raduta [47], two relativistic mean field models ("STOS" and "HS(DD2)") [48,49] with the crust obtained from the model in Hempel and Schaffner-Bielich [50], supposing a temperature of 0.1 keV. The former one contains non-linear interactions whereas the latter one is constructed with density-dependent couplings. One model with hyperons ("SFhoY") [51] completes our list of EoS models. It is a nonlinear model where the crust is again obtained from Ref. [50]. Hyperonic interactions have been chosen to correctly reproduce hypernuclear data and a neutron star maximum mass above current observational limits. Some nuclear and neutron star properties of the different EoS models are listed in Table I. All EoS data are available from the on-line database CompOSE [52]. 1 The model calculation exist only for homogeneous matter and a crust has been added, see the CompOSE entry for details.

FIG. 1. Magnetic field lines in the (x, z)-plane for a MG = 2 M⊙ neutron star model endowed with a magnetic field which central value is bc = 5 × 10 17 G and using the SLy203a EoS of Tab. I. Thick line denotes the surface of the star and the magnetic moment is along the z-axis.

III. GENERIC MAGNETIC FIELD PROFILE

The numerical models of neutron stars endowed with a magnetic field described in Sec. II A consider two components (B r and B θ ) of the magnetic field vector, as measured by the Eulerian observer (see Bocquet et al. [16] for details). In the case of non-rotating stars considered here, this magnetic field is the same as that measured in the fluid rest-frame, denoted as b r and b θ in Chatterjee et al. [15]. As an example, the magnetic field distribution of a full neutron star model is displayed in Fig. 1, for a central value of the magnetic field b c = 5 × 10 17 G. The surface of the star (thick line) does not exhibit any significant deviation from spherical shape, but it is clear that the magnetic field distribution is dominated by the dipolar structure and cannot be accurately described by any spherically-symmetric model.

When trying to parameterize the magnetic field profile, the simplest approach is to consider the norm of the magnetic field, namely

b = g rr (b r ) 2 + g θθ (b θ ) 2 = A (b r ) 2 + (b θ ) 2 , (8)
where the relativistic gravitational potential A(r, θ) has been defined in Eq. [START_REF] Gomes | [END_REF]. Note that b is the quantity that enters the EoSs which take into account magnetization, as explained e.g. in [15]. The central value of this magnetic field norm is denoted as b c = b(r = 0) (independent of θ). In the rest of this work, we will consider this field as the main object of our study. As stated in the introduction, several authors have considered a parameterization of the magnetic field norm by the baryon density n B . In Fig. 2 we have plotted, for the same neutron star model of M G = 2 M ⊙ and b c = 5 × 10 17 G as in Fig. 1, the norm of the magnetic field as a function of baryon number density n B , along two radial directions: for θ = 0 (passing through the pole) and for θ = π/2 (passing through the equator). As these two curves show noticeable differences, including close to the center of the star (n B ∼ 0.8 fm -3 ), it seems that this type of parameterization can induce some inconsistency when describing magnetic field in a neutron star. We therefore try to improve it and adopt a different approach, taking a multipolar expansion of the magnetic field norm (Y m ℓ (θ, ϕ) being the spherical harmonic functions): In Fig. 3 we have plotted the first four non-zero terms of this multipolar decomposition as functions of the coordinate radius r. Note that, because of the symmetry with respect to the equatorial plane, odd-ℓ terms in the decomposition ( 9) are all zero. It appears that, at least in the high-density central regions of the star, the monopolar term b 0 (r), which is spherically symmetric, is dominant over the others. It is important to stress here that, contrary to the magnetic (vector) field, which has no monopolar part in terms of vector spherical harmonics, the norm of the vector field considered here is a scalar field which can possess a monopolar component.

b(r, θ) ≃ Lmax ℓ=0 b ℓ (r) × Y 0 ℓ (θ). (9)
We then look at the behavior of the radial profile of b 0 when varying the neutron star model in Fig. 4. On the left panel, we vary the gravitational mass of the star (either 1.6 M ⊙ or 2 M ⊙ ), as well as the amplitude of the magnetic field central value b c (10 15 G, 10 17 G and 10 18 G). On the right panel of Fig. 4, we vary the EoS used in the stellar model, keeping the gravitational mass to be 2 M ⊙ and the central magnetic field b c = 5 × 10 17 G. These profiles are no longer displayed as functions of the quasi-isotropic coordinate radius r, defined by the line element [START_REF] Gomes | [END_REF], but in view of the application to TOVsystems in Sec. IV, we consider here the Schwarzschild coordinate radius r, defined by the line element (4). The gauge transformation is obtained numerically and profiles are displayed as functions of this radius divided by the star's mean radius r mean which is such that the integrated (coordinate-independent) surface of the star reads A = 4πr 2 mean . Indeed, when the star gets distorted because of the magnetic field, it is difficult to define uniquely a relevant radius. In that sense, r mean is directly connected to the star's surface and some of its emission properties. It is remarkable that, although all possible parameters defining a magnetized stellar model (mass, central magnetic field, EoS) have been varied, all profiles are quite similar and deviate one from another only by a few percent. The only case where a noticeable difference appears is when using quark matter EoS. Therefore, we make the following conjecture: the monopolar part of the norm of the magnetic field follows a universal profile, up to minor variations, when considering different neutron star models with realistic hadronic EoSs. This "universal" profile has been fitted using a simple polynomial:

b 0 (x) = b c × 1 -1.6x 2 -x 4 + 4.2x 6 -2.4x 8 , (10)
where x = r/r mean is the ratio between the radius r in Schwarzschild coordinates (4) and the star's mean (or areal) radius. Let us stress that the aim of the present investigation is to obtain a universal profile for realistic EoSs and that we have therefore excluded polytropic EoSs. A preliminary calculation showed that the general parameterization applicable to the family of realistic EoSs is not applicable directly to the case of polytropes without specific fine tuning.

Going further, we display in Fig. 5 profiles for the dipolar part b 2 (r), defined in Eq. (9). A larger dispersion is visible in these plots, in particular in terms of EoSs (right panel of Fig. 5) and for the largest value of magnetic field (b c = 10 18 G, left panel of Fig. 5). This last point can be understood from the large deformation undergone by the star at this value of central magnetic field, when the contribution from higher-order multipoles starts to become important. Together with Fig. 3, these curves show that spherical symmetry is not a good approximation for strongly magnetized neutron stars, where higher multipoles can have a non-negligible effect. However, the relative robustness of the profiles shown in Fig. 5 indicates that, for not too large magnetic fields, the dipolar correction may be included in an EoS-independent way and some perturbative approach to spherical symmetry may be devised, in a way similar to Konno et al. [18]. We leave these investigations to further studies.

IV. APPLICATION TO A TOV-LIKE SYSTEM

As discussed in the introduction, it is fundamentally inconsistent to solve spherically symmetric equations for magnetized neutron star models since it completely neglects the star's deformation due to the electromagnetic field. It is, however, tempting, to have a simple approach at hand which allows at least to qualitatively reproduce the effects of the magnetic field on (some) neutron star properties performing calculations only slightly more complicated than solving TOV equations. To that end, we modify the TOV system by adding the contribution from the magnetic field to the energy density and a Lorentz force term to the equilibrium equation:

dm dr = 4πr 2 (ε + b 2 µ 0 ) dΦ dr = 1 - 2Gm rc 2 -1 Gm r2 + 4πG p c 2 r dp dr = -ε + b 2 µ 0 + p c 2 dΦ dr -L(r) . (11)
L(r) denotes here the Lorentz force contribution, which is noted dM/dr in Bonazzola et al. [37] (see this reference for more details). Note that these equations [START_REF] Dexheimer | [END_REF] are not derived from any first principle, but only motivated by Eq.( 5), in which we have replaced the diverging term by some phenomenological term, supposed to better take into account the magnetic pressure and the Lorentz force acting on the fluid.

Similar to the magnetic field norm b, we found from the full numerical calculations the following parametric form

L(r) = 10 -41 × b 2 c -3.8 x + 8.1 x 3 -1.6 x 5 -2.3 x 7 , (12)
where x = r/r mean and the central magnetic field, b c , is given in units of G. For the magnetic field b in Eqs. [START_REF] Dexheimer | [END_REF], the profile (10) is applied.

In order to get an idea of the quality of this "TOV-like" approach, we show in Figs. [START_REF] Gomes | 14th International Workshop on Hadron Physics (Hadron Physics[END_REF][START_REF] Gomes | [END_REF] a comparison between stellar sequences (described in mass vs. radius diagrams) obtained with the TOV-like approach in spherical symmetry and the full numerical solution in axial symmetry. In Fig. [START_REF] Gomes | 14th International Workshop on Hadron Physics (Hadron Physics[END_REF] the gravitational mass vs. the mean radius is displayed for a central magnetic field of b c = 10 17 G (l.h.s.) and b c = 10 18 G (r.h.s). As expected, deviations become larger at smaller masses since the ratio of magnetic to matter pressure increases and thus the stars are more strongly deformed and the relation between the mean radius and the radius of a spherically symmetric configuration is no longer obvious. The same is true if the central magnetic field is increased. For further comparisons, we show in addition the solutions obtained with TOV applying the "averaged" shear stress tensor [24], mentioned in Sec. I, and those obtained without any magnetic field (B = 0), too. We display here results from the former approach since it has been used several times in the literature [START_REF] Dexheimer | [END_REF]24,25] for studying properties of magnetized stars. Let us emphasize, however, that the procedure to obtain the "averaged" stress-energy tensor is mathematically ill-defined, as it is not possible to directionally "average" elements of a tensor. Moreover, looking at Fig. 6 results with this approach for the massradius relation most strongly differ from the full solution. Surprisingly, the usual TOV solution with no magnetic field (B = 0), on the contrary, very well reproduces the full solution. Thus, although the star becomes strongly deformed, the mean radius is only marginally influenced by the magnetic field.

Masses should be less sensitive to the deformation. As can be seen in Fig. [START_REF] Gomes | [END_REF], indeed the TOV-like solution for the maximum gravitational mass as well as the gravitational mass for fixed baryon mass as function of the central magnetic field show the correct qualitative behavior. Both increase with b c and the TOV-like approach overestimates the masses up to a factor two in the correction (with respect to the non-magnetized case) at central fields of b c = 10 18 G. This difference shows the difficulty of "spherically-symmetric" approaches to model magnetized neutron stars. The limits of the TOV-like approach can be clearly seen in the determination of effects due to the magnetic field. Finally, note that the "averaged" TOV approach gives a wrong qualitative behavior for the dependence of the maximum mass on the central magnetic field value. Thus, although our investigations can serve as a guideline and reproduce at least for gravitational mass as function of magnetic field the correct qualitative tendency, it should be stressed that it is strongly recommended to use some consistent axisymmetric or three-dimensional approach (e.g. employing publicly available software), to determine properties of magnetized neutron stars or to draw any quantitative conclusion. V. CONCLUSIONS Many attempts can be found in the literature trying to study strongly magnetized neutron stars and to include magnetic field effects on the matter properties. As mentioned in the introduction, most of these investigations suffer from different assumptions and approximations motivated by the complexity of the full system of equations. First, in order to avoid solving Maxwell's equations in addition to equilibrium and Einstein equations, often an ad hoc profile for the magnetic field is assumed, which has no physical motivation. Second, spherical symmetry is assumed for modelling the star.

In this work, we tackle the first point: we proposed a "universal" parameterization of the magnetic field profile (Eq. 10) as a function of dimensionless stellar radius, ob-tained from a full numerical calculation of the magnetic field distribution. We tested this profile against several realistic hadronic EoSs, based on completely different approaches, and with different magnetic field strengths in order to confirm its universality. For the case of quark matter EoSs, preliminary investigations showed that although MIT bag models conform to the universality, other quark matter EoSs may not necessarily do so. The profile is intended to serve as a tool for nuclear physicists for practical purposes, namely to obtain an estimate of the maximum field strength as a function of radial depth, within the error bars observed in our study, e.g. in Fig. (4), in order to deduce the composition and related properties.

We applied the proposed magnetic field profile in a modified TOV-like system of equations, that include the contribution of magnetic field to the energy density and pressure, and account for the anisotropy by introducing a Lorentz force term. Compared with full numerical structure calculations, we find that qualitatively the correct tendency is reproduced and quantitatively the agreement is acceptable for large masses and small magnetic fields (b c 10 17 G). However, we find that the standard TOV system with no magnetic field reproduces much better mass-radius relations, even for strong magnetic fields, than any modified TOV system, with poorly defined magnetic corrections. This is mostly due to the fact that the mean radius is only marginally changed by the magnetic field. We thus think that future studies should employ the profile proposed here to conclude about the importance of magnetic field effects on matter properties, and use TOV system at B = 0 for calculating mass-radius diagrams. For any other property of magnetized stars, we can only recommend the use of a full axisymmetric numerical solution for modelling magnetized neutron stars.

2 FIG. 2 .

 22 FIG. 2. Magnetic field norm b (8) as a function of the baryon density for two angular directions for the same stellar model as in Fig. 1.

6 FIG. 3 .

 63 FIG. 3. Radial profiles of the first four even multipoles, b ℓ , ℓ = 0, 2, 4 and ℓ = 6, see definition (9), of the magnetic field norm b(r, θ) computed for the stellar model described in Fig. 1. From symmetry arguments odd multipoles are all zero.

FIG. 4 .

 4 FIG. 4. Monopolar part of the magnetic field profile b0(r), normalized to its central value for different magnetized neutron star models, as functions of the radius expressed in Schwarzschild coordinates r (see expression (4)) divided by the star's mean radius (see text). Left panel (a): all models are using SLy230a EoS (see Tab. I) but have different masses (1.6 M⊙ or 2 M⊙) and different central magnetic fields bc (10 15 G, 10 17 G or 10 18 G). Right panel (b): all are 2 M⊙ models, with a central magnetic field bc = 5 × 10 17 G, but with different EoSs, see Tab. I for details.

FIG. 5 .

 5 FIG.5. Same as Fig.4but for the dipole term b2(r) of the multipolar expansion (9) of the magnetic field norm.

FIG. 6 .

 6 FIG.6. (color online) Gravitational mass vs. mean radius stellar sequences for the TOV-like solution, TOV solution using a directionally-averaged energy momentum tensor (see text) and the exact (axisymmetric) calculation for a central magnetic field of bc = 10 17 G (left) and bc = 10 18 G (right) employing the SLy230a EoS. Note that, in both figures, the thin-dotted lines (B = 0) correspond to solutions of the TOV system with no magnetic field at all.

ACKNOWLEDGMENTS

This work was financially supported by the action "Gravitation et physique fondamentale" of Paris Observatory. DC acknowledges financial support from CNRS and technical support from LPC Caen and the Paris Observatory in Meudon where part of this study was completed.