Yves Henchoz

Nicola Soldini

Nicolas Peyrot

Davide Malatesta
email: davide.malatesta@unil.ch

Guido Ferretti

Energetics and mechanics of walking in patients with chronic low back pain and healthy matched controls

come

Introduction

Chronic low back pain (cLBP) is a common problem throughout the world, which has a huge impact at the individual, family, community, and societal levels [START_REF] Hoy | The Epidemiology of low back pain[END_REF]. The typical pattern is to avoid movements or activities that exacerbate pain [START_REF] Leeuw | The fear-avoidance model of musculoskeletal pain: current state of scientific evidence[END_REF]. Although studies comparing the overall activity level between patients with cLBP and controls are inconsistent, the distribution of physical activity over the course of a day seems to be different (Griffin e t al. 2012). T he r ecent literature suggests that individuals develop motor control strategies at multiple levels of the motor system that have short-term benefit, but with potential harmful consequences if maintained in the long term [START_REF] Hodges | Pain and motor control: from the laboratory to rehabilitation[END_REF].

Walking is an essential element in everyday activities, and constitutes a major component in current physical activity guidelines for patients with LBP [START_REF] Simmonds | Lower Back Pain Syndrome. ACSM's exercise management for persons with chronic diseases and disabilities, 3rd edn[END_REF]. Compared to asymptomatic individuals, patients with cLBP spontaneously adopt a decreased walking speed [START_REF] Lamoth | Effects of attention on the control of locomotion in individuals with chronic low back pain[END_REF][START_REF] Simmonds | The influence of pain distribution on walking velocity and horizontal ground reaction forces in patients with low back pain[END_REF][START_REF] Lee | Influence of pain distribution on gait characteristics in patients with low back pain: part 1: vertical ground reaction force[END_REF]). More specifically, s everal b iomechanical changes have been reported. Lamoth et al. (2006b) observed a decreased variability in the coordination between trunk and pelvic rotations in the transverse plane, but an increased variability in the frontal plane. The same group reported a decreased ability of patients to adapt trunk-pelvis coordination to changes in velocity (Lamoth et al. 2006a), which was aggravated under the influence o f a n a ttention-demanding task [START_REF] Lamoth | Effects of attention on the control of locomotion in individuals with chronic low back pain[END_REF], suggesting a stronger cognitive regulation of gait. From the study of [START_REF] Crosbie | Coordination of spinal motion in the transverse and frontal planes during walking in people with and without recurrent low back pain[END_REF], it appears that alterations in spinal movement and coordination persist, yet in a small extent, in patients with recurrent low back pain who are pain free at the time of testing. These adaptations have been suggested to result from increased trunk stiffness as a protective strategy to prevent further injury or pain [START_REF] Van Den Hoorn | Mechanical coupling between transverse plane pelvis and thorax rotations during gait is higher in people with low back pain[END_REF].

Muscle activation patterns confirm the hypothesis of a guarding mechanism during walking in patients with cLBP. Increased lumbar muscle activation has been reported during the total stride (van der Hulst et al. 2010b;[START_REF] Vogt | Neuromuscular control of walking with chronic low-back pain[END_REF], including both swing (Lamoth et al. 2006b;van der Hulst et al. 2010b) and double support periods (van der Hulst et al. 2010b). Cocontraction of the erector spinae and rectus abdominis is also increased in patients with cLBP compared to asymptomatic controls (van der Hulst et al. 2010a). In healthy participants, experimentally induced pain or fear of pain during walking had only subtle effects on trunk coordination and erector spinae activity, suggesting that the aforementioned adaptations are long-term rather than short-term consequences of pain [START_REF] Lamoth | Effects of experimentally induced pain and fear of pain on trunk coordination and back muscle activity during walking[END_REF].

In normal subjects, walking is characterized by exchanges between potential and kinetic energy, as illustrated by the inverted pendulum model developed by [START_REF] Cavagna | The sources of external work in level walking and running[END_REF] This transfer of energy minimizes the work to be produced by muscles, and hence the energy cost of walking. In several pathologic conditions, gait pattern alterations are associated with an increased energetic cost of walking [START_REF] Mahaudens | Gait in adolescent idiopathic scoliosis: energy cost analysis[END_REF][START_REF] Bernardi | Cost of walking and locomotor impairment[END_REF]). In patients with cLBP, the aforementioned adaptations (altered trunk-pelvis coordination, increased muscle activation and stronger cognitive regulation of gait) may also cause an augmented production of mechanical work and an increased metabolic demand. To our knowledge, no study has been undertaken on this issue.

The present study aimed to compare the preferred walking speed, the biomechanical and the energetic parameters of walking at different speeds between patients with cLBP and healthy men. It was hypothesized that the preferred walking speed would be lower, and that the mechanical and metabolic demand at fixed/standard speed of walking would be higher in patients with cLBP compared to healthy participants.

Methods

Participants

The present study included 13 men with cLBP (cLBP group) recruited among patients on the waiting list for functional multidisciplinary rehabilitation at the Service of Rheumatology of Lausanne University Hospital (Switzerland). Thirteen healthy control men (CTR group) were recruited through posted announcements, and individually matched to patients for age, body mass and height. For patients, inclusion criteria were age between 30 and 60 years and non-specific low back pain for at least 12 weeks. Exclusion criteria were irritative neurological deficit in progress, sciatica, acute inflammatory rheumatic disease, non-osteoarticular thoracic pain, spinal fracture within the last 3 months, osteoporosis, tumour, severe heart failure or respiratory failure, obesity (body mass index >30), entitled to a total disability pension, active drug addiction, current involvement in litigation related to low back pain, and active psychiatric pathology.

Participant recruitment and data collection were conducted from July to November 2011. Each participant was informed in detail about the study, and provided written consent. They received a financial compensation of CHF 50. Lausanne University Medical School's Clinical Research Ethics Committee approved the study.

Experimental design

Participants attended a single 2-h testing session. After 10 min of treadmill accommodation across experimental walking speeds and a brief rest period, the preferred walking speed was determined. Anthropometric measurements were then taken and subjects were asked to complete some questionnaires. Successively, participants were equipped with a breath-by-breath gas analyser and two inertial sensors fitted with a triaxial accelerometer and gyroscope and the treadmill protocol started.

Preferred Walking Speed (PWS)

After treadmill familiarization across experimental walking speeds (0.83, 1.11, 1.38, 1.67 m s -1 ) and a brief rest period, each subject's PWS was determined according to the methodology proposed by [START_REF] Martin | Effects of age and physical activity status on the speed-aerobic demand relationship of walking[END_REF] Briefly, each subject began treadmill walking at the lowest experimental speed (0.83 m s -1 ), which was then slowly increased until the subject subjectively identified his PWS. This speed was maintained for a minute and was then eventually slightly modified according to subject directives. This procedure was repeated starting with the highest experimental speed (1.67 m s -1 ) and gradually reducing to the preferred speed. The final PWS was considered to be the mean of the two speeds selected by the subject during both the increasing and decreasing speed trials.

Treadmill protocol

The treadmill protocol was determined according to the procedure proposed by [START_REF] Malatesta | Energy cost of walking and gait instability in healthy 65-and 80-year-olds[END_REF] After a 4-min measure of the standing rate of oxygen consumption, participants were asked to complete 5-min walking trials at five different speeds (0.83, 1.11, 1.38, 1.67 m s -1 and PWS) in randomized order, separated by 5-min resting periods. Participants were instructed to walk as naturally as possible in the middle of the belt and without using handrail support.

Assessments

Anthropometric measurements

Each subject was weighed with a standard balance. Stature was measured to the nearest 0.5 cm using a standardized wall-mounted height board, and body mass index (BMI) was calculated as body mass divided by height squared. The leg length was determined with the subject in standing position as the distance between the great trochanter and the ground for the left leg.

Questionnaires

The physical activity level of each participant was estimated using the short form of the International Physical Activity Questionnaire (IPAQ) [START_REF] Gauthier | Psychometric properties of the IPAQ: a validation study in a sample of northern Franco-Ontarians[END_REF], which covers activities performed at work, at home, when on the move and during leisure time. Algorithms that take into account the frequency, volume and intensity of the reported physical activities classify participants as 'high', 'moderate' or 'low' physical activity level. The disability level of cLBP group participants was measured with the Oswestry Disability Index (ODI), validated French version [START_REF] Vogler | Cross-cultural validation of the Oswestry disability index in French[END_REF]. This 10-item scale provides a score ranging from 0 to 100 %. A high score indicates a high degree of disability. A value of 12 was reported as a cut-off which separates cLBP patients with vs without disability [START_REF] Tonosu | The normative score and the cut-off value of the Oswestry Disability Index (ODI)[END_REF]. Fear of movement was assessed in cLBP participants using the Tampa scale for kinesiophobia (TSK), a 17-item scale answered using a 4-point Likert scale (from 'strongly disagree' to 'strongly agree') [START_REF] Miller | The Tampa Scale: a measure of kinesiophobia[END_REF]. A total score ranging from 17 to 68 was calculated after inversion of the scores on items 4, 8, 12 and 16. A high score indicates a high degree of fear of movement. A value over 37 has been reported to differentiate low-fear and high-fear subjects [START_REF] Vlaeyen | Fear of movement/(re)injury in chronic low back pain and its relation to behavioral performance[END_REF].

Pain and perceived exertion

Low back pain intensity for the cLBP group was measured with a visual analogue scale (VAS) [START_REF] Price | The validation of visual analogue scales as ratio scale measures for chronic and experimental pain[END_REF], before the treadmill protocol, after each walking trial and at the end of the treadmill protocol. The Borg scale was used to evaluate the perceived exertion after each walking trial and at the end of the protocol for both cLBP and CTR group participants [START_REF] Borg | Psychophysical bases of perceived exertion[END_REF].

Sitting and standing oxygen consumption

Oxygen uptake ( V O 2 ), carbon dioxide output ( V CO 2 ) and ventilation were measured using a breath-by-breath gas analyser (Metalyzer, Cortex Medical, Leipzig, Germany). The calibration of the Metalizer (pressure, volume and standard gases) was realized before each session. Sitting V O 2 and V CO 2 were measured during 10 min at a steady state, after the participant had maintained a sitting position for 20 min. Standing V O 2 and V CO 2 were measured for 4 min at the beginning of the treadmill protocol. Average V O 2 and V CO 2 were calculated over the last minute. Meta- bolic rate (W kg -1 ) was calculated based on the energetic equivalent of O 2 [START_REF] Åstrand | Textbook of work physiology : physiological bases of exercise. McGraw-Hill series in health education, physical education, and recreation[END_REF].

Energetics of walking and heart rate

Average V O 2 , V CO 2 and heart rate (HR; Polar S810, Kem- pele, Finland) were calculated over the last minute of each 5-min walking trial. Gross metabolic rate (W kg -1 ) for each speed was calculated based on the energetic equivalent of O 2 [START_REF] Åstrand | Textbook of work physiology : physiological bases of exercise. McGraw-Hill series in health education, physical education, and recreation[END_REF]. Gross energy cost of walking (gross C w , J kg -1 m -1 ) was then calculated as the gross metabolic rate (W kg -1 ) divided by the walking speed. A second-order least squares regression was used to model the U-shaped relationship between gross C w and walking speed and calculate the optimal walking speed (OWS: the most economical speed) for each subject. Net energy cost of walking (net C w , J kg -1 m -1 ) was calculated by subtracting standing metabolic rate (W kg -1 ) from gross metabolic rate and dividing by the walking speed.

Mechanics of walking

Two inertial sensors equipped with a triaxial accelerometer and gyroscope (MTx, Xsens, Enschede, The Netherlands) were used to calculate the three-dimensional accelerations of walking as previously described by [START_REF] Peyrot | Do mechanical gait parameters explain the higher metabolic cost of walking in obese adolescents?[END_REF]. An inertial/gyroscope sensor was taped and secured directly to the skin on the lower part of the back, at the L3 level (close to the centre of body mass, COM) using an adhesive strap. This sensor was used to measure the COM accelerations, under the assumption that changes in the relative positions of the sensor and COM over walking time may be neglected. Then, a second sensor was also taped and secured on the instep of subjects' left foot to detect foot events.

Data were processed according to the methodology previously described by [START_REF] Peyrot | Do mechanical gait parameters explain the higher metabolic cost of walking in obese adolescents?[END_REF] Briefly, the inertial sensors provided orientation data in the earth reference system in the form of Euler angles (roll, pitch and heading), which represent rotations of the sensor system into the earth reference system. Since the magnetic north corresponded to the anteroposterior axis in the present study, the three-dimensional accelerations of the two sensors were repositioned in the earth reference system using rotation matrices. All data were sampled at 100 Hz and low-pass filtered at 30 Hz (fourth-order, zero-lag, low-pass Butterworth).

Spatiotemporal parameters

Foot angular velocities were used to determine heel strike and toe-off [START_REF] Jasiewicz | Gait event detection using linear accelerometers or angular velocity transducers in able-bodied and spinal-cord injured individuals[END_REF]. Stride duration, stance duration and single support duration of the contralateral limb were computed as heel strike to consecutive heel strike, heel strike to consecutive toe-off, and toe-off to consecutive heel strike, respectively. Stance and single support durations were expressed relatively to stride duration (%). Stride frequency (Hz) was calculated as the inverse of stride duration. Stride length (m) was calculated by multiplying walking speed and by stride duration.

Mechanical work

Mechanical analyses were performed over ten consecutive strides taken during the last minute of each walking trial. A stride was defined as the period between two consecutive left heels. Mechanical parameters were computed for each stride and then averaged to describe a mean typical stride. According to [START_REF] Cavagna | The sources of external work in level walking and running[END_REF], vertical, forward and medio-lateral (M-L) COM velocities and positions were obtained by integrating twice the three-dimensional accelerations of the mean stride.

The total instantaneous potential (E p ; J) and kinetic (E k ; J) energies of the COM were calculated as follows:

where m is the body mass (kg); g is the gravitational constant; h is the vertical position of the COM; E kf , E kv and E kl are the forward, vertical and lateral kinetic energies, respectively; and V x , V y and V z are the forward, vertical and lateral velocity components (m s -1 ), respectively. The total mechanical energy of the COM (E tot ; J) was calculated as the sum of E p and E k and the external mechanical work (W ext ; J kg -1 m -1 ) as the sum of the positive increments in E tot divided by stride length and by body mass.

The inverted pendulum recovery of mechanical energy of the COM (Recovery; %) was calculated as follows:

where ∆E p (J) and ∆E k (J) are the fluctuations in potential and kinetic energy of the COM.

Mechanical internal work (W int , J kg -1 m -1 ) was calculated with the equation of [START_REF] Nardello | Measured and predicted mechanical internal work in human locomotion[END_REF] as follows:

where f is the stride frequency (Hz); v is the walking speed (m s -1 ); and d is the duty factor. The total mechanical work (W tot , J kg -1 m -1 ) was calculated as the sum of W ext and W int .

(1)

E p = m × g × h (2) E k = E kf + E kv + E kl = 0.5 × m × (V x ) 2 + V y 2 + (V z ) 2
(3)

Recovery = E p + E k -W ext E p + E k × 100 (4) W int = f × v × 1 + d 1 -d 2 × 0.08

Medio-lateral COM displacement

The medio-lateral COM displacement was equal to the total amplitude (from left to right) of the medio-lateral COM position calculated by integration of the medio-lateral velocity.

Statistical analysis

Data are expressed as means (standard deviation) for all variables. The sample size was based on a previous study [START_REF] Mahaudens | Gait in adolescent idiopathic scoliosis: energy cost analysis[END_REF]) that reported a net C w of 2.4 (0.4) J kg -1 m -1 and 1.8 (0.3) J kg -1 m -1 in adolescents with or without idiopathic scoliosis, respectively. To detect a 0.6 J kg -1 m -1 difference with a standard deviation of 0.4 J kg -1 m -1 , 90 % power and at the 5 % significance level, a minimum of 11 participants in both study groups was required. A t test was used to test differences between anthropometric characteristics of the two groups. A χ 2 was used to compare the distribution of the two groups in the three physical activity levels of IPAQ. A two-way repeatedmeasures mixed design ANOVA [walking speed (n = 4; 0.83-1.67 m s -1 ) × group (cLBP vs CTR)] followed by contrasts was used to determine the effect of pathology and speed on energetic and mechanical variables of walking. As it is well established that speed influences these variables, main effects of speed are not reported. A t test was used to determine difference in PWS and in energetic and mechanical at PWS between groups. Difference between PWS and OWS for the two groups was tested with a t test. The level of significance was set at P ≤ 0.05. The analyses were conducted using SPSS 21 software (IBM, Armonk, NY).

Results

Participants' characteristics

The characteristics of the study participants are presented in Table 1. There were no significant differences in body mass index (P = 0.64), lower limb length (P = 0.80) and individuals distribution in the number of children (P = 1.00), smoking status (P = 0.18), origin (P = 0.69) and the three physical activity levels of IPAQ (P = 0.14) between the two groups. Patients with cLPB were significantly more often manual workers than CTR participants (P = 0.05), which may explain the tendency towards a higher level of physical activity in the cLBP compared to CTR group. For cLPB participants, mean duration of symptoms was 5.6 (7.3) years, ODI and TSK scores were 29.2 (10.2) and 46.8 (6.2), respectively. These values indicate moderate levels of disability and fear of movement. The minimum and maximum ODI scores were, respectively, 14 and 50, indicating that all participants from the cLBP group were above the cut-off value of 12 defining patients with disability. The minimum and maximum TSK scores were, respectively, 36 and 56. A total of 12 patients had a TSK score >37 defining a high fear.

PWS (Table 2) and OWS [cLBP: 1.33 (0.09) m s -1 and CTR: 1.41 (0.09) m s -1 ] were significantly slower in cLBP than in CTR group (P = 0.002 and P = 0.03, respectively). OWS was significantly higher than PWS in cLBP group (P = 0.003), but not in CTR group (Fig. 1; P = 0.08).

Rest metabolic rate and energetics of walking

There were no significant differences in sitting [cLBP: 1.18 (0.22) W kg -1 ; CTR: 1.32 (0.19) W kg -1 ] and standing [cLBP: 1.38 (0.20) W kg -1 ; CTR: 1.46 (0.20) W kg -1 ] metabolic rate between the two groups (P = 0.09 and P = 0.32, respectively). Gross C w and net C w were not significantly different between the two groups at the four standard walking speeds (P = 0.88 and P = 0.36, respectively) and at PWS (Fig. 1; Table 2; P = 0.65 and P = 0.88, respectively). The group × walking speed interactions were not significant (P ≥ 0.05).

Mechanics of walking

Spatiotemporal parameters

There were no significant differences in stride frequency (P = 0.91), stride length (P = 0.99) and medio-lateral displacement (P = 0.16) at the four fixed walking speeds between the two groups (Fig. 2). In addition, for these speeds, the single and double support stride durations (%) were similar in the two groups (data not shown; P = 0.80 for both). At PWS, stride length and frequency were significantly lower (P = 0.004 and P < 0.001, respectively) and medio-lateral displacement was significantly higher (P = 0.03) in cLBP compared with CTR group (Table 2).

In contrast, at this speed, there were no significant differences in the single and double support stride durations (%) between the two groups (Table 2; P = 0.14 for both). The group × walking speed interactions were not significant (P ≥ 0.11).

Mechanical works and recovery

W tot , W ext , W int and recovery at the four standard walking speeds were not significantly different between the two groups (Fig. 3; P ≥ 0.44). At PWS, W tot , W ext , W int were significantly lower in cLBP than in CTR group (Table 2; P = 0.007, P = 0.04 and P < 0.001, respectively). There was no significant difference in recovery at PWS between the two groups (Table 2; P = 0.06). The group × walking speed interactions were not significant (P ≥ 0.10).

Heart rate, pain and perceived exertion

There was no significant difference in HR at the four fixed walking speeds between the two groups (Fig. 4a; P = 0.58). Perceived exertion showed a significant main group effect (P = 0.02) with Borg scale scores significantly higher at 1.11 and 1.38 m s -1 (P = 0.01 and P = 0.02, respectively) and tends to be higher at 0.83 and 1.67 m s -1 in cLBP compared with CTR group (Fig. 4b; P = 0.06 and P = 0.09, respectively). For both variables, the group × walking speed interactions were not significant (P ≥ 0.26). At PWS, perceived exertion and HR were similar in both groups (Table 2; P = 0.16 and P = 0.68, respectively). At the end of the experimental trail, Borg scale score was significantly higher in cLBP [11.7 (2.4)] than in CTR group [9.9 (1.1); P = 0.01]. For cLBP individuals, during experimental session, there was no significant time evolution in low back pain intensity assessed with VAS (Fig. 4c; P = 0.21).

Discussion

The aim of this study was to compare the preferred walking speed, the biomechanical and the energetic parameters of walking between patients with cLBP and matched healthy controls. Preferred walking speed was lower in patients with cLBP, but no significant difference was observed between groups in mechanical work, spatiotemporal parameters and energy cost of walking. These results do not support the hypothesis of a less efficient walking pattern in patients with cLBP.

Considering the multiple gait adaptations previously reported in patients with cLBP (i.e. altered trunk-pelvis coordination, increased muscle activation and stronger cognitive regulation of gait) [START_REF] Vogt | Neuromuscular control of walking with chronic low-back pain[END_REF][START_REF] Lamoth | Effects of experimentally induced pain and fear of pain on trunk coordination and back muscle activity during walking[END_REF]Lamoth et al. , 2006aLamoth et al. , b, 2008;;van der Hulst et al. 2010avan der Hulst et al. , b, 2012;;[START_REF] Crosbie | Coordination of spinal motion in the transverse and frontal planes during walking in people with and without recurrent low back pain[END_REF], an augmented production of mechanical work and an increased metabolic demand were expected in patients with cLBP compared to control participants. The absence of any difference in the energetics and mechanics of walking between groups may be explained by several reasons. The aforementioned motor control changes may not be sufficiently large to cause a less efficient walking pattern. In fact, it has previously been shown that the peripheral musculoskeletal disorders, inducing lowest-level gait disturbances, did not affect W ext and R, whereas central nervous system pathologies, causing middle-level gait disorders, altered the mechanics of walking impairing the pendular energy transfer and increasing W ext compared with normal gait [START_REF] Detrembleur | Motion of the body centre of gravity as a summary indicator of the mechanics of human pathological gait[END_REF]. This suggests that W ext and recovery may be "a sensitive indicator to distinguish between central and peripheral neurological disorders" [START_REF] Detrembleur | Motion of the body centre of gravity as a summary indicator of the mechanics of human pathological gait[END_REF]. Alternatively, these changes may occur without compromising the walking pattern efficiency, or in conjunction with compensating adaptations as a strategy to maintain the energetic cost of walking at a low level. As previously shown, unilateral amputee patients walk with C w , W ext and recovery ranging within normal values because they are able to compensate the lower W ext and higher recovery during prosthetic limb step with greater W ext and lower recovery during normal limb step [START_REF] Tesio | The 3-D motion of the centre of gravity of the human body during level walking II Lower limb amputees[END_REF]). However, with our study design, we cannot determine the potential mechanisms involved in this possible compensation. Thus, future studies, which would include kinematic and EMG data in addition to the parameters already assessed in the present study, are needed to further investigate these potential compensating mechanisms in patients with cLBP.

The lower preferred walking speed in patients with cLBP, and consequently the shorter stride length and lower stride frequency at preferred walking speed, is consistent with previous studies [START_REF] Lamoth | Effects of attention on the control of locomotion in individuals with chronic low back pain[END_REF][START_REF] Simmonds | The influence of pain distribution on walking velocity and horizontal ground reaction forces in patients with low back pain[END_REF][START_REF] Lee | Influence of pain distribution on gait characteristics in patients with low back pain: part 1: vertical ground reaction force[END_REF]). These changes are similar to those induced by healthy ageing after age 50 [for review see [START_REF] Beijersbergen | The biomechanical mechanism of how strength and power training improves walking speed in old adults remains unknown[END_REF]]. After this age, PWS is an indicator of general physical health and predicts daily function, independent living and many other clinical conditions [START_REF] Beijersbergen | The biomechanical mechanism of how strength and power training improves walking speed in old adults remains unknown[END_REF]. These earlier gait modifications in our younger patients may thus limit the ability to perform everyday tasks and increase sedentary behaviour in this population. Several hypotheses can be put forward to explain lower PWS in patients with cLBP. First, it may result from a strategy aiming to limit pelvic and thoracic rotations and decrease the expected loads on spinal tissues. In healthy individuals, pelvic and thoracic transverse rotations typically evolve from in-phase (pelvis and thorax move in the same direction) towards antiphase rotations (pelvis and thorax move in opposite directions) as walking speed increases [START_REF] Lamoth | Pelvis-thorax coordination in the transverse plane during walking in persons with nonspecific low back pain[END_REF]. Patients with cLBP exhibit less antiphase coordination [START_REF] Lamoth | Pelvis-thorax coordination in the transverse plane during walking in persons with nonspecific low back pain[END_REF]2006b), which may explain why they prefer lower walking speeds. Second, a lower walking speed may be used to cope with postural instability. Given the decreased ability of patients with cLBP to provide adequate spinal stability [START_REF] Panjabi | Clinical spinal instability and low back pain[END_REF], walking slowly may allow them to adapt to unexpected perturbations during walking with a greater margin of safety [START_REF] Taylor | The effect of walking faster on people with acute low back pain[END_REF]. Besides, at PWS, the medio-lateral displacements of COM were higher in cLBP than in control group (+20 %). This is indirect evidence that cLBP patient takes wider steps as compared to control individuals and may represent an active strategy to increase dynamic balance during walking as previously shown in pathological gait [START_REF] Kuo | Dynamic principles of gait and their clinical implications[END_REF]. Third, patients may spontaneously adopt a decreased speed to set intensity at a lower percentage of their maximal aerobic capacity and thus decrease the physiological relative effort at PWS. However, this would be an inefficient strategy since the energetic cost of walking was minimal at a speed (OWS) that was higher than their preferred walking speed. Finally, fear of pain is an important factor in the development and persistence of cLBP [START_REF] Vlaeyen | Fear-avoidance model of chronic musculoskeletal pain: 12 years on[END_REF]. The belief that pain would be exacerbated if walking speed gets higher may have rendered patients comfortable at a lower speed. In a preliminary analysis of the present data, there was no significant correlation between PWS and the TSK score (Pearson r = -0.136; P = 0.66). This association was also not significant in the study of Lamoth et al. (2006b). However, the number of participants in both studies precludes from drawing firm conclusions.

Patients with cLBP did not report higher levels of pain at fast compared to low speeds. Similarly, Lamoth et al. (2006b) did not find a significant increase in pain intensity before and after a walking protocol. A decrease in pain intensity was even reported by other authors in patients with recurrent LBP [START_REF] Lee | Influence of pain distribution on gait characteristics in patients with low back pain: part 1: vertical ground reaction force[END_REF]). In the present study, whereas pain VAS remained around 30 at all walking speeds, the task was perceived as more difficult by patients with cLBP than controls, as indicated by a higher level of perceived exertion, but this difference was not intensified with increasing speed. In a previous study, a back fatiguing task was subjectively perceived as more strenuous by patients with cLBP than healthy controls, whereas objective EMG recordings showed no significant difference between groups [START_REF] Lariviere | Specificity of a back muscle exercise machine in healthy and low back pain subjects[END_REF]. Ratings of perceived exertion may be influenced by factors such as pain and fear of movement [START_REF] Barker | Perceived and measured levels of exertion of patients with chronic back pain exercising in a hydrotherapy pool[END_REF][START_REF] Wallbom | Concordance between rating of perceived exertion and function in persons with chronic, disabling back pain[END_REF].

A major strength of the present study is that matching between patients and controls prevented from a confounding effect of age, body mass and height on energetics and mechanics of walking. Some limitations must, however, be mentioned. First, only males were included to limit the required sample size and to avoid an effect of sex on the results. It has been shown that the energy cost of walking is 10 % higher in females than males [START_REF] Browning | Effects of obesity and sex on the energetic cost and preferred speed of walking[END_REF]. The present findings must therefore be interpreted with caution in females. However, considering the very small differences observed between groups in mechanical work, spatiotemporal parameters and energy cost of walking, and under the assumption that sex accounts for a different energy cost of walking in both patients with chronic LBP and healthy individuals, the probability that different conclusions would be drawn in females appears to be small. Second, a wide range of walking speeds was required to model the U-shaped relationship between gross C w and walking speed and to calculate OWS, which precluded from recruiting highly disabled patients. Although the study sample is representative of the target population in terms of duration of symptoms, pain severity and fear of movement, the results cannot be generalized to highly disabled patients, which represent 15-20 % of patients with cLBP [START_REF] Payares | Validation in Colombia of the Oswestry disability questionnaire in patients with low back pain[END_REF][START_REF] Osthus | Cross-cultural adaptation of a German version of the Oswestry Disability Index and evaluation of its measurement properties[END_REF]. Third, although inertial sensors are sufficiently accurate to measure the mechanics and the spatiotemporal parameters of walking [START_REF] Meichtry | Criterion validity of 3D trunk accelerations to assess external work and power in able-bodied gait[END_REF][START_REF] Peyrot | Do mechanical gait parameters explain the higher metabolic cost of walking in obese adolescents?[END_REF]), approximations may have arisen from skin movements or changes in the relative positions of the sensor and COM over walking time. Although these factors cause an overestimation of the mechanical work and power [START_REF] Meichtry | Criterion validity of 3D trunk accelerations to assess external work and power in able-bodied gait[END_REF], such errors cannot have substantially affected the comparisons between patients with chronic LBP and controls because the same methodology was used in both groups.

In conclusion, patients with cLBP have a lower preferred walking speed than matched control participants, but no significant difference in mechanical work and energy cost of walking was observed between groups. These findings do not argue in favour of prescribing training modalities aimed at improving gait efficiency. At the same time, an increase in energy cost of walking does not appear to be a factor that discourages patients from engaging in physical activity. The moderate levels of perceived exertion, and the absence of pain exacerbation over the study protocol, imply that high walking speeds are well tolerated by patients with moderately disabling cLBP, and may be used effectively to increase the preferred walking speed, and to improve aerobic capacity in deconditioned patients.

Fig. 1

 1 Fig. 1 Energy cost of walking for the two groups at standard walking speeds. a gross energy cost of walking (gross C W ); b net energy cost of walking (net C W ). cLBP chronic low back pain, CTR control, filled square preferred walking speed (PWS) for cLBP group, hollow square PWS for CTR group, filled triangle optimal walking speed (OWS) for cLBP group, hollow triangle OWS for CTR group. Values are means (standard deviation)

Fig. 2

 2 Fig. 2 Spatiotemporal parameters for the two groups at standard walking speeds. a stride frequency, b stride length, c medio-lateral displacement. cLBP chronic low back pain, CTR control. Values are means (standard deviation)

Fig. 3 Fig. 4

 34 Fig. 3 Mechanical works and recovery for the two groups at standard walking speeds. a total mechanical work (W tot ), b mechanical external work (W ext ); c mechanical internal work (W int ); d recovery. cLBP chronic low back pain, CTR control. Values are means (standard deviation)

Table 1

 1 Participant characteristics

	Mean values (standard deviation)
	cLBP chronic low back pain group and CTR control group, n number
	of participants, BMI body mass index, IPAQ International Physical
	Activity Questionnaire
	* Significant difference between cLBP and CTR groups (P < 0.05)

Table 2

 2 Energetic and mechanical variables at preferred walking speed (PWS) for the two groups

	Variable	cLBP (n = 13)	CTR (n = 13)
	PWS (m s -1 )	1.17 (0.13)*	1.33 (0.11)
	Energy cost
	Gross C w (J kg -1 m -1 )	3.22 (0.28)	3.16 (0.36)
	Net C w (J kg -1 m -1 )	2.06 (0.25)	2.07 (0.21)
	Spatiotemporal parameters
	Stride frequency (Hz)	0.97 (0.06)*	1.29 (0.15)
	Stride length (m)	0.90 (0.04)*	1.06 (0.15)
	Medio-lateral displacement (m) 0.05 (0.01)*	0.04 (0.01)
	Single support duration (%)	78.2 (2.8)	79.8 (2.5)
	Double support duration (%)	21.8 (2.8)	20.2 (2.5)
	Mechanical works and recovery
	W tot (J kg -1 m -1 )	0.69 (0.12)*	0.82 (0.11)
	W ext (J kg -1 m -1 )	0.40 (0.10)*	0.48 (0.09)
	W int (J kg -1 m -1 )	0.29 (0.03)*	0.34 (0.03)
	Recovery (%)	64.9 (8.1)	59.2 (6.2)
	Heart rate (bpm)	84.7 (11.1)	86.5 (10.3)
	Perceived exertion, score	10.6 (1.9)	9.5 (1.4)

Mean values (standard deviation) cLBP chronic low back pain group and CTR control group, n number of participants, Gross C w gross energy cost of walking, Net C w net energy cost of walking, W tot total mechanical work, W ext mechanical external work, W int mechanical internal work, bpm beats per minute * Significant difference between cLBP and CTR groups (P < 0.05)

Acknowledgments Financial compensation to participants was funded by the Service of Rheumatology. No other funding was received for this work.

Compliance with the ethical standards

Conflicts of interest

The authors have nothing to disclose.


Jonathan Viquerat

Stéphane Lanteri

Simulation of near-field plasmonic interactions with a local approximation order discontinuous Galerkin time-domain method

Keywords: computational electromagnetics, nanophotonics, plasmonics, time-domain Maxwell equations, discontinuous Galerkin method, local adaptivity

During the last ten years, the discontinuous Galerkin time-domain (DGTD) method has progressively emerged as a viable alternative to well established finite-di↵erence time-domain (FDTD) and finite-element time-domain (FETD) methods for the numerical simulation of electromagnetic wave propagation problems in the time-domain. The method is now actively studied for various application contexts including those requiring to model light/matter interactions on the nanoscale. In this paper we further demonstrate the capabilities of the method for the simulation of near-field plasmonic interactions by considering more particularly the possibility of combining the use of a locally refined conforming tetrahedral mesh with a local adaptation of the approximation order.

Introduction

Generalities about the DGTD method

The DGTD method can be considered as a finite element method where the continuity constraint at an element interface is released. While it keeps almost all the advantages of the finite element method (large spectrum of applications, complex geometries, etc.), the DGTD method has other nice properties which explain the renewed interest it gains in various domains in scientific computing:

-It is naturally adapted to a high order approximation of the unknown field. Moreover, one may increase the degree of the approximation in the whole mesh as easily as for spectral methods but, with a DGTD method, this can also be done locally i.e. at the mesh cell level. In most cases, the approximation relies on a polynomial interpolation method but the method also o↵ers the flexibility of applying local approximation strategies that best fit to the intrinsic features of the modeled physical phenomena.

-When the discretization in space is coupled to an explicit time integration method, the DG method leads to a block diagonal mass matrix independently of the form of the local approximation (e.g the type of polynomial interpolation). This is a striking di↵erence with classical, continuous FETD formulations. Moreover, the mass matrix is diagonal if an orthogonal basis is chosen.

Email address: Stephane.Lanteri@inria.fr (Stéphane Lanteri)

-It easily handles complex meshes. The grid may be a classical conforming finite element mesh, a nonconforming one or even a hybrid mesh made of various elements (tetrahedra, prisms, hexahedra, etc.).

The DGTD method has been proven to work well with highly locally refined meshes. This property makes the DGTD method more suitable to the design of a hp-adaptive solution strategy (i.e. where the characteristic mesh size h and the interpolation degree p changes locally wherever it is needed).

-It is flexible with regards to the choice of the time stepping scheme. One may combine the discontinuous Galerkin spatial discretization with any global or local explicit time integration scheme, or even implicit, provided the resulting scheme is stable.

-It is naturally adapted to parallel computing. As long as an explicit time integration scheme is used, the DGTD method is easily parallelized. Moreover, the compact nature of method is in favor of high computation to communication ratio especially when the interpolation order is increased.

As in a classical finite element framework, a discontinuous Galerkin formulation relies on a weak form of the continuous problem at hand. However, due to the discontinuity of the global approximation, this variational formulation has to be defined at the element level. Then, a degree of freedom in the design of a discontinuous Galerkin scheme stems from the approximation of the boundary integral term resulting from the application of an integration by parts to the element-wise variational form. In the spirit of finite volume methods, the approximation of this boundary integral term calls for a numerical flux function which can be based on either a centered scheme or an upwind scheme, or a blend of these two schemes.

DGTD methods for nanophotonics/plasmonics

Numerical modeling of electromagnetic wave propagation in interaction with metallic nanostructures at optical frequencies requires to solve the system of Maxwell equations coupled to appropriate models of physical dispersion in the metal. In general, the Drude and Drude-Lorentz models are adopted although there are practical situations for which these models can fail to describe correcly the behavior of some materials (e.g. transition metals [START_REF] Ulbricht | Carrier dynamics in semiconductors studied with time-resolved terahertz spectroscopy[END_REF]- [START_REF] Wol↵ | Simple magnetooptic transition metal models for timedomain simulations[END_REF] and graphene [START_REF] De Abajo | Graphene nanophotonics[END_REF]). Furthermore at some scales, non-local e↵ects starts to play an important role [START_REF] Moreau | Impact of nonlocal response on metallodielectric multilayers and optical patch antennas[END_REF]. The FDTD [START_REF] Taflove | Computational electrodynamics: the finite-di↵erence time-domain method -3rd[END_REF] method is a widely used approach for solving the systems of partial di↵erential equations modeling nanophotonic applications. In this method, the whole computational domain is discretized using a structured (cartesian) grid. In spite of its flexibility and second-order accuracy in a homogeneous medium, the Yee scheme su↵ers from serious accuracy degradation when used to model curved objects or when treating material interfaces. Indeed, the socalled stair-casing approximation may lead to local zerothorder and at most first-order accuracy; it may also produce locally non-convergent results [START_REF] Ditkowski | Convergent cartesian grid methods for Maxwell's equations in complex geometries[END_REF]. Furthermore, for Maxwell's equations with discontinuous coe cients, the Yee scheme might not be able to capture the possible discontinuity of the solution across the interfaces [START_REF] Ditkowski | Convergent cartesian grid methods for Maxwell's equations in complex geometries[END_REF].

Thus, with all its above-listed features, the DGTD method seems to be well suited to the numerical simulation of complex time-domain ectromagntic wave propagation problems. As a matter of fact, the DGTD method for solving the time domain Maxwell equations is increasingly adopted by several physics communities. Concerning nanophotonics, unstructured mesh based DGTD methods have been developed and have demonstrated their potentialities for being considered as viable alternatives to the FDTD method. The most remarkable achievements in the nanophotonics domain since 2009 are due to Busch et al. Busch [START_REF] Niegemann | Higher-order time-domain methods for the analysis of nano-photonic systems[END_REF]- [START_REF] Stannigel | Discontinuous Galerkin time-domain computations of metallic nanostructures[END_REF]- [START_REF] Busch | Discontinuous Galerkin methods in nanophotonics[END_REF] has been at the origin of seminal works on the development and application of the DGTD method in this domain. These works not only deal with the extension of the DGTD method with regards to the complex material models and source settings required by applications relevant to nanophotonics and plasmonics [START_REF] König | The Discontinuous Galerkin Time-Domain method for Maxwells equations with anisotropic materials[END_REF]- [START_REF] Matyssek | Computing electron energy loss spectra with the Discontinuous Galerkin Time-Domain method[END_REF]- [START_REF] Wol↵ | Simple magnetooptic transition metal models for timedomain simulations[END_REF], but also to core contributions aiming at improving the accuracy and the e ciency of the proposed DGTD solvers [START_REF] Niegemann | Using curved elements in the discontinuous Galerkin time-domain approach[END_REF]- [START_REF] Diehl | Comparison of low-storage Runge-Kutta schemes for discontinuous Galerkin time-domain simulations of Maxwell's equations[END_REF]- [START_REF] Niegemann | E cient lowstorage Runge-Kutta schemes with optimized stability regions[END_REF]- [START_REF] Demirel | E cient multiple time-stepping algorithms of higher order[END_REF].

Noteworthy, all these studies adopt a di↵usive DGTD formulation based on upwind numerical fluxes. We have recently adapted the non-dissipative DGTD method initially introduced in [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF] to deal with various local dispersion models for metallic nanostructures. An ADE formulation has been adopted. The resulting ADE-based DGTD method is detailed in [START_REF] Viquerat | Theoretical and numerical analysis of local dispersion models coupled to a discontinuous galerkin time-domain method for Maxwell's equations[END_REF] where we also study the stability and a priori convergence of the method. We first considered the case of Drude and Drude-Lorentz models and, further extend the proposed ADE-based DGTD method to be able to deal with a generalized dispersion model in which we make use of a Padé approximant to fit an experimental permittivity function. The numerical treatment of such a generalized dispersion model is also presented in [START_REF] Viquerat | Theoretical and numerical analysis of local dispersion models coupled to a discontinuous galerkin time-domain method for Maxwell's equations[END_REF]. The possibility of using non-conforming hybrid structured/unstructured (cubic/tetrahedral) meshes in the nanophotonics context has been considered in [START_REF] Léger | A parallel non-conforming multi-element DGTD method for the simulation of electromagnetic wave interaction with metallic nanoparticles[END_REF]. More recently, a high order di↵usive DGTD method formulated on curvilinear tetrahedral meshes has been studied in [START_REF] Viquerat | A 3D curvilinear discontinuous Galerkin time-domain solver for nanoscale light-matter interactions[END_REF].

Objectives of the present study

In most of the existing works on the development of high order DGTD methods for the numerical modeling of light/matter interactions on the nanoscale, the formulation of the method is derived assuming a uniform distribution of the polynomial order to the cells of the underlying mesh. However, in the case of a mesh showing large variations in cell size, the time step imposed by the smallest cells can be a serious hindrance when trying to exploit high approximation orders. Indeed, a potentially large part of the CPU time is spent in the update of the physical field inside small cells where high polynomial orders might not be essential, while they are necessary in the larger cells.

To overcome this limitation, several strategies can be considered. The first one consists in replacing the explicit scheme by an implicit scheme. This, however, is at the expense of the resolution of a large linear system at each time step [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF]. In [START_REF] Moya | Locally implicit time integration strategies in a discontinuous Galerkin method for Maxwell's equations[END_REF], the authors consider a hybrid formulation, where only the smallest cells are treated via an implicit scheme, while keeping an explicit time integration for the rest of the tesselation, thus limiting the time step constraint. Although this strategy provides very interesting results in terms of CPU speedup, maintaining high order time integration is di cult. Another possiblity is to exploit local time stepping combined to an explicit scheme: based on the size of the elements, the mesh is divided in regions, each of which being assigned an appropriate timestep for an explicit time integration (see [START_REF] Piperno | Symplectic local time stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF]- [START_REF] Diaz | Energy conserving explicit local time stepping for second-order wave equations[END_REF]- [START_REF] Grote | Explicit local time-stepping methods for Maxwells equations[END_REF]- [START_REF] Grote | Runge-Kutta-based explicit local time-stepping methods for wave propagation[END_REF]- [START_REF] Demirel | E cient multiple time-stepping algorithms of higher order[END_REF]). As shown for example in [START_REF] Diaz | Energy conserving explicit local time stepping for second-order wave equations[END_REF], high order convergence in time can be preserved with such explicit local time stepping strategies. However, it seems di cult to ensure a good load balance in the case of a parallel implementation, given that the natural repartition is based on a time step criterion, instead of a parallel-related one.

A complementary strategy which is considered here relies on the use of non-uniform distribution of the polynomial order in the framework of a global time step DGDT method. By imposing low orders in small cells and high orders in large cells, it is possible to significantly alleviate both the global number of degrees of freedom and the time step restriction with a minimal impact on the method accuracy. Strategies exploiting locally adaptive (LA) formulations usually combine both hand p-adaptivity in or-der to concentrate the computational e↵ort in the areas of high field variations. Here, the adopted point of view is quite di↵erent: starting from a given mesh and an uniform distribution of the polynomial order P k , the LA strategy exploits all the polynomial orders P p with p  k to obtain a solution of similar accuracy with a reduced computational cost.

In section 2, a DGTD formulation is presented that account for a variable polynomial order. In section 3, the resulting DGTD is tested for convergence, and an extended performance study is provided, both for sequential and parallel executions. Finally, in section 6, the CPU gains are evaluated on two configurations relevant to nanophotonics.

DGTD formulation

Problem statement

For the sake of simplicity, we consider the non-dispersive time-domain Maxwell equations for this presentation. The extension of the formulation presented below to the case of a local dispersion model is straightforward. The continuous, normalized Maxwell system is

µ r @H @t = r ⇥ E, " r @E @t = r ⇥ H J. (1)

Notations

Let ⌦ ⇢ R 3 be a bounded convex domain, and n the unitary outward normal to its boundary @⌦. Let ⌦ h be a discretization of ⌦, relying on a quasi-uniform triangulation T h verifying T h = S N i=1 T i , where N 2 N ⇤ is the number of mesh elements, and (T i ) i2J1,N K the set of simplices. The internal faces of the discretization are denoted a ik = T i T T k if T i and T k are adjacent cells, and n ik is defined as the unit normal vector to the face a ik , oriented from

T i toward T k . For each cell T i , V i is the set of in- dices {k 2 J1, NK | T i T T k 6 = ;}. We define the following approximation space V h V h = n v 2 L 2 (⌦ h ) 3 , v Ti 2 (P p (T i )) 3 8 T i 2 T h o ,
where P p (T i ) is the space of polynomials of maximum degree p on T i . The semi-discrete fields, seeked in space V h , are hereafter denoted (H h , E h , J h ), and on each cell

T i the restrictions (H i , E i , J i ) = (H h Ti , E h Ti , J h

One now seeks the approximations E h and H h of E and H in space V

h . The contribution of each cell is therefore defined as

E i = E h Ti .
Here, one must notice that, for a 3D system, E i is actually a vector that has 3 components

E i = 2 4 E x i E y i E z i 3 5 ,
each of which is locally expanded on the chosen set of basis functions

E v i = di X j=1 E v ij ij , v 2 {x, y, z} . (2)
Therefore, for practical purpose, one defines three vectors of d i components

E v i = 2 6 4 E v i1 . . . E v idi 3 7 5 , v 2 {x, y, z} ,
as well as the following 3d i components vector

E i = 2 6 4 E x ij 16j6di E y ij 16j6di E z ij 16j6di 3 7 5 .

Discretization in space

We consider the system of time-domain Maxwell equations [START_REF] Ulbricht | Carrier dynamics in semiconductors studied with time-resolved terahertz spectroscopy[END_REF] with no current source term. The starting point is the following weak formulation

Z Ti µ r @H i @t • v ik + Z Ti E i • r ⇥ v ik = X l2Vi Z a il (E ⇤ ⇥ n il ) • v ik , Z Ti " r @E i @t • v ik Z Ti H i • r ⇥ v ik = X l2Vi Z a il (H ⇤ ⇥ n il ) • v ik .
We recall that the i subscript refers to the cell T i of the underlying tetrahedral mesh where a P p polynomial approximation of the field components is used. Exploiting the local field expansions from (2), one can cast the first volumic integral of the above system as the following matrix-vector products

Z Ti " r @E i @t • x ik = ✓ M "r i @ E x i @t ◆ k ,
where M "r i is the mass matrix, of dimension

d i ⇥ d i (M "r i ) jk = Z Ti " r ij ik , with (j, k) 2 J1, d i K 2 .

3

In a similar fashion, the curl integral can be cast as

Z Ti H i • r ⇥ x ik = K i ⇥ H i x k .
Here, the three sti↵ness matrices were introduced

(K v i ) jk = Z Ti ij @ ik @v for v 2 {x, y, z} , with (j, k) 2 J1, d i K 2 .
From the latter definition, we define the general 3d i ⇥ d i sti↵ness matrix that will be used in the final system

K i = 2 4 K x i K y i K z i 3 5 .
The flux integral, however, must undergo a di↵erent treatment. For the sake of simplicity, the centered flux is considered, the generalization to other fluxes being straightforward. Consider a neighbor cell T l of T i with a P m polynomial expansion. The flux integral on their common face from the T i side for the x component is

Z a il (H ⇤ ⇥ n il ) • x ik = Z a il (H y ⇤ n z il H z ⇤ n y il ) ik = Z a il ✓ H y i + H y l 2 n z il H z i + H z l 2 n y il ◆ ik . (3)
Consider the following expansions for H y i and H y l , and their counterparts for H z i and H z l

H y i = p X q=1 H y iq iq and H y l = m X r=1 H y lr lr . (4)
Plugging ( 4) in (3) leads to

Z a il (H ⇤ ⇥ n il ) • x ik = Z a il n z il 2 p X q=1 H y iq iq ik + m X r=1 H y lr lr ik ! Z a il n y il 2 p X q=1 H z iq iq ik + m X r=1 H z lr lr ik ! = 1 2 p X q=1 H y iq n z il H z iq n y il Z a il iq ik + 1 2 m X r=1 (H y lr n z il H z lr n y il ) Z a il lr ik = 1 2 p X q=1 H iq ⇥ n il Z a il iq ik + 1 2 m X r=1 H lr ⇥ n il Z a il lr ik = S il H ⇤,i ⇥ n il x k + S ⇤ il H ⇤,l ⇥ n il x k . (5)
We introduce the flux or interface matrices

(S il ) jk = Z a il ij ik and (S ⇤ il ) rk = Z a il lr ik . (6)
The derivation (3) easily extends to the other components, as well as other flux choices. We recall one of the most spread flux formulation

E ⇤ = 1 Y i + Y l ({Y E} il + ↵n ⇥ JHK il ) , H ⇤ = 1 Z i + Z l ({ZH} il ↵n ⇥ JEK il ) , (7)
where

{A} il = A i + A l is twice the mean value of A across the interface, JAK il = A l A
i is the jump of A across the interface, and Y and Z are respectively the inductance and impedance of the considered materials. ↵ 2 [0, 1] is a tunable parameter that allows to vary between the centered flux for ↵ = 0, to a fully upwind flux for ↵ = 1.

To summarize, the flux integral is cut in two parts: (i) the part corresponding to local information, which is integrated via the regular flux matrix S il , and (ii) the part corresponding to the neighbor information, which is integrated via the S ⇤ il matrix. In the case of a conformous interface (i.e. p = m), S ⇤ il = S il . However, for a

P p P m interface, S ⇤ il is a non-conforming rectangular matrix of size s(p) ⇥ s(m) (we recall that s(p) = (p+1)(p+2) 2
). From this point, the remaining of the derivation is similar to the standard case (see for example [START_REF] Viquerat | A 3D curvilinear discontinuous Galerkin time-domain solver for nanoscale light-matter interactions[END_REF]). The final semidiscrete scheme is

M µr i @ H i @t = K i ⇥ E i + X l2Vi S il E ⇤,i ⇥ n il + X l2Vi S ⇤ il E ⇤,l ⇥ n il , M "r i @ E i @t = K i ⇥ H i X l2Vi S il H ⇤,i ⇥ n il X l2Vi S ⇤ il H ⇤,l ⇥ n il , (8)
with the following definitions

(S il ) jk = Z a il ij ik (S ⇤ il ) rk = Z a il lr ik . (9)

Time integration

In this work, we use the optimized 14 stages fourthorder low-storage Runge-Kutta scheme developped by Niegemann et al. [START_REF] Niegemann | E cient lowstorage Runge-Kutta schemes with optimized stability regions[END_REF], hereafter designated as LSRK4 scheme. Its stability region is optimized for the upwind DG eigenspectrum, and is thus an e cient time-integration scheme for our needs. The reader can refer to the aforementionned reference for more details.

Numerical study of the h-convergence

In this section, we perform a preliminary validation of the LA-DGTD (Locally Adaptive DGTD) method by studying its numerical convergence. For this purpose, we consider the test problem of the propagation of an eigenmode in a cubic PEC cavity. Four meshes of increasing refinement, M1 to M4, are used. Unless stated otherwise, a fully upwind flux is used, coupled to the LSRK4 time scheme. In the present case, the used meshes are uniform, and the mesh cells all have the same size. To begin, we separate the computational domain into two vertical stripes of equal sizes, with di↵erent values of the polynomial order (see figure 1(a)). The obtained results are given in table 1. As expected for a P k P l configuration, the asymptotic h-convergence order is min(k, l) + 1.

P l P k (a) Test configuration for h-convergence P k P l P k P l (b)
Influence of the number of non-conformeous interfaces

P 3 P 4 P 1 P 2
(c) Influence of the type of non-conformeous interfaces Figure 1: Order distribution for the h-convergence validation. To test convergence, the domain is arbitrarily cut in two halves, each part receiving a di↵erent order (figure 1(a)). In order to assess the impact of non-conformeous interfaces, the domain is cut in four equal stripes, thus doubling the number of P k P l faces while keeping the same number of tetrahedron per order (figure 1(a)). Another test is conducted by setting a di↵erent order in each quarter (figure 1(c)).

In order to evaluate the impact of non-conforming interfaces, the computational domain is now distributed into four equal size stripes (see figure 1(b)). Hence, the total number of tetrahedra for each value of the interpolation order remains the same as in figure 1(a), but the number of interfaces is twice as high. As can be seen in table 1, for the coarsest meshes the higher amount of non-conforming faces yields a slightly higher (but still acceptable) L 2 error. However, this error overhead vanishes for refined meshes. In a last numerical test, another four stripes division of the computational domain is used, where each part receives a di↵erent order (see figure 1(c)). As in the previous cases, and since all the cells are of same size, the numerical error is driven by the lowest approximation order.

Order distribution strategy

Starting with a given unstructured mesh, it is clear that the distribution of the interpolation order to the mesh cells will have a major impact on the obtained accuracy, as well as on the time required to obtain the numerical solution. Let us assume that the solution is obtained on the given mesh with a homogeneous polynomial order P k . The point is here to see how, with a good distribution of polynomial orders P l with l  k, a solution of similar accuracy can be obtained at a lower computational e↵ort. At first glance, it seems that configurations including small geometrical details, or small gaps between two structures, could benefit from such a strategy. For this reason, for any given mesh, we define the following quantity

h = h max h min ,
which represents the heterogeneity in terms of cell size of the mesh. In the following, t i represents the time step corresponding to the cell T i , computed following the formula ..., while t p i represents the e↵ective time step obtained if cell T i is discretized with a P p polynomial expansion t p i = CFL(p) t i , where t i is the maximal acceptable timestep for cell T i , computed as

t i = c p V Ti A Ti , (10)
where V Ti and A Ti are respectively the volume and the area of cell T i , and c p is an order-dependant constant. The normalized time step includes a rough estimate of the computational charge induced by the polynomial order, and is defined as

t p i = t p i n(p) = CFL(p) t i n(p) ,
where n(p) is the number of degrees of freedom inside a P p cell. Finally, we define p min and p max the minimal and maximal user-authorized orders in the mesh. We also add the following constraint: non-conforming faces cannot connect cells with an order jump higher than one (the allowed configurations are presented on figure 2). Indeed, it is preferable to restrain the number and size of matrices in memory in order to improve data locality. Additionaly, it leads to a robust distribution strategy, as will be shown hereafter.

The first step of the algorithm consists in computing the local time steps t i and sorting them by ascending order. The cell with the smaller t receives order p min and we compute its normalized time step t pmin 1 . Two temporary variables, p loc and t loc , respectively store the current order assigned to the cells and the current restrictive normalized time step. For a given cell, the normalized time step for increased order p loc + 1 is compared to the current limiting normalized time step t loc . In the case where the first is higher than the second, switching to the higher order is assumed to have a limited impact on the final time step. Hence, p loc is increased by one, and t loc is updated. The procedure is summarized in algorithm 1, whose performances are assessed in next section.

Order distribution strategy

Single-core speedup

To evaluate the gains provided by the LA-DGDT strategy, we consider the three meshes M1, M2 and M3 shown on figure 3. These meshes are obtained by the tesselation of a cubic PEC cavity corresponding to the test problem considered previously, a local refinement being imposed on one side of the box. The characteristics of these meshes are summarized in table 2. For each mesh, the cavity mode is computed sequentially for 5 periods. As before, CPU time and maximum L 2 error are stored. The results obtained for homogeneous and mixed orders are presented in table 3. First, we note that the memory occupation has not been optimized in the present implementation of the LA-DGTD strategy. Hence, we can expect that mixed order computations will require the same memory occupation than homogeneous order ones. For mixed order solutions, the speedup of a P k P l computation is obtained by comparing its CPU time with that of a full P max(k,l) computation. For the three considered meshes, mixing two polynomial orders leads to speedups ranging from 1.5 to 2.2. One can remark that there exists a relation between r and the obtained speedup. An interesting point is that the obtained L 2 -errors are less than 1% higher than those of the homogeneous polynomial order computations. Mixing three orders does not provide any improvement for the M1 mesh, i.e. the distribution algorithm did not map the highest polynomial order to a cell of the mesh. This can be easily understood by looking at the compared t distribution of the three meshes (see figure 4). Since the algorithm imposes the lowest order for the cell of smallest t, a certain range for the t is required to exploit the highest polynomial orders. For example, the P 1 P 4 distribution is shown on the same figure for mesh M3. For meshes M2 and M3, however, very interesting gains are obtained, with speedups ranging from 3 to 4.5. In this case, it seems that a higher r implies a higher benefit from the LA-DGTD strategy based on the distribution algorithm 1. Finally, mixing polynomial orders 1 to 4 roughly provides a speedup of 6, while increasing the global L 2 error by less than 1%.

Parallel load balancing

Parallel computing is a mandatory path for reducing the cost of simulations of complex electromagnetic wave propagation problems such as those pertaining to nanophotonics and plasmonics. In this section, we present the results obtained when trying to balance the computational load for a parallel implementation of the LA-DGTD method. The MeTiS [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] graph partitioning tool is used to split the computational domain in subdomains, each of which is then mapped on a core of parallel system. The communications between the cores are handled via the MPI standard. At the end of each computation, each CPU core returns its own CPU time, excluding the time spent in the MPI communication routines. Hence, a good load balance between the cores will manifest as nearly-identical CPU times for all cores. To reach this result, the partitioning algorithms of the MeTiS package can be provided with a weight w i for each cell T i . During the partitioning, this weight is taken into account so that the total weight of the various subdomains are as close as possible. Here, the following definition of the weight is used

w i = n(p i ) + X l2Vi max(s(p i ), s(p l )) (11)
This weighted partitioning is applied to the test problem considered in section 5.1 using the M3 mesh and running 100 time iterations. Given the small number of tetrahedra in the mesh, the study is limited to 4 and 8 subdomains partitions. For a simpler comparison between the two partitionings, the relative deviation (in %) to the mean value is computed for each core. First, the e↵ect of the weighting is assessed by comparing relative deviations obtained from weighted and non-weighted partitions. The results for P 1 P 4 approximation with 4 and 8 CPUs are shown on figure 5. As can be seen, the use of a weighting pattern is mandatory to preserve good parallel performances. For both 4 and 8 cores, applying the weighted partitionning results in a maximal relative imbalance lower than 5%. To further explore the behavior of the algorithm in parallel, the CPU load balances with 4 and 8 cores are re-plotted with matching scales. As expected, increasing the number of cores also increases the maximal imbalance between cores, though in reasonable bounds.

Plasmonic simulations

Plasmonic nanolens

This section presents the computation of the field enhancement obtained in a plamonic nanolens device. To overcome the limitation of the di↵raction limit, it is possible to exploit the focusing e↵ect provided by coupled surface plasmons [START_REF] Dai | Electrodynamics e↵ects in plasmonic nanolenses[END_REF]. A typical nanolens is composed of a chain of metallic nanoparticles (nanospheres being the most common) of decreasing size, aligned with the polarization direction of the incident field (see figure 6). When the nanospheres are of significantly di↵erent sizes, the local field enhancement of the first particle is not perturbed by the second one because of its small relative size. As a result, the locally enhanced field of the first particle acts as an incident field for the second particle, resulting in a second enhancement, and so on. Eventually, the strongest enhancement is obtained in the gap between the two smaller particles [START_REF] Li | Self-similar chains of metal nanospheres as an e cient nanolens[END_REF].

Here, we consider a nanolens made of three gold spheres. The geometry is taken from [START_REF] Li | Self-similar chains of metal nanospheres as an e cient nanolens[END_REF]. The respective radii are 45, 15 and 4.5 nanometers, while the spacings between the sphere surfaces are respectively 4.5 and 1.5 nanometers. Gold is described by a Drude model of parameters " 1 = 3.7362, ! d = 1.387 ⇥ 10 7 GHz and d = 4.515 ⇥ 10 4 GHz. The Drude model is a basic model used to describe the response of free electrons in metal in the THz regime [START_REF] Drude | Zur elektronentheorie der metalle[END_REF]. The nanolens is illuminated via a total-field/scattered-field (TF/SF) interface, with a wideband plane wave of central frequency 700 THz, whose polarization is aligned with the natural axis of the lens (here, the x+ direction). The domain is terminated with a complex frequency-shifted perfectly-matched layer (CFS-PML) [START_REF] Bérenger | Perfectly matched layer for computational electromagnetics, 1st Edition[END_REF]. Finally, a probe point is set at half distance between each pair of spheres. At these positions, the discrete Fourier transform of the field is computed, and the field enhancement g = | Êx| | Êi| is deduced. To obtain a proper resolution of the geometry, very small elements must be used on the surface of the smallest sphere, as well as in the smallest gap, while the rest of the geometry (largest sphere, vacuum and PMLs) are meshed with much coarser elements (see mesh on figure 7). As a result, the r factor is here superior to 800.

To obtain convergence with an homogeneous order over the whole mesh, P 3 approximation is required. The corresponding numerical solution is obtained in 49 hours 48 minutes on 16 cores, and is taken as a reference. To exploit the LA-DGTD method, P 1 to P 3 approximations are used. The order distribution with respect to time step is shown on figure 8(a). To further understand the behavior of the repartition algorithm, a visual representation is added on (%) As expected, first order polynomials are assigned to the elements lying on the surface of the smallest sphere, while second order approximation is used in its vicinity and for the surface of the second sphere. The rest of the mesh is discretized with third order polynomials.

The computed field enhancements are presented on figure 9. As stated in the litterature, particularly intense fields are obtained between the two smallest sphere, where enhancements up to 700 are observed. The P 1 P 3 solution is obtained in 19 hours and 17 minutes on 16 cores, hence providing a speedup of 2.6 over the full P 3 solution. The observed error over the frequency range of interest is less than 1 %. Field maps of the E y component are presented on figure 10 for P 1 , P 1 P 3 and P 3 polynomial approximations. As expected, the field intensity between the spheres is underestimated for the full P 1 computation, while they are almost identical for P 1 P

3 and P 3 approximations.

Bowtie nanoantenna

As a second application, we consider here the computation of the extinction cross section of a metallic bowtie nanoantenna. These structures are actively studied for the very strong field enhancement they provide between the tips of the two triangular nanoparticles (see figure 11), which is known to be inversely proportionnal to the size of the gap. Hence, bowtie nanoantennas are good candidates for surface-enhanced Raman spectroscopy (SERS) applications [START_REF] Hatab | Free-standing optical gold bowtie nanoantenna with variable gap size for enhanced Raman spectroscopy[END_REF]. Recent advances in lithography techniques allowed the creation of structures with gaps as small as 3 nm [START_REF] Zhang | Bowtie nanoantenna with single-digit nanometer gap for surface-enhanced Raman scattering (SERS)[END_REF], while the typical size of the full structure can get close to 200 nm. Additionally, realistic geometries of such antennas include small roundings at the edges and tips, whose typical size is between a few to a few tens of nanometers [START_REF] Grzeskiewicz | Near and farfield properties of nanoprisms with rounded edges[END_REF].

In the present case, we consider a pair of 10 nm thick, equilateral prisms of edge length 100 nm, with a spacing gap of 3 nm. The rounding radius is 2 nm, and is uniform As can be seen on figure 12, the typical setup for such computations requires very small elements (geometrical details of the nanoantenna) as well as very large ones (vacuum and PML cells), and could therefore make good use of the LA-DGTD formulation presented before. To compute the exctinction cross-section, the bowtie is illuminated from above with a wide-band plane wave of central frequency 500 THz, polarized along the major axis of the antenna. With an homogeneous polynomial order on the whole mesh, convergence is obtained for a P 3 approximation, requiring 30 hours and 37 minutes on 16 cores. This is taken as a reference for the LA-DGTD strategy, for which P 1 to P 3 approximations are used. The repartition of orders with regards to the time-step is presented on figure 13, along with a visual representation of the order selection in the mesh. As expected, the first order is attributed to the small cells of the edges and tips, which are then enclosed into a second layer of P 2 elements. All the remaining cells (not represented) receive a third order interpolation.

The computed exctinction cross sections are presented on figure 14. A single very large resonance is observed around 418 THz. As can be seen, the P for the P 1 P 3 case. The white elements correspond to P 1 approximation, while the purple cells are second-order. The remaining cells (not represented) receive third than the full P 1 solution is also plotted. In terms of performance, the P 1 P 3 solution is obtained in 13 hours and 59 minutes on 16 cores, hence yielding a 2.2 speedup factor, which is lower than what was observed in section 6.1. The di↵erence can be attributed to the lower proportion of high-order (P 3 ) elements compared to low-order ones (less than 20 % of P 3 elements here, against 40 % for the nanolens case). However, this remains an appreciable gain for a solution of similar accuracy. As an illustration, a field map of |E| is plotted on figure 15 for P 1 P 3 approximation.

Conclusion

In this paper, the use of local polynomial approximation in the DGTD method is presented. The convergence of the algorithm is demonstrated on a standard PEC cavity case. Then, an order repartition algorithm is proposed that proved to be e cient, both for textbook and realistic nanophotonics-related cases. Although the obtained speedups are lower for realistic cases (between 2 and 2.6 for P 1 P 3 ) than for academic cases (up to 4.5 for P 1 P 3 and 6.15 for P 1 P 4 ), the LA-DGTD strategy represents an interesting gain in speed for long time simulations. However, the repartition algorithm being based on the time step, it also implicitely relies on a basic knowledge of the physical behavior of the computed system (the preliminary grasp Figure 14: Extinction cross-section of the bowtie nanoantenna obtained with P 1 , P 3 and P 1 P 3 approximations. Less than 2 % of relative error is observed between full P 3 and P 1 P 3 computations, for a speedup factor superior to 2.

of the positions of intense fields, basically). A good remedy would consist in a coupling of the algorithm with an a posteriori error estimate, in order to dynamically adapt the polynomial order in the mesh.

Ti

  ) are defined. A set of scalar basis functions ( ik ) 16k6di is defined for each T i , where d i is the number of degrees of freedom (d.o.f.) per dimension. Additionally, to each scalar basis function, the three vectors v ik are associated

Figure 2 :

 2 Figure 2: Authorized interfaces in the local order of approximation implementation. Order jumps are limited to one, yielding three types of interfaces for (p min , pmax) = (1, 4).

Figure 3 :

 3 Figure 3: Meshes for the cubic cavity mode with local refinement on one side.

  Order distribution for mesh M3Compared M1, M2 and M3 and order distribution for the P 1 P 4 case on M3 (right).

8

 8

Figure 5 :

 5 Figure5: Weighted against non-weighted parallel load balance on 4 and 8 cores with M3 mesh. Each bar corresponds to the relative imbalance of a single CPU to the average value computed over all processors. As can be seen, the weighting restores a good balance of the CPU loads, with no relative imbalance exceeding 5%.

Figure 6 :

 6 Figure 6: Nanolens composed of three metallic spheres.

Figure 7 :

 7 Figure7: Mesh setup for a metallic nanolens. The gray cells correspond to the metallic spheres, the blue cells to vacuum, while the red cells constitute the PML region. For this mesh, the ratio r is above 400.

Figure 8 :Figure 9 : 1 (b) P 1 P 3 (

 8913 Figure 8: Polynomial order repartition for the nanolens mesh with respect to time-step (top), and geometrical repartition (bottom)for the P 1 P 3 case. The red elements correspond to P 1 approximation, the green ones to P 2 , and the gray ones to P 3 .

Figure 10 :

 10 Figure 10: Ey field map in the nanolens device at t = 10 fs.For the three views, the field values are scaled to[START_REF] Demirel | E cient multiple time-stepping algorithms of higher order[END_REF][START_REF] Demirel | E cient multiple time-stepping algorithms of higher order[END_REF].

Figure 11 :

 11 Figure 11: Bowtie nanoantenna with rounded edges.

Figure 12 :

 12 Figure12: Mesh setup for a bowtie nanoantenna. The gray cells correspond to the nanoantenna, the blue cells to vacuum, while the red cells constitute the PML region. For this mesh, the ratio r is close to 275.

1 P 3 Figure 13 :

 1313 Figure 13: Polynomial order repartition for the bowtie mesh with respect to time-step (top), and geometrical repartition (bottom) for the P 1 P 3 case. The white elements correspond to P 1 approximation, while the purple cells are second-order. The remaining cells (not represented) receive third

Figure 15 :

 15 Figure 15: |E| field map in the nanolens device at t = 12.3 fs, obtained with a P 1 P 3 approximation. The field values are scaled to [0, 10].

Table 1 :

 1 Error levels and convergence rates of the cubic cavity case for mixed orders of approximation on meshes of increasing refinement. In the case of mixed orders P

						k	P	l , 1 refers to a domain
	cut in two halves (see figure 1(a)), and 2 to a domain cut in four
	stripes (see figure 1(b)). The case P 1 P 4 corresponds to a domain
	cut in four quarters, as depicted on figure 1(c). All simulations were
	run with upwind fluxes and LSRK4 time integration.
			M1		M2
			kE E	h k	r	kE E	h k	r
	P 1	-	2.87 ⇥ 10 1	-	6.05 ⇥ 10 2	2.25
	P 2	-	1.47 ⇥ 10 2	-	1.36 ⇥ 10 3	3.43
	P 3	-	9.24 ⇥ 10 4	-	5.87 ⇥ 10 5	3.98
	P 4	-	9.45 ⇥ 10 5	-	3.11 ⇥ 10 6	4.92
	P 1 P 2	1 2	2.65 ⇥ 10 1 2.07 ⇥ 10 1	--	3.38 ⇥ 10 2 3.55 ⇥ 10 2	2.97 2.54
	P 2 P 3	1 2	8.50 ⇥ 10 3 8.70 ⇥ 10 3	--	9.21 ⇥ 10 4 9.16 ⇥ 10 4	3.21 3.25
	P 3 P 4	1 2	6.73 ⇥ 10 4 6.81 ⇥ 10 4	--	4.41 ⇥ 10 5 4.47 ⇥ 10 5	3.93 3.93
	P 1 P 4	-	2.65 ⇥ 10 1	-	3.38 ⇥ 10 2	2.97
			M3		M4
			kE E	h k	r	kE E	h k	r
	P 1	-	8.66 ⇥ 10 3	2.80	1.46 ⇥ 10 3	2.57
	P 2	-	1.75 ⇥ 10 4	2.96	2.19 ⇥ 10 5	3.00
	P 3	-	3.72 ⇥ 10 6	3.98	2.33 ⇥ 10 7	4.00
	P 4	-	1.98 ⇥ 10 7	3.97	1.15 ⇥ 10 8	4.11
	P 1 P 2	1 2	6.15 ⇥ 10 3 6.30 ⇥ 10 3	2.46 2.49	1.42 ⇥ 10 3 1.43 ⇥ 10 3	2.12 2.14
	P 2 P 3	1 2	1.18 ⇥ 10 4 1.17 ⇥ 10 4	2.96 2.97	1.48 ⇥ 10 5 1.48 ⇥ 10 5	3.00 2.98
	P 3 P 4	1 2	2.80 ⇥ 10 6 2.85 ⇥ 10 6	3.98 3.97	1.76 ⇥ 10 7 1.79 ⇥ 10 7	4.00 3.99
	P 1 P 4	-	6.15 ⇥ 10 3	2.46	1.42 ⇥ 10 3	2.12

Table 2 :

 2 Characteristics of the locally refined cubic cavity meshes. ns is the number of vertices, nt the number of tetrahedrons and r is the ratio between the largest and the smallest cells in the mesh.

		M1	M2	M3
	ns	427	1429	11975
	nt	1513	5042	42652
	r	10.2	100.9	1000.8

Table 3 :

 3 CPU times, memory consumption and error levels for mixed orders of approximation on locally refined meshes. The order repartition was obtained via the procedure described in section 4. All simulations were run with upwind fluxes and LSRK4 time integration.

				M1	M2	M3
		CPU (s)	1.77	48.7	6 0 3 0
	P 1	Mem. (MB) kE E h k	14.1 3.38 ⇥ 10 2	2 2 .1 3.59 ⇥ 10 2	1 0 7 .3 3.87 ⇥ 10 2
				6.50	180	22460
	P 2			17.8 3.68 ⇥ 10 3	3 4 .4 3.71 ⇥ 10 3	2 1 1 3.82 ⇥ 10 3
				21.8	6 1 1	9 0 0 2 0
	P 3			23.9 2.14 ⇥ 10 4	5 4 .6 2.39 ⇥ 10 4	3 8 2 2.53 ⇥ 10 4
				73.8	2 1 0 6	2 2 8 2 2 0
	P 4			32.9 1.42 ⇥ 10 5	8 4 .7 1.76 ⇥ 10 5	6 3 5 1.97 ⇥ 10 5
		CPU (s)	3.84	101	15000
		Speedup	1.69	1.78	1.50
	P 1 P 2	Mem. (MB)	17.8	3 4 .4	2 1 1
		Tet. ratios	0.18, 0.82	0.29, 0.71	0.33, 0.67
		kE E	h k	3.84 ⇥ 10 3	3.71 ⇥ 10 3	3.82 ⇥ 10 3
				14.3	3 7 2	4 0 4 7 0
				1.52	1.64	2.22
	P 2 P 3			23.9	5 4 .6	3 8 2
				0.25, 0.75	0.38, 0.62	0.43, 0.57
				2.38 ⇥ 10 4	2.39 ⇥ 10 4	2.53 ⇥ 10 4
				49.2	1 3 9 0	1 3 0 7 3 0
				1.5	1 .51	1.74
	P 3 P 4			33.0	8 4 .8	6 3 5
				0.17, 0.83	0.28, 0.72	0.31, 0.69
				1.42 ⇥ 10 5	1.76 ⇥ 10 5	1.97 ⇥ 10 5
					180	19890
					3.39	4.53
	P 1 P 2 P 3				54.7	3 9 2
					0.29, 0.31, 0.40	0.33, 0.38, 0.29
					2.39 ⇥ 10 4	2.53 ⇥ 10 4
					695	64820
					3.03	3.52
	P 2 P 3 P 4				84.8	6 3 5
					0.38, 0.22, 0.40	0.43, 0.27, 0.30
					1.76 ⇥ 10 5	1.97 ⇥ 10 5
					347	37130
					6.07	6.15
	P 1 P 2 P 3 P 4				84.9	6 3 6
					0.29, 0.32, 0.13, 0.26	0.33, 0.38, 0.15, 0.14
					1.80 ⇥ 10 5	1.99 ⇥ 10 5

Acknowledgements. The authors gratefully acknowledge support from the Direction Générale de l'Armement (DGA) which partially supports the doctoral thesis of Jonathan Viquerat.