Aleksei S Dmitrichev

Vladimir I Nekorkin

Rachid Behdad

Stéphane Binczak

Jean-Marie Bilbault

Alexey S Dmitrichev

Stephane Binczak

Anti-phase wave patterns in a ring of electrically coupled oscillatory neurons

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

Introduction

Having a clear view of the neuronal system dynamics is significant for an understanding of their functions and it is also irreplaceable when one endeavors to devise mathematical models of the neurons or neuronal networks. The most popular and practically important branch in such kind of researches is an investigation of small ensembles of neurons, where both neurons and their mutual coupling are described by sufficiently realistic mathematical models, based on the experimental data [START_REF] Ermentrout | Mathematical Foundations of Neuroscience[END_REF]. But in spite of simplicity of such models, they possess a great number of collective behavior regimes on account of ultrahigh multistability. Certainly, some of these regimes have no definite physical interpretation, but nevertheless, it can be of service to understand the basic information transmission processes in the brain.

The functioning of such complex neuronal systems as retina and other sensory networks is one of the crucial problems of neuroscience nowadays [START_REF] Kaplan | [END_REF][3][START_REF] Field | Information Processing in the Primate Retina: Circuitry and Coding[END_REF]. Retinal structure is adequately high-ordered, which would have simplified to devise mathematical models and to simulate the networks numerically. But at the same time this results in general view of connectivity patterns and brain principles of operation at large. At present the new mathematical models of retinal-like structures are suggested and some progress in understanding basic properties of such systems is achieved. But we take the view, that the very neuron models mentioned above can be useful in investigating the ordered neuronal structures.

In Fig. 1 the most typical impulse responses of retinal ganglion cells are shown. It is important that the surround cell response is out of phase with the central one. This makes possible to apply an anti-phase synchronization model of two electric-coupled Fig. 1. The first order responses of a P-type primate ganglion cell center and surround, adopted from [START_REF] Kaplan | [END_REF].

neurons to account for the fact of such complicated oscillatory structures existence that will be discussed later.

The paper is organized as follows. In Section 2 we describe the model studied and briefly discuss its basic dynamical properties. A special class of localized collective activity obtained in the model which was called "anti-phase wave patterns" is presented in Section 3. We show that the formation of anti-phase wave patterns is related to establishing of a certain phase relations between oscillations of neighboring elements. In Section 4 we use weakly coupled oscillators theory to study the regions of the existence of such patterns. Finally, the main conclusions are summarized in Section 5.

The dynamical model

We consider the following network composed of identical neural-like units or cells arranged on a one-dimensional ring and coupled locally in a nearest-neighbor fashion:

     dv i dt = -g L (v i -v L ) -g Ca M ∞ (v i )(v i -v Ca ) -. . . . . . -g K n i (v i -v K ) + I + d(v j-1 -2v j + v j+1 ) dni dt = ε n∞(vi)-ni τn(vi) i = 1, 2, . . . N, v j+N (t) ≡ v j (t).
(

) 1
The variable v i describes the evolution of the membrane potential of a neuron located at the jst site of the ring, n i describes the dynamics of the activation level of the potassium ion channels. The terms on the right-hand side of the first equation of (1) accounts for the currents flowing through the cell membrane of the neuron. The first three terms correspond to the leakage, the calcium and the potassium ionic currents, respectively. The last ones signify the external (applied) and the electrical synaptic currents. The parameters g L , g Ca , and g K are the maximum conductances, v L , v Ca , v K are the equilibrium potentials (reverse potentials) for the corresponding ion channels. The functions M ∞ (v) and n ∞ (v) obeying the following laws:

M ∞ (v) = 0.5 [ 1 + tanh ( v-v1 v2 )] , n ∞ (v) = 0.5 [ 1 + tanh ( v-v3 v4 )] , (2)
define, respectively, the stationary activation levels for calcium and potassium ion channels. Function τ n (v) which obeys the law:

τ n (v) = [ cosh ( v-v3 2v4 )] -1 (3)
and parameter ε (ε > 0) both control the characteristic relaxation time for potassium ion channels. Finally the parameter d characterizes the strength of electrical inter-neuron coupling. The action of electrical coupling is based on the existence of so-called gap junction between neighboring neurons [START_REF] Kandel | Principles of neural science[END_REF]. Such junction leads to direct connection of the cytoplasm of the neurons and allows various molecules and ions to pass freely between them. Hereupon the current passing through the gap junction obeys Ohm's law and depends directly on gap conduction (or inversely on gap resistance) and directly on difference in electric potential between neurons. Thus the electrical coupling is well-approximated by linear resistive coupling. Therefore in (1) the coupling terms are expressed in the form of standard one-dimensional discrete Laplacian.

Local cell dynamics

Let us first consider the dynamics of a single unit. It is described by the dimensionless version of the well known Morris-Lecar system [START_REF] Morris | [END_REF]:

{ dv dt = -g L (v -v L ) -g Ca M ∞ (v)(v -v Ca ) -g K n(v -v K ) + I dn dt = ε n∞(v)-n τn(v) (4)
To derive it from the original system we rescaled the time properly and used the following transformations of variables and parameters:

v = v v0 , v m = vm v0 , g m = g m g0 , v i = vi v0 , v 0 = I Ig0 , g 0 = φC ε m = L, Ca, K, i = 1, 4 (5)
where C is the membrane capacitance and φ is the maximum relaxation time for potassium ion channels. The original Morris-Lecar system is known to provide a wealth of dynamical behaviours. Detailed bifurcation analysis of its dynamics can be found, for example, in [7]. Here we briefly discuss only the dynamical regimes essential for the formation of anti-phase wave patterns. For our purposes we fix the parameters: and take ε and I as a control parameters. Figure 2 illustrates the partition of a part of (ε, I)-plane interested (highlighted by shaded area) onto regions with different dynamical behaviour of the system (4). According to the partition the system yields three different dynamical regimes which are depicted in its phase planes in Fig. 3.

g L = 2, v L = -60, g Ca = 4, v Ca = 120, g K = 8, v K = -80, v 1 = -1.2, v 2 = 18, v 3 = 12, v 4 = 17.4, (6)
The main features of the regimes are existence of stable limit cycle C s (see Fig. 3) which corresponds to periodic oscillations of the membrane potential of a neuron. Depending on the value of parameter I there are also either one, or three steady states. Respectively, for I < I * , that is in the region 1 (see Fig. see Fig. 3(b)). The curves H, Q and segment of the line I = I * below the point P correspond to bifurcations that leads to disappearance of the stable limit cycle C s . In particular, for points on Q the unstable and stable limit cycles C u and C s merge together and tangent bifurcation takes place. For points on H the limit cycle is glued to the separatrix loop (homoclinic orbit) formed by separatrices w u + and w s + of the saddle O 2 . Finally for points on the line I = I * below the point P the limit cycle is glued to the homoclinic orbit of saddle-node steady state formed by merging of the states O 1 and O 2 .

Anti-phase traveling patterns

Typical wave patterns in neuronal systems include travelling excitation pulses (action potentials), pulse trains (packets of action potentials) and wave fronts (moving interfaces between different levels of neuronal activities). The important property of such waves is that their principal characteristics (profile, velocity) are defined only by intrinsic properties of the medium and do not depend on initial conditions or perturbations within certain limits. But this is the case when the coupling between the elements is strong enough. When the coupling is weak the discrete nature of the system has an effect on its collective behavior. In particular, propagation of pulses, pulse trains and wave fronts becomes impossible (so-called propagation failure phenomenon) due to the existence of a large number of spatially inhomogeneous stationary states in the system [START_REF] Nekorkin | Sinergetic phenomena in active lattices[END_REF][START_REF] Keener | [END_REF][10][11].

Numerical study of the system [START_REF] Ermentrout | Mathematical Foundations of Neuroscience[END_REF] showed that it yields an interesting kind of localized activity which we called "anti-phase wave patterns". As opposed to aforementioned wave patterns, the latter ones exist even in the case of weak coupling. Figure 4(a) illustrates the evolution of typical anti-phase wave pattern and Fig. 4(b) depicts its snapshot. One can see that the pattern shares some features with envelope solitons found in the nonlinear Schrödinger equation [12,[START_REF] Malomed | Encyclopedia of Nonlinear Science[END_REF] and envelope quasisolitons found recently in reaction-cross-diffusion systems [START_REF] Biktashev | [END_REF]. It has the form of spatiotemporal oscillations with a smooth localized envelope which propagates along the system preserving its shape and velocity. Note, however, that spatial oscillations due to discreteness are not monotonic but represent instantaneous jumps back and forth between the top and the bottom part of the envelope. It is found that system supports formation of anti-phase wave patterns whose envelope has one (Fig. 4(b)) or more humps (Fig. 5). The patterns with even (odd) number of humps exist in systems with even(odd) number of elements only. Moreover the number of the patterns that exist in a system is proportional to its size. Thus the systems with a large number of elements will show high multistability of the anti-phase wave patterns.

Let us consider the dynamical mechanisms underlying the existence of the antiphase wave patterns. Figure 6 illustrates the projections of phase trajectory of the system (1) corresponding to one-hump anti-phase traveling pattern (Fig. 4) on the planes of membrane potentials of the first two neighboring (v 2 , v 1 ) and the first two next-nearest neurons (v 3 , v 1 ) and evolution of membrane potentials of these neurons. One can see that oscillations of the elements have the same periods. Moreover, oscillations of the neighboring elements are almost anti-phase synchronized, that is the phase shift between them is close (but not equal) to π, while oscillations of the nextnearest elements are almost complete synchronized (i.e. the phase shift is close to 0). The existence of such phase relations is caused by the existence of branches of slow and fast movements [15] on the local limit cycles C s of the elements. By-turn, it can be due to the proximity of the parameters of the elements to the saddle or saddle-node separatrix loop bifurcations, or it can be caused by their relaxational behaviour for ε ≪ 1.

The shape of the pattern considered suggests that phase relations analogous to aforementioned ones take place between oscillations of any three neighboring elements. Let us introduce the phase of oscillations of the elements and analyze the phase relations more carefully. Generally speaking, the phase can be rigorously defined for

t t n-1 t n t n+1 tv i v 1 v 3 v 2 (b)
Fig. 6. The projections of phase trajectory of the system (1) corresponding to the pattern shown in Fig. 4 on the planes of membrane potentials of the first three neighboring neurons (a) and evolution of their membrane potentials (b).

harmonic oscillations only. Still, for nonharmonic periodic or even chaotic oscillations, a quantity that possesses some properties of the classical phase can be introduced [START_REF] Anishchenko | Nonlinear Dynamics of Chaotic and Stochastic Systems[END_REF], and this quantity is also called "the phase". In this paper we define the phase of oscillations as follows. For definiteness, let us consider the i-th element. We introduce the sequence {t i n } of instants of the time at which the trajectory of the system crosses the line v i = 0 in the direction corresponding to the increase of variable v i (see Fig. 6(b)):

t i n = {t : v i (t) = 0, v ′ i (t) > 0}.
Then, the phase of oscillations of the i-th element at instant of the time t ′ is determined by the expression

φ i = 2π t ′ -t i n T i n for t ′ ∈ [t i n , t i n+1 ] (7)
where T i n = t i nt i n-1 is the "period" of oscillations. Note that such definition of the phase is meaningful only if T i n is constant or sufficiently close to each other. In the former case oscillations of the i-th element are regular and the phase increases steadily from 0 to 2π all the time. In the latter case oscillations of the element are phasecoherent chaotic and phase piecewise linearly depends on time. Notice that the phase introduced describes the dynamics of an individual element of the system (1) separatelly from the other ones. To describe the dynamics of the system as the whole it is more convenient to consider the quantity: ϕ j i = φ iφ j , that is the phase of oscillation of i-th element calculated relatively to oscillations of the j-th element.

We studied the behaviours of T i n and, for definiteness, ϕ 1 i for the patterns presented. It is found that in all cases neither T i n nor ϕ 1 i vary with time and T i n is constant over the entire network. In Fig. 7 the stationary phase distributions corresponding to the patterns are shown. Thus an anti-phase wave pattern represents a certain phase locked solution of the system [START_REF] Ermentrout | Mathematical Foundations of Neuroscience[END_REF]. In this solution each pair of neighboring elements is almost anti-phase synchronized and there is a constant phase shift between pairs that provides the wave nature of the pattern.

Summing up the results of numerical experiments one can obtain the following dependence of the phase on the unit number:

ϕ 1 i = (π + ∆ϕ)(i -1), (8)
where ϕ i is the phase of oscillation of the i-th element, i = 1, N , and ∆ϕ is a constant phase shift between every pair of elements. Note that since only anti-phase wave patterns with integral number of humps can fit into the full length of the system, the constant phase shift ∆ϕ can possess only discrete values:

∆ϕ = πk N (9)
where k = 0, ±1, ±2, .. (k ≪ N ) is an anti-phase wave mode number, its absolute value coincides with the number of its humps. Here k = 0 corresponds to a standing anti-phase wave. The positive values of k correspond to anti-phase wave patterns propagating to the left, while the negative ones correspond to the patterns propagating to the right. Moreover from ( 8) and ( 9) it follows that ϕ 1 N +1 = π(N + k). It implies that even modes including mode k = 0 can exist in the system with even number of elements, while odd modes can exist in the system with odd number of elements. Notice that it fully coincides with numerical findings presented.

The regions of existence of anti-phase wave patterns

Let us now study the regions of existence of anti-phase wave patterns. Assume that coupling between elements of the system (1) is weak enough (d → 0). In this case the theory of weakly coupled oscillators [START_REF] Malkin | Metodi Puankare i Liapunova v teorii nelineinix kolebanii[END_REF][START_REF] Malkin | Nekotorye zadachi teorii nelineinix kolebanii[END_REF][START_REF] Neu | [END_REF][20][START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF][START_REF] Hoppensteadt | Weakly Connected Neural Networks[END_REF][START_REF] Lewis | [END_REF] can be applied to reduce oscillatory dynamics of the system to the dynamics of a set of first order differential equations describing the evolution of phases of oscillations of its elements. We rewrite the system (1) in the general form:

dXi dt = F(X i ) + dS(X i , X i-1 ) + dS(X i , X i+1 ); ( 10
)
where

X i = [v i , n i ] T , F(X i
) is the local cell vector field defined by right-hand sides of (4) and S(X j ,

X k ) = [v k -v j , 0]
T . Recall that, a local element of the system has the stable limit cycle C s (see Fig. 3). Let T be a period of the cycle and G(t) = [v 0 (t), n 0 (t)] T a trajectory corresponding to the cycle. Following, Kuramoto [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF], we make the transformation of coordinates, ψ i = Θ(X i ). Without going into details, notice, that Θ(X) is an asymptotic phase function, which maps the points from some neighborhood U of the cycle (within its basin of attraction) to a unique phase. And it should satisfy the equation dΘ(X) dt = 1 for all X ∈ U when d = 0. Therefore, if we differentiate ψ i with respect to time using the chain rule for vector functions, we obtain the equation:

dψi dt = ∇ X Θ(X i ) dXi dt = = ∇ X Θ(X i )F(X i ) + d∇ X Θ(X i ) [S(X i , X i-1 ) + S(X i , X i+1 )] = = 1 + d∇ X Θ(X i ) [S(X i , X i-1 ) + S(X i , X i+1 )] (11)
One can see that the equation ( 11) is exact but it is not closed with respect to the phase variables ψ i . Hoverer, since d is weak, then X i (t) is close to G(t), which is close to G(ψ i ). Thus (11) becomes:

dψ i dt = 1 + d∇ X Θ (G(ψ i )) [S(G(ψ i ), G(ψ i-1 )) + S(G(ψ i ), G(ψ i+1 ))] .
By introducing the variables, φ i = -t + ψ i , we have:

dφi dt = d∇ X Θ (G(t + φ i )) [S(G(t + φ i ) , G(t + φ i-1 ))+ +S(G(t + φ i ), G(t + φ i+1 ))] (12)
Note that φ i evolves on a much slower timescale than T . Since the functions in the equation ( 12) are smooth and G(t) is T -periodic we can eliminate the explicit time dependence from its right-hand side by "averaging" over the period T [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF][START_REF] Ermentrout | [END_REF]. Then, by exploiting the fact, that the second component of coupling function S is zero, we arrive at:

dφi dt = d[H(φ i-1 -φ i ) + H(φ i+1 -φ i )], (13)
where

H(φ) = 1 T T ∫ 0 ∂Θ(G(t)) ∂v [ v 0 (t + φ) -v 0 (t) ] dt, H(φ) = H(φ + T ). ( 14
)
Finally to complete the reduction we must compute the partial derivative of the asymptotic phase function Θ(X) with respect to v. Unfortunately, in practice the straightforward calculation of Θ(X) for arbitrary X is not a trivial task. However, there are several ways of direct computation of the partial derivatives of Θ(X) evaluated at the limit cycle solution G(t) which have shown to be very efficient. Here we use the adjoint method based on the results of application of the singular perturbation approach to a weakly coupled oscillators [START_REF] Malkin | Metodi Puankare i Liapunova v teorii nelineinix kolebanii[END_REF][START_REF] Malkin | Nekotorye zadachi teorii nelineinix kolebanii[END_REF][START_REF] Neu | [END_REF][20]. They state that the partial derivatives of Θ(G(t)) are components of T -period solution of the adjoint equation for the local cell linearized around the limit cycle:

dZ dt = -DF(G(t)) T Z (15)
subject to the normalization constraint

Z(0)G ′ (0) = 1.
In practice, the solution to ( 15) is found by integrating the equation backward in time [26] since the system has the opposite stability of the original system.

In previous section we showed that an anti-phase wave pattern represents phase locked solution of the system (1) with the form:

φ i (t) = Ωt + ϕ i , ϕ i -ϕ i-1 = π + ∆ϕ. (16)
By substituting ( 16) into ( 13) one can obtain the following relation:

Ω = H(-0.5 -σ) + H(0.5 + σ). (17)
where σ = k 2N is a wave number and H(φ) = H(T φ). In fact the equation ( 17) defines the dependence of an anti-phase wave pattern local oscillations frequency Ω on its wave number σ. The example of relation Ω(σ) and the functions it depends on calculated for ε = 3.28, I = 48.3 are shown in Fig. 8.

Let us now analyze the linear stability of phase locked solution [START_REF] Anishchenko | Nonlinear Dynamics of Chaotic and Stochastic Systems[END_REF]. We suppose that ϕ i are not constant but depend on the time and introduce the variables, ρ i (t) = φ i (t)φ i-1 (t). It can be readily shown that the dynamics of ρ i satisfies:

dρi dt = d[-H(-ρ i-1 ) + H(-ρ i ) -H(ρ i ) + H(ρ i+1 )]. ( 18
)
Note that anti-phase wave patterns correspond to spatially-homogeneous stationary states of equation [START_REF] Malkin | Nekotorye zadachi teorii nelineinix kolebanii[END_REF]. Thus by analyzing stability of the stationary states of [START_REF] Malkin | Nekotorye zadachi teorii nelineinix kolebanii[END_REF] one can know stability of anti-phase wave patterns. Linearizing [START_REF] Malkin | Nekotorye zadachi teorii nelineinix kolebanii[END_REF] about the state ρ i = ρ 0 = 0.5 + σ we have

dξi dt = d[aξ i-1 -(a + b)ξ i-1 + bξ i+1 ], (19)
where ξ i = ρ iρ 0 are infinitesimal deviations from the stationary state and a = H ′ (-0.5σ), b = H ′ (0.5 + σ). Then, taking into account the periodic boundary conditions, we can represent ξ i as discrete Fourier series [START_REF] Nekorkin | Sinergetic phenomena in active lattices[END_REF] and derive the spectra of Lyapunov characteristic eigenvalues associated with the stationary state ρ 0 :

λ i = -(b + a) + (b + a) cos [ (i -1) 2π N ] + I(b -a) sin [ (i -1) 2π N ] , ( 20
)
where I is the imaginary unit. By analyzing the spectra (20) one can obtain the following stability condition for spatially-homogeneous stationary states of ( 18) and so for anti-phase wave patterns:

H ′ (-0.5 -σ) + H ′ (0.5 + σ) > 0. ( 21
)
The fig. 9(a) illustrates in (ε, I)-plane the region of existence of stable anti-phase wave patterns with wave numbers belonging in the interval σ ∈ [0, 0.05]. We see that for the values of parameter I < I * (see section 2.1) the region coincides with the domain of existence of oscillatory behaviour in a local cell. On the other hand, for I > I * the region is localized in neighborhood of upper boundary of oscillatory cell behaviour, that is curve Q. Though the boundary of the region of existence is sharp enough its structure is by no means trivial. Indeed, in Fig. 9(b),(c) two typical pictures of the structure of the boundary far to the left and to the right from the point P , respectively, are shown. If we cross the boundary to the left of P (Fig. 9(b)) then the smaller an anti-phase a wave pattern has wave number the earlier the pattern becomes unstable. On the other hand if we cross the boundary to the right of P (Fig. 9(c)) then the larger the pattern has a wave number the earlier it becomes unstable. Note that the point P is a point of intersection of the boundaries of existence regions of anti-phase wave patterns with σ = 0 and σ = 0.05. In its neighborhood the structure of the boundary is even more complex. Here the boundary of existence of each anti-phase wave pattern intersects the boundaries of other anti-phase wave patterns with smaller wave numbers forming fractal-like parameters set.

Conclusions

Spatio-temporal dynamics of a ring system modeling collective behavior of electrically coupled oscillatory neurons have been investigated. It has been found that for small values of coupling strength, when the discreteness affects its collective behavior, the system is capable to produce a special class of localized patterns so called "anti-phase wave patterns". The patterns look like the envelope solitons. Indeed they have the form of spatiotemporal oscillations with a smooth localized envelope which propagates along the system preserving its shape and velocity. However, due to discreteness, spatiotemporal oscillations represent instantaneous jumps back and forth between the top and the bottom part of the envelope. It has been shown that the patterns with different number of humps can be formed under the same values of parameters, indicating high multistability. Moreover the patterns with even number of humps can be formed only in the system with even number of elements, while the patterns with odd number of humps exist only in the system with odd number of elements.

The mechanisms of the formation of anti-phase wave patterns has been discussed. In particular, it has been shown that formation of the patterns is related to establishing of a certain phase relations between oscillations of neighboring elements. Indeed each pair of neighboring elements in the system is almost anti-phase synchronized and there is a constant phase shift between pairs that provides the wave nature of the pattern.

By using the weakly coupled oscillators theory the regions of the existence of anti-phase wave patterns are studied. It has been shown that in the domain where a cell has three stationary states the region of existence of patterns coincides with the domain of oscillatory cell behaviour, if not it is localized in a neighborhood of upper boundary of oscillatory cell behaviour. More over it is found that its boundary has a complex, even fractal-like, structure. Note that the anti-phase wave patterns can exist also in systems with non-ring structure. For example, Figure 10 illustrates the typical anti-phase wave pattern found in a chain of locally coupled neural-like units (4) with zero flux boundaries.

This work has been supported by the Russian Foundation for Basic Research (Grants 12-02-00526, 12-02-31825, 13-02-00858) and the Federal Target Program "Scientific and scientific-pedagogical personnel of the innovative Russia" in 2009-2013 (Contracts 8497, 8205).

10 Fig. 2 .Fig. 3 .

 1023 Fig. 2. Partition of a part of (ε, I)-plane on regions corresponding to different dynamical regimes of the system (4) for parameters given by (6), where I * = 39.69 and P (3.8, 39.69) is intersection point of the curve H, and line I = I * .

1 and O 2 Fig. 4 .

 124 Fig. 4. Space-time plot (a) and snapshot (b) of typical one-hump anti-phase wave pattern obtained in the ring of 101 units. Parameter values: ε = 3.28, I = 48.3, d = 0.01.

Fig. 5 .

 5 Fig. 5. Examples of multi-hump anti-phase wave patterns obtained in the rings of (a) 100 and (b) 101 units. Parameter values: ε = 3.28, I = 48.3, d = 0.01.

Fig. 7 .

 7 Fig. 7. Phase distributions corresponding to the anti-phase wave patterns shown, respectively, in Figures 4 and 5.

Fig. 8 .

 8 Fig.8. The images from top to bottom are waveforms of membrane potential v0 corresponding to the local limit cycle C s , v-component Z v of a solution to the adjoint problem(15), interaction function H defined by (14), and the relation Ω(σ) for solutions of the form[START_REF] Anishchenko | Nonlinear Dynamics of Chaotic and Stochastic Systems[END_REF], respectively. Parameter values: ε = 3.28, I = 48.3.

Fig. 9 .

 9 Fig. 9. The region (a) of existence of anti-phase wave patterns with wave numbers σ ∈ [0, 0.05] on (ε, I)-plane and enlarged fragments (a), (b) of its bottom boundary near the points (2.6661, 60.09), (0.9774, 110.06), respectively, in coordinate system (ε ′ , I ′ ) rotated tangent to the boundary. The point R(1.4233, 98.88) is a point of intersection of the bottom boundaries of existence regions of anti-phase wave patterns with σ = 0 and σ = 0.05.

Fig. 10 .

 10 Fig. 10. Space-time plot (a), snapshot (b) and corresponding phase distribution (c) of typical anti-phase wave pattern obtained in the network of 355 units with non-periodic (zero flux) boundaries. Parameter values: ε = 4, I = 38.2, d = 0.01.


Mohamed Jebahi

Damien André

Frédéric Dau

Jean-Luc Charles

Ivan Iordanoff

Simulation of Vickers indentation of silica glass

Keywords: Brittle, Amorphous, Glassy silica, Densification, Fracture, Indentation, Discrete elements, discrete-continuum coupling

The indentation response of glasses can be classified under three headings: normal, anomalous and intermediate, depending on the deformation mechanism and the cracking response. Silica glass, as a typical anomalous glass, deforms primarily by densification and has a strong tendency to form cone cracks that can accompany median, radial and lateral cracks when indented with a Vickers tip. This is due to its propensity to deform elastically by resisting plastic flow. Several investigations of this anomalous behavior can be found in the literature. The present paper serves to corroborate these results numerically using the discrete element method. A new pressuredensification model is developed in this work that allows for a quantitative estimate of the densification under very high pressure. This model is applied to simulate the Vickers indentation response of silica glass under various indentation forces using the discrete element method first, and then a discrete-continuum coupling method with large simulation domains to suppress the side effects and reduce the computational time. This coupling involves the discrete element method (DEM) and the constrained natural element method (CNEM). The numerical results obtained in this work compare favorably with past experimental results.

Introduction

Because of their excellent physical (i.e. mechanical, thermal and optical) properties in addition to excellent chemical resistance, several types of glasses have become the main focus of several researchers. According to the dependence of various properties such as density, hardness, refractive index, Young's and shear moduli with the pressure and fictive temperature (defined in a simple manner as the temperature at which the glass would be in thermal equilibrium), glasses can be classified under three headings: normal, intermediate and anomalous. Normal glasses, such as soda-lime silicate glass, exhibit properties decreasing with increasing fictive temperature and pressure. Anomalous glasses, such as germania, exhibit the opposite trend in these properties. Intermediate between normal and anomalous glasses are expected to have some properties that are independent of these parameters (i.e. pressure and fictive temperature). Of particular interest is silica glass which is known to exhibit anomalous behavior in its thermal and mechanical properties [1,2]. Furthermore, certain properties of this type of glass such as the Young's and shear moduli and the density show an anomalous dependence on the fictive temperature. When indented with a Vickers tip, silica glass exhibits a cracking pattern that is more complex than in a typical normal glass, such as soda-lime silicate [3,4]. Specifically, silica glass has a strong tendency to form ring cracks at the periphery of the indenter impression, but other crack systems typical of normal glasses, i.e. median, radial and lateral cracks, can be introduced [5,6]. The anomalous cracking response in silica glass is due to its propensity to deform elastically by resisting inelastic flow. In the region of ordinary hydrostatic pressures (up to approximately 8 GPa), silica glass behaves in a perfectly elastic manner. However, and contrary to crystalline solids, when the hydrostatic pressure is in the 8 GPa range or beyond, silica glass begins to exhibit signs of permanent deformation. The permanent deformation of silica glass, which is an amorphous solid, is very different from the plastic flow found in crystalline solids. Indeed, the plastic flow in crystalline solids is volume-conservative and initiates under shear stress only. The spherical part of the stress tensor has no influence on the plastic flow. However, the permanent deformation of silica glass is always accompanied by a volume change. In addition, as found in several previous works [7,[START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF][START_REF] Ji | [END_REF]10], it can initiate under hydrostatic pressure. There are several works on this subject in the literature, some of which are controversial. It began with the works of Bridgman and Šimon [7]. They demonstrated that amorphous silica exhibits permanent densification which translates to a permanent volume decrease of approximately 7.5 % at approximately 20 GPa. They also demonstrated that this propensity decreases when other substances are added, such as sodium oxide, Na 2 O. However, at that time, technical difficulties prevented the generation of strictly hydrostatic conditions at pressures above 5 GPa [7]. Uniaxial pressures were applied on thin disks to study the densification of glasses. Unfortunately, these tests induced additional shear stresses, which explain the systematic fracturing of the samples. Cohen and Roy [11] studied the densification of silica powder under uniaxial pressures. They found that the densification begins at approximately 2 GPa (10 GPa in Bridgman and Šimon's work [7]) and the volume reduction can far exceed 7.5 %. According to Reference [12], the difference between the previous results is due to the shearing effects, which are most important in the case of silica powder and allows a greater densification level to be reached. This explanation was belied later by Cohen and Roy [13], who excluded the shearing effects on the densification level. For more details, the reader can refer to [START_REF] Keryvin | Contribution à l'étude des mécanismes de déformation et de fissuration des verres[END_REF]. Several other works treating the shearing effects can be found in the literature [START_REF] Rouxel | [END_REF]16,17]. At present, this subject has not been resolved, and the role of the shearing stresses in densification remains ambiguous. In the present work, the densification of silica glass is studied numerically using a discrete element method developed by André et al. [18,[START_REF] André | Modélisation par éléments discrets des phases d'ébauchage et de doucissage de la silice[END_REF] for brittle continuous media. A new densification model based on the compressive stress in the cohesive bonds between the neighboring discrete elements is proposed. As will be shown later, this model provides a relatively good quantitative estimate of silica glass densification under high hydrostatic pressure (in comparison to previous numerical studies [START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF]). Subsequently, the influence of the densification on the cracking response of silica glass indented with a Vickers tip is investigated using a new fracture model recently developed by André et al. [20]. This model, based on the notion of virial stress [21,[START_REF] Zhou | Proc. R. Soc. Lond. A[END_REF], accurately simulates the cracking mechanisms at both the macroscopic and microscopic scales. This paper is organized as follows: in Section 2, the elastic model of silica glass is briefly reviewed. Subsequently, the densification model developed in this work is presented. Then, the calibration procedure using hydrostatic compression on thick samples (spheres) is detailed. Section 3 reviews the fracture model used to model the cracking behavior of silica glass indented with a Vickers tip. In Section 4, the Vickers indentation process in silica glass is simulated using the discrete element method first, and then the discrete-continuum coupling method, between the discrete element method (DEM) and the constrained natural element method (CNEM) [START_REF] Jebahi | [END_REF], is introduced to reduce the side effects and the computational time. Low indentation forces are applied, initially, to study the densification behavior in the absence of fracturing. Next, high indentation forces are applied to study the influence of densification on the cracking behavior. Section 5 presents the conclusions.

Silica behavior modeling

Elasticity

The variation of the discrete element method (DEM) used here is that developed by André et al. [18] and implemented in the GranOO workbench (software) 1 . The elastic behavior of silica glass is modeled using cylindrical elastic "cohesive beams" having geometric (length l µ and dimensionless radius rµ , defined as the ratio of the beam radius to the mean particle radius) and mechanical properties (Young's modulus E µ and Poisson's ratio ν µ ) between the adjacent discrete elements (Fig. 1). The Euler-Bernoulli analytical model of a beam [START_REF] Timoshenko | History of Strength of Materials: With a Brief Account of the History of Theory of Elasticity and Theory of Structure[END_REF] is used to compute the beam forces and torques on the discrete elements connected to this beam ((1), ( 2), ( 3) and ( 4)). Hereafter, the subscript "µ" denotes microscopic quantities.

F DE1 µ = +E µ S µ ∆l µ l µ x + 6E µ I µ l 2 µ (-(θ 2z + θ 1z )y + (θ 2y + θ 1y )z) (1)
F DE2 µ = -E µ S µ ∆l µ l µ x - 6E µ I µ l 2 µ (-(θ 2z + θ 1z )y + (θ 2y + θ 1y )z) (2)
T DE1 µ = + G µ Io µ l µ (θ 2x + θ 1x )x - 2E µ I µ l µ ((θ 2y + 2θ 1y )y + (θ 2z + 2θ 1z )z) (3)
T DE2 µ = - G µ Io µ l µ (θ 2x + θ 1x )x - 2E µ I µ l µ ((2θ 2y + θ 1y )y + (2θ 2z + θ 1z )z) (4)
Where:

• R(x, y, z) is the beam local frame where x is the beam axis.

• F DE1 µ and F DE2 µ are the beam force reactions acting on discrete elements 1 and 2 (connected to this beam bond).

• T DE1 µ and T DE2 µ are the beam torque reactions acting on discrete elements 1 and 2.

• l µ and ∆l µ are the initial beam length and the longitudinal deflection.

• θ 1 (θ 1x , θ 1y , θ 1z ) and θ 2 (θ 2x , θ 2y , θ 2z ) are the rotations of the beam cross sections expressed in the beam local frame.

• S µ , Io µ and I µ are the beam cross-sectional area, the polar moment of inertia and the moment of inertia with respect to the y and z axes, respectively.

• E µ and G µ are the Young's and shear moduli, respectively.

Table 1 presents the microscopic elastic cohesive beam properties that produce the expected macroscopic elastic behavior of silica glass. These properties have been determined by calibration tests as explained in Reference [18].

Densification

Densification model

This subsection details how the densification behavior of silica glass is modeled. As explained in Section 1, several studies examined the densification of glasses, in particular silica glass. Nearly all previous studies agreed that this behavior takes place under high hydrostatic pressures. However, the effects of shearing on this phenomenon remain a central issue of several scientific studies on glasses. In this work, it is assumed that the effects of shearing on the densification are insignificant and these are dominated by the hydrostatic pressure effects. Therefore, densification of silica glass will be supposed to occur only under high hydrostatic pressures. There are many physical and mechanical reasons that support this assumption:

• According to the characteristics of materials, two main mechanisms can explain a permanent deformation due to shear stresses: (i) In the case of crystalline solids, plasticity is due to dislocation slip, defined as crystallographic defects within the crystalline structure. Silica glass which is an amorphous material lacks the long-range order characteristic of a crystal. Therefore, it is very difficult for dislocation to develop and move within its irregular structure [START_REF] Kanel | [END_REF]. (ii) In the case of brittle materials, such as rock and concrete, permanent deformation is caused predominantly by slip at micro-cracks. This generally weakens the material in question and degrades its mechanical properties [START_REF] Lubliner | Plasticity theory[END_REF]. That is not the case for silica glass. In fact, permanent deformation in silica glass enhances the mechanical properties [START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF][START_REF] Keryvin | Contribution à l'étude des mécanismes de déformation et de fissuration des verres[END_REF].

• Because silica glass is an isotropic material, permanent deformation due to shear stress is volume-conservative. Therefore, even if a shear stress can cause silica glass plasticity, this does not affect the densification level (volume change at given pressure). For an infinitesimal strain, the volume change can be expressed as:

dV V 0 = J -1 trace(ε)
where dV is the volume change, V 0 is the initial volume, J is the determinant of the Jacobian matrix and ε is the strain tensor. If only shear stresses are applied, trace(ε) = 0, and therefore dV = 0.

Based on the previous assumption (the shearing effects on the densification are insignificant compared with those of hydrostatic pressure), a densification model adapted for discrete element models is developed in this work. A "spring-slider" system is connected in series with each cohesive beam linking two particles, as shown in Figure 2.

When the compression stress in the associated beam σ c beam reaches the microscopic densification pressure σ d µ , the slider activates and the densification mechanism takes place. The densification level beyond σ d µ is controlled by the microscopic tangential modulus E t µ , which is adjustable (Fig. 3). Finally, to model the saturation stage, the slider stops slipping above a certain limit controlled by the maximum microscopic permanent deformation of cohesive bonds ε

p max µ .
As shown in the literature [START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF][START_REF] Keryvin | Contribution à l'étude des mécanismes de déformation et de fissuration des verres[END_REF][START_REF] Zha | [END_REF], the macroscopic mechanical properties of silica glass increase with densification. According to Ji [START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF], Young's modulus (E) and Poisson's ratio (ν) increase by 46% and 56%, respectively, at the end of the densification process (Fig. 4). For the sake of simplicity, in the present work, the variations in these two parameters with hydrostatic pressure (Fig. 4) are approximated by linear piecewise functions. The mechanical properties remain constant until the densification pressure P d is reached, then increase linearly between P d and P s (the saturation pressure) and remain constant above P s . To take into account the macroscopic mechanical properties dependence with pressure in the densification model, it is necessary to find the variations in the microscopic cohesive beam properties, which implies such macroscopic variations. In other words, it is necessary to find relationships between the macroscopic and microscopic (of cohesive beams) mechanical properties. Based on the work of André et al. [18], the following conclusions can be drawn:

• The microscopic Poisson's ratio of cohesive beams ν µ has no influence on the macroscopic elastic behavior.

The latter depends only on the dimensionless cohesive beam radius rµ and the microscopic Young's modulus E µ . • The macroscopic Poisson's ratio ν depends only on rµ . The function relating these two parameters may be approximated by a second-order polynomial.

• The macroscopic Young's modulus E is a function of both rµ and E µ . For a given value of rµ , E depends linearly on E µ , whereas it is a second-order function of rµ for a constant E µ .

According to these facts, the relationships between the microscopic and macroscopic parameters can be written as in ( 5) and ( 6):

rµ = a 0 + a 1 ν + a 2 ν 2 (5)
           E(E µ , rµ ) = b 0 (r µ ) + b 1 (r µ ) E µ b 0 (r µ ) = b 00 + b 01 rµ + b 02 r2 µ b 1 (r µ ) = b 10 + b 11 rµ + b 12 r2 µ (6)
The various parameters of Equations ( 5) and ( 6) were determined by calibration tests. Knowing of these relationships (( 5) and ( 6)), the variations of the microscopic properties of silica glass with pressure can be determined and taken into account in the densification model. Figure 5 presents the cohesive bond behavior taking into account the variation of the microscopic mechanical properties with densification. As can be shown in this Figure, the relaxation slope (current microscopic Young's modulus) increases with the densification state.

Calibration of the microscopic densification parameters

The densification model involves three microscopic parameters which are the microscopic densification pressure σ d µ , the microscopic tangential modulus E t µ and the maximum permanent deformation of the beam bonds ε p max µ . The macroscopic densification parameters are the densification pressure P d , the saturation pressure P s and the volume change ratio ∆V V 0 (or density change ratio ∆ρ ρ 0 ). Several studies have measured the densification of silica glass, and different values of the macroscopic densification parameters can be found in the literature (Table 2). In the present work, Ji's results [START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF], which were obtained using thick samples, are used (Table 2) (because Ji conducted a complete experimental study of densification and Vickers indentation for several types of glasses). In this section, the microscopic densification parameters will be calibrated numerically to find the desired macroscopic densification behavior. A 100 mm-diameter sphere model (thick model) subjected to hydrostatic compression is used for the calibration study. One major problem of the discrete methods is the dependence of results on the discrete element number (N p ) used to discretize the numerical sample. Therefore, the influence of this parameter on the densification response is first studied. The sphere model was discretized using various numbers of discrete elements N p (from 100 to 20 000). The microscopic densification parameters were fixed as follows: the microscopic densification pressure σ d µ = 5 GPa, the microscopic tangential modulus E t µ = 0 GPa and the maximum permanent deformation of the beam bonds ε p max µ = 0.05. Figure 6 presents the results of this study. For a small number of discrete elements (N p ), the densification response fluctuates greatly with this parameter (Fig. 6a). These fluctuations can be reduced by increasing N p . Beyond N p = 10 000, the densification response becomes very weakly affected by the discrete element number (Fig. 6b).

In the remainder of this section, N p = 10 000 discrete elements will be used to discretize the numerical samples.

To simplify the calibration process, the influence of the microscopic densification parameters on the macroscopic silica behavior is studied. Figure 7a shows that the macroscopic densification pressure P d depends linearly on the microscopic one σ d µ . In addition, as observed in Figures 7b and7c, this parameter is effectively independent of the other microscopic properties, i.e. E t µ and ε p max µ . Figure 8 shows the influence of the various microscopic parameters on the macroscopic saturation pressure P s . This parameter depends linearly on the microscopic densification pressure σ d µ , and it depends on the microscopic tangential modulus E t µ (around 12% variation for E t µ up to 100 GPa ). However, it is effectively independent of the maximum permanent deformation of the beam bonds ε p max µ (Fig. 8c), only 5% P s variation for ε p max µ up to 0.12. Concerning the third macroscopic parameter ∆V V 0 , Figure 9 shows that this parameter depends only on E t µ and ε p max µ and it is effectively independent of σ d µ . After studying the influence of the microscopic parameters on the macroscopic densification behavior, calibration becomes fairly easy. Because the macroscopic densification pressure P d depends only on the microscopic one σ d µ , this parameter must be calibrated first. Subsequently, the microscopic tangent modulus E t µ can be calibrated to ensure the expected value of the macroscopic saturation pressure. Finally, the maximum permanent deformation of the beam bonds ε p max µ can be set to find the appropriate global volume change. Table 3 gives the calibrated microscopic densification parameters that must be used to ensure the correct macroscopic densification behavior of silica glass.

Validation

This subsection attempts to validate the microscopic densification parameters obtained by calibration (Table 3). The same geometric model as for the calibration study is used in this section. The various microscopic densification parameters of cohesive beams are set as indicated in Table 3. This model is subjected to various hydrostatic Hydrostatic compression, P (GP a) pressures. The associated permanent deformations are measured and compared to Ji's [START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF] experimental results (Fig. 10). As can be observed in Figure 10, the numerical results are in agreement with the experimental ones.

Fracture

The approaches most commonly used to model fracture in discrete methods, e.g. DEM [18,[START_REF] André | Modélisation par éléments discrets des phases d'ébauchage et de doucissage de la silice[END_REF], are based on the computation of bond strain [28] or stress [29]. The bond breaks if the associated strain (or stress) exceeds a given value, generally determined by calibration tests. On the macroscopic scale, this criterion is moderately satisfactory [START_REF] André | Modélisation par éléments discrets des phases d'ébauchage et de doucissage de la silice[END_REF]. However, on the microscopic scale, it cannot reproduce correctly the cracking mechanisms [START_REF] André | [END_REF]. Recently, another criterion based on the computation of an equivalent of the stress tensor at each particle has been proposed by André et al. [20]. This criterion is detailed in the remainder of this section. The equivalent Cauchy stress is computed with the help of the notion of virial stress which is the most commonly used definition of stress in discrete particle systems, particularly in Molecular Dynamics (MD) [21]. This "stress" includes two parts [21]: the first part depends on the mass and the velocity (or, in some versions, the fluctuation of the velocity) of particles; the second part depends on inter-particle forces and particle positions, providing a continuum measure for the internal mechanical interactions between particles. In this definition, the average virial stress over a volume Ω around a particle i is given by ( 7)

Π = 1 Ω (-m i ui ⊗ ui + 1 2 j i r i j ⊗ f i j ) ( 7
)
where m i and ui are, respectively, the mass and the velocity of the particle i, r i j = r j -r i , r i is the position vector of particle i, f i j is the inter-particle force acting on particle i exerted by particle j and ⊗ denotes the tensor product. The virial stress ( 7) is widely used to compute an equivalent of Cauchy stress in discrete systems. Recently, Zhou [START_REF] Zhou | Proc. R. Soc. Lond. A[END_REF] demonstrated that this quantity is not a measure for mechanical forces between material points and cannot be regarded as a measure of mechanical stress. Zhou proposed another formulation to compute the average mechanical stress, including only the second part of (7). For a region Ω around a particle i, the average stress is:

σ = 1 2Ω i∈Ω j i r i j ⊗ f i j (8)
The expressions ( 7) and ( 8) are developed for Molecular Dynamics (MD), where the inter-particle forces are a function of the Lennard-Jones potential Φ(r i j ) and the distance between the particles. In this instance, the two expressions of stress ( 7) and ( 8) lead to symmetric tensors. However, this cannot be generalized to all other discrete models. The variation of the DEM used here, where the particles are linked by cohesive beams, does not ensure the symmetry condition of the resulting tensors. André et al. [START_REF] André | [END_REF] proposed a slight modification of ( 8) to restore the symmetry condition of the stress tensor as follows:

σ = 1 2Ω i∈Ω j i 1 2 (r i j ⊗ f i j + f i j ⊗ r i j ) (9)
An equivalent Cauchy stress tensor for each discrete element i is computed using [START_REF] Ji | [END_REF], where Ω = Ω i , the volume of i (10).

σi = 1 2Ω i j i 1 2 (r i j ⊗ f i j + f i j ⊗ r i j ) ( 10
)
According to Griffith theory [30,31], the fracture of brittle solids initiates under tensile stress. Therefore, the new criterion postulates that a discrete element i is released from its neighbors when the hydrostatic stress ( σhyd

i = 1 3 trace( σi )
) is positive (tension state in i) and exceeds the microscopic fracture strength (σ hyd f ). In this case, all the beam bonds linking this particle to its neighbors break and do not intervene any more to compute the inter-particle forces and torques in the next time step. The microscopic fracture strength of silica glass, determined by calibration tests, is σ hyd f = 67 MPa.

Vickers indentation

As shown previously, silica glass experiences permanent densification under high compressive pressure. Comparable hydrostatic stresses can be reached beneath a sharp indenter during indentation loading. This section studies Vickers indentations in silica glass using various levels of indentation forces. First, low indentation forces are used to study the densification behavior at micro-scales (without fracture). Subsequently, high indentation forces are used to study to densification influence on the cracking behavior of silica glass. A hemispherical model with a radius of 12 µm is used for the Vickers indentation studies (Fig. 11). This model, whose the spherical part (dome) is fixed, is discretized using N p = 100 000 discrete elements (i.e., the average radius of the discrete elements is around 0.2 µm) to capture the microscopic mechanisms of the densification. It is indented with a square-based pyramidal indenter having included face angles of 136 • (Vickers indenter) and 2 µm in height (Fig. 11). This indenter undergoes a moderate displacement in the z direction until the reacting force reaches the desired indentation force. Table 4 presents a comparison between the experimental [START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF]32] and the present numerical indentation results for low indentation forces (F indent = 0.1 N and F indent = 0.5 N). For F indent = 0.1 N, the DEM simulation yielded a relatively good estimate of the volume of the indentation print V -and the volume of the piled-up material around the indenter V + . The densified volume V d which is a function of V -and V + compare favorably with the experimental results [32,[START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF] (in comparison to previous numerical studies [START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF]). Furthermore, a good numerical result was obtained for the diagonal length of the indentation print (L) which is used to compute the hardness. However, the indentation depth D was low compared with the experimental results [32,[START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF]. This can be explained by the fact that in the experimental tests, micro-cracks can initiate under the indenter tip (zone of high stress concentration), causing measurement errors in D. The experimental value of D can includes not only the depth of the indentation but also the depth of the crack. For F indent = 0.5 N, the densification state becomes more significant and a far greater densification print is obtained (Fig. 13). Compared with the experimental results, the DEM numerical simulation provided a relatively good estimate of the volume of the piled-up material around the indenter V + and the diagonal length L of the indentation print. Concerning D, the analysis of the F indent = 0.1 N case remains valid. However, this time, the volume of indentation print V -is slightly high. This can be explained by the side-effects. In the case of F indent = 0.5 N, the densified zone is relatively important compared with the global geometric model. Consequently, the effects of the overstated fixation conditions applied to the dome (the spherical part of the geometric model) become significant. This can increase the region where the densification pressure P d can be reached and then V -. One way to reduce the side effects is to use larger domains. However, such a solution can amplify considerably the computational time, because a very fine discretization must be applied to capture the microscopic densification effects. To circumvent the computational time problem, the discrete-continuum coupling model is used, a hybrid between the discrete element method and the constrained natural element method (CNEM-DEM coupling) [START_REF] Jebahi | [END_REF]. The discrete element model is applied in the indentation region where the densification pressure (P d ) can be reached, whereas the continuum model (CNEM) is applied over the remainder of the body being modeled. Figures 12 presents the coupling model which is a box of 30 µm on each side. The discrete subdomain where the DEM is applied is the same as in the previous study (DEM computation). The remainder of the box is modeled using CNEM and discretized with 10 000 nodes. A 2 µm-thick bridging zone is used between these two models to ensure correct wave propagation between them. Because the CNEM domain is far from the indentation region, it will experience only elastic deformations. Therefore, only the elastic behavior is modeled in this region. As shown in Table 4, in the case of F indent = 0.1 N, the DEM and CNEM-DEM coupling results are practically the same. However, in the case of F indent = 0.5 N where the side-effects became significant, the CNEM-DEM coupling results are better than those obtained using the discrete model (DEM) only. Furthermore, the CNEM-DEM coupling improved the shape of the densified zone as shown in Figure 14 which presents a cut-away view of the indented model in the case of F indent = 0.5 N. It can be observed that the densified zone obtained from the coupling method (Fig. 14c) is closer to Ji's experimental results (Fig. 14a) [START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF] than the densified zone obtained from the DEM computations (Fig. 14b). Indeed, it is less deep and broader in the former (coupling result) than in the latter (DEM result). The previous paragraph showed that the densification model, developed in this work, gives a relatively good quantitative estimate of silica glass densification at micro-scales when low indentation forces are applied (in comparison to previous numerical studies [START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF]). Next, a high indentation force (F indent = 30 N) is applied to study the influence of densification on silica glass cracking behavior. V -and V + are, respectively, the volumes of the indentation print and of the piled-up material around the indent, V d is the densified volume, D is the indentation depth and L is the diagonal length of the indentation print.

A model similar to the one presented in Figure 11 is used for this study. Because the indentation force is high, the region where the densification pressure P d can be reached becomes large. Hence, the radius of the model is increased to 120 µm. As shown in previous studies [33,34], brittle materials indented with a spherical indenter produce a cone crack. However, when a Vickers tip is used, normal glasses produce median, radial and lateral cracks. Instead, silica glass has a strong tendency to form a cone crack, even when indented with a Vickers tip. Because of its important densification behavior, a spherical densified zone is formed beneath the Vickers indenter which in turn operates as a spherical indenter, so that a cone crack is set up that can accompany median, radial and lateral cracks. Figure 15 presents the numerical cracking response of silica glass indented with a Vickers tip at the beginning of fracture (before the cracking becomes unstable). In this figure, only the discrete elements where the fracture criterion is reached are shown. Qualitatively, the result is in good agreement with those obtained in other studies [6,5]. Moreover, the radius of the cone crack obtained numerically (r cone num = 41.14 µm) is in good agreement with the experimental works of Arora et al. [6] (r cone exp ≈ 40 µm).

Conclusion

This paper attempted to reproduce numerically the anomalous behavior of silica glass and increase the comprehension of this behavior. Understanding permanent deformations in this type of glass represents a considerable challenge. An assumption has been made that the densification takes place only under high compressive pressure and that no shearing stresses influence the final permanent deformation levels. Based on this assumption, a densification model, adapted for discrete element methods, is developed in this work. The densification behavior is modeled by a "spring-slider" system incorporated between beam-linked particles at the micro scale. This system activates only under compression stress. The densification model involves three microscopic parameters which can be determined by calibration tests. As shown previously, the calibration process of these parameters is fairly easy. The macroscopic densification pressure P d is a linear function of microscopic densification pressure σ d µ only. Therefore, σ d µ can be determined independently from the other microscopic parameters. Subsequently, the microscopic tangential modulus can be calibrated taking into account the macroscopic saturation pressure P s . Finally, the maximum permanent deformation can be set to adjust the macroscopic volume change ∆V V 0 (or the density change ∆ρ ρ 0 ). It is shown that, after calibration, this model gives good macroscopic densification results compared with those obtained experimentally (Fig. 10). To validate this model at microscopic scales, simulations of the Vickers indentation process were performed with three different indentation forces. First, low forces (0.1 N and 0.5 N) were used to validate microscopically the densification model in the absence of fracture. A comparison of these results with experimental ones [START_REF] Ji | Mécanique et physique de l'indentation du verre[END_REF] shows that this model is also validated at microscopic scales (when microscopic loading is applied) and represents an advancement in the comprehension of permanent deformation in silica glass. The influence of the densification behavior on the cracking response was studied by applying a high indentation force (30 N). As can be observed in Figure (Fig. 15), the silica glass forms primarily a cone crack. This is due to the densification behavior, which leads to the formation of a spherical densified region. This, in turn, behaves as a spherical indenter that will induce in most brittle materials a cone crack at the beginning of fracture. When the indentation force further increases, other types of cracking appear and induce the complete destruction of the sample. The developments in this work have been implemented in the GranOO workbench, which now provides a robust numerical tool to simulate brittle continuous media. For more details, the reader is referred to http://www.granoo.org.

Figure 1 :

 1 Figure 1: Elastic beam model between two adjacent discrete elements

Figure 2 :

 2 Figure 2: Densification model between two adjacent discrete elements

Figure 3 :Figure 4 :

 34 Figure 3: Equivalent cohesive bond behavior

Figure 5 :

 5 Figure 5: Behavior of equivalent cohesive bond, taking into account the variation of the mechanical properties with densification

Table 2 :Figure 6 :

 26 Figure 6: Sensitivity of the macroscopic densification behavior to the number of discrete elements (N p )

(a) E t µ = 0

 0 GPa, ε p µ = 0.05 (b) σ d µ = 5 GPa, ε p µ = 0.05 (c) σ d µ = 5 GPa, E t µ = 0 GPa

Figure 8 :

 8 Figure 8: Variation of macroscopic saturation pressure (P s ) with microscopic densification parameters: microscopic densification pressure (σ d µ ), microscopic tangential modulus (E t µ ) and maximum permanent deformation of the beam bonds (ε p max µ )

Figure 10 :

 10 Figure 10: Density change with pressure: comparison between experimental [8, 9] and numerical results.

Figure 11 :

 11 Figure 11: DEM indentation model

Figure 12 :

 12 Figure 12: CNEM-DEM coupling indentation model

Figure 13 :

 13 Figure 13: Top view of indentation print (blue: negative volume change) and piled-up material around indenter (red: positive volume change). Left: indentation force F indent = 0.1 N; right: indentation force F indent = 0.5 N.

Figure 14 :

 14 Figure 14: Cut-away view of the indented region; F indent = 0.5 N; top: experimental results from Ji [8], middle: DEM numerical results; bottom: CNEM-DEM coupling numerical result

Figure 15 :

 15 Figure 15: Cracking pattern of silica glass indented with a Vickers tip; F indent = 30 N; only the discrete elements where the fracture criterion is reached are shown in (b) and (c).

Table 1 :

 1 The microscopic elastic parameters of silica glass (obtained by calibration)

	Young's Modulus Poisson's ratio Dimensionless radius
	E µ = 414 GPa	ν µ = 0.3	rµ = 0.57

Table 3 :

 3 The microscopic densification parameters of silica glass (obtained by calibration)

	Micro. dens. pressure Micro. tangential modulus Max permanent def in the bonds
	σ d µ = 7.13 GPa		E t µ = 98.82 GPa			ε	p max µ	= 0.067
	0	5	10	15	20	25	30	35	40

Table 4 :

 4 Vickers indentation results: low indentation forcesV -(µm 3 ) V + (µm 3 ) V d (µm 3 ) D (µm) L (µm)

		Num (DEM): Present	1.304	0.136	1.168	0.159	4.758
		Num (CNEM-DEM): Present	1.302	0.131	1.171	0.159	4.758
	F indent = 0.1 N	Num [8]	--	--	0.058	--	--
		Exp [8]	0.727	0.108	0.605	0.377	4.707
		Exp [32]	1.034	0.018	0.952	0.421	5.157
		Numerical (DEM)	14.534	1.174	13.360	0.397	9.318
		Numerical (CNEM-DEM)	13.628	0.920	12.708	0.409	9.957
	F indent = 0.5 N	Num [8]	--	--	0.61	--	--
		Exp [8]	9.971	1.138	8.805	0.908 10.371
	/	Experimental error [8]	±4%	±7%	±10%	±6%	±2%

http://www.granoo.org


K Saleh

A Fernandez

O Llopis

G Cibiel

Fiber ring resonators with Q factors in excess of 10 10 for Time and Frequency Applications

Keywords: Microwave-optics, optoelectronic oscillator, optical resonator, optical fiber, phase noise I

The design of ultra high Q factor fiber ring resonators is described. The resonators are characterized with a microwave frequency modulation technique, which allows a precise measurement of their performance. Optoelectronic oscillators are then realized. The phase noise performance is compared to the results obtained with lower Q factor resonators, and the relatively high intra-cavity bandwidth noise is discussed.

INTRODUCTION

In the field of low phase noise microwave sources, the opto-electronic oscillator (OEO) is now very popular. Since the first delay line feedback OEOs in the 1990s [START_REF] Yao | High frequency optical subcarrier generator[END_REF][START_REF] Yao | Progress in the optoelectronic oscillator -a ten year anniversary review[END_REF], based on optical fiber spools of a few kilometers, the design of these systems has evolved towards more compact oscillators based on optical resonators. These resonators can be whispering gallery mode resonators, such as mini disk monocrystaline resonators [START_REF] Savchenkov | Phase noise of whispering gallery photonic hyper-parametric microwave oscillators[END_REF][START_REF] Ilchenko | Miniature Oscillators Based on Optical Whispering Gallery Mode Resonators[END_REF][START_REF] Merrer | Characterization technique of optical whispering gallery mode resonators in the microwave frequency domain for optoelectronic oscillators[END_REF], passive fiber loop resonators [START_REF] Merrer | High-Q Optical Resonators for Stabilization of High Spectral Purity Microwave Oscillators[END_REF][START_REF] Saleh | Optoelectronic oscillator based on fiber ring resonator: overall system optimization and phase noise reduction[END_REF] or active fiber loops resonators, i.e. coupled optical and electronic oscillator (COEO) [START_REF] Yao | Coupled optoelectronic oscillators for generating both RF signal and optical pulses[END_REF][START_REF] Yu | Ultralow-noise mode-locked laser with coupled optoelectronic oscillator configuration[END_REF].

Recent work [START_REF] Saleh | Optoelectronic oscillator based on fiber ring resonator: overall system optimization and phase noise reduction[END_REF] has shown the capability of fiber ring resonators (FRRs) of featuring very high optical quality factors (Q opt ). These FRRs are compact and easy to use. However, their high quality factor, combined with significant intra-cavity power enhancement factor, lead to the generation of nonlinear optical effects inside the resonator, especially Rayleigh and Brillouin scattering [START_REF] Saleh | Optical Scattering Induced Noise in Fiber Ring Resonators and Optoelectronic Oscillators[END_REF]. When the FRR is used to stabilize the frequency of an optoelectronic oscillator (OEO), these nonlinear optical effects degrade the OEO phase noise. In order to eliminate a great part of these effects, an optical isolator has been inserted inside a 100 m long FRR that featured a Q opt of 2 10 9 . The use of this FRR to stabilize the frequency of a 10 GHz OEO has resulted in a phase noise level as low as -52 dBc/Hz at 10 Hz offset frequency [START_REF] Saleh | Optoelectronic oscillator based on fiber ring resonator: overall system optimization and phase noise reduction[END_REF]. The investigations described in this paper tried to extend the performance of such a system to higher Q factors, in order to reach a lower phase noise level for the OEO. However, even if we have been successful in increasing the Q factor of our resonators above 10 10 , the resulting phase noise is higher than the one obtained with the resonator with Q opt of 2 10 9 . This problem is discussed and some hypothesis on the origin of this noise are given.

II. OEO BASED ON FIBER RING RESONATORS

The fiber ring resonator structure is depicted in Figure 1. It uses one or two optical couplers, preferably two if a filtered output is needed (resonator used in transmission mode). The parameters which can be optimized to get the desired Q factor and insertion losses values are the coupling coefficients (and couplers residual losses), the residual losses of all the elements of the loop (including fiber splices) and the fiber type. Of course, a trade-off has to be found between high Q factor and low insertion losses. The previously realized ring resonators were in the range of 20 m to 100 m length. Because it is the optical resonator bandwidth which is preserved in the microwave application, the microwave Q factor of the OEO is degraded by the ratio of the microwave to optical carrier frequency, which is about 20 000 for a 10 GHz application (optical carrier at 194 THz, or 1550 nm) . Therefore Q opt must be higher than 10 9 in order for this technology to be competitive with the best microwave resonators, which means a fiber ring larger than 10 m and a low coupling factor, or a fiber ring of 100 m if lossy elements are included in the loop. With this approach, most of our designs were based on 20 m fiber rings, featuring Q factors of about 3 10 9 , using either SMF fiber or PM fiber. A 100 m ring has also been designed, using PM fiber and including an isolator in order to prevent the growing up of Stimulated Brillouin Scattering (SBS), which appears at relatively low input power because of the large intra cavity power enhancement factor inside the resonant loop.

These resonators have been used to stabilize various OEOs. All the OEOs have the same topology, which is described in previously published papers [START_REF] Saleh | Optoelectronic oscillator based on fiber ring resonator: overall system optimization and phase noise reduction[END_REF][START_REF] Saleh | Optical Scattering Induced Noise in Fiber Ring Resonators and Optoelectronic Oscillators[END_REF] and depicted in Figure 2. The laser is first locked to the resonator frequency thanks to a Pound Drever Hall (PDH) feedback loop ; the oscillation is started thanks to a second loop which is a microwave loop including one or two amplifiers, a phase shifter and a dielectric resonator filter centered at 10.2 GHz. The phase noise results obtained on 10.2 GHz OEOs realized using this system are depicted in Figure 3. In each case, the oscillator system is the same and only the optical resonator is removed and replaced by another. Three types of resonators are used in the OEOs described in this The best result is obtained with the 100 m FRR OEO, featuring -52 dBc/Hz at 10 Hz offset, and -128 dBc/Hz at 10 kHz. With this resonator, we have been able to remove the main nonlinear optical effects limiting the performance of the 20 m FRR OEO. The phase noise of this oscillator is only limited by the RF amplifiers phase noise, and the measured spectrum can be described by a Leeson type approach. III. ULTRA HIGH Q RESONATORS DESIGN Our goal was then to increase again the Q factor of our FRR, above 10 10 (which means an equivalent Q factor at 10 GHz above 500 000). To this purpose, two ways can be followed : realizing a 100 m long FRR without isolator in the loop, or realizing a much longer FRR, such as a 1000 m FRR, and keeping the isolator which prevents some of the optical nonlinear effects (particularly the growing up of SBS).

In such a design, a trade-off has to be found between Q factor enhancement and the optical losses through the resonator (high Q means low coupling, which induces losses). Optical losses higher than 5 dB means RF losses higher than 10 dB, which are added to the losses of the microwave to optical and optical to microwave conversion which are generally high. We have thus chosen to limit optical IL to about 5 dB. The resonator is modeled either using an equivalent circuit approach on ADS software, or using a dedicated Matlab approach (both give the same results). The resonator performance is simulated varying the coupling factors of the optical couplers (the other parameters being the residual losses of the optical devices provided by the devices manufacturers).

For the 100 m resonator without isolator, a coupling factor close to 1 % meets our requirements in Q and insertion losses (IL). In the case of the 1000 m resonator, a coupling factor of about 10 % is necessary to maintain relatively low IL. This coupling factor also corresponds to the maximum intra cavity power enhancement factor.

The resonators have then been realized, and included in metallic boxes with passive thermal isolation. The devices characterization is performed still using the PDH stabilization loop : the laser is stabilized on one resonance, and the lateral resonances are measured using a microwave vector network analyser (VNA), i.e. the green loop in Figure 2 is replaced by the VNA.

Figure 4 depicts the transmission curve (S 21 ) obtained on the 100 m resonator (actually, more precisely, it is a 120 m resonator). A Q opt of 1.25 10 10 has been measured, with optical IL of 7.5 dB (a bit degraded compared to our simulation).

Figure 5 depicts the same measurement performed on the 1000 m resonator. In this case, a Q opt of 1.7 10 10 has been reached, with optical IL of 4.1 dB, which was perfectly in agreement with our simulated data. The measured S 21 phase curve is also depicted in this figure (it is this parameter which stabilizes the oscillator), and the high slope looks very clean and can be easily measured using our microwave modulation technique.

IV. OEOS WITH Q IN EXCESS OF 10 10

These two resonators have then been included in the OEO system described in Figure 2. A 10.2 GHz oscillation has been obtained, and its phase noise has been measured. Frequency (GHz)

Figure 5 : Measurement of the transmission response of the 1000 m resonator, using the PDH lock and a frequency modulation approach (vector network analyser). The measured Q factor and IL are in this case identical to the ones simulated. The S21 phase response is also depicted in this figure .   The phase noise of the 120 FRR OEO is depicted in Figure 6. As shown in the figure, a sudden increase of the noise is observed inside the resonator half bandwidth. We were expecting such a phenomenon, because this resonator is not protected against optical nonlinear effects, and the intra cavity power enhancement factor is huge. We have attempted to reduce its influence by reducing the resonator optical input power, but the best phase noise level we could get with the PDH locking still going on was -40 dBc/Hz @ 10 Hz, which has to be compared to the -52 dBc/Hz @ 10 Hz obtained with the preceding 100 m FRR with a Q factor 6 times lower (but made of PM fiber and including an isolator).

We were more confident in the 1000 m FRR OEO, because of the isolator in the loop, but once again a sudden increase of the noise is observed inside the resonator half bandwidth, has shown in Figure 7. Also shown in this figure is the effect of a reduction of the resonator input power : the changes observed in the carrier to noise ratio outside the resonator bandwidth are proportional to this power, but the phase noise improvement inside the cavity is not (but there is effectively an improvement at low power). The FSR for this resonator is also very small : 200 kHz, which means that the dielectric resonator filter is no more sufficient to filter only one resonant mode. However, at low optical power, we still were able to get a single oscillation, even if the lateral modes can be seen on the spectra every 200 kHz offset.

For this 1000 m FRR OEO, the best phase noise results is -27 dBc/Hz @ 10 Hz and -105 dBc/Hz @ 10 kHz, which is still far from the best result described in section II. The reason may lie in two points : either the 1000 m section between two isolators is too large to prevent the built up of Rayleigh scattering noise, or a polarization noise is predominant because conventional SMF fiber has been used in this design instead of PM fiber (because of the cost of 1000 m PM fiber). ). Again, a large increase of the noise is clearly seen inside the resonator half bandwidth, in spite of an attempt of reducing the intra cavity power (the three spectra presented correspond to different optical power levels at the resonator input).

V. CONCLUSION

Optical resonators with Q factor in excess of 10 10 have been successfully designed and realized using the fiber ring resonator approach. The stabilization of a laser on one resonance has been demonstrated, using the Pound Drever Hall technique. The measurement of the resonators characteristics have thus been performed thanks to a microwave modulation and a vector network analyser. The measured data (Q factor, insertion losses, out of band rejection…) are very close to the simulated data for these ultra high Q resonators.

OEOs have then been realized with these two resonators and an approach based on the PDH laser stabilization on the resonator. The measured phase noise spectra show a sudden increase of the noise inside the resonator half bandwidth, which is the signature of a noise induced by the resonator itself.

Further studies will try to cancel this parasitic noise contribution, and also to improve the performance of the best OEO realized up to now, which was based on a 2 10 9 Q factor resonator.

Figure 1 :

 1 Figure 1 : Fiber ring Resonator topology

FMFigure 2 :

 2 photodiode

Figure 3 :

 3 Figure 3 : Phase noise of previously realized OEOs, stabilized with the same PDH circuit, including one OEO based on a CaF2 mini-disk (Q = 10 8 ), one OEO realized with a 20 m FRR (Q = 3.5 10 9 ) and one OEO stabilized with a 100 m FRR which includes an isolator (Q = 2 10 9 ).

Figure 4 :

 4 Figure 4 : Measurement of the transmission response near 10.2 GHz of the 120 m resonator, using the PDH lock and a frequency modulation approach (microwave network analyser).

  RF = 8.2 dB Δν -3dB =12.3 kHz Q opt =1.68 10 10

Figure 6 :Figure 7 :

 67 Figure 6 : Phase noise of the 10.2 GHz OEO based on the 120 m resonator (Q = 1.25 10 10 ). The two curves correspond to two different loop parameters. A large increase of the noise is clearly seen inside the resonator half bandwidth.

ACKNOWLEDGMENT

This study is sponsored by the French National Space Centre (CNES) and the "Direction Générale de l'Armement" (DGA), French MoD procurement agency, under the contract ANR-11-ASTR-0029 MINOTOR.