V T Tikhonchuk
email: tikhonchuk@u-bordeaux.fr

M Bailly-Grandvaux

J J Santos

A Poyé

Quasi-stationary magnetic elds generation with a laser-driven capacitor-coil assembly

Recent experiments are showing possibilities to generate strong magnetic elds on the excess of 500 T with high energy nanosecond laser pulses in a compact setup of a capacitor connected to a single turn coil. Hot electrons ejected from the capacitor plate (cathode) are collected at the other plate (anode) thus providing the source of a current in the coil. However, the physical processes leading to generation of currents exceeding hundreds of kiloamperes in such a laser-driven diode are not suciently understood. Here we present a critical analysis of previous results and propose a self-consistent model for the high current generation in a laser-driven capacitor-coil assembly. It accounts for three major eects controlling the diode current: the space charge neutralization, the plasma magnetization between the capacitor plates, and the Ohmic heating of the external circuit the coil-shaped connecting wire. The model provides the conditions necessary for transporting strongly super-Alfvenic currents through the diode on the time scale of a few nanoseconds. The model validity is conrmed by a comparison with the available experimental data.

I. INTRODUCTION

Schemes of strong magnetic eld generation with high energy laser pulses provide a compact setup for plasma magneto-hydro-dynamic experiments in an open geometry accessible for various diagnostics. In dierence from pulsed power sources of strong magnetic elds [START_REF] Pollock | High magnetic eld generation for laser-plasma experiments[END_REF][START_REF] Albertazzi | Production of large volume, strongly magnetized laser-produced plasmas by use of pulsed external magnetic elds[END_REF] evolving on a microsecond time scale and maintaining magnetic elds of an order of a few tens of Tesla, the laser driven sources operate on a nanosecond time scale and may create elds in vacuum exceeding hundreds of Tesla. Korobkin and Motylev [3] were the rst who suggested to use a laser for driving a current in a coil. With a relatively low energy (1 J) an long (20 ns) laser pulse they measured an electric current of 2 kA, producing in the center of a 1.4 mm loop a magnetic eld of 2 T. A possibility to upscale such a scheme was shown by Seely [START_REF] Seely | Pulsed megagauss elds produced by laser-driven coils[END_REF]. The rst large energy experiment was reported by Daido et al. [START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF]. The authors used a single loop, 1 mm diameter copper coil connected to two capacitor plates as schematically shown in Fig. 1. One plate (cathode) is irradiated by a laser pulse passing through a hole drilled in another plate (anode). In the experiment [START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF], the cathode was driven by a 1 ns 100 J pulse of a CO 2 laser at the wavelength of 10.6 µm. By using a loop antenna (B-dot probe) and extrapolating the data to the coil center, the authors reported a magnetic eld of 40 -60 T. They explained qualitatively that the current in the coil, attaining the value of 60 kA, is maintained by a high voltage of more than 200 kV generated in the capacitor by a thermo-electron emission from the laser heated cathode. That is, the capacitor operates as a laser-driven diode.

While no further details were provided concerning the target operation, it was mentioned that the maximum distance between the capacitor plates and thus the magnetic pulse duration are limited by the magnetization of the diode current. Adaptation of this scheme for the Nd laser at the wavelength of 1.06 µm was reported by Courtois et al. [START_REF] Courtois | Creation of a uniform high magnetic-eld strength environment for laser-driven experiments[END_REF]. The authors were using the Vulcan laser system with the pulse energy of 300 J and duration of 1 ns, and succeeded to generate a smaller magnetic eld (up to 28 T) in a copper loop of a diameter of 2 mm with the laser intensity of 4 × 10 16 W/cm 2 . They studied the magnetic eld dependence on the laser energy and pulse duration and compared single and double coil geometries, the latter allowing a better magnetic eld homogeneity. The authors have discussed also diculties in measuring the magnetic elds with coil probes positioned at a distance of a few tens of centimeters from the coil. First, extrapolation of the measured magnetic eld and current to the coil position could be a source of a large error. Indeed, the dipole magnetic eld decreases as a cube of the distance, and a dierence between the coil size and the probe position of a factor of 100, implies a dierence between the measured and estimated magnetic elds of 6 orders of magnitude. So even a small error in the measured eld may result in the error of a factor of a few times in the coil. Moreover, the authors noticed that the current in the diode is also a source of an electromagnetic eld, which could be of a compatible level and duration with the eld generated in the coil. Because of proximity of the capacitor and the coil, it is dicult to make a clear distinction between the signals coming from these two sources.

The detected magnetic eld presents in fact a pulse in the gigahertz domain, which has to be separated from other parasite signals generated in such a harsh laser plasma interaction environment.

A major breakthrough in the development of capacitor-coil targets has been presented by Fujioka et al. [START_REF] Fujioka | Kilotesla magnetic eld due to a capacitor-coil target driven by high power laser[END_REF]. With a 1 kJ laser pulse of the Gekko laser system at the wavelengths of 1.06 and 0.53 µm and 1 mm diameter coils, the authors reported magnetic elds exceeding 1.5 kT, that is, by raising the laser pulse energy by a factor of 3, they achieved an increase of the magnetic induction by a factor of 50, compared to Ref. [START_REF] Courtois | Creation of a uniform high magnetic-eld strength environment for laser-driven experiments[END_REF]. Apparently, there was a systematic error in the measurements, as that value corresponds to the energy stored in the magnetic eld comparable to the laser pulse energy. These data were acquired with two diagnostics: the pickup magnetic coil and the optical Faraday rotation in a piece of glass placed near the target. Very likely, the measurement error is related to a poor knowledge of the Verdet constant of the glass sample irradiated by a X-ray ash from the laser spot. These data were corrected in the subsequent experiment, carried out on the LULI2000 facility by Santos et al. [START_REF] Santos | Laser-driven platform for generation and characterization of strong quasi-static magnetic elds[END_REF],

where magnetic elds up to 600 T were reported for a laser energy of 500 J, but for a smaller coil, of a diameter of 0.5 mm. There, the authors used laser intensities above 10 17 W/cm 2 at the 1.06 µm wavelength, and magnetic B-dot probe measurements were compared with the data obtained from proton deectometry. This was the rst unambiguous measurement of the magnetic eld in the coil, clearly separating it from the capacitor. However, the proton data were consistent with the B-dot measurements only at the raising part of the current, up to 350 ps, yielding a magnetic eld of 100 T. The proton deections have decreased at later times, which was explained by electrostatic screening the coil magnetic eld by a plasma cloud emerging from the capacitor. This hypothesis was conrmed in a more recent experiment of the same group [START_REF] Law | Direct measurement of kilo-tesla level magnetic eld generated with laser-driven capacitor-coil target by proton deectometry[END_REF], where the proton radiography was obtained in the moment of the magnetic eld maximum thanks for shielding of the coil from the diode. Figure 2 shows the magnetic elds measured with dierent diagnostic methods on the laser installations LULI2000 and GEKKO-XII. All these data have however been collected in experiments with laser pulses of the energy less than 1 kJ and at wavelengths larger than 0.53 µm. Experiments in the domain of higher laser energies and shorter wavelengths of 0.35 µm were reported in Refs. [START_REF] Gao | Ultrafast proton radiography of the magnetic elds generated by a laser-driven coil current[END_REF][START_REF] Goyon | Ultrafast probing of magnetic eld growth inside a laser-driven solenoid[END_REF]. The authors of Ref. [START_REF] Gao | Ultrafast proton radiography of the magnetic elds generated by a laser-driven coil current[END_REF] used two Omega EP laser beams with the total energy of 2.5 kJ in a 1 ns pulse and with the intensity 3 × 10 16 W/cm 2 . The maximum magnetic eld measured in the center of an U-shaped loop with a curvature radius of 0.3 mm was at the level of 40 T. The measurements were accomplished with two diagnostics: the B-dot probes and proton deectometry. Moreover, by comparing the experiments with single and double loop coils the authors observed that the current in each coil remains the same.

That is the signature of the laser diode operating as a voltage source, that is, the current in the circuit is proportional to the external resistance, which has to be much larger than the internal resistivity of the diode itself. The eciency of converting the laser pulse energy into the magnetic eld in this experiment was about 0.01%, which is 10 times lower than the number reported in Ref. [START_REF] Courtois | Creation of a uniform high magnetic-eld strength environment for laser-driven experiments[END_REF] and more than 100 times less than in Ref. [START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF]. This fact indicates a strong dependence of the coil-capacitor target performance on the laser wavelength. More ecient generation of magnetic eld of 200 T on the same installation Omega EP and with a similar laser and target setup was reported recently [START_REF] Goyon | Ultrafast probing of magnetic eld growth inside a laser-driven solenoid[END_REF].

A large scatter of the magnetic eld measurements along with very basic qualitative estimates proposed in Refs. [START_REF] Seely | Pulsed megagauss elds produced by laser-driven coils[END_REF][START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF],

gave rise to discussions concerning the quality of the collected data and the value of maximum current that can be delivered by a laser driven diode [START_REF] Fiksel | A simple model for estimating a magnetic eld in laser-driven coils[END_REF]. In this paper we propose a more rigorous analysis of the capacitor-coil assembly operation and, in particular, an analysis of the performance of the laser driven diode. In Sec. II we recall the initial model [START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF] and further improvements proposed in Refs. [START_REF] Goyon | Ultrafast probing of magnetic eld growth inside a laser-driven solenoid[END_REF][START_REF] Fiksel | A simple model for estimating a magnetic eld in laser-driven coils[END_REF] and show their contradictions and limits. Then in Sec. III propose a more physically justied model that accounts for major processes taking place during the laser pulse: the space charge neutralization, the plasma magnetization and Ohmic heating of the coil wire. The major experimental results are recalled in Sec. IV and compared with the model. Section V contains one conclusions.

II. MODEL OF A VOLTAGE SOURCE

A. Simple model of the capacitor-coil system A model of operation of the laser driven capacitor-coil target was proposed in Refs. [START_REF] Seely | Pulsed megagauss elds produced by laser-driven coils[END_REF][START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF]. A simplied scheme with the corresponding notations is shown in Fig. 1. The tension V c at the coil edges is assumed to be given and the current I c in the external circuit is described by the equation

V c = R c I c + L c dI c /dt (1)
where R c is the coil resistance and L c is its inductance. The magnetic eld in the center of the coil B 0 is related to the current by the formula of the magnetic dipole, B 0 = µ 0 I c /2a, where a is the coil radius and µ 0 is the vacuum magnetic permeability. For the coil radius a = 0.5 mm, the current I c = 100 kA corresponds to the magnetic eld of 125 T in the coil center. These numbers provide the reference point for our further estimates.

In order to estimate the current and tension, let us consider the following representative example of a copper target.

The coil resistance R c = ηl w /s w and inductance L c µ 0 a ln(a/b) depend on the wire length l w 2πa, the cross section s w = πb 2 and the wire radius b. The copper resistivity at the room temperature is η Cu 2 × 10 -8 Ω•m. For b 30 µm and a = 1 mm, the coil resistance is about 0.03 Ω and the inductance is ∼ 3 nH. The corresponding relaxation time τ c = L c /R c ∼ 100 ns is larger than the typical laser pulse duration, t las ∼ 1 ns. That is, during the laser pulse, the inductance dominates and the coil current grows linearly with time, assuming a constant tension.

Later in time when the capacitor is short cut and the laser is switched o, the current decreases exponentially with the relaxation time τ c . This behavior corresponds qualitatively to the observations presented in Ref. [START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF], although the relaxation time has to be corrected accounting for the wire Ohmic heating. The value of the maximum tension generated in the diode is therefore dened by the laser pulse duration, V c L c I c /t las ∼ 300 kV.

The main problem is to explain the generation of such a high tension in the diode at the nanosecond time scale.

The system is described in [START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF][START_REF] Courtois | Creation of a uniform high magnetic-eld strength environment for laser-driven experiments[END_REF] as an empty capacitor with the capacitance C d = 0 S d /d, where 0 is the vacuum dielectric permittivity, S d ∼ 10 mm 2 the surface of the capacitor plate and d 1 mm is the distance between the plates. The diode capacitance is very small, C d 0.1 pF, so it can be rapidly charged to a high voltage even by a small current. Indeed, the tension of 300 kV in such a capacitor corresponds to the charge Q d = C d V c ∼ 30 nC. It could be charged with a 30 A current during the laser pulse time of 1 ns.

The charging current in the diode, I d , is provided by emission of hot electrons from the laser heated cathode.

According to Ref. [START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF], these electrons are generated due to the resonance absorption of the laser pulse energy at the target surface. For the hot electron density n h and the distribution function f h , the charging diode current writes:

I d = eπr 2 h n h ∆Ω ∞ eVc dε v f h (ε), (2)
where r h is the radius of the emission zone, which can be assimilated to the radius of the laser spot, ∆Ω is the solid angle of the hot electron emission, v and ε are the electron velocity and kinetic energy. The low limit in the integral denes the minimum electron energy needed to access the anode having the potential -V c . (We assume the cathode potential to be zero.) Then the diode tension is dened by the following relation:

C d V c = t dt I d (t ). (3)
By solving this relation one nds the tension V c (t), which then can be injected in Eq. ( 1) for calculation of the external current I c . Similarly to the standard voltage source, the diode tension is independent on the external current and is dened by the internal process in the diode. Assuming a Maxwellian distribution of hot electrons with a temperature T h , electron emission in the solid angle ∆Ω ∼ π and replacing the time integration in Eq. ( 3) with a product of the characteristic current by the laser pulse duration, one obtains the following estimate for the diode potential:

eV c T h ln ω 2 ph t las d v h r 2 h S d (4)
where ω ph = (e 2 n h / 0 m e ) 1/2 is the plasma frequency associated with the hot electron density n h and v h = (T h /m e ) 1/2 is the hot electron thermal velocity. The term inside the logarithm can be easily of the order of 10 6 or more, which corresponds to the tension of 15-20 times the hot electron temperature, T h . So, for a modest hot electron temperature of 15 -20 keV, one obtains a tension of the order of 200 -300 kV in good agreement with the observations. However, these rapid estimates are oversimplied and contain a caveat, which has been overseen by the authors of Refs. [START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF][START_REF] Courtois | Creation of a uniform high magnetic-eld strength environment for laser-driven experiments[END_REF]. Indeed, the high ratio eV c /T h 15 implies a low current 100 A and a low charge in the capacitor, Q d 100 nC. In contrast, the high discharge current I c ∼ 100 kA transports during the same time a charge thousand times larger. This strong charge imbalance is conrmed also by energy considerations: the estimated energy stored in the induction magnetic eld in the coil is more than thousand times larger than the electrostatic energy stored in the capacitor. A large energy needed to feed the coil current cannot be stored in the empty capacitor, and the separate analysis of the internal and external circuits leads to inconsistent results. Two solutions for this problem have been proposed: either suppose that the system operates in a quasi-steady regime, where a continuous charge ow through the diode is maintained by the laser pulse energy deposition, or assume that the diode capacitance is much larger than estimated above. The rst hypothesis proposed in Ref. [START_REF] Fiksel | A simple model for estimating a magnetic eld in laser-driven coils[END_REF] implies that the diode capacitance is small and the internal current is equal to the external current, I c = I d . The second hypothesis proposed in Ref. [START_REF] Goyon | Ultrafast probing of magnetic eld growth inside a laser-driven solenoid[END_REF] implies that the capacitor is lled with a suciently dense plasma screening the electrostatic eld.

According to Ref. [START_REF] Goyon | Ultrafast probing of magnetic eld growth inside a laser-driven solenoid[END_REF], the capacitor is lled with a plasma produced by laser ablation from the cathode. Then the distance between the plates d in the expression for the capacitance C d = 0 S/d has to be replaced by the plasma Debye length λ D = ( 0 T e /e 2 n e ) 1/2 , where T e and n e are the plasma temperature and density. The needed capacitance can be obtained assuming the ratio d/λ D to be of the order of 10 4 . This hypothesis unfortunately is not self-consistent:

First, the needed plasma density is of the order of a few percent of the critical density, so it would take a time of a few hundred picoseconds to ll the total capacitor volume, while the current is generated in a shorter time scale.

Second, such a dense plasma would absorb incident laser radiation thus reducing the energy ux to the cathode and thus the available current. Third, a large diode capacitance implies a large charging current, which cannot be larger than the discharge current. Then one should suppose that the charging and discharging currents are approximately equal, which brings us back to the rst hypothesis. Fourth, such a dense plasma is a good conductor and it would short cut the circuit by providing the return current through the diode. As it is shown in the next section, the current continuity condition, I c = I d , can be satised without requiring a large diode capacitance.

B. Plasma expansion from the cathode

The current continuity is an important propriety of the diode operation. Another important feature, already mentioned in the earlier paper [START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF], is the plasma expansion from the cathode under the energy deposited by the laser.

The plasma isothermal expansion is described by the one dimensional self-similar model [1315]. The plasma density has an exponential spatial prole n h (z) = n h0 exp(-z/c s t) extending from the initial position z = 0 to the edge z f (t).

Here, n h0 is the initial density of hot electrons, c s = (ZT h /m i ) 1/2 is the ion acoustic velocity, Z and m i are the ion charge state and mass. Position z f of the edge of expanding plasma is dened by the condition of quasi-neutrality, that is, the plasma local scale length z f has to be equal to the local Debye length λ D (z f ) = [ 0 T h /e 2 n h (z f )] 1/2 . Solving the relation λ D (z f ) z f , one nds z f (t) and the plasma expansion velocity, v f = żf : z f 2c s t ln(z f /λ Dh ) and v f 2c s ln(2c s t/λ Dh ).

(

) 5
where λ Dh = ( 0 T h /e 2 n h0 ) 1/2 is the hot electron Debye radius. The expanding plasma is carrying an electric eld eE z T h /(c s t). Multiplying this value by the plasma thickness z f , one obtains the following expression for the potential drop in the diode, eV f -2T h ln(z f /λ Dh ). That is, the potential retaining the electrons at the cathode increases logarithmically with time thus limiting the current owing through the diode.

The authors of Ref. [START_REF] Fiksel | A simple model for estimating a magnetic eld in laser-driven coils[END_REF] proposed, however, another estimate for the potential drop in the diode eV ∼ -T h d/(c s t). This estimate could be obtained by multiplying the electric eld E z by the distance between diode plates d. However, this cannot be the case as the one-dimensional self-similar model is limited to short expansion times and to distances much smaller than the distance between the diode plates. As soon as the plasma size z f v f t becomes comparable with the laser spot radius, r h , the non-stationary planar plasma expansion transforms into a stationary spherical expansion, and the logarithmic time-dependent term in Eq. ( 5) has to be replaced by a constant value ln(r h /λ Dh ) [START_REF] Robson | Scaling of proton acceleration driven by petawatt-laser-plasma interactions[END_REF]. Therefore, at the distance z ∼ r h from the cathode the plasma expansion velocity and the potential jump are stabilized at the levels v m 2c s ln(r h /λ Dh ) and eV m -2T h ln(r h /λ Dh ).

For the typical hot electron temperature T h ∼ 20 -40 keV and the hot electron density n h0 ∼ 10 21 cm -3 , the typical value of the ion acoustic velocity is c s 1 mm/ns and the Debye length λ Dh is 0.1 µm. Then for the typical value of the laser focal spot r h 50 µm, the one-dimensional plasma expansion is limited to a narrow zone of less than 100 µm near the cathode, and the non-stationary expansion ends in less than 10 ps.

The important (but intermediary) conclusion following from that estimate is that ions at the edge of the planar expansion zone have already acquired a high velocity v f (10 -15)c s 10 -15 µm/ps. That corresponds to the ratio of the ion energy to the hot electron temperature ε i /ZT h 2 ln 2 (r h /λ Dh ) which is of the order of 25 -50. Assuming that the ion acceleration is stopped at the distance r h from the cathode, we conclude that the hot plasma arrives at the anode at the time t s d/v f ∼ 100 ps, and after that moment the diode is operating in an approximately stationary regime throughout the laser pulse duration, which is of the order of 1 ns. Thus the diode is indeed lled with a tenuous plasma, while supplying the current in the external circuit.

III. MODEL OF THE LASER DRIVEN DIODE

In this section we propose a model for the laser-driven diode operation that describes the electrostatic potential distribution inside the capacitor and the current that can be supplied to the external circuit. We start with the characterization of the laser-driven electron population, describe the potential distribution assuming a stationary operation and then account for two complementary eects: the change of the parameters of the external circuit due to the Ohmic heating and magnetization of the current inside the diode.

A. Hot electron characterization

The laser radiation interacting with the cathode plate generates hot electrons. In the intensity domain of interest, I las λ 2 las 10 14 Wµm 2 /cm 2 , the major source of hot electrons is the laser resonance absorption. The theoretical analysis, numerical simulations and comparison with experiments [START_REF] Forslund | Theory of hot-electron spectra at high laser intensity[END_REF] suggest the following expressions for the density of hot electrons and their temperature:

n h0 /n c 0.2(I las λ 2 las /T c ) 0.5 , T h 9(I las λ 2 las ) 0.25 keV. (7)
Here, n c = 1.1 × 10 21 λ -2 las cm -3 is the electron critical density, T c is the temperature of the cold background plasma expressed in keV, the laser intensity I las is expressed in the units of PW/cm 2 and the laser wavelength λ las is in microns. These expressions are valid however for the laser irradiances I las λ 2 las 10 16 Wµm 2 /cm 2 . According to the numerical simulations [START_REF] Estabrook | Properties of resonantly heated electron distributions[END_REF] and experimental measurements [START_REF] Priedhorsky | hard X-ray measurements of 10.6 µm laser-irradiated targets[END_REF], the hot electron temperature scales for higher laser irradiances as T h 12(I las λ 2 las ) 0.42 keV,

where the units are the same as in Eq. [START_REF] Fujioka | Kilotesla magnetic eld due to a capacitor-coil target driven by high power laser[END_REF]. The correspondent hot electron density can be estimated by assuming a power balance, hI las n h0 v h T h , where v h = (T h /m e ) 1/2 is the hot electron average velocity and h is the fraction of laser energy transferred to electrons. According to estimates presented in Ref. [START_REF] Goyon | Ultrafast probing of magnetic eld growth inside a laser-driven solenoid[END_REF] it varies in the range from 0.01 to 0.1.

B. Potential distribution in the diode

As it is explained in the previous section, the space inside the diode can be divided in two parts: (i) the plasma acceleration zone corresponding to a quasi-planar plasma expansion near the cathode, where the density decreases exponentially, and (ii) the transport zone corresponding to a quasi-spherical expansion at distances from the cathode larger than the laser focal spot size. A transition between these zones takes place at the distance z r h , where the plasma density can be estimated as n f n h0 (λ Dh /r h ) 2 . For the typical ratio r h /λ Dh ∼ 10 2 , the plasma density at the edge of the acceleration zone is already quite low, n f 10 -4 n h0 and the potential drop given by Eq. ( 6) exceeds the hot electron temperature already by a factor 10 or more. Assuming the cathode potential to be equal to zero and the hot electrons temperature T h ∼ 20 keV, the potential at the edge of acceleration zone V m would be of the order of -200 kV or more. The corresponding potential and density distributions are shown schematically in Fig. 3 by the dashed lines.

Knowing the plasma density and electron temperature, one can estimate also the value of the current in the diode.

As it operates in a quasi-stationary regime, the current I d is the same at any position, and it can be evaluated at the edge of the zone of planar expansion I d πr 2 h en f v h πλ 2 Dh en h0 v h . Thus, the electron current depends only on the hot electron temperature Estimating the diode impedance Z d = V f /I d , one nds quite a large value:

I d π 0 T h v h /e. (9
Z d 2 π 0 v h ln r h λ Dh = 240 Ω × c v T h ln r h λ Dh . ( 10
)
It is evident that the diode impedance is two-three orders of magnitude larger than the impedance of the external circuit, and the laser-driven diode current is strongly limited by its internal impedance.

Such a conguration, however, cannot be stationary. Indeed, the potential at the anode V c is dened by the impedance of the external circuit, Z c R c + L c /t las ∼ 1 Ω. According to the hypothesis of a quasi-stationary diode operation, I d = I c , the absolute value of the anode tension V c = Z c I d is much smaller than |V m |. Therefore, the potential inside the diode has a deep minimum at z r h . It decreases from zero value at the cathode to V m and then increases back to V c in the transport zone z > r h (see dashed curves in Fig. 3). Consequently, the electric eld in the diode changes sign and the ions accelerated in the zone of planar expansion would be decelerated in the transport zone. However, the plasma accumulation in the transport zone would then aect the potential distribution. It will evolve so to reduce the ion over-acceleration and make their distribution smoother. The only stationary conguration possible is the one corresponding to a monotonous prole of the electric potential and the ion velocity. It is shown in Fig. 3 by the solid lines: the potential drop in the plasma acceleration zone is connected to the anode with a plateau in the transport zone. Such a potential distribution is in agreement with the general theory of a plasma lled diode originally proposed by Pierce [START_REF] Pierce | Limiting stable current in electrons beams in the presence of ions[END_REF] and further developed in Refs. [2123]. This conguration is the major result of our analysis: the potential at the end of the acceleration zone has to be equal to the anode potential, -V c = V m . The value of the potential jump -V c in the plasma acceleration zone is dened by the plasma distribution in the transport zone. The hot electron density decrease in the acceleration zone is limited by the value n hc = n h0 exp(-eV c /T h ) at the edge of the transport zone z = z m , which is needed to maintain that potential jump. Then the electron current that is able to pass through that potential barrier can be estimated

as I d eπr 2 h n hc v h = I 0 exp[-eV c (t)/T h ], where I 0 eπr 2 h n h0 v h (11)
denes the maximum current that can be extracted from the diode. According to the hot electron scaling [START_REF] Fujioka | Kilotesla magnetic eld due to a capacitor-coil target driven by high power laser[END_REF], that value could be easily 10 MA or more. However, the current in the circuit is much smaller. It depends on the laser pulse duration and the impedance of the external circuit Z c . The current-voltage characteristic of the diode is dened by relation [START_REF] Goyon | Ultrafast probing of magnetic eld growth inside a laser-driven solenoid[END_REF] along with the condition of the current continuity,

I d = I c , I c (t) = I 0 exp[-eV c (t)/T h ] (12)
For a crude estimate of the diode current one can substitute the tension in the external circuit by a simplied relation, V c ∼ Z c I c . Solving then Eq. ( 12) in the limit I 0 T h /eZ c one nds:

I c (T h /eZ c ) ln(eZ c I 0 /T h ). (13)
Therefore, the dimensionless parameter controlling the diode current is eZ c I 0 /T h . It can be easily of the order of 10 2 -10 3 . That is, for the hot electron temperature T h ∼ 30 -40 keV and the impedance of Z c ∼ 1 Ω, the tension in the diode could be of the order of 7 T h /e 200 kV and the current of the order of 200 kA. That crude estimate shows a possibility to generate suciently high currents in a diode lled with a tenuous plasma. The potential distribution in our model is quite similar to the one proposed by the authors of Ref. [START_REF] Goyon | Ultrafast probing of magnetic eld growth inside a laser-driven solenoid[END_REF], but the required plasma density in the diode is much lower, and the plasma lls not the whole volume between the anode and cathode, but only the zone where current is transported. However, such a simplied model is not yet sucient for applications. It does not account for such physical eects as the temporal evolution of the parameters of the external circuit and the magnetic eld created in the diode.

C. Time dependent characteristics of the external circuit

Large currents owing in the circuit induce the time variation of the coil parameters. The magnetic stress force I c Ba ∼ 10 5 -10 6 N acting on the coil with a mass m w 0.1 mg corresponds to the displacement velocity ∼ 1 -10 µm/ns, which is insucient to change the coil shape. Therefore, the coil inductance can be considered as a constant. This estimate is in agreement with the conclusion of Ref. [START_REF] Goyon | Ultrafast probing of magnetic eld growth inside a laser-driven solenoid[END_REF] based on the numerical modeling of the coil dynamics.

The Ohmic heating of the coil with the current of the order of 100 kA corresponds to the energy release of the order of 1 J, assuming the wire resistance of 0.1 Ω. In order to estimate the wire temperature one should know the heat capacity of the copper C v = 0.385 J/gK and the latent heat of the phase transitions, melting and vaporization at the temperatures ∼ 1380 and 2860 K, respectively. The melting enthalpy H m 200 J/g is relatively small, and it is sucient to deliver 20 mJ to melt a wire with a mass of 0.1 mg. Heating to the vaporization temperature would require about 100 mJ. Both numbers are small compared to the expected energy release of 1 J. In contrast, the vaporization enthalpy H v 4600 J/g is large, and the vaporization energy is comparable with the delivered energy.

Therefore, one may expect the wire temperature to increase to the level of the order of vaporization temperature, which is approximately 3000 K, that is, ten times the room temperature. The wire resistivity η being linear function of the temperature is then increasing by a factor of 10 with respect to the estimate made in Sec. II A for a cold material. Similar arguments apply also to other materials, gold and nickel, used in the experiments [START_REF] Santos | Laser-driven platform for generation and characterization of strong quasi-static magnetic elds[END_REF][START_REF] Law | Direct measurement of kilo-tesla level magnetic eld generated with laser-driven capacitor-coil target by proton deectometry[END_REF][START_REF] Goyon | Ultrafast probing of magnetic eld growth inside a laser-driven solenoid[END_REF]. Optical shadowgraphy (at 532 nm) of the coil expansion. On the left, 0.65 ns-gated images corresponding to a cold coil (top, t = 0) and two dierent timings after laser-driving of the capacitor-coil target. On the right, the ns-scale expansion is resolved in a single shot by imaging the coil diameter into the slit of a streak camera. The average expansion velocity of the coil rod was measured to be 11 ± 3 µm/ns (as indicated by the dashed orange line). Data acquired on the experimental campaign [START_REF] Santos | Laser-driven platform for generation and characterization of strong quasi-static magnetic elds[END_REF].

The wire radial expansion also contributes to the wire temperature limitation. Shadowgraphic imaging of the coil during and after the laser pulse performed in the LULI2000 experiment [START_REF] Santos | Laser-driven platform for generation and characterization of strong quasi-static magnetic elds[END_REF] and shown in Fig. 4 indicates that the wire expands due to the Ohmic heating induced by the current. The wire radius increases with the velocity ∼ 11 µm/ns.

The energy transmitted to the wire expansion can be estimated assuming that the velocity is increasing linearly with the radius. Than the wire kinetic energy can be written as m w v 2 w /8, which is ∼ 1 J for the radial velocity v w 10 µm/ns. This estimate provides an additional argument for considering the gure of 3000 K as a reasonable estimate of the wire temperature.

Another important issue is the skin eect. The nanosecond time scale of current variation is too short for a homogeneous distribution of the current in the wire. The current is localized in the skin layer of a thickness l s (ηt las /µ 0 ) 1/2 , which is of the order of 4 µm for a cold copper. This thickness is about 10 times less than the wire diameter. The current is thus localized at the wire surface leading to faster heating and surface ablation. However, as soon as the wire is heated to the temperature about 3000 K, the skin layer thickness is increasing by a factor of 3 and becomes comparable with the wire thickness. Therefore, for large currents corresponding to the maximum of magnetic eld generation, the current distribution can be considered as suciently homogeneous in the wire. Nevertheless, the wire core remains cold thus assuring the mechanical stability of the coil during the discharge time.

D. Super-Alfvenic currents in the laser driven diode

In the experiments with the laser driven diode, the size of the laser spot at the cathode is rather small of the order of 30 -50 µm. This is necessary for obtaining a hot plasma and producing a large amount of energetic electrons. However, this size is more than 20 -30 times smaller than the typical distance anode-cathode. Thus, the eects of transverse forces on the plasma dynamics and electric current have to be evaluated.

Let us consider the plasma in the transport zone inside the diode as a column with a given mass ow and a given electric current. Its radius r p is of the order of the laser spot radius at the cathode, and it increases with the distance. According to the stationary model of the plasma expansion presented in Section III B, at the distance z ∼ r h from the cathode, where the spherical expansion is formed, the hot electron density decreases to the level n hc n h0 exp(-eV c /T h ) and the ions are accelerated in the expanding plasma to the velocity v i c s eV c /T h , which is of the order of 3 -5 mm/ns. The corresponding ion kinetic energy

ε i = 1 2 m i v 2 i 1 2 ZT h eV c T h 2 (14)
could attain the value of a few hundred keV, about 10 times larger than the hot electron temperature. We assume in what follows a constant axial ion velocity v i and a constant mass ow πr 2 p n hc m i v i .

As plasma propagates from the cathode to the anode its radius changes under the eect of two forces: the hot electron pressure gradient, dP e /dr, and the Lorentz force, j × B. Assuming a constant hot electron temperature and the maximum plasma density at the axis, the pressure gradient is directed toward the axis, thus creating a radial electric eld, E r T h /er p . This eld induces a plasma expansion. Another contribution to the radial electric eld is related to the electric current. The azimuthal magnetic eld at the plasma edge can be estimated as B θ µ 0 I d /2πr p . The radial Lorentz force ev h B θ acts on electrons pushing them towards the axis. That also corresponds to a radial electric eld directed to the axis. The total radial electric eld acting on ions can be estimated as a sum of these two contributions,

E r -v h B θ + T h /er p .
One can estimate the current needed to fully compensate the hot electron pressure from the force balance, ev h B θ T h /r p . These two forces are equal, if

I d I A v h /2c
, where I A = 4πm e c/eµ 0 = 17 kA is the Alfven current. This condition does not depend on the plasma column radius. Consequently, super-Alfvenic currents are pinching the plasma column by pushing the ions towards the axis.

The evolution of the plasma radius with the distance inside diode can be estimated by assuming a small radial plasma velocity v r = dr p /dt v i . Combining then the equation for the ion acceleration m i dv r /dt = ZeE r with the relation for the ion axial displacement dz = v i dt, one obtains an equation for the radius of plasma column in the following form:

d 2 r p dz 2 - g r p where g = T h eV c 2 2 I d c I A v h -1 . (15)
Two terms in the parenthesis in the expression for g account for the Lorentz force and the electron pressure, respectively. The ion energy is related to the anode potential by Eq. ( 14). This equation can be integrated with the initial condition dr p /dz = α 0 near the cathode, z z m d, where r p = r h . Then we nd a relation between the plasma divergence angle and the radius:

dr p dz 2 = α 2 0 -2g ln r p r h . (16)
According to this equation, the plasma radius increases with the distance from cathode from the initial value r h at z r h to the maximum value r max = r h exp(α 2 0 /2g) and then decreases back to zero. The value r p = 0 corresponds to the collapse of the plasma column. The distance where the collapse takes place corresponds to the maximum distance of the current transport. The collapse length z coll is dened by the integral of Eq. ( 16):

z coll = r h exp(α 2 0 /2g) 0 dx α 2 0 -2g ln x .
Equating the collapse distance to the distance between the diode plates, we nd a condition for the maximum current that can be transported through the diode:

2g/π = (r h /d) exp(α 2 0 /2g).

(

) 17
This equation relates the coecient g [START_REF] Yu | Ion acceleration in expanding multispecies plasmas[END_REF] to the geometric characteristics of the diode: the plasma ow divergence α 0 and the ratio d/r h . In the domain parameters of interest α 0 2 and d/r h 20 this function is shown in Fig. 5.

It increases with the plasma divergence and decreases with the diode width remaining below 1 in the whole domain.

Considering now the denition of g in Eq. ( 15), we nd then the expression for the maximum diode current:

I p 0.5I A v h c 1 + g eV c T h 2 . ( 18
)
The maximum diode current is thus a quadratic function of the tension. In the domain of low tensions, eV c T h , it is limited by the Alfven current. For high tensions, eV c > T h the maximum current can be increased above this value if the coecient g is suciently large. By increasing the initial plasma divergence α 0 one may signicantly increase the maximum current. For example, for d/r h = 10 the coecient g increases from 0.03 for α 0 = 0 to 1 for α 0 = 2.

Thus, while the current magnetization imposes a limitation on the maximum current, it may exceed several times the Alfven current in the considered geometry. The crucial parameters in that relation are the ratio of the electrostatic potential to the hot electron temperature and the plasma stream divergence. By increasing these parameters one enforces the rigidity of the plasma ow and correspondingly its capacity to transport larger electron currents.

IV. DISCUSSION

A. Current-voltage characteristic

Combining the current-voltage characteristics related to the space charge limit (12) I sc (V c ) and the plasma magnetization [START_REF] Estabrook | Properties of resonantly heated electron distributions[END_REF] I pm (V c ), we can dene the diode current-voltage characteristic as a minimum between these two Typically, I * is of the order of the Alfven current and V max is a few times larger than the hot electron temperature. For I c > I * the diode current is limited by the cathode potential jump. Three panels in Fig. 6 show the variation of the current generation conditions in function of the laser intensity and wavelength. In the experiment [START_REF] Daido | Generation of a strong magnetic eld by an intense CO2 laser pulse[END_REF] (panel a) a CO 2 laser at the wavelength λ las = 10.6 µm was delivering to a copper target a 100 J pulse of 1 ns duration to a focal spot of r h = 120 µm. The intensity of 1.3 × 10 14 W/cm 2 corresponds to the hot electron temperature of 40 keV and the hot electron density n h0 /n c 0.1. The model shows that the maximum tension is expected to be about 4T h , that is 160 kV, and the maximum current could attain 6I A , that is 100 kA. However, because of the laser pulse duration shorter than the relaxation time τ c ∼ 10 ns, the current is limited to about 40 -50 kA. These numbers are in qualitative agreement with the measurements.

In the experiment [START_REF] Courtois | Creation of a uniform high magnetic-eld strength environment for laser-driven experiments[END_REF] (panel b) a Nd glass laser at the wavelength λ las = 1.05 µm was delivering to a copper target a 300 J pulse of 1 ns duration to a focal spot of r h = 30 µm. The intensity of ∼ 10 16 W/cm 2 corresponds to the hot electron temperature of 30 keV and the hot electron density n h0 /n c 0.07. The model shows that the maximum tension is expected to be about 5T h , that is 150 kV, and the maximum current could attain 6I A , that is 100 kA. However, because of a short laser pulse duration and a relatively large coil inductance of ∼ 4 nH, the current is limited to about 1.3I A , that is 20 kA. These numbers are also in qualitative agreement with the measurements.

In the experiment [START_REF] Gao | Ultrafast proton radiography of the magnetic elds generated by a laser-driven coil current[END_REF] (panel c) two 1 ns laser beams at the third harmonic of a Nd glass laser at the wavelength λ las = 0.35 µm were delivering 2500 J to a focal spot of r h = 60 µm. The intensity of ∼ 2 × 10 16 W/cm 2 corresponds to the hot electron temperature of 18 keV and the hot electron density n h0 /n c 0.045. The model shows that the maximum tension is expected to be about 7T h , that is 120 kV, and the maximum current could attain 6I A , that is 100 kA. However, because of a short laser pulse duration and a relatively large coil inductance of ∼ 2 nH, the current is limited to about 1.6I A , that is 30 kA. These numbers are also in qualitative agreement with the measurements.

A quantitative comparison would require to account for the temporal laser pulse prole, the coil heating, resistivity evolution and other factors. The present model is aiming at a qualitative evaluation of the major physical processes in the diode and a possible optimization of its performance. Figure 6 demonstrates the validity of our model for a large domain of laser wavelengths and intensities.

B. Temporal prole of the diode current

Knowing the current-voltage characteristic one can solve the external circuit equation ( 1) and nd the temporal prole of the current that can be achieved in the coil. Figure 7 shows the temporal proles of the current in the coil and the tension for the parameters of the experiment by Santos et al. [START_REF] Santos | Laser-driven platform for generation and characterization of strong quasi-static magnetic elds[END_REF]. Here, the laser pulse at the wavelength λ las = 1.05 µm of energy 500 J and duration 1 ns was delivered to a nickel target at the focal spot of a radius a = 15 µm achieving intensity of ∼ 10 17 W/cm 2 . The coil inductance is 2 nH, the resistance of a cold wire is 0.1 Ω and the distance between diode plates is 0.85 mm. The corresponding hot electron temperature is 84 keV and n h0 /n c 0.22. In that simulation the coil resistance was increased to R c = 1 Ω accounting for the wire heating as discussed in Section III C. The corresponding relaxation time τ c = L c /R c 2 ns is twice longer than the laser pulse duration. The current shown in Fig. 7a continuously grows with time attaining the value of 110 kA at the end of the laser pulse.

It is smaller than the value of 250 -300 kA claimed in [START_REF] Santos | Laser-driven platform for generation and characterization of strong quasi-static magnetic elds[END_REF], but that dierence is acceptable accounting to the model simplications. The tension stays approximately constant at the level ∼ 300 kV, which corresponds to the regime of current limitation by the cathode potential jump. The circuit impedance is decreasing with time approaching the value of the coil resistance, which was maintained constant in the present model. According to Eq. ( 13), the dimensionless parameter controlling the current saturation in this regime is the ratio of the hot electron temperature to the product of the coil impedance to the maximum diode current, T h /eZ c I 0 . So the current can be further increased by increasing the hot electron temperature, that is, the laser irradiance, I las λ 2 las , and the impedance of the external circuit. The former implies increasing the laser pulse energy as the intensity is already dened.

V. CONCLUSIONS

The analysis of capacitor-coil system presented in this paper allows to identify its mode of operation and the major parameters controlling the current intensity and duration. Considering this system as an electric circuit, we conclude that the optimal regime is achieved if the internal and external impedances are comparable, and with a nanosecond laser pulse the diode operates in the quasi-stationary regime because of its small capacitance. The diode operation time can be optimized by adjusting three characteristic times: the laser pulse duration, the response time of the external circuit (corrected by the eect of Ohmic heating) and the time of plasma propagation through the diode.

The latter one depends on the ratio of the size of the laser focal spot and the distance between the capacitor plates.

The important point is the time delay in the discharge current, which is equal to the time needed for fast ions to propagate through the diode. It is demonstrated that the charge accumulation in the Debye sheath limits the vacuum diode current to a low value, and only neutral plasma is capable to maintain large electron currents for a suciently long time. The maximum diode tension for the diode current smaller than the Alfven current is limited by the current magnetization. The higher diode currents are limited by the internal impedance of the diode. The present model neglects physical eects related to the cold plasma expansion from the cathode and anode and the return current in the diode if the plasma density increases. Analysis of these eects is out of the scope of the paper. They appear on a time scale of a few nanoseconds and may eventually short cut the current in the external circuit.

The maximum current delivered by the diode is proportional to the hot electron temperature. The major control parameter is the laser irradiance I las λ 2 las . The use of high laser intensities (tight focusing) and long wavelengths may improve the diode performance providing more energetic electrons capable to maintain a high current at large potential jumps. The capacitor-coil setup is limited to the currents in the range up to a few hundred kiloamperes and to durations of a few nanoseconds. A systematic experimental study of the system performance in function of the laser intensity and capacitor geometry is desirable for comparison with the theoretical model. It would be also desirable to measure the plasma characteristics inside the diode, to verify the plasma neutrality, to measure its density, the magnetic elds and the energies of ions accelerated in the cathode sheath.

=b 2 Figure 1 .

 21 Figure 1. Scheme of the capacitor coil target with notations used in the paper. The materials commonly used for target fabrication are aluminum, copper, nickel and gold.

Figure 2 .

 2 Figure 2. Measurements of magnetic elds at the coil center of nickel targets with dierent diagnostics as reported in Refs. [8] (a) and [9] (b).

Figure 3 .

 3 Figure 3. Potential (a) and density (b) distribution in the diode. Stationary distribution is shown with solid lines, dashed lines show the transient distribution before the potential equilibration.

Figure 4 .

 4 Figure 4.Optical shadowgraphy (at 532 nm) of the coil expansion. On the left, 0.65 ns-gated images corresponding to a cold coil (top, t = 0) and two dierent timings after laser-driving of the capacitor-coil target. On the right, the ns-scale expansion is resolved in a single shot by imaging the coil diameter into the slit of a streak camera. The average expansion velocity of the coil rod was measured to be 11 ± 3 µm/ns (as indicated by the dashed orange line). Data acquired on the experimental campaign[START_REF] Santos | Laser-driven platform for generation and characterization of strong quasi-static magnetic elds[END_REF].

Figure 5 .

 5 Figure 5. Dependence of the parameter g on the plasma ow divergence α0 and the ratio d/r h according to Eq.[START_REF] Forslund | Theory of hot-electron spectra at high laser intensity[END_REF]. The values of α0 are shown with numbers near the correspondent curves.

/Figure 6 .

 6 Figure 6. Current-voltage characteristic of the diode Vc(I) accounting for the eects of the space charge and magnetization. The chosen set of parameters corresponds to the experiments [5] a), [6] b) and [10] c). The distance between the diode plates is set to 1 mm and the resistance Rc = 1 Ω

Figure 7 .

 7 Figure 7. Temporal dependence of the current Ic (a) and tension (b) in the diode-coil circuit (a) and the circuit impedance Zc = Vc/Ic (c). The parameters correspond to the experiment [8].

ACKNOWLEDGMENTS

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reect those of the European Commission. We also acknowledge nancial support from the French National Research Agency (ANR) in the frame work of `'the investments for the future Programme IdEx Bordeaux-LAPHIA (ANR-10-IDEX-03-02).


J J Moreau

NUMERICAL INVESTIGATION OF SHEAR ZONES IN GRANULAR MATERIALS

The Contact Dynamics method is a numerical technique for computing the dynamical motion of collections of bodies submitted to the unilateral constraints of mutual non-interpenetrability, taking into account dry friction in the case of contact and also the velocity jumps which may arise from collisions. The essential nonsmoothness of the problem is faced without resorting to any regularizing approximation such as artificial repulsion or artificial viscosity. The method is illustrated here by the simulation of a rectilinear shear zone occurring in a slowly strained granulate. The concept of Internal Moment is introduced in order to analyse the transmission of forces.

Introduction 1.Computation in Multibody Dynamics

After the pioneering work of P.A. Cundall [START_REF] Cundall | A computer model for simulating progressive large scale movements of blocky rock systems[END_REF], a number of authors have developed codes for computing the motion or the equilibrium of collections of bodies through the approximate replacement of the non-interpenetrability constraints by repulsion laws which enter into action when a pair of bodies come sufficiently close to each other. This is the mechanical analogue of the penalty methods commonly used in Constrained Optimization. Some of these codes are directy inspired by the numerical techniques of Molecular Dynamics.

Significant simulations of geological structures, of masonery works or of granular materials have been obtained in this way. The main drawback is that the need of precision calls for very stiff repulsion laws, so that the time-discretization algorithms (usually of the explicit type), applied in integrating the equations of Dynamics, have to resort to fine time-steps and possibly also to artificial damping in order to secure numerical stability. In some instances the effects of artificial elastic repulsion and of artificial damping may become difficult to distinguish from those of the mechanically intentional elasticity or damping.

A different class of numerical methods are qualified as event-driven. They consist in applying to the unilateral constraints of non-interpenetrability (usually subject to friction) the theoretical procedures of Rational Mechanics. The second order differential equations of Dynamics are exploited with account of the contact constraints, if some of them are effective during the concerned time interval, and of the associated contact forces. Standard numerical techniques are applied to these equations, allowing one to determine instants at which the status of constraints changes. Such a change may arise from a collision suddenly introducing new constraints, in which case some empirical collision laws are invoked in order to determine the new velocities from which the subsequent motion have to be computed. The status of constraints may also change from the cessation of contacts, at instants which have to be detected by calculating the contact forces in every contact point. If at some instant a contact force is found with a direction incompatible with unilaterality, one concludes that the further motion has to be calculated by assuming a different set of effective contacts. But it has long been known [START_REF] Delassus | Mémoire sur la théorie des liaisons finies unilatérales[END_REF] that the contacts which cease at the detected instant are not necessarily those for which an unfeasible reaction force was calculated. The heart of the matter is a Complementarity Problem of the sort commonly met in Convex Optimization [START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF].

For the application of the above approach to the Dynamics of Machines, which involves a relatively small number of bodies, the reader may refer to [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF]. In what concerns Granular Materials, event-driven methods are primarily applied to the dynamics of dilute states, so that the interaction of grains solely rests on collisions. A trick to approximate durable contact by motions of this sort consists in introducing some auxiliary agitation. Anyway, event-driven computation becomes unfeasible when the average time between collision is too small.

The Contact Dynamics Method

This method, abbreviatively labelled as 'the CD method' originated in [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], a text in which the unilateral mechanical constraint of (possibly frictional) noninterpenetrability receives a formulation, resting on a theoretical background of elementary Convex Analysis, which proves suitable for numerical computation. This formulation also applies to deformable bodies, spatially discretized by Finite Elements [START_REF] Jean | Frictional contact in collections of rigid or deformable bodies: numerical simulation of geomaterials[END_REF], with applications in such domains as Geodynamics or the industrial processes of metal forming.

Only mechanical systems with finite number of degrees of freedom will be considered in this lecture. An exposition of the method in this case may be found in [START_REF] Moreau | Some numerical methods in multibody dynamics : application to granular materials[END_REF]. The possibility of collisions makes one expect velocities to be discontinuous functions of time. Actually, in the Dynamics of systems affected with dry friction, velocity jumps may even occur in the absence of any collision as the result of some sort of dynamic locking [START_REF] Lecornu | Sur la loi de Coulomb[END_REF].

The mathematical formulation of the corresponding evolution problems rests on the assumption (for which some mechanical justification exists [START_REF] Moreau | An expression of classical dynamics[END_REF]) that velocities are vector functions of time with locally bounded variation [START_REF] Moreau | Bounded variation in time[END_REF]. Then the role of the accelerations is played by vector measures on the considered time interval. In the traditional smooth case, these measures happen to possess density functions relatively to the Lebesgue measure dt, which are nothing else than the accelerations, as functions of time. Possible isolated velocity jumps manifest themselves in this formalism by 'atoms' of the above measures, i.e. some vector-valued Dirac measures localized at the corresponding instants. The system dynamics is expressed in general as a measure differential inclusion instead of the conventional differential equation of 'smooth' Dynamics.

The analysis of the solutions, i.e. the theoretical study of their existence and of their approximations, is a rather technical matter. The interested reader may refer to [START_REF] Monteiro Marques | Differential Inclusions in Nonsmooth Mechanical Problems: Shocks and Dry Friction[END_REF] for these mathematical aspects, but the principal merit of the formulation lies in its ability to generate time-discretization algorithms which, at every step, are ready to treat velocity jumps on the same footing as smooth evolutions. The discretized equation merely expresses the balance of momentum on the corresponding time-interval, should the involved mechanical actions be ordinary forces or percussions. The nature of the calculations to be effected at each time-step makes the resulting algorithms belong to the implicit type, so that they are never hampered by numerical instability.

The physical information concerning the contact interactions is fed into the calculation in a standardized form called a complete contact law, devised in order to take care automatically of the geometric conditions of non-interpenetrability. Of course, these contact laws have to encapsulate not only the description of unilaterality and friction in persistent contacts but also the model accepted for predicting the outcome of possible collisions.

At each time-step, an iterative procedure is applied in order to calculate simultaneously all contact impulsions, collisional or not, and the resulting changes in velocities. The computation cost of this procedure naturally increases with the number of contacts to be handled. Parallelisation is currently investigated but, up to now, only microcomputers or workstations have been employed. For a dense pack of 4000 round grains, the computation time is typically 15 s per step in two-dimensionsal situations, 30 s per step in three-dimensional ones, on a Power Macintosh 8500 microcomputer. The detection of contacts is implemented in such a way that it takes only a fraction of a second in the above case. In near-equilibrium situations, computation becomes much quicker because the iterative procedure uses as first guess the contact forces calculated at the preceding step (in particular, if an equilibrium state is attained at a specified precision level, each step only consists of a precision check). Of course, in the collision-dominated dynamics of rapid granular flows, the latter trick has no effect, but in this case speed is favoured by the reduced number of contacts. The choice of the step-length, in order to reach a given precision depends, for compact states on the magnitude of the acting forces and, in rapid flows, on the magnitude of the collision velocities. Of course, quantitative rules could be formulated in the framework of a dimensional inspection. For instance, reducing the acting forces is equivalent to increasing the mass of each grain.

Collisions

Physically, a collision is a 'short' episode during which the contact forces between the involved bodies assume 'large' values. The phenomenon cannot be properly analyzed without taking body deformations into account, not only in what concerns the vicinity of the impact locus but also in discussing to what extent some wave-like propagation may transfer momentum to other parts of the investigated system. However, from the beginnings of Mechanics, it has been attempted to incorporate collisions into the dynamical study of bodies though observed at such a scale that microscopic deformations are invisible. This is called the theory of Rigid Body Collisions, a recent survey of which may be found in [START_REF] Brogliato | Nonsmooth Impact Mechanics. Models, Dynamics and Control[END_REF]. The objective is to calculate the velocities of the system parts after the collision from the values these velocities had before it. The 'collision operators' commonly used to this end can only provide rather crude predictions, as long as a detailed analysis of the collision process is not made, usually carried out through multiple scaling techniques involving a 'microtime'. Furthermore, one should keep in mind that if a body engaged in a collision is part of cluster of previously contacting objects, all contact points in the cluster are liable to transmit percussions. In such cases of multiple collisions the traditional concept of restitution coefficient may produce grossly incorrect results.

In CD algorithms, the following improvement on the restitution concept provides a treatment of multiple collisions which at least has the merit of consistency. Contact laws, similar to those applied to persistent contacts, are invoked, in order to connect the (unknown) contact percussions with some formal velocities. The latter are constructed from the (known) pre-collision velocities and from the (unknown) post-collision ones, by a weighted averaging operation. If all the weight is placed on post-collision velocities, the consequent prediction turns out to be that of a collision of the fully inelastic sort. On the contrary, by placing equal weights on both velocities, one generates a model of fully elastic collision. At each contact point, different weights may be assigned to the tangential and to the normal components of the relative velocities.

The invoked contact laws also include friction. In the simple case of a binary collision, this approach results in a collision operator involving three parameters, identifiable as the friction coefficient, the normal restitution coefficient and the tangential restitution coefficient. For spherical bodies, the resulting formula comes out identical to what has been proposed, from an entirely different background in [START_REF] Walton | Numerical simulation of inelastic, frictional particle-particle interactions[END_REF][13], fairly agreeing with the experiments of [ [START_REF] Foerster | Measurements of the collision properties of small spheres[END_REF]].

The computation of the shaking of dense granulates related in [START_REF] Moreau | Numerical experiments in granular dynamics: vibrationinduced size segregation[END_REF] provides typical examples of multiple collisions. Physically, one expects that the transmission of collisional impulses throughout a pack of contacting grains involve the elastic deformation of these grains in a process quite similar to the propagation of sound. Since it has been decided here to model grains as perfectly undeformable, such propagative effects are ignored. The success of the method in the investigated instances should be ascribed to the fact that, physically the time taken by sound-like disturbances to travel the extent of the experimental cell is definitely shorter than the period of shake.

In contrast, multiple collisions play a minor role in the dynamics of dilute granulates. The greater efficiency of the CD method, if compared with eventdriven ones, mainly stems from that, once a time-discretization mesh has been chosen, all the collisions which are detected as occurring on a given time-step are treated together. This of course entails some trade-offs which have to be assessed. It seems of little consequence that the ordering of the collisions, which mechanically should be successive, is only internal to the algorithm. In fact, the problems in view are physically undeterministic, since a slight change in the initial conditions is enough to produce, after a short time, a completely different sequence of positions and collisions. A more critical observation is that each collision is treated as occuring only once between the involved bodies in the considered time-step. An accumulation of bounces, such as those of a ping-pong ball coming to rest, is thus viewed as a single collision as soon as the expected successive bounces are all comprised in the same time-step. This could result in underestimating the energy loss arising from inelastic collisions. This source of error may be checked by repeating the computation with reduced time step : if no change is found in the energy vs time curve, one may conclude that no harm was done.

Previous Applications of the Method

Calculations of shaken granular materials, in two or three dimensions, have been the first published illustrations of the CD method. They gave the occasion of very convincing comparisons [START_REF] Moreau | Some numerical methods in multibody dynamics : application to granular materials[END_REF] with the physical experiments related in [START_REF] Clément | Experimental study of heaping in a two-dimensional "sandpile[END_REF], evidencing the effects of boundary friction : boundary vortices and heaping.

Concerning the popular phenomenon of size-segregation the role, in certain settings, of the friction-induced boundary currents and the consequent convection effects has been revealed [START_REF] Moreau | New computation methods in granular dynamics[END_REF], independently of the experiments of [START_REF] Knight | Vibration-induced size separation in granular media : the convection connection[END_REF]. Furthermore, CD simulation has disclosed how such boundary currents are generated. [START_REF] Moreau | Numerical experiments in granular dynamics: vibrationinduced size segregation[END_REF] Above applications pertained to the rapid dynamics of dense granular materials. Actually, some simulations of slow deformation have confirmed that, even if the external driving which compels a granular sample to deform is so slow that, at the level of the overall motion, acceleration effects seem negligible, this is not usually the case microscopically. [START_REF] Jean | Numerical treatment of contact and friction: the Contact Dynamics method[END_REF] It is well known (see, in particular [START_REF] Meftah | Evidence of local 'seisms' of microscopic and macroscopic stress fluctuations during the deformation of packings of grains[END_REF]) that the deformation involves a succession of crises which essentially are dynamical processes. Consequently, in numerical simulations, the grain masses must be entered with their true physical values, without attempting to improve computation speed by artificially increasing them.

Provided dynamical effects are taken into account, the average stress due to contact in a portion of deforming granulate remains a meaningful element of the analysis. We therefore have to revisit this concept before finally presenting in Sec.3 below a detailed example of shear.

The Definition of Stress

A number of authors have proposed definitions of the stress in a portion of granular material in equilibrium or quasi-equilibrium, all essentially equivalent [see, among others, [START_REF] Weber | Recherches concernant les contraintes intergranulaires dans les milieux pulverulents[END_REF] [28] [START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF]]. If dynamical effects are also to be taken into account, one should first recall that the convection of momentum significantly contributes to the transmission of macroscopic efforts in the standard kinetic theory of fluids as well as in the modelling of dilute granular flows [see e.g. [START_REF] Savage | The stress tensor in a granular flow at high shear rate[END_REF]]. But, in the evolution of dense granulates which is the subject of this lecture, the mean free path appears negligible in regard to particle size. Consequently, a concept pertaining only to contact interactions, possibly collisional, is all what is needed. Anyway, the following [START_REF] Moreau | Internal moment and contact forces in fragmented materials, oral lecture at EUROMECH Colloquium 351: Systems with Coulomb Friction[END_REF] remains consistent in all situations. By a fragmented system, we mean a collection of (deformable or rigid) bodies which interact only through contact. Examples: a granular material, a piece of masonry. In order to define, for such a system, some concept able to play, as far as possible, the role that the Cauchy stress tensor holds in conventional continuous media, we are to introduce the 'internal moment tensor' of an arbitrary material system. This does not require of the number of constituents to be large and the concept makes sense even for a single granule. In particular, when constructing averages over a granular sample, one may also investigate the statistical dispersion of the individual terms involved.

The Virtual Power Formalism

According to this standard formalism of Classical Mechanics, [START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: microstructures[END_REF] any effort (this is not necessarily a force in the traditional sense) experienced by a material system is defined by the expression of the power that this effort would develop if the system elements were affected by a class of imagined velocity fields, called virtual velocities or test fields. The rule of the game is that this class of vector fields constitutes a linear space and that the functional 'power' is linear on this space.

Let S be a material system occupying a bounded portion of space and let ϕ denote an affine test field, with Cartesian components

ϕ i (x) = ϕ i (o) + b ij x j . (1)
Since the power P int (S, ϕ) of the internal efforts of S when ϕ is taken as test field depends linearly on ϕ, i.e. linearly on the parameters ϕ i (o) and b ij , there exist R i and M ij such that

P int (S, ϕ) = R i ϕ i (o) + M ij b ji .
The virtual power formalism generalizes Newton's third law in the form:

The internal efforts of any material system have zero power whenever ϕ equals the velocity field of a rigid motion, i.e. a field as in Eq.( 1) with b ji = -b ij (provided Cartesian coordinates are orthonormal).

This readily implies that R i vanishes and that M ij constitutes the components of a symmetric tensor M of rank 2, independent of the choice of the origin and of the (orthonormal) axes in use. Let us call this tensor the Internal Moment of S.

Example. Let S consist of a portion of a classical continuous medium occupying a bounded domain Ω of space, with σ as Cauchy stress field. Then classically, for every continuously differentiable ϕ, P int (S, ϕ) = Ω σ ji ϕ i,j dω.

(

) 2
This actually constitutes the very definition of the Cauchy stress in a synthetic construction of Continuum Mechanics [START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: microstructures[END_REF].

Here we comply with the common usage in Geomechanics and Civil Engineering of counting a stress as positive when it is directed as a pressure. In Mechanical Engineering and in general Continuum Mechanics, the reverse convention is applied, which would require a minus sign in the right-hand side of Eq. [START_REF] Delassus | Mémoire sur la théorie des liaisons finies unilatérales[END_REF].

By taking as ϕ the same affine vector field as before, one obtains

M ij (S) = Ω σ ij dω. (3)

The Static Case

Assume in this Subsec. that the considered system is in equilibrium. Then the Principle of Virtual Power states that the total power of all the efforts it experiences, namely the internal and the external ones, equals zero whatever is ϕ, so P int = -P ext . If, for instance, the external efforts experienced by S simply consist of a finite collection of forces f α exerted at points x α , one may take ϕ as in Eq.1 and find

M ij = - α x α i f α j (4)
i.e. in the event of equilibrium, the internal moment equals the opposite of the moment tensor of the external forces. The latter is by itself independent of the choice of the origin and symmetric because the equations of Statics make the resultant of external forces vanish, as well as their skew moment. All the above writing readily extends to continuously distributed external efforts, such as gravity forces, and also to external efforts of higher order. The latter may include space distributions of magnetic actions or the boundary couples commonly introduced in the phenomenological description of the resistance to rolling, but this does not prevent M from being a symmetric tensor of rank 2.

Additivity

Suppose that S equals a collection of subsystems S 1 , S 2 ,. . . which interact only through contact actions. The internal efforts of S consist of the respective internal efforts of the subsystems and of the mutual efforts that these subsystems may exert upon each other. Since the efforts exerted, say by S 1 upon S 2 and by S 2 upon S 1 are contact actions (possibly involving couples of resistance to rolling), they are located at the same points of space and, by Newton's third law, have opposite values. Therefore, the total contribution of these mutual efforts in the internal moment of S vanish, so one has M (S) = M (S 1 ) + M (S 2 ) + . . . More generally, by assigning to every finite collection of subsystems of S the corresponding internal moment, one defines a finitely additive tensor-valued measure.

The occurence of an integral in Eq.( 3) makes additivity evident in this special case and shows that, in a classical continuous medium, the Cauchy stress field is nothing else than the density function of the internal moment measure relatively to the volume (i.e. Lebesgue) measure. The same equation expresses that the volume average of the Cauchy stress field over Ω equals M /vol(Ω).

Therefore, if a portion S of granular material is viewed as 'enveloped' in a bounded domain Ω, it is natural to say that the tensor M (S)/vol(Ω) constitutes the average stress of S. The uncertainty in defining such an enveloping domain has little consequence if S is sufficiently dense and comprises sufficiently many grains. It is only at this stage that some assumption of 'large number' has to be made. It will be shown on the example of Sec.3 that the tensor calculated in that way does possess the main property expected of the Cauchy stress, namely it allows one to predict how the force transmitted across an ideal cut depends on the direction of this cut.

The above construction, based on the Virtual Power formalism, has the advantage of possible extension to many other models of continuous media. For instance the internal moment measure of a string or a chain admits a density function relatively to the arc-length measure. The visualisation of force chains in the granular sample of Sec.3 through the internal moments of individual grains suggests that the construction of models of granular media could rest on formulations in which the internal moment measure would not be represented by a density with respect to volume.

The Dynamical Terms

The existence of the tensor M (S) and its additivity property for a fragmented material system are inherent in the principles of Classical Mechanics. There only remains to know how to evaluate this tensor in specific situations.

Equation ( 4) pertains to the special case of equilibrium. If, more generally, the motion of S relative to some inertial reference frame, under some specified mechanical actions, has been calculated, one has to express that, for every test field ϕ, the virtual power of all the efforts experienced by the system equals that of the vector measure a(x) dµ(x), where a denotes the acceleration field and µ the mass mesure. By using as ϕ the same affine vector field as before, one obtains, for every fragment s of the system, in coordinate-free notation,

M (s) = -tens.mom(ext.efforts) + s x a(x) dµ(x). (5)
This applies, in particular, to each element of a granulate material. The total expression on the right-hand side, like M (s) itself, does not depend on the choice of the origin. When s is an indeformable grain, it proves convenient to take as origin the center of mass of this rigid body, then express a(x) by the standard formulas of rigid body kinematics and use the inertia tensor of s in order to calculate the integral in closed form. The result involves in general the acceleration of the center of mass of s, the spin vector ω and its time-derivative ω'.

Let us only report the simplification which occurs in two-dimensional situations, when the rigid grain possesses rotational symmetry. Then Eq.( 5) is found to reduce to M (s) = -sym.mom(ext.efforts) -1 2

ω 2 I 1. (6)
where I denotes the moment of inertia of s about its center and 1 the rank 2 unit tensor. The symmetric moment tensor of external efforts is evaluated about the center of s. Interestingly, ω' does not appear in this special case. What precedes readily adapts to instants at which efforts are of percussional nature. Recall in this connection that CD computation consists in solving, on each time-step, the balance of momentum of the system, in which forces manifest themselves through impulsions. Proper forces are retrieved after division by the step-length.

Slow Deformations of a Dense Granular Material

Due to space limitation, we restrict ourselves in this paper to the description of a single numerical experiment, concerning a compact two-dimensional assembly of strictly indeformable circular objects. Of course, the CD method could as well be applied to some models approaching natural grain shape more closely [START_REF] Laalai | Micromechanical modelling of ballasted tracks[END_REF].

The choice of the circular shape is not uniquely motivated by computational ease, but by the want of avoiding the mix of distinct effects. In [START_REF] Oda | Micro-fabric and couple-stress in shear bands of granular materials[END_REF] physical experiments on the slow shear of assemblies of elongated particles are reported.

There is shown that shear may induce orientation effects, making of a 'shear band' a structurally individualized part of the granulate, strongly affecting in particular the pattern of force chains. Similarly in [START_REF] Rothenburg | Micromechanical features of granular assemblies with planar elliptical particles[END_REF], the numerical simulation of assemblies of elliptic particles (performed by a method of elastic repulsion) reveals some structural effects of particle elongation during shear. Let a two-dimensional pack of circular grains be squeezed between two parallel platens with locations y = y inf , y = y sup . In laboratory, the two lateral free boundaries of such a pack would be sustained by the pressure of a liquid, separated from the granulate by a membrane. We rather use here an ideal containing action, easier to analyze. The grains are sorted according to y, so as to divide the height y sup -y inf into horizontal bands of uniform width w. In each band the five ending grains on the left are submitted, in proportion to their respective masses, to horizontal forces totalizing wP , where P denotes a fixed pressure. The five ending grains on the right are similarly submitted to horizontal forces totalizing -wP . The width w, i.e. the number of bands, is adjusted in order that five grains make sufficient barriers at both ends of each band.

Platen friction is chosen greater than interparticle friction. If platens were forced to move vertically toward each other, this would generate the same classical X-shaped shear as in the simulations related in [START_REF] Bardet | A numerical investigation of the structure of persistent shear bands in granular media[END_REF]. Here we instead produce a single straight shear zone which proves more informative. This is done by maintaining the lower platen fixed, while the upper one is given a 45 o downward motion. A short time after start, the velocity vectors of the grain centers assume the distribution shown on Fig. 1. Nearly undeforming packs look attached to the respective platens, while an oblique shear zone lies in-between.

On Fig. 2, a grid of 16 strips, parallel to the general direction of shear, has been superimposed to the zone. For each strip, the average of the longitudinal components of velocities is calculated. The result plotted in the histogram reveals a rather smooth transition from the zero velocity of the lower part to the constant value corresponding to the forced motion of the upper platen. This takes place on a length approximately equal to 25 times the average diameter of grains, but no clear limit can be ascribed to the shear zone.

Incidentally, one can also detect in all regions some grains with aberrant velocities. Through zooming, these grains would be seen flying inside persistent cavities. This is the very effect commonly observed when watching shear-box experiments performed on Schneebeli materials, i.e. piles of cylindrical rods stacked parallel. Such grains originate from a partial collapse of the cavity boundary and their velocity of expulsion may be much larger than the velocity of the overall motion. Of course, a free grain ends by hitting the cavity boundary, possibly triggering some rearrangement in the granulate.

Rapidly changing ondulations affect the velocity pattern in the shear zone. Sometimes, one may distinguish assemblies of a few grains whose velocity distribution is approximately that of a solid rotation. This reminds one of the persistent rigid clusters observed in real sand experiments. [START_REF] Drescher | Photoelastic verification of a mechanical model for the flow of a granular material[END_REF] That the effect here is much less pronounced and very fugacious probably comes from the grains being circular. Anyway such swirles have nothing to do with the vortices generated when two fluid flows merge. In the latter case, vortices last sufficiently long for their convection by the fluid to be observed. Furthermore, the velocity distribution in a fluid vortex is quite different, since velocities become larger toward the center.

In Fig. 1, three rectangular probing areas have been placed on the field. At the center of each one, a cross represents the principal directions and principal values of the average stress tensor, as defined in Sec.2. For each area, the histogram of contact normals is drawn, constructed by dividing the plane into 24 sectors. Each sector is filled in gray up to a radius proportional to the number of the contact points at which the normal n to the contacting grains (in either direction, so that the histogram is symmetric) belongs to the sector. The elongated shapes of these histograms reveal a strong anisotropy in the distribution of contact orientations. The dominant direction may be estimated by constructing the fabric tensor relative to the probing area, i.e. the symmetric tensor whose components equal the averages of n i n j and determining the major principal direction of this tensor. The table at the top shows that this major principal direction is, in each of the three cases, close to that of the average stress tensor.

The dynamical terms, as expressed in Subsec.2.4, have been taken into account when computing the components of average stress tensors. Actually, they have been found to contribute by less that 0.1 % in the result. The point is that the reported numerical experiment exhibits little dynamical effects. This should be ascribed in part to the relative magnitudes of the upper plate velocity, of the lateral pressure and of the grain masses. The present setting also differs from that of a shear-box by the presence of the lateral free boundaries. Exerting a distributed contention pressure upon these boundaries generates less frustration than acting on the granulate by means of a rigid wall. The global effect of dynamical crises may be assessed by calculating the normal components of the resultant forces exerted by the grains upon the two platens. Since the lateral contention forces are horizontal a difference between these two components purely reflects dynamical effects. In the present experiment the difference remains inferior to 1%, while in the analogous setting of a granulate confined between parallel platens under periodic lateral conditions (which are more frustrating than a pressure) the difference may sometimes attain 50%.

The main observation to be made here is that the central oblique probing area, which covers the shear zone and the two other ones covering the vicinity of the platens do not reveal significant differences in average stress and in fabric anisotropy.

On Fig. 2, the force chains are visualized by drawing the crosses representing the internal moments of all grains individually. Another visualisation of these chains is provided on Fig. 3 by representing each pair of intergranular forces by a line segment centered at the contact point. Both representations make clear that the pattern of force chains is completely indifferent to the presence of the shear zone.

A circular probing area has been placed in the upper part of the sample on Fig. 3. The program has calculated, for a series of diameters, the resultant force exerted by the grains with centers on one side upon those on the other side. The corresponding plot allows one to compare the values with the result of a calculation based on the average stress tensor of the area. In spite of the relatively small number of grains involved for each diameter and of the inhomogeneity due to force chains, the two plots are remarkably concordant.