Marco Stucchi

Fabrizio Pittorino
email: fabrizio.pittorino@unibocconi.it

Matteo Di Volo

Alessandro Vezzani

Matteo Di Volo

Raffaella Burioni

Order symmetry breaking and broad distribution of events in spiking neural networks with continuous membrane potential

Keywords: neuron dynamics, synchronization, oscillations, criticality, chaos, symmetry breaking 2010 MSC: 00-01, 99-00

Introduction

Networks of living neurons exhibit complex collective patterns of activity, ranging from asynchronous to synchronous firing even in absence of stimuli [START_REF] Van Vreeswijk | Partial synchronization in populations of pulse-coupled oscillators[END_REF][START_REF] Fox | Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging[END_REF]. The different degree of synchrony associated with these regimes is a fundamental property of the network dynamics and it is at the basis of the collective oscillatory rhythms observed in neural activity in the cortex [START_REF] Buzsáki | Mechanisms of gamma oscillations[END_REF]. A particularly interesting dynamical regime on the border between synchronous and asynchronous dynamics is represented by the outbursts of activity called neural avalanches [START_REF] Beggs | Neuronal avalanches in neocortical circuits[END_REF][START_REF] Petermann | Spontaneous cortical activity in awake monkeys composed of neuronal avalanches[END_REF][START_REF] Bellay | Irregular spiking of pyramidal neurons organizes as scale-invariant neuronal avalanches in the awake state[END_REF]. These are cascades of activity clustered in time and separated by periods of quiescence and they have been observed both in vitro and in vivo in different cortical areas and across many species [START_REF] Beggs | Neuronal avalanches in neocortical circuits[END_REF][START_REF] Persi | Modeling of synchronized bursting events: the importance of inhomogeneity[END_REF][START_REF] Bellay | Irregular spiking of pyramidal neurons organizes as scale-invariant neuronal avalanches in the awake state[END_REF][START_REF] Shriki | Neuronal avalanches in the resting meg of the human brain[END_REF]. Neuronal avalanches represent a peculiar mode of activity in the cortex, whose function and origin is still much debated [START_REF] Buzsaki | Rhythms of the Brain[END_REF]. Their name stems from the fact that the distribution of the number of spiking neurons, i.e the size of the outburst, and the distribution of their durations lack a characteristic scale and are compatible with power laws, whose exponents fulfill scaling relations [START_REF] Friedman | Universal critical dynamics in high resolution neuronal avalanche data[END_REF][START_REF] Muñoz | Avalanche and spreading exponents in systems with absorbing states[END_REF]. The dynamical activity in the avalanches regime also exhibits peculiar correlation properties and optimal computational performance [START_REF] Mora | Are biological systems poised at criticality?[END_REF][START_REF] Muñoz | Colloquium: Criticality and dynamical scaling in living systems[END_REF][START_REF] Wilting | Operating in a reverberating regime enables rapid tuning of network states to task requirements[END_REF]. These findings suggest that neural networks might operate close to a critical point, whose type and origin is still debated [START_REF] Chialvo | Emergent complex neural dynamics[END_REF] and could be associated to a hybrid synchronization regime [START_REF] Di Santo | Landau-ginzburg theory of cortex dynamics: Scale-free avalanches emerge at the edge of synchronization[END_REF][START_REF] Buendía | Hybrid-type synchronization transitions: where marginal coherence, scalefree avalanches, and bistability live together[END_REF].

Another open point is what are the minimal microscopic ingredients at the level of neuron dynamics that can account for such an emergent dynamical regime. In a recent paper [START_REF] Pittorino | Chaos and correlated avalanches in excitatory neural networks with synaptic plasticity[END_REF] we have shown that an avalanches dynamical regime, with peculiar synchronization and correlation properties, can emerge in a sparse heterogeneous (in the sense of neuronal connectivity) network of leaky-integrate and fire (LIF) neurons [START_REF] Brunel | Quantitative investigations of electrical nerve excitation treated as polarization[END_REF][START_REF] Brunel | Lapicque's 1907 paper: from frogs to integrate-and-fire[END_REF] with short-term synaptic plasticity modeled through the Tsodyks-Uziel-Markram mechanism [START_REF] Tsodyks | The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability[END_REF][START_REF] Tsodyks | Synchrony generation in recurrent networks with frequency-dependent synapses[END_REF]. Interestingly, this regime emerges from a transition to chaos observed in the corresponding homogeneous network (all-to-all connections), which turns into an avalanches regime in the presence of heterogeneous connectivities, suggesting a fundamental role of neuronal heterogeneity. While the primary role of excitatory connections and a regulatory mechanism, such as synaptic plasticity, have been shown to be key ingredients, it is not understood to what extent this scenario depends on the structural organization of neuron connections and on the choice of single neuron intrinsic dynamics.

In particular, even if the LIF model has the great advantage of being analytically integrable, it suffers from a discontinuity at the spiking time that may have strong effects on a dynamical regime like the one we are studying, where neurons have a high synchronous activity during the burst. In order to investigate the relevance of single neuron time scales, we introduce a new LIF model which is analytically integrable and removes the hard discontinuity in voltage dynamics at the spiking event, thus mimicking the continuous dynamic of the action potential (see [START_REF] Izhikevich | Dynamical systems in neuroscience[END_REF] for a discussion of models with continuous membrane potential). This requires an additional variable which, using the physics jargon, plays the role of a mass in the context of a second order differential equation. According to its continuous features, we call such model the c-LIF (continuous LIF).

We study a homogeneous fully connected network of excitatory c-LIF neurons with short-term synaptic plasticity and we show that it exhibits a complex bursting dynamical regime with characteristic scale free dynamics, even without heterogeneity in the network. The microscopic inertial timescale of the c-LIF has a desynchronizing effect leading to avalanches of activity when acting in the region of parameters where the homogeneous network has a chaotic regime. By using a mean field approach, in this simple model we can trace back the mechanism for the onset of the broad activity distributions and show that the avalanche dynamical regime results from an order symmetry breaking, i.e. from neurons crossings in the sequence of firing times of the network. The possibility to break this symmetry has already been investigated in [START_REF] Kirst | Sequential desynchronization in networks of spiking neurons with partial reset[END_REF][START_REF] Kirst | Partial reset in pulse-coupled oscillators[END_REF] in the presence of partial resets. Here we show that order symmetry breaking is a crucial ingredient to produce bursts with an internal structure and scale free avalanches distributions in a homogeneous network, with peculiar syncronization properties. This also clarifies that networks with homogeneous connections can give rise to an avalanche regime, even if neurons are all identical, both in their intrinsic dynamics and in the way they are connected to the surrounding network.

The article is organized as follows. In Sec. 2 we briefly recall the LIF equation and we extend it to the c-LIF single neuron model. We then consider N globally coupled excitatory c-LIF neurons, with the dynamics of their interactions following the TUM model for short-term synaptic plasticity, accounting for the synaptic transmitter dynamics. In Sec. 3 we recall the dynamical features of networks of LIF neurons with short-term synaptic plasticity and in Sec. [START_REF] Beggs | Neuronal avalanches in neocortical circuits[END_REF] we investigate the dynamics of the c-LIF fully connected network in the parameter space, observing quasi-synchronous, asynchronous and bursty chaotic dynamics. In Sec. 5 we link the complex collective dynamics to a symmetry breaking of the spiking order of neurons. In Sec. 6 we study the synchronization properties of the model. Finally, Sec. 7 is reserved for discussions and future perspectives.

The c-LIF model

The LIF neuron model [START_REF] Dayan | Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems[END_REF] retains the basic features of a biological neuron, while being easy to analyse math-ematically and computationally. The dynamical equation is:

τ 1 Vi = E c -V i + R in GI i , (1)
where τ 1 is the membrane time constant, R in is the membrane resistance, GI i is the synaptic current received by neuron i from all its pre-synaptic neurons (G is a parameter regulating its strength) and E c is the contribution of a constant external current (multiplied by R in ). Whenever the potential V i (t) reaches the threshold value V th it is reset to V r and a spike is sent to all post-synaptic neurons. It is useful to rescale time with the membrane time constant τ 1 and to introduce the dimensionless quantities

a = Ec-Vr V th -Vr , g = RinG V th -Vr , v = V -Vr V th -Vr
, so that v r = 0 and v th = 1. With this rescaling, the difference between the potential at the threshold and at the reset is set to 1 and the unit of v is of the order of 10mV (considering V th = -55mV and V r = -65mV ). The membrane potential v i of a neuron i with an interaction current I i (t) evolves in time according to the differential equation:

τ 1 vi (t) = a -v i (t) + gI i (t) . ( 2
)
The potential is reset to 0 at times t i (m) when it reaches the threshold v i (t i (m)) = 1.

The LIF model describes with an instantaneous reset of the voltage V a process that in reality determines the continuous and non-instantaneous rise and the decay of the membrane voltage V during the emission of a spike. This continuous dynamics of the voltage V at the spike emission has been explained by the seminal work of Hodgkin and Huxley [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], through a four dimensional set of differential equations accounting for the kinetics of potassium and sodium channels in neurons' membranes. In the direction of simplifying this description, but still accounting for its fundamental features, many models have been proposed. Some of them were able to reduce the dimensionality of the system and to keep track of the voltage dynamics at the spike emission, as in the Morris Lecar [START_REF] Morris | Voltage oscillations in the barnacle giant muscle fiber[END_REF] or FitzHugh-Nagumo models [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF]. These models are nevertheless not integrable, at variance with the LIF model. In this work we explicitly include a decay time scale at the spiking event and we focus on the effects of this time scale on the collective dynamics of a neural network. Instead of modelling potassium and sodium channels, which would add further nonlinearities in the model, we limit ourselves to consider the simplest modification of a LIF model able to keep track of a continuous evolution of the membrane voltage during the spike emission, the c-LIF. This model has also the big advantage to be integrable. Interestingly, dynamical regimes similar to those observed in our model have recently been observed in a non-integrable LIF neural model with an additional non-linear term that describe the dynamics of the voltage at the spike emission [START_REF] Wang | Collective firing patterns of neuronal networks with short-term synaptic plasticity[END_REF].

The c-LIF model evolves according to an extension of the LIF equation: We call the term on the left "inertial term", in analogy with classical mechanics. Eq. (3) reduces to LIF Eq. ( 2) in the limit τ m → 0.

τ 2 m vi = -τ 1 vi + a -v i + gI i (t). (3)
With the introduction of the time scale τ m , the equation for the membrane potential is turned into a second order differential equation: as a consequence the state of neuron i is characterized both by the value of its membrane potential v i and by its first derivative vi . At the spiking times t i (n), when v i (t i (n)) reaches the threshold v i (t i (n)) = 1, we impose to the membrane potential v i (t i (n)) to be continuous and we set its first derivative to the negative value:

vi

(t i (n)) = - τ 1 τ 2 m . (4)
With this choice, the membrane potential of an isolated c-LIF neuron decays at the firing event in a time of order

τ 2 m
τ1 , as shown in Fig. 1. In particular, for τ m → 0 we recover the behavior of the LIF model. This is clarified in the inset (b), where we show that for τ m → 0 the minimum of the continuous potential tends to zero (i.e. the reset potential of the LIF model) and the time difference δ t between the minimum and the spiking event vanishes. In analogy with the standard LIF model, we call

ISI i (n) = t i (n+1)-t i (n) the n-th inter-spike interval of the c-LIF neuron i.
Notice that this rule is not equivalent to a manual reset of the voltage to a value different from zero, known as partial reset rule [START_REF] Bugmann | Role of temporal integration and fluctuation detection in the highly irregular firing of a leaky integrator neuron model with partial reset[END_REF][START_REF] Leng | Common stochastic inputs induce neuronal transient synchronization with partial reset[END_REF]. Indeed, in our model we are dealing with interactions that are continuous in time, so that neurons continuously cross the threshold of firing and, as a result, the partial reset rule would have not effects. Therefore the mechanism that we are describing here is actually only due to the decay time and not to a time varying reset rule.

Short-term plasticity tunes the dynamics of neurotransmitter, temporarily affecting neural responses. The widely used TUM model, which we consider here, has been shown to correctly reproduce experimental results for neurotransmitter dynamics in different conditions [START_REF] Tsodyks | Synchrony generation in recurrent networks with frequency-dependent synapses[END_REF]. According to the TUM dynamics, the amount of neurotransmitter at each synapse is finite and the state of the i-th synapse is determined by three variables: y i (t), z i (t) and x i (t) = 1-y i (t)-z i (t), representing the fraction of neurotransmitters respectively in the available, active or inactive states. The synaptic currents y i evolve in time according to the following dynamical equations:

ẏi (t) = - y i (t) τ in + u(1 -y i (t) -z i (t))S i (t) (5) żi (t) = y i (t) τ in - z i (t) τ r , (6)
where S i (t) = m δ(t -t i (m)) represents the spike train of neuron i (t i (m) is the time when neuron i fires its m-th spike). The parameter u regulates the percentage of available resources released at every spike and τ in and τ r are the characteristic times of decay of active resources y i and recovery of inactive ones z i , respectively. The synapses are depressed (i.e. the effect of the presynaptic neuron on the postsynaptic one is reduced) if the firing rate is high, since neurotransmitters need time to recover from the inactive to the available state. We assume that all parameters appearing in equations ( 5) and ( 6) are independent of neuron index and that all efferent synapses of a given neuron follow the same evolution.

We consider an all-to-all network of N c-LIF neurons, where the input current received by neuron i reads:

gI i (t) = gY (t) = g N j y j , (7)
i.e. every neuron i feels the same mean field Y . Eq. ( 7) can be modified and extended to take into account different topologies, including the inhomogeneities of the system [START_REF] Burioni | Average synaptic activity and neural networks topology: a global inverse problem[END_REF] and inhibition [START_REF] Bertolotti | Synchronization and long-time memory in neural networks with inhibitory hubs and synaptic plasticity[END_REF]. Eq.s (3-7) define our network model which we will study numerically by an event-driven approach [START_REF] Brette | Exact Simulation of Integrate-and-Fire Models with Exponential Currents[END_REF][START_REF] Di Volo | Synchronous dynamics in the presence of short-term plasticity[END_REF] (see Appendix for details). To summarize, the model is characterized by four characteristic time scales: τ 1 sets the typical oscillation time of an isolated single neuron, τ 2 m /τ 1 is the decay time of the action potential, τ in rules the decay time of synaptic resources and τ r their typical recovery time. The constant a determines the shape of the neural oscillation and g is the interaction strength. Setting τ 1 = 1, our time unit is set to a typical value of the order of 10ms. Moreover, we expect that τ in and τ 2 m /τ 1 are much shorter than τ 1 , while τ r has a larger value. In particular, we set the parameters to values that are typically employed in the literature [START_REF] Pittorino | Chaos and correlated avalanches in excitatory neural networks with synaptic plasticity[END_REF]: τ in = 10 -3 , τ r = 10, a = 1.3 and u = 0.5 unless otherwise stated. We vary the values of τ m and g to analyse the different dynamical regimes of the network. We also set g > 0 since we are only considering excitatory neurons. Notice that we choose these values for the parameters a and τ in because they allow to clearly observe and characterise the dynamical effects of the decay timescale over a broad range of values of τ m and g. Nevertheless, the dynamical regimes that we describe hereafter can be observed even for different values of the parameters. In particular, the avalanche dynamics is observed also at biologically realistic physical values of decay time τ 2 m /τ 1 (i.e. comparable with the duration of the action potential, see Section 6).

3. The dynamics in the limit τ m → 0

In the limit τ m → 0, the c-LIF reduces to the classic LIF model with synaptic plasticity, that we recently investigated in [START_REF] Pittorino | Chaos and correlated avalanches in excitatory neural networks with synaptic plasticity[END_REF]. In this limit the dynamics of the allto-all network is perfectly synchronous [START_REF] Di Volo | Synchronous dynamics in the presence of short-term plasticity[END_REF] and can be reduced to an effective single neuron equation. Despite its simplicity, the mean field LIF with TUM synaptic plasticity features an interesting variety of dynamical regimes, depending on the parameters. In particular in Fig. 2 we show the bifurcation diagram of the interspike intervals for different values of the coupling g. The bifurcation diagram shows a chaotic regime for intermediate values of g. This regime emerges because of the competition between the slow time scale τ 1 , which describes the intrinsic dynamics of the neuron, and the fast interaction time scale τ in of the synaptic resources. In particular, an heuristic estimate of the effective time scales and interaction strength derived in [START_REF] Pittorino | Chaos and correlated avalanches in excitatory neural networks with synaptic plasticity[END_REF] shows that for τ in τ 1 the dynamics should display a chaotic behavior if K 1 τ r /τ in < g < K 2 (τ r τ 1 )/τ 2 in where the numerical prefactors K 1 and K 2 are expected to barely depend on the system parameters. In this intermediate regime, indeed, in Eq. 2 the fast interaction term gI i (t) = gY (t) and the slow evolving terms a -v(t) have the same order of magnitude and one cannot treat one of them perturbatively. As a consequence, the simultaneous presence of two different timescales gives rise to chaos via a period doubling mechanism.

When introducing disorder (e.g. in the coupling g or, equivalently in the topology of the network), the synchronous chaotic regime is desynchronized and turned into a richer regime with power-law distributed avalanches [START_REF] Pittorino | Chaos and correlated avalanches in excitatory neural networks with synaptic plasticity[END_REF]. Neurons with heterogenous interactions have different intrinsic speed and this unpacks the synchronous events observed in the all-to-all model. As a result, the underlying chaotic dynamics in the presence of heterogeneity leads to the emergence of population avalanches. In the rest of the paper we fix τ in = 10 -3 τ 1 = 1, in such a way to span the synchronous, chaotic and periodic regimes by varying g. We will show that the introduction of a finite τ m > 0 gives rise to a desynchronizing effect that produces power-law distributed avalanches even in the homogeneous all-to-all network. g

ISI(n) ISI(n)

The role of the time scale τ m in collective dynamics

The timescale τ m has strong effects on the collective dynamics of the homogeneous all-to-all case, which we now investigate. We focus on the conservation of the order of neurons spiking events [START_REF] Kirst | Sequential desynchronization in networks of spiking neurons with partial reset[END_REF][START_REF] Kirst | Partial reset in pulse-coupled oscillators[END_REF], showing that this observable is tightly related to the emergence of bursty dynamics, characterized by avalanches of activity with heavy tailed distribution.

Let us first note that the mean field LIF model, beside the already mentioned discontinuity of the membrane potential, displays another unphysical feature, that is the conservation of fire ordering. Indeed, at τ m = 0 the firing order is trivially maintained by the dynamical evolution. If neuron j + 1 fires after neuron j, we can define ∆v j (t ) = v j (t ) -v j+1 (t ), which is positive at the firing time t of neuron j + 1. From the mean field equations with τ m = 0, we obtain for the dynamics between two firing events of neurons j+1 and j, that τ 1 ∆v j (t) = -∆v j (t). Hence, v j (t) remains larger than v j+1 (t) until the next firing event where the neuron j is again going to fire before neuron j + 1, preserving therefore the firing order during the evolution. On the other hand, for τ m > 0 we have τ 2 m ∆v j (t) = -τ 1 ∆v j (t) -∆v j (t), therefore the sign of ∆v j (t) is not preserved between firing events and the firing ordering can be violated.

Hereafter, hence, in the stationary regime we will order neurons according to their firing events and we will study if their firing order is preserved in time. Raster plots are plotted according to this labelling of neurons. In what follows, we study (otherwise stated) an all-to-all network of N = 500 neurons, discarding the first 5•10 6 firing events in order to reach the stationary regime.

Small coupling. We first study the case of small coupling g, meaning g < K 1 τ r /τ in where the all-to-all network with τ m = 0 displays a fully synchronous non-chaotic collective dynamics. Fig. 3 shows that population spiking events are now not completely synchronous as soon as τ m > 0. In particular, in the left column we show the raster plot (for each firing event we plot a dot at the index of the corresponding firing neuron), and in the right column we plot the return map ISI i (n + 1) vs. ISI i (n) for the neuron i. Let us remark that even if the interspike interval explicitly depends on the neuron, the curve ISI i (n + 1) vs. ISI i (n) does not because the system is homogeneous. Hereafter, therefore, we will drop the index i. For very small τ m (upper panels) the spikes are almost synchronous, but the ISI n+1 vs. ISI n plot shows the presence of a very small closed attractor, accounting for the presence of a quasi-periodic behavior with small oscillations of the ISI [START_REF] Van Vreeswijk | Partial synchronization in populations of pulse-coupled oscillators[END_REF][START_REF] Burioni | Microscopic mechanism for self-organized quasiperiodicity in random networks of nonlinear oscillators[END_REF]. For intermediate values of τ m , middle panels, the quasi-periodic behavior becomes more pronounced and it can be directly observed in the raster plot where firing now is clearly asynchronous. The effect of increasing τ m is to lead to a stronger and stronger desynchronization, up to higher values of τ m where the dynamics falls in the well known splay state (lower panels), an extremely regular asynchronous regime where neurons firing times are equally spaced [START_REF] Van Vreeswijk | Partial synchronization in populations of pulse-coupled oscillators[END_REF]. Remarkably, the firing order is conserved in all these regimes, independently of the value of τ m , as it is shown by the raster plots in Fig. 3 (we also explicitly check that all pairs of consecutive firing neurons do not exchange their mutual firing order along the dynamics).

Intermediate coupling. The τ m = 0 case for intermediate values of g is characterized by synchronous chaotic dynam- As shown by this latter panel the dynamics is chaotic, but there is no firing order symmetry breaking. Lower panel c): in correspondence to this dynamical phase characterized by the preservation of neurons firing order, the distribution of events is not broad (at increasing values of N = {500, 2000, 8000}).

ics, and also in this regime a time constant τ m > 0 induces a desynchronization. At very small τ m (see Fig. 4 upper panels), the synchronous burst structure is only slightly modified. One can indeed observe the presence of clearly separated almost synchronous bursts where each neuron fires exactly once. The firing order in this regime is preserved during the evolution. By increasing τ m (see Fig. 5 upper panels) bursts are so broadened that they overlap. The ISI(n + 1) vs. ISI(n) plot shows the presence of a non-trivial attractor, pointing at an underlying chaotic dynamics, where the firing order is not preserved. At larger τ m the system enters into a quasi-periodic regime, clearly visible from the one-dimensional closed attractor in the ISI(n + 1) vs. ISI(n) in Fig. 6). Here, the neurons firing order is again preserved. Finally, at even larger τ m the dynamics becomes asynchronous and the system falls in the splay state (see Fig. 7).

Large coupling. In the case of large coupling, i.e. for g > K 2 (τ r τ 1 )/τ 2 in the dynamics at τ m = 0 is again synchronous and periodic. However now a positive value of τ m > 0 fully destabilizes the synchronous dynamics and we enter in an asynchronous regime. Indeed, by perturbing the synchronous state at τ m = 0 the dynamics reaches the asynchronous splay state for any value of τ m > 0; in particular, the smaller τ 2 m , the longer it takes to the system to desynchronize. Clearly, in the splay state, the firing order is conserved for any values of τ m .

Order symmetry breaking and broad distributions of events

As already observed, in the regime described in Fig. 5 fire ordering is not conserved by the dynamical evolution. Table 1: Measured values of the minimum ISI value ISI min at τm = 0 compared with the value τ * = τ 2 m /τ 1 at which firing order begins to be violated in the model with τm > 0, for three different values of the coupling g.

Two conditions seem to be necessary to violate the ordering: first the fully synchronized dynamics at τ m = 0 should be chaotic, second a desynchronizing mechanism is needed allowing for the overlap of different bursts, and this can be caused by disorder in the couplings [START_REF] Pittorino | Chaos and correlated avalanches in excitatory neural networks with synaptic plasticity[END_REF] or by the inertial time scale τ m . In the latter case, the reset time of a neuron τ 2 m /τ 1 provides an estimate of the duration of a broadened burst and, as far as this time scale is much smaller than the minimal time between two consecutive bursts, no overtaking is possible. In this perspective, we verify that the minimal ISI min measured at τ m = 0 (i.e. the minimal ISI values appearing in the vertical bars in the bifurcation diagram in Fig. 2) provides the order of magnitude of τ 2 m /τ 1 = τ * at which we start to observe a violation of firing order in the presence of τ m . For different values of g in the chaotic regime, these values coincide with what is observed in Fig. 2. Measured values of ISI min and τ * as a function of g are reported in Table 1.

In the dynamical regime where fire ordering is broken, we now measure the size and duration distribution of avalanches. We apply the avalanche definition adopted in experimental settings [START_REF] Beggs | Neuronal avalanches in neocortical circuits[END_REF]. An event or avalanche is defined τs-1 , as expected for a system near criticality.

as a set of consecutive neuronal spikes such that all interspike intervals t(spike n+1 ) -t(spike n ) are smaller than a threshold δ [START_REF] Beggs | Neuronal avalanches in neocortical circuits[END_REF], setting δ of the same order of magnitude of the average interspike interval of the network. The size s of an event is then defined as the number of spikes it contains. Analogously the duration of an event is the time elapsed between the first and the last spikes of the event. We then consider the distributions of sizes h(s), durations h(T ) and mean avalanche size s as a function of the duration T , as shown in Fig. 8.

Figures 45678summarize our results. A broad distribution in the event size is observed only in the regime where the firing order is not conserved. For small values of τ 2 m (Fig. 4) indeed we observe an irregular dynamics but the neurons fire in quasi-synchronous bursts and typically one observes only large events while small events are rare. In the regular quasi-periodic regime (Fig. 6) and in the splay state (Fig. 7) no large events are present and a cut-off in the characteristic size h(s) is observed. Instead, the overlap of the different bursts and the overtaking in firing order that characterizes the regime in Fig 5 gives rise to a broad distribution in h(s). In particular, such distribution is compatible with a power law h(s) ∼ s -τs with τ s 1.5. The distribution does not appear to depend significantly on the threshold δ and it scales with N to higher values of s as N is increased, since a larger N allows the presence of larger avalanches. Analogously, as shown in Fig. 8, the distribution of avalanches durations T shows the same power law behavior, with exponent τ T = 1.9. Interestingly, the exponents satisfy the scaling relation τ = τ T -1 τs-1 , where τ = 1.8 is the exponent for the average size as a function of fixed duration, as expected for a system near criticality [START_REF] Muñoz | Avalanche and spreading exponents in systems with absorbing states[END_REF].

Synchronization and the Kuramoto parameter

In this section we consider the Kuramoto parameter [START_REF] Kuramoto | Self-entrainment of a population of coupled nonlinear oscillators[END_REF] as a measure of the synchronicity level of the neurons, to characterize the different dynamical regimes of the model. The Kuramoto parameter reads:

R(t) = 1 N N i=1 e iφi(t)
where φ i (t) is the phase of neuron i at time t:

φ i (t) = 2π t -t i (m) t i (m + 1) -t i (m)
and t i (m) is the m-th spike of neuron i and t ∈ [t i (m), t i (m+ 1)]. The Kuramoto parameter R and its average in time R take values in the interval [0, 1], ranging from an asynchronous system at R = 0 to a perfectly synchronous one at R = 1.

Figure 9 shows a summary of the different behaviors. Continuous lines corresponds to the Kuramoto parameter and dashed vertical lines delimit the region where overtaking in neuron dynamics are observed and broad distributions of avalanches are present. At small g = 1.5 • 10 4 no overtaking is present. In the quasi-periodic regime the average Kuramoto parameter R is positive, then a transition to the splay state asyncronous regime with R = 0 is observed at larger τ m . For the intermediate values of g at large τ m we also observe a transition from the quasi periodic regime to the splay state where R = 0. At lower τ m vertical dashed lines delimit the region where exchanges in neuron firing order occur and avalanche sizes and durations are in the scaling regime. Interestingly in this regime the average Kuramoto parameter is quite large .918 R .995 signaling the presence of almost synchronous events even if the sizes of such events fluctuate. Then, for even smaller τ m , the dynamics is chaotic and almost synchronous as indicated by the Kuramoto parameter ( R ≈ 0.999). Finally for larger g = 5 • 10 6 the average Kuramoto parameter R vanishes for any value of τ m > 0: this means that the system is in an asynchronous splay state.

We show the standard deviation ∆R of the Kuramoto parameter in the inset of Fig. 9. We observe as expected ∆R ≈ 0 in the quasi-periodic and splay states, while we measure non-negligible values .025

∆R

.12, indicating strong fluctuations of the Kuramoto parameter, in correspondence to the region where exchanges in neuron firing and broad distributions of avalanches are present. Strong fluctuations may be present also at higher values of τ 2 m , after the firing crossing region, where a quasi-periodic phase with a complex two-dimensional orbit is observed. This phase is characterised by lower synchronization, R 0.8. A similar scenario concerning the fluctuations of this order parameter has been reported in [START_REF] Pittorino | Chaos and correlated avalanches in excitatory neural networks with synaptic plasticity[END_REF]. for different values of the coupling parameter g spanning the whole quasi-synchronous to asynchronous transition (each curve corresponds to a value of g according to the legend). Here a = 1.3 and τ in = 10 -3 . Regions of the parameter space characterized by exchanges in neuron firing order and broad distribution of events are delimited by vertical lines (corresponding to the values of g according to their color). Lower panel: Kuramoto parameter and its standard deviation (inset) for a = 1.01, τ in = 10 -2 , g = 4 • 10 2 . Notice that for these value of the parameters, we observe the bursty phase (delimited by the dashed vertical lines) up to a value of τ 2 m = 4 • 10 -2 , which corresponds to a decay time of 0.4ms in physical units (if we consider τ 1 = 10ms).

The dynamical regime that we observe and characterize is expected to be present also for more realistic value of the involved timescales. The lower panel of Figure 9 shows that, for different values of the parameters, the bursty dynamical regime with avalanches is observed even at larger decay times τ 2 m /τ 1 , which are comparable with the realistic duration of the action potential. In particular, the bursty phase is present up to values of τ 2 m = 4 • 10 -2 , which corresponds to a decay time of 0.4ms in physical units if the membrane time constant is set to τ 1 = 10ms.

Conclusions

We have investigated the synchronization dynamics of a fully connected neural network with short-term synaptic plasticity with a new membrane potential time scale in the equations governing the single neuron dynamics. The c-LIF model that we study here avoids unphysical discontinuities in the membrane potential characterizing the usual LIF neurons, while still being integrable. We have shown that if the time constant τ m related to the new term is not too large the system maintains its variety of dynamical phases. For small values of the coupling g the network is in a quasi-periodic, quasi-synchronous state, while for large values the system is completely asynchronous. Interestingly, the desynchronizing effect of this time scale induces a more complex regime for intermediate values of g, in the presence of a chaotic dynamics. The bursting regime is characterized by a broad and robust distribution of event sizes and durations, even in a pure mean field model with homogenous couplings. We have numerically measured the exponents of these distributions and shown that they satisfy a scaling relation. The onset of this regime is related to the order symmetry breaking of the spike times of single neurons, which produces a non trivial dynamical regime even in the absence of disorder or heterogeneity. Notice that in the context of phase oscillators the presence of a mass term has been found to give rise to interesting dynamical regimes with breaking of symmetry in an homogeneous system, called chimeras [START_REF] Olmi | Intermittent chaotic chimeras for coupled rotators[END_REF]. In our neural model the mass term emerges naturally from the introduction of a continuous dynamics of the membrane voltage and is the essential factor giving rise to chaotic dynamics and avalanches as observed in experimental data of neuronal network.

Appendix A. Event-driven map

It is convenient for numerical simulations to transform the set of equations (3), ( 5), ( 6) and ( 7) into an eventdriven map [START_REF] Brette | Exact Simulation of Integrate-and-Fire Models with Exponential Currents[END_REF][START_REF] Di Volo | Synchronous dynamics in the presence of short-term plasticity[END_REF], which does not require the numerical integration of the neural equations. Indeed, it is possible to find the exact analytical solution of the LIF and c-LIF equations as they are linear between two consecutive spikes of the network. The idea is to find the analytical solution to these equations from the time immediately after the last spike of the network, t n , to the time immediately before the next spike of the network, t n+1 . As long as t n < t < t n+1 the last term on the right in eq. ( 5) does not contribute.

In order to make the equations readable, we define the following coefficients: v i,n = v i (t n ), vi,n = vi (t n ), Y n = g kY (t n ), ∆ = t-t n , α = τ

y i (t) = y i (t n )e -∆ τ in (A.2) z i (t) = z i (t n )e -∆ τr + τ r τ r -τ in y i (t n )(e -∆ τr -e -∆ τ in ) (A.
3)

The time of the next spike of the network is t n+1 = min i {t i } where the set {t i } is determined by imposing v i (t i ) = 1 in eq. (A.1). Once that t n+1 is known it is possible to update the state variables of the i-th neuron: v i (t n ) ← v i (t n+1 ), vi (t n ) ← vi (t n+1 ), y i (t n ) ← y i (t n+1 ), z i (t n ) ← z i (t n+1 ).

Moreover, for the particular set of neurons r spiking at time t n+1 the value of v has to be reset to -τ1 and the activated resources ux need to be added to y and consequently to the global field Y . See also [START_REF] Di Volo | Synchronous dynamics in the presence of short-term plasticity[END_REF][START_REF] Pittorino | Chaos and correlated avalanches in excitatory neural networks with synaptic plasticity[END_REF] for details on the event-driven simulation of the LIF model. By repeating these steps we obtain the full evolution of the system.

Figure 1 :

 1 Figure 1: Time evolution of the rescaled membrane potential v of a free (g = 0) c-LIF neuron (Eq. 3) compared to the usual LIF model (Eq. 2). Parameters: a = 1.3, τ 1 = 1 for both models and τ 2 m = 10 -2 for the c-LIF one. Inset (a): zoom on the membrane potential decay region of the c-LIF model to be compared with the LIF discontinuity. Vertical dashed lines highlight the decay time τ 2 m τ 1 of the isolated c-LIF neuron. Inset (b): in the c-LIF model, plot of the minimum Vm of the potential V and of the time difference δt between the spiking event and the subsequent minimum of V as a function of τm.

Figure 2 :

 2 Figure 2: Bifurcation diagram for the all-to-all LIF model with TUM short-term synaptic plasticity and τ in = 10 -3 , τm = 0. The attractor for the interspike interval of the network (ISI) is shown in function of the bifurcation parameter (the coupling constant g). Upper panel: full bifurcation diagram. Lower-left panel: zoom on the first transition, revealing a period-doubling transition to chaos. Lower-right panel: zoom on the second transition. Red boxes are zooming regions.

Figure 3 :

 3 Figure 3: Dynamical phases of an all-to-all network of N = 500 c-LIF neurons with TUM synaptic plasticity at a small value of the coupling, g = 1.5 • 10 4 . Left column: raster plots (we have discarded a transient of t = 7 • 10 3 ). Right column: plot of the (n + 1)-th interspike interval ISI(n + 1) as a function of the n-th interspike interval ISI(n) of a single neuron inside the network. Each row in the plot shows the dynamics for a different value of τ 2 m : τ 2 m = 6.9 • 10 -4 for panels a) and b) (quasi-periodic phase); τ 2 m = 6.8•10 -2 for panels c) and d) (quasi-periodic phase); τ 2 m = 1.3 • 10 -1 for panels e) and f) (splay state).

Figure 4 :

 4 Figure 4: Dynamical phase and distribution of events of an all-toall network of N = 500 c-LIF neurons with TUM synaptic plasticity at an intermediate value of the coupling, g = 1.0 • 10 5 and for a low value of τ 2 m = 1.4 • 10 -5 . Upper left, panel a): raster plot (we have discarded a transient of t = 2 • 10 3 ). Upper right, panel b): plot of the (n + 1)-th interspike interval ISI(n + 1) as a function of the n-th interspike interval ISI(n) of a single neuron in the network.As shown by this latter panel the dynamics is chaotic, but there is no firing order symmetry breaking. Lower panel c): in correspondence to this dynamical phase characterized by the preservation of neurons firing order, the distribution of events is not broad (at increasing values of N = {500, 2000, 8000}).

Figure 5 :

 5 Figure 5: Dynamical phase and distribution of events of an all-toall network of N = 500 c-LIF neurons with TUM synaptic plasticity at an intermediate value of the coupling, g = 1.0 • 10 5 and for an intermediate value of τ 2 m = 7.0 • 10 -4 . Upper left, panel a): raster plot (we have discarded a transient of t = 2 • 10 3 ). Upper right, panel b): plot of the (n + 1)-th interspike interval ISI(n + 1) as a function of the n-th interspike interval ISI(n) of a single neuron inside the network. As shown by this latter panel the dynamics is chaotic, and in this phase there is firing order symmetry breaking. Lower panel c): in correspondance to this dynamical phase characterized by the breaking of neurons firing order, there is a broad distribution of events (at increasing values of N = {500, 2000, 8000}).

Figure 6 :

 6 Figure 6: Dynamical phase and distribution of events of an all-toall network of N = 500 c-LIF neurons with TUM synaptic plasticity at an intermediate value of the coupling, g = 1.0 • 10 5 and for τ 2 m = 9.5 • 10 -3 . Upper left, panel a): raster plot (we have discarded a transient of t = 2 • 10 3 ). Upper right, panel b): plot of the (n + 1)-th interspike interval ISI(n + 1) as a function of the n-th interspike interval ISI(n) of a single neuron inside the network. As shown by this latter panel the dynamics is quasi-periodic, and the firing order becomes preserved again. Lower panel c): in correspondance to this dynamical phase characterized by the preservation of neurons firing order, the distribution of events is not broad (at increasing values of N = {500, 2000, 8000}).

Figure 7 :

 7 Figure 7: Dynamical phase and distribution of events of an all-toall network of N = 500 c-LIF neurons with TUM synaptic plasticity at an intermediate value of the coupling, g = 1.0 • 10 5 and for a high value of τ 2 m = 1.8 • 10 -2 . Upper left, panel a): raster plot (we have discarded a transient of t = 2 • 10 3 ). Upper right, panel b): plot of the (n + 1)-th interspike interval ISI(n + 1) as a function of the n-th interspike interval ISI(n) of a single neuron inside the network. As shown by this two panels the dynamics is a splay state, in which the firing order is preserved. Lower panel c): in the splay state, the distribution of events is not broad (at increasing values of N = {500, 2000, 8000}).

7 •

 7 10 -5 1.5 • 10 -4 4.6 • 10 5 2.6 • 10 -5 3.1 • 10 -5 10 6 1.4 • 10 -5 1.0 • 10 -5

Figure 8 :

 8 Figure 8: (color online) Size (s) and duration (T) distributions of all-to-all networks of cLIF neurons endowed with TUM synaptic plasticity for g = 1.0 • 10 5 and τ 2 m = 7.0 • 10 -4 . Upper panel: the size distribution is a power-law with exponent τs = 1.5. Middle panel: duration distribution, showing a power-law with τ T = 1.9. Lower panel: average size for fixed duration, showing a power law with exponent τ = 1.8. In all cases the power-law scaling holds for increasing values of N = {4000, 8000, 16000}). The values of the exponents are fitted over the range of the corresponding dashed line and satisfy the relation τ = τ T -1τs-1 , as expected for a system near criticality.

RFigure 9 :function of τ 2 m τ 1

 91 Figure 9: (color online) Upper panel: the time-averaged value R of the Kuramoto parameter and of its standard deviation ∆R as a function of τ 2 m τ 1

2 1 -4τ 2 m , β = τ1+α 2τ 2 m, γ = τ1-α 2τ 2 m,T 1 2 m τ 2 in α 3 -2Y n τ 2 in α 3 -

 22212323 = τ 2 m + τ in (τ in -τ 1 ), and we obtain the following form for the event-driven map of the c-LIF model:v i (t) =a + 2τ (τ 1 τ in -2τ 2 m ) in α + e -β∆ -aT 1 (τ 1 -α) + Y n τ in 2τ 2 m + τ in (-τ 1 + α) + T 1 2τ 2 m vi,n + v i,n (τ 1 -α) + e -γ∆ aT 1 (τ 1 + α) + Y n τ in -2τ 2 m + τ in (τ 1 + α) -T 1 2τ 2 m vi,n + v i,n (τ 1 + α) (τ 1 τ in -2τ 2 m ) 2 α e -∆ τ in 2Y n τ in α -e -β∆ [-2aT 1 + Y n τ in (τ 1 -2τ in + α) + T 1 (2v i,n + vi,n (τ 1 + α))] -e -γ∆ [2aT 1 + Y n τ in (-τ 1 + 2τ in + α) -T 1 (2v i,n + vi,n (τ 1 -α))] (A.1)


A Simon

R Pascaud

L Liard

O Pascal

Experimental Comparison of Wideband and Narrowband Plasma-based Microstrip Power Limiters

Keywords: plasma devices, gas discharge devices, low temperature plasma, microwave breakdown, nonlinear wave propagation

This article presents the experimental work performed as part of the use of the non linear plasma-wave interactions in planar circuits to design microwave power limiters. A first non-resonant plasma-based microwave power limiter composed of a pre-ionized plasma discharge integrated into a 50 Ω microstrip line has been experimentally characterized with steady-state and transient measurements. These results have been then confronted at each step to a resonant one using the same plasma discharge. This study highlights the effect of the resonance on the non linear plasma-wave interactions and its consequences on the behaviour and the performances of a plasma-based microwave power limiter.

I. INTRODUCTION

When a high power microwave (HPM) threat impinges on typical receiver systems, they can be strongly deteriorated [START_REF] Backstrom | Susceptibility of electronic systems to high-power microwaves: Summary of test experience[END_REF]. Microwave power limiters are commonly involved to protect them from these threats. For low HPM signals, active lumped elements such as PIN diodes or MEMS are used to design limiting functions [START_REF] Yang | A novel analysis of a Ku-band planar p-i-n diode limiter[END_REF], [START_REF] Shojaei-Asanjan | Tunable RF MEMS-based frequency-dependent power limiter[END_REF]. However, these integrated components remain, for most of them, fragile when the incident microwave power becomes more important.

Beyond a few tens of watts, the use of plasma discharges has proved its effectiveness in achieving microwave power limitation. In plasma-based microwave power limiters, plasma discharge arises when the induced electric field by the HPM is sufficiently intense to cause gas breakdown. The resulting plasma then interacts with the incident wave which leads to the limitation of the output power. This mechanism has been widely used in T/R (Transmit/Receive) tubes in order to protect RADAR receivers from high power microwave signals [START_REF] Smullin | Microwave Duplexers[END_REF]. However, T/R tubes can not be easily integrated into compact planar circuits.

Several studies have then focused on the integration of plasma discharges into planar circuits for power limitation [START_REF] Patel | Microstrip plasma limiter[END_REF]. However, the high limiting threshold power was really challenging to design operational systems. Consequently, resonant circuits [START_REF] Parsons | Microwave plasma formation within a 2d photonic crystal[END_REF], pre-ionization techniques, or both [START_REF] Pascaud | Low insertion loss microplasma-based limiter integrated into a microstrip bandpass filter[END_REF]- [START_REF] Semnani | A high-power widely tunable limiter utilizing an evanescent-mode cavity resonator loaded with a gas discharge tube[END_REF] were used to lower and/or control this threshold.

This work was supported by the Université Fédérale Toulouse-Midi-Pyrénées (UFTMiP), the Région Occitanie, and the Direction Générale de l'Armement (DGA) under Grant ANR-15-ASTR-0020-02, DIOMEDE project, within the framework of the 2015 ASTRID program. In this paper, we experimentally compare the behaviour and performance between non-resonant and resonant, that is to say wideband and narrowband, plasma-based microstrip power limiters.

II. EXPERIMENTAL SETUP

A dedicated experimental setup has been developed to perform steady-state and transient power measurements as shown in Fig. 1. A microwave unit allows to set up a power ramp and measure the input, output, and reflected powers by means of power sensors for the steady-state measurements. Regarding the transient measurements, a microwave switch is first inserted between the microwave generator and the power amplifier. Then, power sensors are replaced by a 4-channels fast oscilloscope. The microwave switch and the oscilloscope are synchronized with a control signal which pulse width and repetition frequency are equal to 1 ms and 100 Hz, respectively. Moreover, the experimental setup includes additional units that control the type of gas, its pressure, and the DC voltage and current for the pre-ionization of the microwave power limiter.

III. CIRCUITS DESIGN A. Non-resonant microstrip power limiter

Firstly, a non-resonant microwave power limiter has been designed according to the schematic presented in Fig. 2. It consists of a 50 Ω microstrip transmission line printed on a 1.524 mm thick Rogers RO4003C substrate with ε r = 3.55 and tan δ = 0.0021. This line is drilled in its center with a Φ ms = 1.5 mm diameter hole and connected to 50 Ω coaxial straight connectors at both ends. This microwave circuit is combined with a pre-ionization system using a DC plasma discharge [START_REF] Pascaud | Low insertion loss microplasma-based limiter integrated into a microstrip bandpass filter[END_REF]. The pre-ionization of the microstrip hole is done by a micro-hollow cathode discharge (MHCD) located under the ground plane. It has been optimized to work at a pressure of 10 Torr (i.e., 1333 Pa) in argon that is to say h M HCD = Φ M HCD = 1.5 mm. Fig. 3 presents the magnitude of the simulated and measured small-signal S-parameters without pre-ionization as a function of the frequency.

B. Resonant microstrip power limiter

A similar resonant microwave power limiter has been designed. As shown in Fig. 4, it consists of an electromagnetically coupled annular ring resonator in microstrip technology printed on the same substrate. The microstrip width and the hole diameter are the same as previously. Note that the hole has been drilled on the location where the electric field is maximum at the resonant frequency (i.e., ϕ = 0 • in Fig. 4). The ring radius r is equal to 11.7 mm while the stub length L s is 15.75 mm. Fig. 5 presents the magnitude of the simulated and measured small-signal S-parameters without pre-ionization as a function of the frequency. The resonant frequency of the ring resonator is equal to 2.5 GHz.

IV. EXPERIMENTAL RESULTS

A. Steady-state measurements

Measurements using a VNA have shown that the additional insertion loss due to pre-ionization remains below 0.05 dB for both circuits when the pre-ionization current I DC is under 5 mA. Consequently, a pre-ionization current of 3 mA is considered here. Such limiting devices then exhibit almost no insertion loss. Fig. 6a presents the measured output power P out as a function of the input power P in at 2, 3, and 4 GHz for the nonresonant circuit. It should be noted that without pre-ionization (I DC = 0 mA) no non-linear interaction is observed for the considered power range, that is to say up to 40 dBm. With a pre-ionization (I DC = 3 mA), a breakdown occurs and a non-linear behaviour is then observed. The pre-ionization system thus generates an initial electron density in the upper hole of the DUT which drastically decreases the required microwave electric field for gas breakdown [START_REF] Bandel | Effect of preionization on microwave breakdown in neon[END_REF]. The measured breakdown threshold, defined here as the 1 dB compression point, are equal to 29.9, 28.9, and 30.7 dBm at 2, 3, and 4 GHz, respectively. Almost no dependence with frequency is noticed here since for such a reduced frequency range, the breakdown electric field does not vary substantially at 10 Torr in argon [START_REF] Roth | Industrial Plasma Engineering[END_REF]. These measurements have been repeated varying the current flowing through the DC plasma discharge from 3 to 40 mA without any noticeable change in the limiting threshold power or dynamic range.

Fig. 6b presents the measured output power P out as a function of the input power P in at 2.5 GHz for the resonant case. As previously, with a DC pre-ionization current of 3 mA, a microwave breakdown occurs in the hole and the output power is then reduced. One can observe that the limiting threshold of 22.8 dBm is lower than the non-resonant case. Due to the microwave resonance of the ring resonator, the electric field is locally increased in the hole which lead to higher electric field for the same input power. Thus, the gas breakdown that happens for a given electric field value can occur for lower input power. Numerical simulations using Ansys HFSS has shown that the steady-state maximum electric field in the hole is 2.1 times larger for the resonant circuit than for the simple transmission line considering the same input power. Thus, the power must be multiplied by 4.2, or 6.2 dB larger, in the non-resonant circuit to obtain the same electric field as in the resonant circuit. This theoretical result is very close to the observed difference between the measured limiting thresholds.

To understand the interaction between the HPM signal and the plasma discharge, it is also interesting to observe the reflected, output, and dissipated power rates relative to the input power that are denoted R ref l , R out , and R dis , respectively. The dissipated power is expressed as P dis = P in -P out -P ref l .

Fig. 7a presents the measured power rates as a function of the input power at 3 GHz for the microstrip line. It is observed that the decrease of the output power rate R out is mostly due to the dissipation of the input power by the plasma. Considering for instance P in = 40 dBm, the measured output power rate R out is then reduced to 6.6 % while the reflected and dissipated power rates reach 18 % and 75.4 %, respectively.

Contrary to the non-resonant power limiter, the output power of the ring resonator is mainly reduced due to reflection as shown in Fig. 7b. One can observe that for P in between 22.8 dBm and 28 dBm the plasma dissipates most of the incoming power, but for P in larger than 28 dBm the reflected power rate increases until reaching 58 % at 40 dBm while the dissipated and output power rates decrease to 39.7 % and 2.3 %, respectively. The plasma, here, acts as a disruptive impedance that tends to mismatch the microstrip device [START_REF] Parsons | Microwave plasma formation within a 2d photonic crystal[END_REF].

B. Transient measurements

Transient high power measurements of both circuits have also been conducted. Fig. 8a shows the measured envelope of the input, reflected, and output microwave signals considering a modulated pulse of 3 GHz with a power of 36 dBm with I DC = 3 mA for the microstrip line. No significant change has been observed on the time response when varying I DC from 3 to 40 mA or the modulation frequency from 2 to 4 GHz. However, the response time of the power limiter clearly depends on the input power of the microwave pulse. The higher the input power, the faster the device responds, thus reacting in 87 µs at 36 dBm and 8.6 µs at 45 dBm for example. This response time has been defined as the time needed for the output voltage to reach to within 5 % of its final steady-state value.

Regarding the transient analysis of the resonant circuit, Fig. 8b shows the measured envelope of the input, reflected, and output microwave signals considering a modulated pulse of 2.5 GHz with a power of 36 dBm and I DC = 3 mA. Here, a 1 µs response time is measured which is 87 times lower than the non-resonant case with the same input power. Actually, the higher electric field due to the resonance leads to lower response time [START_REF] Löfgren | Breakdown phenomena in microwave transmit/receive switches[END_REF].

V. CONCLUSION

In this article, we experimentally compared the behaviour of non-resonant and resonant plasma-based microstrip power limiters. Between 2 and 4 GHz, the non-resonant power limiter dissipates most of the incoming power from around µs. The comparison with the resonant power limiter clearly highlights the impact of the microwave resonance on the power limitation. Thus, the resonant microwave power limiter has a lower limiting threshold of 22.8 dBm, reflects most of the incoming power, and reacts 87 times faster compared to the non-resonant case with the same input power. Such a study may be promising to understand the non linear plasma-wave interaction and design optimized plasma-based microwave power limiters depending on the electromagnetic environment.

Fig. 1 .

 1 Fig. 1. Experimental setup for steady-state (blue configuration) and transient (red configuration) high-power microwave measurements.

Fig. 2 .

 2 Fig. 2. Cut view of the non-resonant 50 Ω microstrip line combined with a pre-ionization system.

Fig. 3 .

 3 Fig. 3. Simulated and measured small-signal S-parameters of the nonresonant transmission line without pre-ionization.

Fig. 4 .

 4 Fig. 4. Front view of the 50 Ω characteristic impedance ring resonator in microstrip technology.

Fig. 5 .

 5 Fig. 5. Simulated and measured small-signal S-parameters of the resonant ring resonator without pre-ionization.

  (a) Microstrip transmission line. (b) Ring resonator.

Fig. 6 .

 6 Fig. 6. Measured output power Pout as a function of the input power P in for (a) the microstrip transmission line and (b) the ring resonator.

  (a) Microstrip transmission line. (b) Ring resonator.

Fig. 7 .

 7 Fig. 7. Output Rout, reflected R ref l and dissipated R diss power rates as a function of the input power P in for (a) the microstrip transmission line at 3 GHz and (b) the ring resonator at 2.5 GHz.


Experimental study of a tunable sub-wavelength plasma-based VHF-UHF resonator

Vincent Laquerbe and Romain Pascaud Thierry Callegari, Laurent Liard, and Olivier Pascal

ISAE-Supaéro LAPLACE, CNRS, UPS, INP Université dé Toulouse

Université dé Toulouse 31055 Toulouse, France 31400 Toulouse, France Email: vincent.laquerbe@isae.fr Abstract-In this paper, a tunable sub-wavelength resonator using a hemispherical plasma discharge is investigated at VHF-UHF frequencies. The overall setup is detailed. The resonance is experimentally investigated with an indirect method which consists in studying the influence of the plasma hemisphere on the transmission between two near-field monopole probes. The design presented here opens new perspectives for the realisation of plasma-based electrically small antennas thanks to a nonintrusive and adequate integration of the whole plasma control system below a ground plane.

Index Terms-Electrostatic resonator, Plasma devices.

I. In t r o d u c t i o n

S UB-WAVELENGTH particles are known to exhibit very peculiar resonant phenomena at optical frequencies [START_REF] Bohren | Absorption and Scattering o f Light by Small Particles[END_REF]. For example, the quasi-static analysis of an electrically small sphere of dielectric constant e imersed in vacuum shows that its polarisability a, which represents its response to an incident static electric field E 0, is given by [START_REF] Sihvola | Character of surface plasmons in layered spherical structures[END_REF] :

6 e -1 a = 4 na6 (1) e + 2 w
where a is the radius of the considered sphere. When e = -2 , the polarisability a becomes theoretically infinite and reveals the electrostatic dipole resonance, also known as the localized surface plasmon resonance (LSPR). At microwave frequencies, it has been suggested to use this property for the realisation of efficient electrically small antennas by coupling such an ENG (for Epsilon NeGative) resonator to a feeding dipole or monopole-like probe [START_REF] Ziolkowski | Metamaterial-based efficient electrically small antennas[END_REF], [START_REF] Stuart | Electrically small antenna elements using negative permittivity resonators[END_REF]. Practical implementations are based on metamaterials that exhibit negative permittivity in the microwave range [START_REF] Smith | Loop-wire medium for investigating plasmons at microwave frequencies[END_REF], [START_REF] Stuart | Subwavelength microwave resonators exhibiting strong coupling to radiation modes[END_REF]. However, it was also suggested in [START_REF] Ziolkowski | Metamaterial-based efficient electrically small antennas[END_REF], [START_REF] Stuart | Electrically small antenna elements using negative permittivity resonators[END_REF] to use plasma discharges to obtain ENG properties at microwave frequencies.

In this paper, a plasma-based electrostatic resonator is presented. Its experimental characterisation is proposed by measuring its effect on the transmission between two monopole probes. Radiation enhancement and frequency tuning are experimentally demonstrated.

(thereafter plasma) is ignited inside a hemispherical bell jar. The plasma relative permittivity follows the Drude model [START_REF] Lieberman | Principles o f Plasma Discharges and Materials Processing[END_REF]:

£p (w ) = 1 - W W (w + iv ) ' with wp = 2 e2 n e £o me ( 2
)
where v, u p, e, n e, m e, and e0 are the electron collision frequency, the plasma angular frequency, the Coulomb charge, the electron density, the electron mass, and the vacuum permittivity, respectively. For a given operating angular frequency u , it is clear from (2) that within a dense plasma (i.e., when u p > u ), the dielectric constant becomes negative and a socalled ENG medium is created. In our study, as shown in Fig. 1, the spherical resonator is replaced by a hemisphere above a ground plane. The jar is composed of a pyrex-like glass whose dielectric constant is equal to 5.5. Its outer diameter 2a is 215 mm (A/4.65 at 300 MHz) and it has a thickness of 5.5 mm. The inner gas composition and its pressure are controlled by a vacuumsealed system, comprising a two stage pump and a gas bottle with controllable flow rate. Ignition power is then brought to the discharge through a Faraday shield thanks to a 4-turns inductive coil. The coil, placed below the Faraday shield and outside the sealed volume, is connected to a matching system and driven by a RF current at 13.56 MHz. The RF matching The experimental setup is presented in Fig. 1. A low pressure non-equilibrium non-magnetized plasma discharge system is composed of two capacitors, mounted in parallel, that can be tuned to maximize the power transfer between the coil and the Inductively Coupled Plasma (ICP) discharge [START_REF] Chen | Capacitor tuning circuits for inductive loads[END_REF]. The Faraday shield has been designed numerically with Ansys HFSS and allows an efficient inductive energy coupling between the coil and the discharge by suppressing the capacitive part, which generates more losses in the discharge [START_REF] Mahoney | Electron-density and energy distributions in a planar inductively coupled discharge[END_REF].

The main advantage of the proposed design is that the vacuum generation system and the plasma ignition source are placed below the ground plane. Therefore, no other element is protruding above the ground plane except the resonator itself.

As can be seen from Fig. 2a, the resulting ICP discharges fills the entire hemispherical volume. It is also clear from Fig. 2a that the plasma is inhomogeneous, which is in fact the case for any discharge. Although homogeneous ENG structure has been considered in [START_REF] Ziolkowski | Metamaterial-based efficient electrically small antennas[END_REF], [START_REF] Stuart | Electrically small antenna elements using negative permittivity resonators[END_REF], it has been recently shown analytically in [START_REF] Laquerbe | Analytical model to study the electrostatic resonance of subwavelength radially inhomogeneous negative permittivity spheres[END_REF] that a radially inhomogeneous spherical medium with a continuous profile might also excite the desired electrostatic resonance. In this case, the electron density at the center of the structure must be larger to balance its fall at outer interfaces.

III. Ex p e r i m e n t a l St u d y

The resonance of the plasma hemisphere is experimentally investigated with two identical monopole probes, placed symmetrically at a distance of 40 mm of the bell jar (see Fig. 2b). In order to minimize side effects, the ground plane has been extended beyond the probes and forms a metallic disk with an outer diameter of approximately 2 m. The probes radius is 1 mm and their length is 4 cm. They are soldered to N-type flange connectors that are directly connected to a Vector Network Analyser. The whole setup without plasma has been simulated with Ansys HFSS. The transmission coefficient between the two probes is presented in Fig. 3 and is in good agreement with measurements. At these frequencies, the transmission is low because the probes are inefficient, since their nominal operating band is rather expected around 1.875 GHz.

When a plasma discharge is ignited with the appropriate electron density, the transmission coefficient is strongly affected in the VHF-uHF, where its permittivity is negative. To quantify the influence of the plasma and its transmission enhancement, it is more interesting to consider the relative | between the two probes when no plasma is ignited. transmission coefficient A S 21, simply calculated as the difference, in decibels, between the amplitude of the transmission coefficients with and without plasma, and defined by :

A S2 i S (p 0 21 dB S. dB (3)
In this study, an argon plasma is considered at the pressure P = 5 mTorr (approximately 0.67 Pa). Fig. 4 represents A S 2i as a function of the frequency for various RF input power Pin injected to the coil. For each value of P in, A S 2i becomes locally positive. This transmission enhancement can be explained by the electromagnetic coupling occuring between the two probes and the plasma discharge, which behaves as a resonator at these frequencies. At this pressure, the largest measured transmission enhancement reaches 5 dB and broadens over 50 MHz. The resonant frequency can now be tuned by simply adjusting the RF power. Indeed, as predicted by theoretical approaches, tuning the RF input power modifies the electron density inside the discharge. This modification of the electron density changes the dielectric constant, which in turn, tunes the resonant frequency. From Fig. 4, it can be seen that a large tuning capability is achieved through a RF power range of 20 W. The measured frequency agility varies roughly linearly with P in. At this pressure, it has been observed that the current device could not be operated for lower values of P in, since the discharge would simply extinguish.

In addition, the influence of the inner pressure is now investigated. The RF input power is set equal to Pin = 138 W and the pressure P varies. The corresponding relative transmission coefficients are shown in Fig. 5. In this case, the transmission between the two probes is again amplified and reflects the excitation of the electrostatic resonance. Variations of inner pressure also provide frequency agility of the resonance. Indeed, when P varies, the number of particles that can be ionized to create electron changes, which allows to control n e. However, this method might not be desired for frequency tuning because the pressure directly affects the The ENG resonance has been demonstrated using a two probes system that measures the influence of the discharge on the transmission coefficient. Despite the non-negligible electrical size (ka at least equal to 0.58) of the structure, transmission amplifications were measured up to 10 dB, which suggests very large amplification can be obtained for smaller ka. The main limitation of the current prototype remains its reliability in the VHF band. So far, no resonance was observed below 260 MHz as no stable discharge could be sustained at both low input power and low inner pressure. Further work is needed to overcome this limitation and create highly miniature resonator at VHF frequencies and beyond.

o IV t v 'M 1 » . f •' ,i! . j ! 1 ' i .i > y ; : ÿjjJh.i , ! ' i f U U- it! : ! I, jil ... i V .'¿ V Y ; / V .. ..-v.-v A - ■ A • : ■ V / V i / v
This original implementation of the LSPR using plasma discharge can be further extended to other frequency band as long as the plasma electron density follows accordingly, as stated by ( 2) and the value of the ratio wp/w. By coupling a simple probe to this resonator, one may create a miniature ENG structure that radiates efficiently, as initially suggested in [START_REF] Ziolkowski | Metamaterial-based efficient electrically small antennas[END_REF], [START_REF] Stuart | Electrically small antenna elements using negative permittivity resonators[END_REF]. The prototype developed in this paper might be well suited for this application, thanks to its placement above a ground plane. A hermetically-sealed coaxial probe, placed at the center of the sphere and terminated by a monopole, could bring the microwave power to be radiated. The symmetry of the structure may ease its mechanical integration and isolate the RF power (used to sustain the discharge) from the signal of interest feeding the monopole itself. Such an electrically small antenna will exhibit interesting electrical properties like furtivity (when no discharge is ignited) and frequency agility (with the plasma on).

V. ACKNOW LEDGEM ENT

The authors would like to thank the French Defence Agency (DGA) for its financial support. electron collision frequency v [START_REF] Lieberman | Principles o f Plasma Discharges and Materials Processing[END_REF], and so the ohmic losses in the discharge. Thus, larger values of P increase losses. As long as the discharge can be sustained for relatively low value of input power Pin, an inner pressure P as low as possible shall be preferred. Similarly, the pressure plays a significant role in the amplitude of the resonance. Fig. 5 indeed highlights a stronger transmission enhancement for the lowest values of pressure (up to 10 dB enhancement at P = 4 mTorr), due to lower ohmic losses in the discharge. The same extinction issue has been observed as before and the inner pressure could not be lowered as the discharge would also shut down.

IV. Co n c l u s i o n a n d D i s c u s s i o n s In this paper, a novel hemispherical electrostatic resonator has been developed in the VHF-uHF band using a plasma discharge. The plasma here acts as an ENG medium and its parameters (RF input power P in and discharge pressure P ) can be adjusted to achieve frequency tuning.

Re f e r e n c e s

  (13,56 M H z) II. De v i c e a n d e x p e r im e n t a l s e t u p

Fig. 1 .

 1 Fig.1. Schematic of the experimental setup comprising gas management and RF plasma ignition systems. Except for gas inlet and outlet, only half of the structure is represented for clarity.

■40Fig. 2 .

 2 Fig. 2. Photography of (a) the inductively coupled plasma discharge in the vacuum-sealed hemisphere and (b) the two monopole probes placed on an extended ground plane.

FrequencyFig. 3 .

 3 Fig. 3. Simulated and measured transmission coefficients || between the two probes when no plasma is ignited.

Fig. 4 .

 4 Fig. 4. Relative transmission coefficients A S 2i measured at P = 5 mTorr for various RF input power Pi n .

Fig. 5 .

 5 Fig. 5. Relative transmission coefficients A S 2i measured at Pin = 138 W for various gas pressure P .