PhD Ann Hallemans
email: ann.hallemans@uantwerpen.be

Evi Verbecque

Raphael Dumas

Laurence Cheze

Angèle Van

Thomas Robert

Universiteitsplein Movant

Angèle Van Hamme

Developmental Changes In Spatial Margin Of Stability In Typically Developing Children Relate To The Mechanics Of Gait

Keywords:

 The spatial margin of stability decreases with increasing age  Development of dynamic balance control results in changes in step-time parameters  Changes in step width and single stance relate to the decrease in margin of stability  These features suggest improvements in dynamic balance control with increasing age

Introduction

At walking onset, between 10 and 15 months of age, infants walk at slow speed with steps shorter than their leg length [START_REF] Hallemans | 3D joint dynamics of walking in toddlers -A crosssectional study spanning the first rapid development phase of walking[END_REF][START_REF] Badaly | Beyond the average: walking infants take steps longer than their leg length[END_REF]. A high variability and low stability of gait is characteristic [START_REF] Bisi | Measures of gait stability: performance on adults and toddlers at the beginning of independent walking[END_REF]. Already a few months after the onset of walking, speed and step length dramatically increase [START_REF] Hallemans | 3D joint dynamics of walking in toddlers -A crosssectional study spanning the first rapid development phase of walking[END_REF][START_REF] Hallemans | Changes in 3D joint dynamics during the first 5 months after the onset of independent walking: A longitudinal follow-up study[END_REF] and occasionally steps larger than leg length can be observed [START_REF] Badaly | Beyond the average: walking infants take steps longer than their leg length[END_REF]. Other features of an immature gait pattern are a wide base of support and a prolonged duration of stance as well as increased double support [START_REF] Sutherland | The Development of Mature Gait[END_REF][START_REF] Sala | Gait component changes observed during independent ambulation in young children[END_REF]. Generally, walking gait is considered to be adult-like around the ages of 7 -8 years old [START_REF] Sutherland | The Development of Mature Gait[END_REF][START_REF] Vaughan | Neuromaturation of human locomotion revealed by non-dimensional scaling[END_REF]. Nevertheless, recent work has shown that immature gait characteristics retain during late childhood and early adolescence when walking at adult-like speeds and this relation is dependent upon leg length [START_REF] Froehle | Age-related changes in spatiotemporal characteristics of gait accompany ongoing lower limb linear growth in late childhood and early adolescence[END_REF]. Parameters such as double support time, single support time and base of support are only mature in females around the age of 14 years and in males age of maturation is as high as 18 years [START_REF] Froehle | Age-related changes in spatiotemporal characteristics of gait accompany ongoing lower limb linear growth in late childhood and early adolescence[END_REF]. Also consistency in variability of gait develops well into adolescence [START_REF] Muller | Intra-individual gait speed variability in healthy children aged 1-15 years[END_REF]. Thus, maturation of walking continues well after childhood. Changes in gait parameters can partially be explained by growth and changes in the dimensions of body segments [START_REF] Vaughan | Neuromaturation of human locomotion revealed by non-dimensional scaling[END_REF][START_REF] Froehle | Age-related changes in spatiotemporal characteristics of gait accompany ongoing lower limb linear growth in late childhood and early adolescence[END_REF]. It is accepted that maturationof balance control is another contributing factor in the maturation of gait [START_REF] Muller | Intra-individual gait speed variability in healthy children aged 1-15 years[END_REF][START_REF] Assaiante | Development of postural control in healthy children: a functional approach[END_REF][START_REF] Ivanenko | Kinematics in newly walking toddlers does not depend upon postural stability[END_REF]. Features of an immature gait, such as the prolonged duration of stance, increased double support time and the wider base of support are considered characteristic for an immature control of balance. So, control of balance is immature at the onset of walking and is a factor that drives changes in gait pattern throughout childhood. The spatial margin of stability is a measure derived from biomechanics characterizing dynamic balance control mechanisms during gait. We hypothesize there will be an agedependent decrease in the spatial margin of stability in children with typical development that can be related to changes in the spatio-temporal characteristics of gait.

Methods

Study design

Data on 3D gait analysis in typically developing children (age 1 -10) were retrospectively included in this cross-sectional study to investigate age-related changes in gait stability using the extrapolated Regarding gait stability, balance control is more critical in the ML direction, where it is ensured by active foot placement strategies that require central nervous system control [START_REF] Hof | Control of lateral balance in walking: Experimental findings in normal subjects and above-knee amputees[END_REF]. Based on the above mentioned balance principles, Hof showed that this foot placement strategy is primarily driven by the distance between the CoP and the XCoM at foot-of [START_REF] Hof | The 'extrapolated center of mass' concept suggests a simple control of balance in walking[END_REF], referred to as "spatial margin of stability (MoS)". If the MoS is too small, there is a risk of crossover steps, if it is too large, it forces the person to have a large step width and/or a small stance time leading to a non-efficient gait. Since, this quantity has proved to be a very interesting window into the balance control of gait [START_REF] Hof | Control of lateral balance in walking: Experimental findings in normal subjects and above-knee amputees[END_REF].

In a static (standing) situation, balance is achieved when the vertical projection of the center of mass falls within the base of support [START_REF] Shumway-Cook | The growth of stability: postural control from a development perspective[END_REF][START_REF] Winter | Human Balance and Posture Control during standing and walking[END_REF][START_REF] Hof | The condition for dynamic stability[END_REF]. However, this condition is insufficient in dynamic situations and the velocity of the center of mass should also be accounted for [START_REF] Hof | The condition for dynamic stability[END_REF]. Based on the inverted pendulum model, Hof [START_REF] Hof | The condition for dynamic stability[END_REF] showed it can be simply done by replacing the CoM by the XCoM, a point defined as projection on the ground of the COM augmented by a quantity proportional to its velocity.

By doing so, similar balance principles apply: the XCoM moves away (diverges exponentially) from the center of pressure and balance can be maintained only if the CoP can be placed outward of the XCoM. As the CoP is constrained to remain within the BoS, there are two possible situations: 1) the XCoM lies within the BoS and it can be captured by shifting the CoP within the BoS, or 2) it lies ouside and a step has to be performed in order to enlarge the BoS before being abble to capture the XCoM.

center of mass and spatial margin of stability concepts. The study protocols have been approved by the local ethical committees. Parents gave written informed consent at the time of study inclusion and were aware that data could be used retrospectively for further research. Data collection took place between October 2002 and November 2012.

Setting

Gait analysis data were collected using an integrated set-up with an optometric movement registration system (6 -8 cameras, Vicon Mcam 460 or T10 series, Vicon ® , Oxford, UK or Motion Analysis ® system with 8 Eagle ® cameras, Motion Analysis Corporation, Santa Rosa, USA) and multiple force plates (2 -4 platforms, AMTI OR6-7 or Bertec ® , Columbus, USA, dimensions 0.4 x 0.5m). Both systems were integrated in the walkway and synchronized at sampling frequencies of 100 Hz and 1000 Hz, respectively. Reflective markers (diameter 14 mm) were attached following the Davis model [START_REF] Davis | A gait analysis data collection and reduction technique[END_REF] or custom [START_REF] Samson | Biomechanical maturation of joint dynamics during early childhood: updated conclusions[END_REF].

Children were introduced to the walkway and were given time to explore and get used to the surroundings. After attaching the markers, again a short habituation period was provided. Data registration was started after the children were no longer aware of the attached markers when moving around. After performing a static anatomical calibration trial in which all markers were visible, dynamic trials were collected. During the dynamic trials one or two caregivers were standing on each end of the walkway to encourage the child to move towards them in a straight line. All children walked barefoot at self-selected speed.

Participants

Out of databases from previous studies all 3D gait analysis of typically developing children that met following inclusion criteria were selected: a total body gait analysis with clear foot strikes on the force plates and full marker visibility for at least two consecutive strides. Details of participant selection in the previous studies is dependent upon the study and can be found in previous publications [START_REF] Samson | Biomechanical maturation of joint dynamics during early childhood: updated conclusions[END_REF][START_REF] Zeni | Two simple methods for determining gait events during treadmill and overground walking using kinematic data[END_REF][20][START_REF] Van Dam | A cross-sectional study about the relationship between morphology and kinematic parameters in children between 15 and 36 months[END_REF][START_REF] Van De Walle | Mechanical energy estimation during walking: Validity and sensitivity in typical gait and in children with cerebral palsy[END_REF].

Variables of interest Anthropometric measurements

For each subject, information was obtained on age (in years), body mass, body length and leg length.

All measures were taken according to standard procedures.

Gait parameters

Walking speed, stride length, step width and duration of swing were used to characterize the temporospatial characteristics of gait. Gait parameters were considered as absolute as well as dimensionless values, i.e. normalized to leg length according to Hof [START_REF] Hof | Scaling gait data to body size[END_REF].

Data measurement

Marker trajectories were labeled in the Vicon Workstation (version 4.6 for Windows) or Nexus (version 1.8.x for Windows) software or Cortex software (Motion Analysis). Based on the force plate data and the ankle marker (malleolus lateralis) trajectories instances of foot strike and foot off were determined [START_REF] Zeni | Two simple methods for determining gait events during treadmill and overground walking using kinematic data[END_REF].

The total body center of mass was calculated using either the standard Vicon clinical model (Plug-In Gait application for Vicon Workstation and Nexus software) or custom (Appendix 1).

The .c3d files were then exported to Matlab and a custom written script was used to calculate the XCoM and MoS according to the formulas described by Hof [START_REF] Hof | The condition for dynamic stability[END_REF][START_REF] Hof | The 'extrapolated center of mass' concept suggests a simple control of balance in walking[END_REF]. Spatio-temporal parameters were calculated from the left and right ankle marker (malleolus lateralis) trajectories.

Bias

All children in this study were volunteers so some form of selection bias cannot be excluded. Children were assumed to show a typical development if no developmental problems were reported by the parents. No additional developmental tests were performed. It is therefore possible that a small number of children might be diagnosed with developmental problems or delay at a later age, especially in the younger age groups.

Study Size

The study is descriptive in nature. Therefore, a sample size calculation cannot be performed. As a rule of thumb, a sample size of 50 children is considered good while sample of 100 children is excellent when providing reference data [START_REF] Terwee | Rating the methodological quality in systematic reviews of studies on measurement properties: a scoring system for the COSMIN checklist. Quality of life research : an international journal of quality of life aspects of treatment, care and rehabilitation[END_REF].

Extrapolated Center of Mass (XCoM)

The "position of the extrapolated center of mass" (XCoM) is the vector sum of the center of mass position and a proportion of its velocity. In human movement, balance is typically related to the vertical projection of the body center of mass on the ground, that should fall within the base of support. This condition is insufficient in dynamic situations, in this case the velocity of the center of mass should also be accounted for. Thus, Hof (2005) defined this new vector quantity (XcoM) where the center of mass position is extrapolated in the direction of its velocity.

Spatial margin of stability (MoS)

The minimum spatial margin of stability was defined as the minimum distance between the center of pressure and the XCoM along the mediolateral axis during the single support phases. The mediolateral axis was defined as the axis in the transversal plane (XY plane) perpendicular to the walking direction, which was derived from the CoM coordinates. The spatial margin of stability is expressed in absolute values (mm, MoS) and normalized relative to leg length (MoS_LL)

Statistical methods

Prior to further analysis, outcome parameters were averaged over the different number of strides that were included for each participant. Outliers were either corrected if possible or removed from the dataset. The Shapiro Wilk test showed that the assumption of normality of the distribution was violated for the outcome variables MoS and MoS_LL. Thus, non-parametric statistics were used. To

Results

INSERT TABLE 1 HERE

Outcome data

Descriptives of gait parameters and spatial margin of stability can be found in Table 2.

INSERT TABLE 2 HERE

Main results

Pearson correlation coefficients

Descriptive data

Out of the available databases, a total number of 447 trials in 93 children were eligible for data analysis. After checking for outliers and missing data, data of 84 children with a mean age of 4.7 ± 2.7 years (age range 1 year -10 years 8 months) were included in the final analysis. Table 1 provides an overview of the number of children in each age range and presents their anthropometric characteristics.

stride length (rho = -0.510, p < 0.001), leg length (rho = -0.44. p < 0.001), age (rho = -.368, p = 0.001) and duration of swing (rho = -.350, p = 0.004). A weak but significant correlation was observed between the MoS and walking speed (rho = -.243, p = 0.026) and step width (rho = .285, p = 0.012).

INSERT TABLE 3 HERE

Table 4 presents bivariate Spearman rank correlation coefficients on the normalized data. Normalized spatial margin of stability MoS_LL shows a significant and strong positive correlation with dimensionless step width (rho = .757. p < 0.001). A significant and strong negative correlation was found with age (rho =-.677, p < 0.001) and duration of swing (rho = -.659, p < 0.001). A weak but significant correlation was found between MoS_LL and dimensionless stride length (rho = .295, p = 0.026). No significant correlation was observed between the normalized spatial margin of stability MoS_LL and dimensionless walking speed.

Stepwise linear regression

The test results of the stepwise linear regression model using the non-normalized data indicated one predictor, explained 18% of the variance in spatial margin of stability MoS (R² = 0.178, F(1,36) = 7.780, p = 0.008). It was found that duration of swing significantly predicted the spatial margin of stability ( = -0.422, p = 0.008). Figure 1 shows the linear relation between spatial margin of stability and duration of swing. Age, leg length, stride length, step width or walking speed did not significantly explain any additional variance in spatial margin of stability.

INSERT FIGURE 1 HERE

The test results of the stepwise linear regression model using the normalized data values, indicated three predictors, explained 66% of the variance in normalized spatial margin of stability MoS_LL (R² = 0.662, F(3,34) = 22.159, p < 0.001). It was found that dimensionless step width ( = 0.755, p < 0.001), dimensionless stride length ( = 0.421, p = 0.005) and duration of swing ( = -0.364, p = 0.014) significantly predicted the normalized spatial margin of stability. So, accounting for the leg length, children who take relatively wider steps, who take relatively longer steps and who show a shorter duration of swing have a larger normalized spatial margin of stability (Figure 2).

INSERT FIGURE 2 HERE

Discussion

Key results

The objective of this study was to investigate age related changes in spatial margin of stability as a model of studying development of dynamic balance control during locomotion in childhood. We hypothesized that the spatial margin of stability would be large at the onset of walking and would decrease with increasing age. These changes would be related to changes in the temporo-spatial characteristics of gait. This hypothesis was confirmed. Changes in the spatial margin of stability as a function of age, could be explained by changes in duration of swing. Furthermore, when normalizing for leg length, relations are found between the normalized spatial margin of stability, dimensionless step width, dimensionless stride length and duration of swing. Thus, the spatial margin of stability is primarily related to gait biomechanics. Off-course, these changes in gait biomechanics are agedependent.

Limitations

Despite its interesting findings, this study has some limitations. All children who participated in this study, did so on a voluntary basis which possibly introduces a selection bias. Furthermore, the lab setting creates an unnatural situation that might introduce small changes in the gait pattern of children compared to a real life situation. The distribution of ages is not homogeneous, there are slightly more younger children in the dataset than older children. However, normalized gait parameters are maturelike after the age of 3. Therefore, a minimal effect on our results is expected.

Data from different research centers and study protocols were pooled. Use of different marker protocols might influence the data. However, differences between the different protocols were statistically investigated using a one-way ANOVA and no differences were found.

Interpretation

Changes in gait pattern are driven by the interactions between growth, maturation and experience [START_REF] Froehle | Age-related changes in spatiotemporal characteristics of gait accompany ongoing lower limb linear growth in late childhood and early adolescence[END_REF][START_REF] Lacquaniti | Patterned control of human locomotion[END_REF]26]. It is difficult to separate these items from each other. Here we observed a significant decrease of spatial margin of stability with age. However, when pooled with spatio-temporal parameters, age did not come out as a primary factor in the stepwise regression model of the spatial margin of stability. Does this mean that dynamic balance control during locomotion does not change throughout childhood? Probably not, since significant relations were found between the spatial margin of stability when normalized to leg length and different gait parameters. Changes in these specific gait parameters, i.e. duration of swing, stride length and step width are known to be age dependent [START_REF] Sutherland | The Development of Mature Gait[END_REF][START_REF] Sala | Gait component changes observed during independent ambulation in young children[END_REF].

Interestingly, relations between spatial margin of stability and gait parameters could be linked to the mechanics of gait control in the medio-lateral direction, i.e. the direction where stability constraints are the most pregnant [START_REF] Hof | Control of lateral balance in walking: Experimental findings in normal subjects and above-knee amputees[END_REF]. In particular, during the single stance phase, and by neglecting the small center of pressure adjustments under the stance foot, it showed that the center of mass diverges proportionally to the margin of stability and exponentially with time. Thus, a larger spatial margin of stability induces a larger lateral divergence of the CoM that could be compensated by a larger and/or a quicker step. Our data are in line with these predictions: the normalized spatial margin of stability showed both a positive correlation with step width and a negative one with duration of swing (Figure 2). Moreover, an immature balance control requires an increase of the MoS to avoid frequent falls.

Larger step width and shorter single stance duration could be an adaptation to his immature balance control.

Values of spatial margin of stability found in this study appeared coherent with previous findings on adults. As an example, Hof [START_REF] Hof | The condition for dynamic stability[END_REF][START_REF] Hof | The 'extrapolated center of mass' concept suggests a simple control of balance in walking[END_REF]

Future research

Further research should focus on how children actually regulate their stability in the medio-lateral direction, and how this regulation evolves with age. In particular, it could be interesting to see if we can observe an evolution in the fine regulations of the center of pressure under the foot during stance phase [START_REF] Hof | The 'extrapolated center of mass' concept suggests a simple control of balance in walking[END_REF] with age, that could be related to the biomechanical maturation of ankle joint dynamics around the age of four [START_REF] Samson | Biomechanical maturation of joint dynamics during early childhood: updated conclusions[END_REF].

Conclusion

The spatial margin of stability, MoS, decreases with increasing age. However, development of dynamic balance control, represented by these changes in MoS, is linked to changes in step-time parameters of gait.

Conflict of Interest

None.

[26] Rodriguez EB, Chagas PSC, Silva PLP, Kirkwood RN, Mancini MC. Impact of leg length and body mass on the stride length and gait speed of infants with normal motor development: A longitudinal study. Brazilian Journal of Physical Therapy / Revista Brasileira de Fisioterapia 2013; 17(2): 163-169.

  describe the outcome parameters, median and interquartile range (IQR) were calculated per biological age. To explore the relation between the spatial margin of stability MoS or MoS_LL, age, leg length and gait parameters Spearman rank correlation coefficients were calculated.Variables showing a significant correlation (p < 0.05) with the spatial margin of stability were afterwards added to a stepwise linear regression model starting with the variable showing the highest coefficient of correlation. Two models were generated, for both (non-)normalized variables.Normalized values were considered to remove the effect of growth. Variables were entered into the model in case the probability of F was < 0.05. Variables were removed from the model in case the probability of F was >0.10. Goodness of fit was investigated by R² values. Residuals were checked for normality using a histogram and normal P -P plot of regression standardized residual. Statistical analysis was performed using IBM® SPSS® statistics software (version 23.0.0 for Windows, IBM corporation).

  reported values during walking in the range of 1.5 -3 % of leg length. If younger children displayed larger normalized MoS values, older children reached this range from the age of 7 years old.

Figure 1 :

 1 Figure 1: Relation between the spatial margin of stability and duration of swing. The linear regression line with 95% confidence interval was fitted using the least squares error method.

Figure 2 :

 2 Figure 2: Relation between the normalized spatial margin of stability (% of leg length, MoS_LL) as dependent variable and dimensionless step width (first panel), dimensionless stride length (second panel) and duration of swing (third panel). The linear regression line with 95% confidence interval was fitted using the least squares error method

Table 3

 3 presents bivariate Spearman rank correlation coefficients based on the non-normalized data.The spatial margin of stability MoS shows a highly significant but moderate negative correlation with

Table 1 :

 1 Descriptive data of the study sample including the age ranges, mean ages and anthropometric characteristics (N = number of children, F = number of females)

											l
											e
											g
											l
											e
	Age Groups			Age (years)		mass (kg)	body height (mm)	BMI		n g t
											h
											(
											m
											m
											)
		N	F	Mean SD	Mean SD	Mean SD	Mean SD	Mean	S D
	1 year -1 year 11 months	17 12 1.4	.2	10.98 1.17 791	43	17.7	2.0 311	3
												5
	2 years -2 years 11 months	12 7	2.3	.2	13.55 1.61 915	46	17.3	2.1 350	5
												3
	3 years -3 years 11 months	8	4	3.4	.3	16.10 2.37 1004	50	17.0	1.9 454	4
												5
	4 years -4 years 11 months	11 8	4.5	.3	17.80 2.64 1065	51	16.9	1.9 514	3
												3
	5 years -5 years 11 months	10 6	5.4	.3	18.67 2.36 1098	51	16.6	1.9 527	2
												8
	6 years -6 years 11 months	10 4	6.3	.3	21.12 3.07 1175	61	16.6	1.8 586	3
												3
	7 years -7 years 11 months	4	2	7.5	.2	24.60 .78	1256	18	16.6	1.8 654	2
												7
	8 years -8 years 11 months	4	2	8.5	.3	25.30 2.95 1319	13	16.4	1.8 691	2
												3
	9 years -9 years 11 months	3	1	9.4	.3	33.77 3.40 1383	23	16.2	1.5 717	1
												5
	10 years -10 years 11 months 5	2	10.5	.2	32.78 2.89 1432	46	16.1	1.5 766	2
												1
	All	84		4.7	2.7 18.29 6.77 1054	198			491	1
												4
												8

Table 2 :

 2 Descriptive data of gait parameters and spatial margin of stability presented according to biological age: both absolute values and relative values normalized to leg length (%) are represented as median (IQR) for each biological age group

Acknowledgements

This work was supported by the Flemish Research Council [grant number AUHA/09/006], the University of Antwerp Research Council [grant number ID33549], the University of Antwerp, the Antwerp University Hospital and the University of Lyon.

 298) 165 ( 33) 120 ( 33) 26 ( 16) 31 ( 5)

Appendix: Methodology used to estimate the whole body center of mass from LBMC's data

The following procedure was used to estimate the whole body center of mass from marker trajectories:

1-Sixteen ( 16) reflective markers were placed on particular anatomical points of the trunk, pelvis and lower limbs: left and right acromions, anterior and posterior iliac spines, medial and lateral femoral epicondyles, medial and lateral maleolli and first metatarsal heads. 3-Mass and position of the center of mass of 7 segments (feet, shanks, thighs and Head-Arm-Trunk) were then estimated from these markers and hip joint centers using regressions from Jensen (1989), considering children's age.

4-

The whole body center of mass was then obtained as the weighted barycenter of segment's centers of mass.

This procedure is summarized in the figure below:

Appendix II: The spatial margin of stability and the normalized spatial margin of stability plotted per age group


J Massoni

L Biamino

G Jourdan

O Igra

L Houas

EXPERIMENTAL AND NUMERICAL INVESTIGATION OF WEAK BLAST WAVE INTERACTION WITH A THREE LEVEL BUILDING

It is well known that strong blast waves are lethal for a living creature, while it is believed that weak blast waves are harmless. This might be the case when there is a head-on collision between an on-coming weak blast wave and the considered body; however, a completely different scenario takes place when the initially weak (safe) blast wave hits a body after multiple reflections from existing walls. Such cases take place when a weak blast wave, resulting from a sudden explosion, hits a residential complex.

In the conducted experiment a very simple model of a three level building is subjected to a weak blast wave. The evolved wave pattern inside the building rooms is recorded by a sequence of schlieren photos and the prevailing pressures are deduced from numerical simulations.

INTRODUCTION

While it is well known that strong blast wave are lethal for a living creature, it is commonly accepted that weak blast waves are harmless. While this might be the case when there is a direct impingement of the on-coming weak blast wave on a body, a completely different scenario takes place when the weak blast wave hits a body after multiple reflections from existing walls. Such a case can be found when a weak blast wave, resulting from a sudden explosion, hits a residential complex. Furthermore, explosion injuries can be caused by direct hit of a blast wave, and/or injury caused by flying debris resulting from the explosion, and/or injuries resulting from being dragged/pushed by the blast wave [START_REF] Debien | Lésions par explosion[END_REF]. The present investigation focuses on potential injuries from collision with a multiple reflection of initially weak blast wave.

In the conducted experiments a very simple model of a three level building is subjected to an initial harmless blast wave. Attention is given to find where the most dangerous places inside the building are, places without having death risk, and where are the safest places for people to be. Information given in Table 1 reveals human response to overpressure; this information is used in defining a place as being safe or non-safe.

Initially a weak, planar shock wave propagates inside the shock tube. Upon reaching the tube exit-open-end it immerges into the open atmosphere as a blast wave. The investigated building model is placed in proximity to the shock tube exit as shown in Fig. 1. Supplementing the experimental investigation, numerical simulation of the generated flow field, caused by the blast wave interaction with the building model, is conducted. In urban life blast waves can be generated from a sudden explosion of gas lines placed underground, from a truck carrying explosive substances in a residential environment, and/or from a terrorist attack in the public domain.

Considering weak blast waves, all these examples have one thing in common, i.e., the fact that generally people are not informed of potential risks to their ears or lungs, not to mention life, if they stand in proximity to a wall or near a corner close to where a blast wave propagates. Subsequently it will be shown where are the dangerous places inside a residential building impacted by a weak blast wave.

From a fluid mechanics point of view, the present work deals with the problem of multiple shock/blast wave reflections. Studies concerning shock wave reflections from rigid boundaries can be found in the open literature; e. g. in Ben-Dor [START_REF] Ben-Dor | Shock Wave Reflection Phenomena[END_REF], Igra et al. [3] and more recently, in Volume 2 of the Handbook of Shock Waves by Ben-Dor, Igra and Elperin [START_REF] Van Leer | On the relation between the upwind-differencing schemes of Godunov, Enguist-Osher and Roe[END_REF]. Table 1: Direct effects of blast wave generated overpressure on humans (from [START_REF] Debien | Lésions par explosion[END_REF], [5] and [START_REF] Van Leer | Towards the ultimate conservative difference scheme V. A Second order sequel to godunov's method[END_REF]). Direct effect means without considering secondary effects such as flying glass and/or brick debris.

EXPERIMENTAL SET-UP

The present experiments were conducted using an open end shock tube and thereby generating a very weak blast wave. The generated blast waves were sufficiently weak to ensure that they pose no life endangerment to people.

The building model was constructed from Plexiglas in order to allow high speed shadowgraph and schlieren visualizations; the facility used is shown in Fig. 1. 7!"# 89 $%@@ A%B %@8C D8@E F"@ A%B %@8C 7G8%H " AIFBIFC

NUMERICAL SIMULATION

All computations were performed for air; due to the relatively low Mach number of the incident shock wave it can be treated as a perfect gas. The conservation equations for ideal compressible fluids (the Euler equations), indicating conservation of mass, momentum and energy were solved. Details regarding the used numerical solution can be found in.

RESULTS AND DISCUSSION

A sequence of schlieren photos showing blast wave interaction with the building model is shown in Fig. 3. At t=33 µs the blast wave just started its interaction with the building model. With progressing time, the blast wave penetrated into the building and at t=200 µs it is clearly seen inside the first floor, approaching the room ceiling. From the color intensity of the blast wave it is apparent that the part which penetrated into the building first floor is weaker than the part propagating in the open atmosphere. At t=300 µs the following waves are visible: the progressing main blast wave starts penetrating into the building second floor while the part penetrated into the first floor is seen now reflected from the first floor ceiling and the blast segment that propagated earlier toward the wall separating between the two rooms of the first floor has now penetrated into the second room. As noticed earlier, the strongest part of the blast wave, at t=300 µs, is the one propagating through the open atmosphere. However, the color intensity of the regular reflection of the blast wave from the ceiling of the first floor second room increases, indicating a strengthening blast wave at this location; see in Fig. 3 at t=300 and 333 µs. Further confirmation of the local pressure increase prevailing around the blast reflection in the considered time is evident in Fig. 4 where computed over-pressure is shown. A significant pressure increase is noticed at t=433 and 466 µs in the corner area between the ceiling of the first floor and the room wall; an overpressure, ∆p of 10 kPa, is seen in Fig. 4. It is apparent from Fig. 3 that the cause for this overpressure increase near the corner is the blast reflection from the room ceiling and wall. This is the strongest overpressure inside the building model witnessed during the experiment. As is evident from Table 1 this overpressure is not fatal, however it can cause damage to human eardrums. A slightly stronger blast wave could result in irreversible damage, i.e., fatal results. At later times, t ≥ 600 µs the reflected blast wave from the corner weakens while it propagates towards the room floor; see Fig. 4 at t = 600 and 666 µs.

With increasing time, 800 µs ≤ t ≤ 1033 µs multiple wave reflection between the room walls, floor and ceiling takes place as is evident in Fig. 3. From the color intensity of these waves as shown in the schlieren photos (Fig. 3), it is clear that these reflecting waves are much weaker than those seen outside the building model.

It is of interest to note from Fig. 4 that at later times, at 400 < t < 500 µs the overpressure experienced by those exposed to the travelling blast wave in the open atmosphere is smaller than that prevailing inside the building near the room corners.

CONCLUSIONS

It is clear from the present investigation that blast wave damages are higher for people standing near a wall, or even higher when near a corner. Furthermore, at a late time (t > 400 µs) the overpressure behind reflected blast wave from a room corner is significantly higher than that experienced in an open space. This is due to the pressure increase behind the reflected blast wave from the corner between the room ceiling and walls. The present work might be useful in designing building structures complying with safety standards, as well as for people who are not familiar with explosion generated blast waves.

Finally, the present work could also be helpful to validate numerical codes able to predict such effects without taking into account broken construction materials and broken glasses.

Figure 1 :

 1 Figure 1: View of the Plexiglas building model and its location close to the shock tube exit.

Figure 2 :

 2 Figure 2: Description and dimensions (in mm) of the tested building model. (a) Schematic of the model; (b) photograph of the Plexiglas building model, and (c) its location (in cm) relative to the shock tube exit.

Figure 3 :

 3 Figure 3: Sequences of schlieren photos showing the evolution and multiple reflections of the blast wave inside the Plexiglas building. The incident shock wave Mach number inside the shock tube is 1.17.

Figure 4 .

 4 Figure 4. Numerical simulations showing the overpressure encountered within the three floor building model for an incident shock wave, inside the shock tube, having a Mach number of 1.17. Experimental results are shown on the left and appropriate simulations appear on the right.


A Inglebert

A Ghizzo

T Reveille

D Del Sarto

P Bertrand

F Califano

Multi-stream Vlasov model for the study of relativistic Weibel-type instabilities

We describe and apply a recently proposed model [ EPS 95, 45002 (2012) ] using a Hamiltonian formalism for the study of Weibel-type instabilities in the relativistic limit. Taking advantage from the invariance of the generalized canonical momentum, it allows for a drastic reduction of the computational time when compared to the full Vlasov-Maxwell system of equations. Analytically, the model is exact and we recover the standard fluid dispersion relations in the case of the Weibel and filamentation instabilities, to be considered as out benchmark physical problems.

Introduction

In these last decades, an intense research activity has been carried out for the study of the Weibel and Current Filamentation Instability (hereafter WI and CFI). Indeed, both instabilities are considered as basic processes in a plasma since associated to the possibility of generating "quasi-static" magnetic fields. The WI, first proposed by Weibel [1] (1959), is driven by the presence in the plasma of a thermal anisotropy, for example in the form of an electron bi-Maxwellian plasma with T ⊥ > T , where ⊥ and represent the perpendicular and parallel direction with respect to the wave vector k ≡ k . In this case, the micro-currents produced by any small perturbation, rearrange through the growing waves into filaments, and separate more and more, thus selfconsistently generating a corresponding magnetic field [2].

A similar process occurs in the presence of an electron beam and its return current produced by the background plasma in order to maintain quasineutrality. Any infinitesimal transversal perturbation separates these currents that repeal each other, thus reinforcing the initial perturbation. As a result, a magnetic field similar to the Weibel case is generated by the CFI [3,4]. This instability plays a key role in the context of laser plasma interaction, in particular concerning the problem of fast ignition [5]. In laser-plasmas, fast electron beams are routinely produced during the propagation of an ultraintense, ultra-fast laser pulse, as shown by means of PIC simulations in the over-dense regime [6]. For an under-dense plasma, the CFI has been invoked in order to to explain the quasi-static magnetic fields observed in PIC simulations in the wake of the laser pulse [7] . The concept of CFI was first mentionned by Fried in Ref. [8] studying the propagation of counter-streaming beams.

The physical mechanisms underlying the development of both the Weibel and CF instabilities, is very similar and has to be found in the free energy stored in the initial phase space velocity anisotropy. In both cases, the initially unmagnetized system transforms the initial thermal, or kinetic, energy into magnetic energy associated to the generation of a dipolar magnetic field per wave-length. Due to the mathematical analogies, it is possible to describe both types of instability by the same analytical model based on the description of the perpendicular momentum space in terms of "streams".

Furthermore it is worth noticing that the relativistic formulation of the Weibel instability remains a complex problem since temperature anisotropy distributions cannot be easily introduced starting from the standard Maxwell-Jüttner distribution. This important point of the relativistic kinetic theory was first pointed out by Yoon in Ref. [9] who proposed a modified Maxwell-Jüttner distribution to include temperature anisotropy. However, a straightforward extension of the non relativistic approach to the relativistic regime led to a substantial discrepancy between theory and simulations, as recently observed in PIC simulations [10]. Thus the choice of an equilibrium distribution function in the relativistic regime is still today an important open problem. However, the main difficulties when studying the non-linear evolution of Weibel-like instabilities arise from the complexity to follow numerically the dynamics in phase space, in particular in the relativistic regime where the standard splitting scheme [11] cannot be used.

Our aim here is to apply a new kinetic model, recently proposed in Ref. [12] making analytical calculations more tractable with respect to the full Vlasov-Maxwell system, but capable of representing the evolution of any anisotropic phase space distribution starting from a general initial condition. The model is here extended and apply to characterize the linear and nonlinear regimes of the instability till saturation. Such a detailed description of the different regimes is afforded by following each stream separately. Furthermore from a numerical point of view, the method also allows for a drastic reduction of the computational time.

We adopt a kinetic description and we compare our results with those obtained by solving the full Vlasov-Maxwell system of equations in the case of the CFI (see [13] and refs. therein) and of the WI [14]. The model, based on a Hamiltonian reduction technique, makes use of a multi-stream description in the limit of a one-dimensional plasma, where exact canonical invariants can be defined in the transverse direction to reduce the dimension of the momentum space (here p ⊥ ). The use of such invariants allows to find out a broad class of exact non linear solutions of the Vlasov-Maxwell system. Our approach is similar to the one presented in Ref. [15] for a drift-kinetic model and in Ref. [16] for the study of magnetic reconnection. Finally, we note that since we are limited to a one-dimensional plasma, the coupling with the purely electrostatic two-stream instability is excluded.

The paper is organized as follows. In Section II we present the multistream model as derived from the Vlasov-Maxwell system. In Section III we recover the dispersion relation in the relativistic, fluid regime. The nonrelativistic limit of the dispersion relation is discussed in Section IV where we specifically address the case of the CFI and of the WI. In Section V we present a numerical comparison of our multi-stream model with the full 1D-2V Vlasov-Maxwell system in the case of the CFI. A non-symmetric case of CFI is also investigated in this section, with particular emphasis on the efficiency of the model to describe the saturation mechanism of CFI in detail. Conclusions are given in Section VI.

The multi-stream Vlasov-Maxwell model

We consider a collisionless plasmas and we adopt a kinetic description since, at least concerning the instability saturation mechanism, kinetic effects dominate and are expected also to play a key role in the further non linear dynamics. We assume, as first, the ions as a fixed neutralizing background. The Vlasov equation, self-consistently coupled to the Maxwell's equations, reads:

∂F ∂t + p x m e γ ∂F ∂q + e E + p × B m e γ ∂F ∂p = 0 (1)
γ = 1 + p 2 m 2 e c 2 (2)
where F (q, p, t) is the electron distribution function and γ is the Lorentz factor. Here we limit our study to wave modes propagating along the x-direction and to a 2V (p x , p y ) momentum space with a magnetic field directed along the z-axis. In this case, the Vlasov equation possess an exact invariant giving the possibility to reduce the dimension of the momentum space. In our model, introduced in its main aspects in Ref. [25], the possibility of describing "exactly" the full dynamics in a reduced phase space is achieved by using a finite number N of "streams" (or bunches) of particles, whose dynamics is described by a kinetic Vlasov-type equation, independent of the transverse kinetic momentum of coordinate p ⊥ . In the following, we discuss in detail the main points of the model . We start by considering the Hamiltonian of one particle,

H = m e c 2 (γ -1) + eφ (x, t) (3)
where φ is the electrostatic potential, m e the electron rest mass and e the electron charge. By using the canonical momentum P c = p + eA ⊥ , where A ⊥ is the vector potential, the Lorentz factor, as defined in Eq. ( 2), can be expressed as

γ = 1 + (P c -eA ⊥ (x, t)) 2 m 2 e c 2 1/2 (4)
Now, along the longitudinal x-direction, the canonical momentum simply reduces to P c,x = p x , while in the perpendicular direction to P c,⊥ = p ⊥ + eA ⊥ (note that the choice A = A ⊥ corresponds to the choice of a Coulomb gauge ∇ • A = 0). Thus the Hamilton equation ∂P c /∂t = -∂H/∂q can be split into the longitudinal and perpendicular direction:

dP cx dt = - ∂H ∂x (5) dP c,⊥ dt = - ∂H ∂q ⊥ = 0 (6)
since the Hamiltonian, as defined in Eq. (3), does not depend on the transverse spatial coordinates.

Owing to the invariance of the generalized canonical momentum, we can represent the plasma evolution by a class of initial conditions invariant under the dynamics. Without loss of generality, we represent the plasma in term of N bunches of particles, each "stream" j, j = 1, ..., N , having the same initial perpendicular momentum P c,⊥ = C j = const.

Using Eq. ( 5), we consider, for each j-th particle population, a distribution function f j (x, p x , t) that must satisfy a reduced Vlasov-type equation. The Hamiltonian of a particle belonging to the j-stream, is given by:

H j = m e c 2 (γ j -1) + eφ (x, t) (7)
From Eq. ( 4), the γ j Lorentz factor for the j-stream now reads:

γ j = 1 + p 2 x m 2 e c 2 + (C j -eA ⊥ (x, t)) 2 m 2 e c 2 1/2 (8)
The j-stream is described by the distribution function f j that must satisfy a "reduced" Vlasov equation:

Df j Dt = ∂f j ∂t + [f j , H j ] = 0 (9)
where [f j , H j ] represents the standard Poisson bracket. By casting Eqs. ( 7)- (9), we obtain the reduced Vlasov-type equation for the j-stream:

∂f j ∂t + p x m e γ j ∂f j ∂x + eE x - 1 2m e γ j ∂ ∂x (C j -eA ⊥ ) 2 ∂f j ∂p x = 0 (10)
These N equations, representing what we defined as the multi-stream model, must be now coupled to the Maxwell's equations. By using the Poisson equation in the form

∂E x ∂x = e ε 0 (n (x, t) -n 0 ) ( 11
)
where n is the density defined by

n (x, t) = N j=1 ˆ+∞ -∞ f j (x, p x , t) dp x ( 12
)
we can cast the Maxwell equations into a single one for the perpendicular potential vector A ⊥ :

∂ 2 A ⊥ ∂t 2 -c 2 ∂ 2 A ⊥ ∂x 2 = 1 ε 0 N j=1 J ⊥,j (x, t) (13)
where J ⊥,j is the perpendicular current associated with the particle population j, i.e. the j-stream, and defined by:

J ⊥,j = e m e (C j -eA ⊥ ) ˆ+∞ -∞ f j γ j dp x (14)
In summary, our multi-stream model is given by N reduced Vlasov equations of type (10) for each j-stream, self-consistently coupled to the electrostatic longitudinal field E x , Eq. ( 11), and to the potential vector A ⊥ , Eq. ( 13). One remark must be pointed out:

The full distribution function F (x, p x , p ⊥ , t) can now be written as a sum of Dirac distributions:

F (x, p x , p ⊥ , t) = N j=1 f j (x, p x , t) δ (p ⊥ -(C j -eA ⊥ (x, t))) (15)
thus, reducing the 3D phase space into a 2D one (plus N values for the corresponding C j ).

In the case of a two-dimensional (x, y) plasma, the z-direction is thus reduced to the p z component perpendicular to the (x, y) plane. Since the corresponding Hamiltonian does not depend on the variable z, as F indeed, it is possible to find a canonical invariant P c,z = p z + eA z with a potential vector of the type A = A z (x, y, t) e z . Thus, the condition we need to make use of the conservation of the transverse canonical momentum to reduce the phase space dimension, is the possibility to separate the electrostatic and electromagnetic contributions in the electric field E = -∂A/∂t -∇φ.

The multi-stream model in the fluid approximation

We now derive the dispersion relation of Weibel-type instabilities in the fluid limit making use of the multi-stream model. We consider in Eq. ( 15), the case N = 1, corresponding to a cold plasma distribution function in the perpendicular direction. This limit has been intensively investigated in the context of laser-plasma interaction, in particular for the study of parametric Raman-type instabilities in the semi-relativistic [17] and strongly relativistic regime [18], or in the study of self-induced transparency in over-dense plasmas [19].

The "multi-fluid" model is obtained by taking the moments of the Vlasov equation, Eq. ( 10). Let's now introduce, for each j-stream, the density n j , the averaged longitudinal momentum u j and a mean Lorentz factor γ j . Without loss of generality, we assume a linearly polarized electromagnetic wave. We can write the potential vector and the transverse canonical momentum in the form:

A ⊥ = A y (x, t) e y ; P c,⊥ = P c,y e y = const = C j e y (16) which corresponds to the most general case due to the x-axis symmetry. Here C j is a constant. For each j-stream, the continuity, motion and Poisson equations can be written as:

∂n j ∂t + ∂ ∂x n j u j m e γ j = 0 ( 17
)
∂u j ∂t + u j m e γ j ∂u j ∂x = eE x - 1 2m e γ j ∂ ∂x (C j -eA y ) 2 (18)
∂E x ∂x = e ε 0 N j=1 n j (x, t) -n 0 (19)
Using Eq. ( 13), we obtain the following equation for the vector potential:

∂ 2 ∂t 2 -c 2 ∂ 2 ∂x 2 + N j=1 e 2 m e ε 0 ρ j A y = e m e ε 0 N j=1 C j ρ j (20)
where ρ j is the "stream" density defined as:

ρ j = ˆ+∞ -∞ f j γ j (x, p x , t) dp x (21)
Finally, the fluid limit of the Lorentz factor is obtained from Eq. ( 8):

γ j = 1 + u 2 j m 2 e c 2 + (C j -eA y (x, t)) 2 m 2 e c 2 1/2 (22)
From Eqs. ( 17) to ( 20), together with the definition of the "stream" charge density and of the "averaged" Lorentz factor, Eqs. ( 21)-( 22), represent a closed system.

We now make a linear, normal modes analysis around an equilibrium characterized by a mean density n 0,j (for the stream j) and a zero mean velocity u j0 . We write the fluctuating quantities, i.e. the density, velocity, vector potential, longitudinal electric field, as δn j , δu j , δA y , δE x . By linearizing Eqs. (17) to (20) and performing a Fourier transform in time (∂/∂t = -iω) and space (∂/∂x = ik), we get the following set of equations:

-iωδn j + ikn 0,j m e Γ 0,j δu j = 0 (23)
-iωδu j = eδE x + iek m e Γ 0,j C j δA y ( 24
)
ikδE x = e ε 0 N j=1 δn j (25)
-ω

+ k 2 c 2 + N j=1 ω 2 p,j Γ 0,j δA y - N j=1 ω 2 p,j C 2 j m 2 e c 2 Γ 3 0,j δA y = e m e ε 0 N j=1 C j δn j Γ 0,j (26)
where ω 2 p,j = n 0,j e 2 /m e ε 0 is the "plasma frequency" of the j-stream and the quantity Γ 0,j corresponds to the zero-order Lorentz factor:

Γ 0,j = 1 + C 2 j m 2 e c 2 (27)
Note that, in Eqs. ( 25) and ( 26), we have imposed, at zero order, the charge neutrality condition and an initial zero net current:

N j=1 n 0,j = n 0 ( 28
) N j=1 C j n 0,j Γ 0,j = 0 (29)
By solving the algebraic system of equations, Eqs. ( 23)-( 26), we finally get the dispersion relation of our multi-streams fluid system:

1 - N j=1 ω 2 p,j ω 2 Γ 0,j -ω 2 + k 2 c 2 + N j=1 ω 2 p,j Γ 3 0,j + k 2 c 2 ω 2 N j=1 ω 2 p,j Γ 3 0,j C 2 j m 2 e c 2 = - k 2 c 2 ω 4 N j=1 ω 2 p,j Γ 2 0,j C j m e c 2 (30)
The analytical formulation of the linear dispersion relation, Eq. ( 30), as a discrete summation over a collection of streams, provides a general and exact approach to take into account any anisotropy of the distribution function.

The multi-stream model in the weak relativistic regime

From a mathematical point of view, the multi-stream model can be seen as a "reduction technique" of the phase space dimension, due to the elimination of the perpendicular momentum p ⊥ that remains hidden in the streams C ⊥, j . Of course, an accurate description of a continuous distribution in p ⊥ , can be obtained only using a large value of N . However, in the case of the Weibel instability, the macroscopic description requires only the first moments (see for instance Ref. [20]), indicating that only a small number of streams is necessary to give a complete description of this instability. For the CFI, only two streams can be used to describe the instability. The use of exact invariants simplify the analysis of a Hamiltonian system such as the Vlasov system. For example, in the (x, p) phase space, the exchange of perpendicular momentum is described by the partial derivatives ∂/∂p ⊥ . From a numerical point of view, the derivative operator ∂/∂p ⊥ has to be approximated by some finite difference scheme. Therefore a minimum size for the mesh in the perpendicular space p ⊥ is required and we are faced with the usual sampling problem; on the other hand, a multi-stream description can be used even with a small number of streams, sometimes N = 2 or N = 3, or even with N = 1, with a significative reduction of the numerical effort. The particular case N = 2 corresponds, as already noted, to the current filamentation instability. In the limit of two initial "streams", N = 2, Eq. ( 30) reduces to the counter-streaming instability in the relativistic regime, in agreement with Pegoraro et al. [21].

The interesting result is the fact that Eq. ( 30) contains only quadratic terms in C j . Furthermore, it is always possible to define the quantity

N j=1 ω 2 p,j C 2 j /Γ 3 0,j m 2 e c 2
that can be view as a "relativistic temperature". We also note that in the non relativistic regime, it is easy to construct the different streams associated to a distribution function characterized by a density n 0 and a perpendicular thermal velocity v th , as we will see, in the form N j=1 ω 2 p,j C 2 j /m 2 e = ω 2 p v 2 th , allowing one, for example, to recover exactly the Weibel dispersion relation. Finally, from Eq. ( 30) we point out that:

(i) It is possible to recover the non relativistic regime by taking the limit Γ 0,j → 1 in Eq. (30). Thus, in a multi-stream framework, the non relativistic dispersion relation of a Weibel-type instability is given by:

k 2 c 2 + ω 2 p -ω 2 + k 2 c 2 ω 2 N j=1 ω 2 p,j C 2 j m 2 e c 2 = 0 (31)
This equation shows that the longitudinal, purely electrostatic, and the transverse, purely electromagnetic, contributions are decoupled.

(ii) We can discretize the streams "coarsely", since there is no exchange of particles between the different populations, as known for other Hamiltonian systems, as for example in the water bag model [22,23] or in the contour dynamics technique [24].

In the following, we consider two basic examples of instability in the non relativistic regime and in the approximation of a cold longitudinal plasma temperature.

The current filamentation instability

We take the limit N = 2. Following Eq. ( 15), it is possible to represent the initial distribution function in the following form:

F (x, p x , p y , t = 0) = n 0 δ (p x ) G (p y ) = n 0 δ (p x ) 2 j=1 F j δ (p y -C j ) (32)
where the quantities F j are just normalization factors. We consider the case of the CFI. This instability is driven by electron momentum anisotropy in velocity space, corresponding to two counter-streaming electron beams of momentum p 0,1 < 0 and p 0,2 > 0 and densities n 0,1 and n 0,2 . Quasi-neutrality imposes that the total net current must be zero, corresponding to the condition n 0,1 p 0,1 + n 0,2 p 0,2 = 0, or

2 j=1 n 0,j C j = 0 (33)
In Fig. 1 we show a schematic representation of the function G (p y ) as represented by the two invariants C 1 = p 0,1 and C 2 = p 0,2 . Since the plasma density is normalized to one, ´F (x, p x , p y , t) dp x dp y = 1, we get:

ˆG (p y ) dp y = 2 j=1 F j = F 1 + F 2 = 1 (34)
By considering the non relativistic limit, Γ 0,j = 1, charge neutrality and zero net current conditions, Eqs. ( 28) and ( 33), reduce to

n 0,1 + n 0,2 = n 0 ; n 0,1 C 1 + n 0,2 C 2 = n 0,1 p 0,1 + n 0,2 p 0,2 = 0 (35)
from which we get

2 j=1 C 2 j n 0,j = -n 0 C 1 C 2 = -n 0 p 0,1 p 0,2 (36)
By replacing this expression in the dispersion relation, Eq. ( 31), we recover the standard dispersion relation for the non relativistic CFI:

ω 4 -k 2 c 2 + ω 2 p ω 2 + ω 2 p k 2 c 2 p 0,1 p 0,2 m 2 e c 2 = 0 (37)
in agreement with Refs. [27,28] where the CFI dispersion relation is obtained from the fluid equations in the relativistic regime.

The case of the Weibel instability

By taking N = 3 streams, it is possible to recover the Weibel dispersion relation driven by a temperature anisotropy in velocity space. In Fig. 2 we show a schematic representation of this system and we define an initial distribution function F in the form of a superposition of three of streams:

F (x, p x , p y , t = 0) = n 0 δ (p x ) G (p y ) = n 0 δ (p x ) N =3 j=1 F j δ (p y -C j ) (38)
Let's now consider a homogeneous Maxwellian distribution function F (p), of corresponding density n 0 = 1 and thermal velocity v th = k B T /m e . The zero order momentum ˜F (p) dp x dp y = n 0 , in a multi-streams approach, is given by:

ˆG (p y ) dp y = 3 j=1 F j = 1 (39)
The first-order p y -momentum ˜py F (p) dp x dp y = 0 and imposing the zero net current condition, Eq. ( 33), reads: 41) and (42), we obtain the following relations:

3 j=1 C j F j = 0 (
3 j=1 C 2 j F j = 2F 1 C 2 1 = m 2 e v 2 th ( 43
)
3 j=1 F j C 4 j = 2F 1 C 4 1 = 3m 4 e v 4 th (44)
By solving Eqs. (39), ( 43) and (44) we get

F 1 = F 3 = 1/6, F 2 = 2/3 and C 3 = -C 1 = √ 3m e v th .
By using these parameters in Eq. ( 31), we finally recover the well-known dispersion relation, first obtained by Weibel [1]:

ω 4 -(ω 2 p + k 2 c 2 ) ω 2 -ω 2 p k 2 v 2 th ω 2 = 0 (45)
5 Numerical comparison between the full kinetic approach and the reduced Vlasov models

The Vlasov simulation (V-model)

We now focus on the numerical study of the current filamentation instability aiming at comparing, in the 1D-2V limit, the full Vlasov -Maxwell system, Eqs. ( 46), ( 47), (48), with the reduced Vlasov-Maxwell model presented in Section 2, Eqs. ( 10), ( 11), (13). In the following, these models will be referred as the V-model and the R-model model, respectively. The V-model is solved by making use of a phase space code based on a semi-lagrangian algorithm [29] that integrates the following set of equations. From Eq. ( 1), the 1D-2V Vlasov equation reads,

∂F ∂t + p x m e γ ∂F ∂x + e E x + p y B z m e γ ∂F ∂p x + e E y - p x B z m e γ ∂F ∂p y = 0 (46)
self-consistently coupled to the Maxwell's equations:

∂B z ∂t = - ∂E y ∂x (47) ∂E y ∂t = -c 2 ∂B z ∂x - J y ε 0 (48)
where the the current, in the transverse direction, is given by

J y (x, t) = -e ¨py m e γ F (x, p x , p y , t) dp x dp y (49)
The longitudinal component of the electric field E x is obtained from the Poisson equation, Eq. ( 11), where the electron density n(x, t) is given by

n(x, t) = ¨F (x, p x , p y , t) dp x dp y (50)
The R-model is integrated by solving the j-Vlasov equations as written in Eq. ( 10), with N = 2, self-consistently coupled with the Maxwell and Poisson equations. In this case, the fields are defined in Eqs. ( 12) and ( 14) for the electron density and the current density, respectively.

A set of numerical simulations has been performed for the case of the CFI in the relativistic regime. Both models are initialized by two counterpropagating streams. The initial distribution function is composed of a bi-Maxwellian distribution with beam momenta centered at p 01 /m e c = -3.0e y and p 02 /m e c = -(n 0,1 /n 0,2 ) p 01 /m e c. The d.f. is perturbed by purely transverse magnetic fluctuations (on the magnetic B z -component) with wave vector k 0 = k 0 e x , perpendicular to the electron beams. For the sake of simplicity, we have chosen a symmetric case of equal beam densities, n 01 = n 02 = 0.5n 0 , corresponding to p 02 /m e c = 3.0e y .

From now on, we use the electron mass and charge, the velocity of the light, the plasma frequency and the electron skin depth as characteristic quantities to normalize the equations. The magnetic and electric fields are normalized to B = Ē = m e cω p /e. The numerical space domain, in dimensionless units, is L x = 2π, such that the largest wave lengths are typically the most unstable and we excite initially only the largest one admitted in the numerical domain, corresponding to k 0 = 1. The initial magnetic perturbation, in dimensionless units, is given by:

B z = a 0 sin(k 0 x) (51)
where a 0 = 10 -4 is the (small) initial amplitude given in dimensionless units.

In the V-model, we take an initial electron plasma temperature T e = 2keV in both the p x and p y directions. In the R-model instead, only the p x direction has a non zero temperature and we keep the same temperature as in the Vmodel. The phase space sampling in the V-model is N x N px N py = 256 × 257 2 , i.e. 1.69 × 10 7 grid points. The phase space sampling in the R-model is N x N px = 512×513 grid points. The time step used in both simulations is ∆t = 0.003. In Fig. 3, we plot the time evolution of the magnetic energy resulting from the numerical solution of the V-model. During the linear phase, 10 < t < 25, the initial perturbation with wave vector k = 1 grows with a growth rate η V -model 0.47, in agreement with the theoretical value, η th = 0.47, obtained by solving the dispersion relation of the CFI, Eq. ( 30). The same curve (not shown here) has been obtained by using the multi-stream R-model initialized with two streams only, i.e. N = 2 in Eq. (30).

The relativistic dispersion relation of the CFI in a collisionless plasma [25] and its role as a driver of non linear dynamics leading to wave-breaking processes [26], has been studied many years ago in the framework of a fluidtype description involving two electron fluid populations. These results show that wave-breaking spontaneously develop during the early non linear phase of the process of generation of quasi-static magnetic fields. However, both processes in a collisionless, plasma, are saturated by particle trapping effects as shown, in the non relativistic limit, in Ref. [28], thus requiring a fully kinetic approach.

Our kinetic V-model approach confirms the occurrence of wave breaking, as shown in Fig. 4, top panel, where we plot the time evolution of the electron density at x = π/2 corresponding to the magnetic peak location. In the bottom panel, we show the same plot for the longitudinal electric field E x . We see that the instability has saturated, but showing the occurrence of several plasma wave breaking events characterized by strong density peaks reaching more than three times the initial density mean value.

The first wave-breaking scenario, characterized by the strongest density peak, is observed during the saturation phase, at t 34.3. Such event can be considered as the final explosive phase of the cavitation process started around t 25.

In Fig. 5, top to bottom, we show the spatial profiles at t = 30 of the longitudinal electric field E x , the electron density n e and the magnetic field B z in normalized units. We see the formation of two spikes (per perturbation wavelength) associated to a wave breaking process able to produce particle acceleration depending on the initial temperature. We underline that the kinetic description adopted here allows to follow in detail the non linear relativistic dynamics associated to the cavitation process. Indeed, it provides the physically correct saturation mechanism (electromagnetic particle trapping) for the instability as well as for wave breaking, thus stopping the local growth of the singularity. The impact of particle trapping on the d.f. shape is shown in Fig. 6 where we draw, at t = 30, the p y (p x ) averaged distribution function in the (x, p x ) ((x, p y )) phase space, top (bottom) panel. Similar to the electrostatic case, the d.f. filaments in the (x, p x ) phase space evolves more and more towards vortex-like configurations. However, here the electromagnetic nature of the trapping process is visible in the x -p y phase space, bottom panel. We see that also in p y the beams, concentrated around the electrostatic peaks, rotate around the maxima of the self-generated magnetic field, in agreement with non relativistic results presented in Ref. [13].

It is worth to underline that, even if the phase space dynamics evolution is quite complex, both beams maintain their initial transverse shape indicating that each beam keeps constant its initial "temperature", as a consequence of the invariance of the perpendicular canonical momentum. It is this Hamiltonian constraint that allows us to built the multi-stream model.

In Fig. 7 we show the electric field, the density and the magnetic field as in Fig. 5, but at a later time, t = 33. We see that the peaking process of the electric field and the density has saturated creating localized wider structures, but with a reduced amplitude, around the regions where the absolute amplitude of the magnetic field is maximum. These structures have the shape of strong density cavities centered at x 1 = π/2 and x 2 = 3π/2. The electric and magnetic field are concentrated inside such cavities where the maximum value of the density at the boundaries becomes several time (more or less three times) the initial density value. In Fig. 8 we show the (x, p x ) and (x, p y ) phase space at t = 33, analogous to the one shown in Fig. 6 . We see that the previous vortex structures, typical of a regime where particle trapping is dominated by the electrostatic longitudinal component, have now been replaced by a stretched pattern driven by the magnetic force. Indeed, the magnetic to the electrostatic ratio, of order one at t = 30, has increased by a factor of five during the cavitation process.

In Fig. 9 we show the electron density and the magnetic field in the x -t plane, still obtained using the V-model. We see, at saturation, t 30, the formation of two main density spikes located in x = 0, π. Then, the density spikes propagate at nearly constant velocity, v ≈ 0.37, leading to the formation of two cavities delimitated by the density spikes themselves. In the cavities we observe the concentration of the z-component of the magnetic field with a local plateau. We note that these processes repeat periodically and occur on a characteristic time scale of the order of the inverse of the plasma frequency. The electron non linear dynamics shown in Fig. 9, is a clear indication of the lack here of a new "meta stable" Vlasov equilibrium, as instead often observed in the electrostatic limit in the presence of particle trapping processes leading to the formation of vortex structures in phase space. In other words, the electron population here do not tend to relax towards any asymptotic Vlasov equilibrium, as for example the so called BGK states typical of a 1D-1V nonrelativistic electrostatic system where the ions, pushed by the "static" dipolar electric field, play a key role in the long time dynamics of the system [30]. Here instead, the system evolves on the electron time scale and it is not significantly affected by the ion's motion (as observed in the same simulation, but allowing the ions to move).

The multi-stream simulation (R-model)

We present now the results of a numerical simulation performed with the Rmodel using the same set of physical parameters as in the previous case of the V-model simulation, Section 5.1. All quantities are in dimensionless units as in the previous V-model case. First of all, we note that in the reduced Vlasov equation, Eq. ( 10), the index j is nothing but a label since no differential operator is carried out on the variable p y . What we actually do is to bunch together particles within the same stream j and let each stream evolve following the reduced Vlasov equation. The different streams are coupled through the Poisson's equation (11) and the Maxwell equations, Eqs. ( 47) and (48). As a result, the multi-stream model, that relies on the invariance of the perpendicular canonical momentum, gives a correct description of the plasma dynamics even using a small number of streams.

In Fig. 10 we show the "averaged" electron distribution function at t = 30, 33, 36 in the (x, p x ) phase space resulting from the contribution of both beams, j=1,2 f j (x, p x , t). These results must be compared with the analogous ones obtained by means of the V-model and presented in Figs. 6 and8, top panels. We see that the evolution of the d.f. iso-lines in phase space are very accurately reproduced by the R-model. However, the use of the R-model allows us to separate the specific contributions of each beam during the development of the CFI. In particular, let's now focus our attention on the dynamics of the second beam. In Fig. 11, we show the phase space distribution function of the second stream, f 2 (x, p x , t), at the same times as in Fig. 10, t = 30, 33, 36. By taking the zero moment of the distribution f 2 (x, p x , t), we plot in Fig. 12 the corresponding electron density, i.e. the density associated to the second electron beam.

In Fig. 11 we observe the initial electron stream acquiring a p x momentum component, comparable to the initial p y one, due to the rotation around the self-generated magnetic field. At the same time, local density (and electric) spikes are produced by large amplitude plasma fluctuations [START_REF] Davidson | Methods in linear and nonlinear plasma theory[END_REF], pumped by the electrostatic electric field component driven non linearly by the CFI (see Ref. [28], Fig. 6 in the non relativistic limit). These spikes last for a few plasma times at a given spatial location (see Fig. 12) before disappearing, but sufficiently long to accelerate bunch of particles at energies larger then the initial ones, eventually leading to the formation of filaments starting from the central bulk of the distribution function and winding up more and more in phase space (see Fig. 11, central and bottom panels). Thus the analysis of a specific beam (in this case the second one) allowed by the R-model shows the mechanism of particles acceleration as a consequence of wave breaking events. The particles then evolve in phase space due to the combined effects of the ballistic motion in the presence of electromagnetic forces.

A non-symmetric case for CFI % added sub-section %

In this section we give a last example concerning a non symmetric case of CFI instability. Identifying the magnetic field amplification mechanism and the resulting coupling with particles in the saturation regime is the key for a better understanding of CFI or Weibel-type instabilities. The small spatiotemporel scales that are accessible to global Vlasov-Maxwell simulations and the constraints on computational time and phase space sampling imply that, such global simulations can often not be applied to examine in detail both aspects of magnetic field amplification and the resulting wave-particle dynamics.

Our multi-stream model can give important insight into kinetic processes leading to the saturation process itself. By selecting appropriate initial particle "bunches" , the phase space structures can be isolated and examined in detail, allowing to separate the different physical processes linked to each beam. We consider a slightly different case where the system is initially non symmetric leading to a different structure of the wave breaking process. As first, using the V-model, we have carried out a numerical simulation with parameters n 0,1 = 0.4 and p 0,1 = 4 for the first beam and n 0,2 = 0.6 and p 0,2 = -n 0,1 p 0,1 /(n 0,2 γ 0,1 ) -0.85 for the second counter-propagating beam. Both streams have the same temperatures, namely T ex = T ey = 2keV . The phase space sampling is N x N px N py = 256 × 257 2 in the full Vlasov simulation (the V-model). The perturbation amplitude for the magnetic field is given by δB 0 = 10 -5 and we excite the same wave vector k 0 = 1.

In Fig 13 we draw the time evolution of the (longitudinal) electrostatic energy in top panel together with the corresponding evolution of the electron density measured at the middle of the box. As expected the electron density in Fig 13 reveals a successive of strong density peaks as the result of several wave breaking processes. The same simulation has been now carried out by means of the R-model and the corresponding time evolution of the electrostatic energy is plotted on top panel in Fig. 14. We observe a very good agreement with the numerical results of Fig. 13 obtained with the V-model. The oscillatory behaviour observed in the electrostatic energy takes place during the saturation phase of CFI. The bottom panel in Fig. 14 shows the standard evolution of the magnetic energy: in the linear phase the magnetic energy increases linearly (in a logarithmic scale) with a numerical growth rate of γ num 0.393 in very good agreement with the theoretical value of γ th = 0.398, followed by the regime of saturation till t ≥ 35. The nature of the wave-breaking changes going from the symmetric to the non-symmetric case. The occurrence of strong wave-breaking associated with the formation of two density spikes and the occurrence of a rotating structure in phase space in the fast stream; finally, a large cavitation in the density takes place for the slow electron population.

In Fig. 15 we focus on the shaded isocontours of the electron distribution function corresponding to the fast beam (i.e. with a transverse momentum of p 0,1 = 4 in normalized units). The phase space structures reveal that the growing magnetic fields modify the electron trajectories of the fast beam which now becomes magnetically trapped. In this last case we have used a finer x -p x phase space mesh grid with N x N px = 513 2 grid points. In Figs. 15 electrons make a global rotation between times t = 33.8 (the vertical arrow indicates here the "initial" position of the filament around the peak of density) and the final time t = 48.8. At that instant, the filament, located inside the magnetically trapping structure, exhibits a similar shape in comparison to the start time t = 33.8. Thus the corresponding magnetic bounce frequency is close, in dimensionless units, to ω B 0.40 (or equivalently to a bounce period of T B = 15.48, i.e. the time during three successive peaks in electron density in bottom panel in Fig. 13). The simulation is thus consistent with the magnetic trapping saturation theory as suggested by Davidson et al in [START_REF] Davidson | [END_REF]. In particular in magnetic trapping theory the instability reaches a steady state when the magnetic bounce frequency ω B becomes comparable to the linear growth rate γ th , which is the case here since γ th = 0.398. The magnetic bounce frequency is here, (see Ref. [33] ) in the relativistic regime of the instability, defined by

ω B ω p 2 = k 0 c ω p p 01 mγ 01 ω ce ω p (52)
where p 01 and γ 01 are respectively the transverse momentum and Lorentz factor of the first stream. Here ω ce = eB max /mγ 01 . Using p 01 = 4, γ 01 = 4.123 and B max = 0.70 we obtain an estimate of ω B 0.406 in good agreement with the numerical value obtained by the R-model. In Fig. 15 the first five curves of the distribution function correspond to the first period of the magnetic bounce motion, while the last curve, at time t = 150, shows the trapping structure after several bounce motions (after eight temporal periods). We observe that such a structure can persist over a long time.

CONCLUSION

A multi-stream model based on a Hamiltonian reduction technique, first developped in [12], has been applied to study CFI both in the linear and nonlinear regimes of the instability till saturation. The multi-stream model is a set of kinetic Vlasov-type equations obtained in a Hamiltonian framework, allowing to reduce the dimension of the phase space, but without loosing the accuracy of the system dynamics. The idea is to suppress the perpendicular momentum variable p ⊥ taking advantage from the invariance of the generalized canonical momentum, i.e. by noticing that, for a given j-population, the quantity P c,⊥ = p ⊥ + eA ⊥ = C j is constant in time when the q ⊥ variable is lacking. By approximating the system as a finite number of "streams" or, in other words, as a summation over an ensemble of Dirac distributions, we have obtained the linear dispersion relation of Weibel-type instabilities. This formulation gives also a clear description of the temperature effects in the non relativistic regime. Furthermore, by a direct comparison with the full Vlasov model, we have shown that our R-model provides a very accurate description of the non linear, relativistic regime in the case of the current filamentation instability. We underline that here, for the sake of simplicity, we have considered the case of two symmetric beams as a "typical" case for the CFI; however, the multi-stream model is not limited to such special configurations.

To conclude, the multi-stream model offers an exact description of the plasma dynamics even with a small number of "streams", providing a relatively simpler analytical framework and bringing the link between different types of instabilities, as for example the CFI and the Weibel. From a physical point of view, interesting results can be obtained even in the limit of a small number of "streams", e.g. N = 1, 2 or 3. On the other hand, for problems where a larger N is needed for an accurate description of a continuous distribution, the numerical effort becomes more and more comparable to that required for a standard Vlasov approach (defined as the V-model). However, we believe that the use of the multi-stream sampling should be more accurate because of the lack of numerical diffusivity in the reduced dimension direction as well as for the less important phase space filamentation effects that typically affect Vlasov Eulerian calculations.

The philosophy used in the multi-stream approach is reminiscent of the reduction procedure met in Hamiltonian theory, as for example the canonical realization of a Lie group, which is a subgroup associated with the Lie algebra composed of the ordinary Poisson bracket (see Ref [24]). The water-bag concept in Refs. [22,23] or the Hamilton-Jacobi model are other examples of Hamiltonian reduction techniques. While the water-bag concept makes full use of the Liouville invariance in phase space (we follow only the evolution of phase space contours using the property df /dt = 0), the Hamilton-Jacobi model is based on the property of existence of adiabatic invariants in action-angle variables when well separated time and space scales can be found. Furthermore, in the "multi-streams" model, there are no constraints on the shape of the distribution function which can be far from a Maxwellian one.

m 2 e v 2 thyFm 4 e v 4 th

 24 40)Assuming a symmetry property, i.e. C 2 = 0 and F 1 = F 3 , we get C 1 = -C 3 . Then, by using higher order momenta, ¨p2 y F (p) dp x dp y = n 0 (p) dp x dp y = 3n 0

Figure 1 :

 1 Figure 1: G vs. the perpendicular momentum p y as represented by Dirac masses in the R-model for the case of the counter filamentation instability, in dotted lines the real shapes of the distibution function.

Figure 2 :

 2 Figure 2: vs. the perpendicular momentum p y as represented by Dirac masses in the R-model for the case of the Weibel instability, we have represented in dotted line the shape of the distribution in the form of a Maxwellian.

Figure 3 :

 3 Figure 3: vs time for k = 1 using the V-model. The numerical growth rate, η V -model 0.470, is in very good agreement with the theoretical value, η th = 0.47. The physical parameters are p 0,1/2 = ±3 and n 0,1/2 = 0.5. All units are dimensionless.

Figure 4 :

 4 Figure 4: n e at x = π/2 vs. time. Bottom panel: same plot for the longitudinal electric field E x .

Figure 5 :

 5 Figure 5: E x , the electron density n e and the z-magnetic field component B z vs. x at time t = 30, top, middle and bottom panel, respectively, obtained with the V-model.

Figure 6 :

 6 Figure 6: (x, p x ) and (x, p y ) phase space (top and bottom panel) averaged in the p y and p x direction, respectively, at t = 30. These results have been obtained by using the 1D-2V Vlasov-Maxwell code (V-model) in the case of the current filamentation instability. The initial distribution condition is symmetric and is composed of two Maxwellian with beams moments centered at p 0,1/2 = ±3 e y .

Figure 7 : 29 Figure 8 :

 7298 Figure 7: E x , the electron density n e and the z-magnetic field component B z vs. x at time t = 33, top, middle and bottom panel, respectively, obtained with the V-model

Figure 9 :

 9 Figure 9: n e (top panel ) and the magnetic field B z (bottom ) in the x -t plane obtained by means of the V-model.

Figure 10 :

 10 Figure10: (x, px ) phase space representation of the quantity j f j (x, p x , t) where the summation is made over the two streams at t = 30, 33, 36. These results have been obtained by using the R-model in the case of the current filamentation instability.

Figure 11 :

 11 Figure 11: x, p x ) phase space of the second beam, j = 2, using the R-model at t = 30, 33, 36.

Figure 12 :

 12 Figure 12: j = 2, at t = 30, 33, 36.

Figure 13 :

 13 Figure 13: The electrostatic energy (top panel) and the corresponding electron density n e at x = π vs. time (on bottom panel) using the V-model. The physical parameters are n 01 = 0.40 and p 01 = 4 for the first stream and n 02 = 0.6 and p 02 = -0.86 for the second. The arrow indicates a magnetic bounce trapping period.

Figure 14 :

 14 Figure 14: The electrostatic energy vs. time on top panel, using the R-model. Bottom panel: the magnetic energy in logarithmic scale vs. time. The physical parameters are identical to those used in Fig. 13.

Figure 15 :

 15 Figure 15: x, p x phase space of the first stream, j = 1 (fast beam), using the R-model. The arrows indicate the position of the filament around the density peak during the rotation for the first magnetic bounce motion. Last curve, at time t=150, shows the magnetic trapping structure after eight magnetic bounce periods.

Acknowledgment

We acknowledge the IDRIS computational center, Orsay, France, for computer time allocation on their computers.