Damien Depannemaecker
email: damien.d@cnrs.fr

Alain Destexhe

Viktor Jirsa

Christophe Bernard

Modeling seizures: from single neurons to networks

Dynamical systems tools offer a complementary approach to detailed biophysical seizure modeling, with high potential for clinical applications. This review describes the theoretical framework, allowing theorizing certain properties of seizures and their classifications according to their dynamics properties at onset and offset. We describe various modeling approaches spanning different scales, from single neurons to large-scale networks. We try to offer a large accessible overview through nonexhaustive examples of recent important works.

Introduction

Epilepsies are defined by the occurrence of spontaneous and recurrent seizures [START_REF] Fisher | ILAE Official Report: A practical clinical definition of epilepsy[END_REF]. Seizures are characterized by electrophysiological (field potential) signatures, which can be different from one patient to the next, and even within a given patient [START_REF] Saggio | A taxonomy of seizure dynamotypes[END_REF].

Given such variability, it is important to rely on a taxonomy of seizures, not only for mechanistic studies (different seizures may be associated to different underlying mechanisms, hence to different pharmacotherapy), but also for diagnostic purposes (in particular for presurgical evaluation in the case of drug-resistant epilepsies). One approach to characterize seizure types is to analyze their dynamical properties [START_REF] Saggio | A taxonomy of seizure dynamotypes[END_REF][START_REF] Jirsa | On the nature of seizure dynamics[END_REF][START_REF] Saggio | Fast-Slow Bursters in the Unfolding of a High Codimension Singularity and the Ultra-slow Transitions of Classes[END_REF].

For this, dynamical systems theory is essential [START_REF] Jirsa | On the nature of seizure dynamics[END_REF][START_REF] Stefanescu | Computational models of epilepsy[END_REF][START_REF] Soltesz | Computational neuroscience in epilepsy[END_REF] . In this journal, Stefanescu, and colleagues [START_REF] Stefanescu | Computational models of epilepsy[END_REF] show the importance and the basis of the methods used for computational modeling of seizures. Building on this work, we propose here to highlight some recent works in this field. Numerous approaches emerged, from phenomenological aspects of dynamics to detailed biophysical models, from the subcellular scale to the large scale of whole brain regions [START_REF] Jirsa | The Virtual Epileptic Patient: Individualized whole-brain models of epilepsy spread[END_REF][START_REF] Depannemaecker | A unified physiological framework of transitions between seizures, status epilepticus and depolarization block at the single neuron level[END_REF][START_REF] Toglia | Mitochondrial dysfunction and role in spreading depolarization and seizure[END_REF][START_REF] Wei | Unification of neuronal spikes, seizures, and spreading depression[END_REF][START_REF] Wei | Oxygen and seizure dynamics: II. Computational modeling[END_REF][START_REF] Depannemaecker | Computational study of responses of spiking networks to incoming seizure-like perturbation[END_REF]. We discuss the interest of each of these different approaches to show how all contribute to a better understanding of the mechanisms underlying seizures.

Computational modeling of seizures.

Brain activity constantly changes from one moment to the next. A key factor (among many) determining the temporal evolution of a neuron's activity is the excitability of the neuronal membrane. Thus, the vast majority of models use this variable, from the cell level, to the behavior of neural networks, to brain's global activity.

At any level of analysis, the type of study that can be carried out from a dynamical system point of view depends upon the number of dimensions of the model (i.e., number of variables). In order to maintain a manageable dimensionality, a model often needs to be reduced before moving to the next level of analysis, as shown schematically in Fig 1. The core of a model is what constitutes the sources of the electrophysiological activity (the readout of the model; here seizures). Some models include other variables that interact with electrophysiological activity, such as metabolic activity, variations in ion or oxygen concentration, or the interaction with glial cells.

Another approach consists of modeling the signal itself (i.e. the readout), which can be recorded by different methods such as the electroencephalography (EEG), magnetoencephalography (MEG), and functional magnetic resonance imagery (fMRI).

Since, in these cases, it is the measurement that is modeled and not the source, it becomes difficult to obtain insight in terms of biological mechanisms. The mapping between source and sensor signals is one the most problematic factors, when interpreting human brain imaging signals. The Virtual Brain (TVB) simulator closes this gap to some degree [START_REF] Ritter | The Virtual Brain Integrates Computational Modeling and Multimodal Neuroimaging[END_REF][START_REF] Sip | Computational modeling of seizure spread on a cortical surface[END_REF], even though principled model non-identifiability cannot be overcome by definition. However, this approach has the advantage to remain as close as possible to the clinical data.

It is however possible to establish a link between the source and the signal [START_REF] Cosandier-Rimélé | Computational modeling of epileptic activity: From cortical sources to EEG signals[END_REF], which requires keeping track of the temporal evolution of certain variables. This is the reason why a large majority of models used in computational neuroscience are based on systems of differential equations. In a first-order differential system, the state of a variable at a time t is calculated based on its variation with respect to time t-∆t. The rate of amplitude changes over time will determine the time scale. Each equation composing the system can vary according to its own time scale. This is an important consideration in epilepsy because spontaneous seizures occur on a much slower scale (hour, day, month) than that of the electrophysiological signature of a seizure itself (seconds). Thus, to model epilepsy, a so-called rapid system is needed, which can switch from the "normal" state to a seizure state, as well as a slower system driving the transition between these states. It is possible to construct a slow-fast system [START_REF] Bertram | Multi-timescale systems and fast-slow analysis[END_REF][START_REF] Izhikevich | Dynamical Systems In Neuroscience : The Geometry of Excitability and Bursting[END_REF], which includes both time scales, and in which an external input can be used to drive the transition between states. Such models can produce the different states, or in other words, dynamical regimes. These different dynamical regimes can correspond to different attractors. In small dimension systems, these attractors become simple, such as fixed-points or limit-cycles. This is one of the advantages provided by the reduced dimensional phenomenological models described in the next section.

Phenomenological models of seizure

The phenomenological approach, reducing the models to a minimum number of variables and parameters, allows an exhaustive study of the dynamics. The Epileptor is one such model [START_REF] Jirsa | On the nature of seizure dynamics[END_REF]. It was constructed to reproduce the most common electrophysiological signature of seizures found in patients and experimental models [START_REF] Jirsa | On the nature of seizure dynamics[END_REF][START_REF] Houssaini | Seizures, refractory status epilepticus, and depolarization block as endogenous brain activities[END_REF][START_REF] Houssaini K El | The epileptor model: A systematic mathematical analysis linked to the dynamics of seizures, refractory status epilepticus, and depolarization block[END_REF]. Interestingly, the Epileptor model also includes in its dynamical regime status epilepticus and depolarization block [START_REF] Houssaini | Seizures, refractory status epilepticus, and depolarization block as endogenous brain activities[END_REF][START_REF] Houssaini K El | The epileptor model: A systematic mathematical analysis linked to the dynamics of seizures, refractory status epilepticus, and depolarization block[END_REF]. Using the same approach and a minimal generic model, it is possible to generate sixteen possible seizure types in terms of dynamics, or dynamotypes [START_REF] Saggio | A taxonomy of seizure dynamotypes[END_REF][START_REF] Saggio | Fast-Slow Bursters in the Unfolding of a High Codimension Singularity and the Ultra-slow Transitions of Classes[END_REF]. Each of these classes differs from the others according to the type of transition between the healthy state and the seizure state. In terms of dynamical systems theory, this corresponds to the types of bifurcation at the onset and the offset of the seizure [START_REF] Jirsa | On the nature of seizure dynamics[END_REF][START_REF] Houssaini K El | The epileptor model: A systematic mathematical analysis linked to the dynamics of seizures, refractory status epilepticus, and depolarization block[END_REF]. Each bifurcation can be displayed on a map of brain dynamics [START_REF] Saggio | A taxonomy of seizure dynamotypes[END_REF][START_REF] Saggio | Fast-Slow Bursters in the Unfolding of a High Codimension Singularity and the Ultra-slow Transitions of Classes[END_REF]. Thus, this theoretical approach not only leads to a taxonomy of clinically observed seizures [START_REF] Saggio | A taxonomy of seizure dynamotypes[END_REF], but also offers a theoretical framework to explain the coexistence of healthy and pathological regimes, and how it is possible to navigate the map brain dynamics, explaining why individual patients can display several dynamotypes [START_REF] Houssaini | Seizures, refractory status epilepticus, and depolarization block as endogenous brain activities[END_REF]. This type of model, called phenomenological, precisely presents the interest of describing only the observed phenomenon, independently of the associated biophysical underpinnings. Although it does not inform directly on the underlying biological mechanisms, it provides a theoretical framework and predictions in terms of dynamics, to which more detailed biophysical models can refer, from the single neuron to the large-scale brain model scales. In the next subsections, we will detail some approaches of interest for modeling electrophysiological seizures, across scales from the lowest (i.e., single neuron), to the highest (i.e., whole brain).

Models at the single neuron level

Single neuron-scale models are very diverse and can include very different levels of detail. Some models are highly detailed, taking into account sub-cellular organelles such as mitochondria [START_REF] Toglia | Mitochondrial dysfunction and role in spreading depolarization and seizure[END_REF] or the detailed morphology of the dendritic and axonal tree [START_REF] Santhakumar | Role of mossy fiber sprouting and mossy cell loss in hyperexcitability: A network model of the dentate gyrus incorporating cell types and axonal topography[END_REF][START_REF] Tejada | Combined Role of Seizure-Induced Dendritic Morphology Alterations and Spine Loss in Newborn Granule Cells with Mossy Fiber Sprouting on the Hyperexcitability of a Computer Model of the Dentate Gyrus[END_REF]. In the latter case, it is necessary to divide the geometry of the neuron into different compartments. As a function of the objectives of the model, it is possible to reduce the general dynamics to that of a single compartment. We then speak of a point neuron model (see Fig. 1).

Many approaches are possible for this type of modeling [START_REF] Izhikevich | Which Model to Use for Cortical Spiking Neurons[END_REF]. We will see how these approaches can be extended to model of the seizures in particular. A first approach was that of Cressman and colleagues [START_REF] Cressman | The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics[END_REF], who proposed a model at the single neuron scale. It is based on the formalism developed by Hodgkin and Huxley [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. It can be made more complex by adding the equations describing the variations in ionic concentrations and the interaction with glial cells, which dynamically play the role of slow variables. This model has the advantage of being detailed from a biophysical point of view. It presents a large number of differential equations limiting the possibility to fully study its dynamics, which requires using a reduced model [START_REF] Cressman | The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics[END_REF].

The reduced model reproduces only a few bifurcation types at the onset and offset of a seizure, but not the most commonly observed in clinical and experimental data (SN/homoclinic). Such an approach has been extended [START_REF] Cressman | The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics[END_REF][START_REF] Ullah | Assimilating seizure dynamics[END_REF], resulting in a unified framework for spike, seizure, and spreading depression [START_REF] Wei | Unification of neuronal spikes, seizures, and spreading depression[END_REF]. In this framework, it is possible to study specific biophysical interactions; for example the relationship between oxygenation and electrophysiological activities [START_REF] Wei | Oxygen and seizure dynamics: II. Computational modeling[END_REF]. At this scale, the link can be made with in vitro experiments [START_REF] Ingram | Oxygen and seizure dynamics: I. experiments[END_REF].

Following the same general strategy, a different reduced model has been proposed [START_REF] Depannemaecker | A unified physiological framework of transitions between seizures, status epilepticus and depolarization block at the single neuron level[END_REF]. It offers a unified framework for spiking, burst, seizure, status epilepticus, and depolarization block. It describes the excitability of the membrane and the interaction with the environment (bath/glial cell) through changes in ionic concentrations with only four differential equations. Compared with the previous one, this model is less detailed biophysically, but as it is dimensionally reduced, it allows a direct correspondence to phenomenological models. The underlying dynamics corresponds to what is observed in the signals recorded experimentally.

A third approach has been proposed by Chizhov and colleagues [START_REF] Chizhov | Minimal model of interictal and ictal discharges "Epileptor-2[END_REF], which also takes into account the dynamic constraints described by phenomenological models.

The model can reproduce the different dynamics of seizures in high potassium conditions. The electrophysiological activities described in this model are not only associated with variations in ionic concentrations but also with synaptic activity.

Although it is at the scale of a single neuron, the model considers the effects of the network through synaptic inputs. The network models are discussed in the next subsection.

Models at the network level

Network models describe many cells interacting with each other. If the synaptic activity is often preponderant in these models, some show that all the other interactions (ephaptic, ionic diffusion, etc.) can be sufficient for the genesis and the propagation of seizures [START_REF] De Almeida | Mechanistic hypotheses for nonsynaptic epileptiform activity induction and its transition from the interictal to ictal state-Computational simulation[END_REF][START_REF] Rodrigues | pH during non-synaptic epileptiform activity-computational simulations[END_REF]. Ephaptic interactions may for instance synchronize atin potential propagation in tightly packed fiber bundles [START_REF] Sheheitli | A mathematical model of ephaptic interactions in neuronal fiber pathways: Could there be more than transmission along the tracts?[END_REF]. Slow oscillations can emerge from non-synaptic interactions, enabling the transitions towards seizure-like states [START_REF] De Almeida | Mechanistic hypotheses for nonsynaptic epileptiform activity induction and its transition from the interictal to ictal state-Computational simulation[END_REF][START_REF] Rodrigues | pH during non-synaptic epileptiform activity-computational simulations[END_REF][START_REF] Bikson | Synaptic transmission modulates while non-synaptic processes govern the transition from pre-ictal to seizure activity in vitro[END_REF]. But most studies done at the network scale take into account the interconnection between neurons through models of synapses [START_REF] Naze | Computational Modeling of Seizure Dynamics Using Coupled Neuronal Networks: Factors Shaping Epileptiform Activity[END_REF][START_REF] Liou | A model for focal seizure onset, propagation, evolution, and progression[END_REF][START_REF] Kim | The influence of depolarization block on seizure-like activity in networks of excitatory and inhibitory neurons[END_REF][START_REF] Jacob | A proposed mechanism for spontaneous transitions between interictal and ictal activity[END_REF][START_REF] Destexhe | Spike-and-wave oscillations based on the properties of GABA(B) receptors[END_REF][START_REF] Destexhe | Can GABA A conductances explain the fast oscillation frequency of absence seizures in rodents?[END_REF][START_REF] Destexhe | Corticothalamic Feedback: A Key to Explain Absence Seizures[END_REF]. In particular, such models can account for the genesis of generalized seizures with spike-and-wave EEG patterns, with mechanisms that critically depend on the balance between GABAergic and glutamate conductances [START_REF] Destexhe | Spike-and-wave oscillations based on the properties of GABA(B) receptors[END_REF][START_REF] Destexhe | Can GABA A conductances explain the fast oscillation frequency of absence seizures in rodents?[END_REF], as well as cortico-thalamic interactions [START_REF] Destexhe | Spike-and-wave oscillations based on the properties of GABA(B) receptors[END_REF].

However, an essential question arises at this scale. Network models are constructed from individual models of interconnected neurons. It is not obvious that because we assemble these bricks that the result of the network obtained is a relevant model of biological networks. One solution is to keep only the emerging behaviors of the network, which displayed dynamics comparable to that observed experimentally. Still, it does not imply that the causes of these emergences are the same. In the case of seizures, transitions to the ictal state do not take place uniformly in the network.

Modeling at this scale makes it possible to observe the collective effects and the propagation of the seizure [START_REF] Depannemaecker | Computational study of responses of spiking networks to incoming seizure-like perturbation[END_REF][START_REF] Liou | A model for focal seizure onset, propagation, evolution, and progression[END_REF][START_REF] Ho | Interaction between synaptic inhibition and glialpotassium dynamics leads to diverse seizure transition modes in biophysical models of human focal seizures[END_REF]. The reverse approach, imposing the phenomenological constraint from the upper scale to build a network with a biophysical description, was implemented by Naze and colleagues [START_REF] Naze | Computational Modeling of Seizure Dynamics Using Coupled Neuronal Networks: Factors Shaping Epileptiform Activity[END_REF]. Even if the number of biophysical interactions described remains limited, they show how dynamical system theory makes it possible to build a biophysical model on a smaller scale from a framework derived from macroscopic observations. The construction of networks involves large-dimensional models, which makes it impossible to apply directly the dynamics analysis used in small-dimensional systems. To overcome this problem, one can consider global measures or reduce the description of the network's activity to a reduced number of variables. Considering the network as a whole, aim to study the activity on larger scales, which is discussed in the next section.

From mean-field to large scale models

Mean-Field models describe the activity of a large group of neurons using average values of their activity, such as the mean membrane potential or, more often, the firing rate of the population. Some studies have derived these descriptors of the activity from single neuron models [START_REF] Carlu | A mean-field approach to the dynamics of networks of complex neurons, from nonlinear Integrate-and-Fire to Hodgkin-Huxley models[END_REF][START_REF] Montbrió | Macroscopic Description for Networks of Spiking Neurons[END_REF]. Another method is to consider that a neuron can represent the mean activity of a population of neurons. It can be used in neuralfield (including space representation) or neural mass model to describe seizures [START_REF] Liou | A model for focal seizure onset, propagation, evolution, and progression[END_REF][START_REF] Buchin | Adaptation and inhibition control pathological synchronization in a model of focal epileptic seizure[END_REF][START_REF] Kim | The influence of depolarization block on seizure-like activity in networks of excitatory and inhibitory neurons[END_REF]. One of the historical models at this scale has been developed by Wilson and Cowan [START_REF] Wilson | Excitatory and Inhibitory Interactions in Localized Populations of Model Neurons[END_REF] and has been used to study focal epilepsy through bifurcation analysis [START_REF] Kim | The influence of depolarization block on seizure-like activity in networks of excitatory and inhibitory neurons[END_REF][START_REF] Meijer | Modeling focal epileptic activity in the Wilson-cowan model with depolarization block[END_REF]. These approaches permit a comprehension of the phenomena observable at a global activity scale [START_REF] Kim | The influence of depolarization block on seizure-like activity in networks of excitatory and inhibitory neurons[END_REF][START_REF] Takeshita | Transitions between multistable states as a model of epileptic seizure dynamics[END_REF]. The reduction of networks to mean-field models necessarily leads to the loss of information on certain internal dynamics and collective effects, to describe only certain global aspects. This may be sufficient to characterize seizures. Reducing the number of variables facilitates a bifurcation analysis. This reduction allows the comparison with phenomenological models.

At such scale of analysis, results are relevant to recordings obtained in humans, using scalp or intracerebral EEGs. In particular, it is possible to consider and study interactions between regions, like thalomo-cortical interactions [START_REF] Marten | Derivation and analysis of an ordinary differential equation mean-field model for studying clinically recorded epilepsy dynamics[END_REF]. It also allows studying particular properties of networks in a phenomenological way [START_REF] Benjamin | A phenomenological model of seizure initiation suggests network structure may explain seizure frequency in idiopathic generalised epilepsy[END_REF][START_REF] Terry | Seizure generation: The role of nodes and networks[END_REF]. At this scale, the temporal dissection exists, and the variation of a parameter or slow variables enables the transitions to the ictal state and its return to a "healthy" regime [START_REF] Marten | Derivation and analysis of an ordinary differential equation mean-field model for studying clinically recorded epilepsy dynamics[END_REF][START_REF] Wendling | Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition[END_REF]. As such models are at the scale of a large region, it is possible to connect them to create a large scale network describing the whole brain [START_REF] Melozzi | The virtual mouse brain: A computational neuroinformatics platform to study whole mouse brain dynamics[END_REF][START_REF] Goldman | Brain-scale emergence of slow-wave synchrony and highly responsive asynchronous states based on biologically realistic population models simulated in The Virtual Brain[END_REF].

A whole brain approach using the phenomenological model described previously has been developed [START_REF] Olmi | Controlling seizure propagation in large-scale brain networks[END_REF]. At this level, such models can have important clinical applications. One example is the prediction of the size of the epileptogenic zone in pharmaco-resistant epilepsies [START_REF] Junges | The role that choice of model plays in predictions for epilepsy surgery[END_REF][START_REF] Kini | Virtual resection predicts surgical outcome for drug-resistant epilepsy[END_REF][START_REF] An | Optimization of surgical intervention outside the epileptogenic zone in the virtual epileptic patient (VEP)[END_REF][START_REF] Jirsa | The Virtual Epileptic Patient: Individualized whole-brain models of epilepsy spread[END_REF][START_REF] Lopes | An optimal strategy for epilepsy surgery: Disruption of the rich-club?[END_REF]. Another possible outcome of such modeling approach is the control of seizure through stimulation [START_REF] Ching | Distributed control in a mean-field cortical network model: Implications for seizure suppression[END_REF][START_REF] Ashourvan | Model-based design for seizure control by stimulation[END_REF]. These computational neuroscience methods have developed very rapidly in recent years; they hold great promises in terms of translation to the clinic.

Discussion

Computational models applied to epilepsy already have a long history [START_REF] Soltesz | Computational neuroscience in epilepsy[END_REF].

Recent progress in seizure modeling has brought a new understanding of this complex phenomenon. At the level of biophysical interactions, such as variations in ionic concentrations [START_REF] Depannemaecker | A unified physiological framework of transitions between seizures, status epilepticus and depolarization block at the single neuron level[END_REF][START_REF] Wei | Unification of neuronal spikes, seizures, and spreading depression[END_REF][START_REF] Cressman | The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics[END_REF], it is possible to create direct correspondences with experimentally measurable physical quantities. It is possible to make predictions design intervention strategies targeting known biophysical pathways. In terms of collective network dynamics, computational modeling shows how the complexity leads to the emergence of a rich repertoire of behaviors, including seizures [START_REF] Kim | The influence of depolarization block on seizure-like activity in networks of excitatory and inhibitory neurons[END_REF]. At this mesoscopic scale, reduced models have been proposed, bringing a description of the relation between electrophysiological activities and recorded signals such as EEG [START_REF] Blenkinsop | The dynamic evolution of focal-onset epilepsies -combining theoretical and clinical observations[END_REF].

Networks of mesoscopic models permit the study of the seizure propagations, building the bridge with macroscopic scales and clinical observations. Thus, the phenomena associated with seizures can be described at different scales. A global understanding of this phenomenon requires this pluralist approach in terms of scale. The complexity of the brain does not allow us to establish a simple relationship between the different scales. However, it seems that fundamental invariant properties are needed such as at least two distinct time scales, one representing the rapid electrophysiological activity or the variations of measured signals, and a second, slower, at the origin of the transitions between the ictal and non-ictal states, which can be linked to a very wide variety of phenomena depending on the scales. This temporal dissection is included in the phenomenological model, which also opens the way for classification according to the characteristics of the onset and the offset of the seizures through the different types of bifurcations found in the biophysical models at different scales. This is consistent with clinical observations where the onset is an important criterion [START_REF] Lagarde | The repertoire of seizure onset patterns in human focal epilepsies: Determinants and prognostic values[END_REF]. This offers a direct link between theoretical computational studies and clinical practice through new approaches that have been born in recent years, such as the virtual epileptic patient [START_REF] Jirsa | The Virtual Epileptic Patient: Individualized whole-brain models of epilepsy spread[END_REF], and tends to contribute to develop seizure forecasting methods and presurgical evaluation [START_REF] Kuhlmann | Role of multiple-scale modeling of epilepsy in seizure forecasting[END_REF].

Figure 1 .

 1 Figure 1. Simplistic diagram of modeling at different scales. Details models at a given scale have a high number of dimensions (i.e., number of variables) and can be reduced. From the reduced model, it is possible to build a correspondence to the next larger scales by building a network.


domain: spi phys

text:
Nurdan Anılmış Baci
email: nurdan.baci@tubitak.gov.tr

Aydın Yeniay

Wave-front Aberrations of Laser Beams Propagating through the Atmosphere

Keywords: atmospheric turbulence, wave-front aberrations, laser beam propagation

Small temperature variations carried by the turbulent velocity field in the Earth's atmosphere produce small phase perturbations in laser beam propagating through it. Perturbations in the phase cause wave-front errors and may significantly decrease the Strehl ratio of the laser beam. For compensation of atmospheric distortions, adaptive optics systems are used. Before the design of any adaptive optics system, the simulation of the atmospherically distorted laser beam is important. In this work, we describe the purpose, theory, implementation, and results of a wave optics propagation simulation code LAtmoSim and we investigate the effect of laser system parameters on the degree of the wave-front error caused by turbulence.

INTRODUCTION

Laser beams that are used in directed energy and free-space optical communication applications can become severely distorted by optical turbulence when propagating through the atmosphere. Those distortions mostly originate in the atmospheric boundary layer, which typically extends up to a height of about 2 kilometers and contains 75% percent of the atmosphere's mass. Effects are beam wander, beam broadening and irradiance fluctuations. The use of adaptive optics techniques is possible to mitigate the effects of turbulence and such systems are well known for astronomical applications but implementation and performance in directed energy applications are still not very well known. Any AO system's goal is to eliminate the distortions in the wave-front of the light that is the result of optical path variations through measuring the errors using a wave-front sensor, calculating an appropriate correction, and applying this correction to a deformable mirror. This feedback loop is carried out several hundred times a second to comply with the temporal bandwidth requirement. To determine the characteristics of an AO system, the effects of the atmospheric turbulence on the wave-front has to be simulated. The estimation of the performance of laser systems depends on the accuracy of the atmospheric models being used in the propagation prediction codes. After decades of research, some analytical theories such as geometrical optics [START_REF] Clarke | Analysis of Laser Beam Propagation in a Turbulent Atmosphere[END_REF] , Rytov Method and Markov approximation [START_REF] Klyatskin | On the Displacement of Spatially-Bounded Light Beams in a Turbulent Medium in the Markovian-Random-Process Approximation[END_REF][START_REF] Kon | Dispersion of Light Beam Displacements in the Atmosphere with Strong Intensity Fluctuations[END_REF][START_REF] Mironov | On the Theory of Spatially Limited Light Beam Displacements in a Randomly Inhomogeneous Medium[END_REF] have been developed to calculate the characteristics of laser beam propagation. But these methods are approximation in certain conditions, therefore their applicability is limited and theoretical calculations of the statistics of scintillation are difficult especially when the intensity fluctuations become large. Therefore, numerical methods are developed to get a more realistic representation of atmospheric turbulence effects on laser beam propagation. These methods are known as beam propagation method [START_REF] Van Roey | Beam-Propagation Method: Analysis and Assessment[END_REF] . Other names for the methods are the split-step Fourier technique [START_REF] Tappert | The Parabolic Approximation Method[END_REF] and the random phase-screen method [START_REF] Fejer | The Diffraction of Waves in Passing through an Irregular Refracting Medium[END_REF][START_REF] Lee | Weak Scattering in Random Media, With Applications to Remote Probing[END_REF] . Here we introduce laser beam propagation code LAtmoSim which enables us to assess the effect of the atmosphere on the laser beam wave-front and to determine the design parameters of the AO system by using methods described above. In this paper, we also present the results from the work done to predict the intensity of the atmospheric turbulence. A quantitative measure of the intensity of optical turbulence is called the refractive index structure parameter, C n 2 . The measurement of C n 2 is difficult; it requires expensive and bulky detectors and optical imaging systems. Estimating boundary layer C n 2 from easily available meteorological parameters, such as temperature and humidity makes to predict the effectiveness of a laser system before its use easy. The remainder of the paper is organized as follows. Section 2 introduces the capabilities of the simulation software LAtmoSim. Section 3 describes the basic mathematical models used for simulation operations such as the phase screen model. Section 4 presents the results of the work done to predict the refractive index structure parameter. In Section 5, models of wave-front aberrations are explained. In Section 6, we conclude with the prospect of how this work might be extended.

LATMOSIM: SIMULATIONS OF LASER BEAM PROPAGATION THROUGH THE ATMOSPHERE

LAtmoSim is developed to do numerical calculations on laser beam propagation through the atmosphere. The software is written in Python and Tkinter is used to develop the user interface of the software. The code takes into account atmospheric effects such as absorption, scattering, turbulence and thermal blooming. Including effects of aerosol and molecular absorption and scattering of radiation is achieved by providing a user interface for parameters in FASCODE.

As an alternative to FASCODE, LAtmoSim calculates atmospheric transmittance by utilizing data fitted analytic models of visibility and rain [START_REF] Marshall | The Distribution of Raindrops with Size[END_REF][START_REF] Kim | Comparison of Laser Beam Propagation at 785 nm and 1550 nm in Fog and Haze for Optical Wireless Communications[END_REF] . For slant path engagement scenarios such as ground to air and air to ground, it uses the Hufnagel-Valley model for vertical C n 2 profiles. LAtmoSim can create multiple numbers of focused or collimated Gaussian, super-Gaussian and Bessel Gaussian laser beams with properties including wavelength, beam size, truncation ratio, beam quality M², total power, and super-Gaussian order and Bessel coefficient. Atmospheric turbulence is modeled by generating random phase screens and filtering with the Kolmogorov spectrum. The loss of low spatial frequencies is handled with the sub-harmonic method. Those methods are described in detail in Section 3. The code calculates ensemble averages with a large number of realizations for the calculation of beam wander, long-term spot size, on-axis scintillation index SI; and compares those simulation results with values obtained from analytical calculations from Rytov theory. The code may be run both wave optics and geometric propagation mode to allow the code to be checked against linear analytical models. The other calculations include power-in-the-bucket, tracked-beam Strehl ratio, RMS wave front error, Peak to Valley (PtV) value of the phase distribution. In addition to atmospheric turbulence, a high-energy beam introduces an entirely new physical effect: Thermal blooming. Thermal blooming is the spreading of a laser beam when the beam heats the medium through which it is propagating. As the laser beam passes through the atmosphere, some of the beam's energy is absorbed, heating the atmosphere. The heated region will be less dense and will, consequently, have a lower index of refraction. Since the hottest regions will normally be in the center of the beam, a negative lens will develop in the atmosphere, and this lens will cause the divergence of the beam. LAtmoSim has both analytical [START_REF] Gebhardt | Self-Induced Thermal Distortion in the Near Field for a Laser Beam in a Moving Medium[END_REF] and numerical steady-state model 12 of thermal blooming. Figure 1 shows the results of two different atmospheric cases. In the first case, we take into consideration only the turbulence. In the second case, we added the effect of thermal blooming to turbulence. The simulation parameters are given in Table 1. The short term spots for each case are seen in Figure 1.

ATMOSPHERIC TURBULENCE MODELING

Turbulence is a complex phenomenon, and characterizing and modeling its effect on optical propagation poses a significant challenge. Both analytical [START_REF] Frisch | Wave propagation in random media[END_REF][START_REF] Ishimaru | Wave propagation and scattering in random media[END_REF][START_REF] Flatte | Sound transmission through a fluctuating ocean[END_REF][START_REF] Rytov | Principles of statistical radio physics: wave propagation through random media[END_REF] and numerical [START_REF] Frehlich | Simulation of laser propagation in a turbulent atmosphere[END_REF][START_REF] Martin | Intensity images and statistics from numerical simulation of wave propagation in 3-d random media[END_REF] methods have been developed. A promising method of numerical simulation uses a split-step phase screen approach. The fluctuations in refractive index are very small and their scale size is very large compared to the wavelength of the light, diffraction and scattering can be separated in numerical calculations. The propagation path is divided into many short intervals so that the propagating wave is hardly disturbed by each interval, but their cumulative effect can be significant. The starting point of numerical models is the parabolic approximation to the wave equation where the continuous random media is treated as a series of random phase screens transverse to the propagation direction of the electric field E. The field is propagated between the phase screens using the numerically efficient fast-Fourier-transform (FFT) algorithm. Random phase screens are generated with a spectral-domain algorithm also based on the FFT. Generation of random phase screens φ(x, y) is accomplished by 17

    ) / / ( 2 exp ) , ( ) , ( ) , ( 0 0 y x n m N lm N jn i m n ibz m n a y l x j x N y N             (1)
Where Δx and Δy are the grid spacing in the x and y coordinates, respectively; N x and N y are the array dimensions; L x = Δx N x and L y = Δy N y are the linear dimensions; and a (n, m) and b (n, m) are zero-mean Gaussian uncorrelated random numbers with ) , ( ) , ( ) , (

2 2 y x y x q m q n q q m n b m n a         (2)
Where <> denotes ensemble average, Φ φ (q x , q y ) is the spatial spectrum of the phase screens as a function of propagation distance z, Δ q x = 2π / (N x Δx) and Δq y = 2π / (N y Δy) are the grid size of the simulation in spectral space. The phase spectrum is

n z k     2 2  (3)
Where k = 2π / λ is the wavenumber of the field, Δz is the separation between phase screens in the propagation direction z, and Φ n is the spectrum of refractive-index fluctuations. The Kolmogorov model is the most frequently used for turbulence description because of its relative mathematical simplicity and good agreement with experiment [START_REF] Kolmogorov | Dissipation of energy in locally isotropic turbulence[END_REF] . Kolmogorov model of the index of refraction power spectrum is

3 / 11 2 33 . 0 ) (    q C q n n (4)
Where q is the magnitude of three-dimensional wave vector and C n 2 is the refractive-index structure constant. Because the phase screens produced by the FFT algorithm are periodic, the simulation of the fields has errors for those statistics that are sensitive to large-scale turbulent fluctuations. These errors will be removed by an improved simulation algorithm for the phase screens. The simulation algorithm does not include the power in the spectral region (-Δq x /2, Δq x /2) and (-Δq y /2, Δq y /2) which produces an error in the large-scale phase fluctuations. Lane et al. [START_REF] Lane | Simulation of a Kolmogorov phase screen[END_REF] proposed an algorithm based on the addition of sub-harmonics with random complex amplitudes to improve the large scale statistics i.e.,

   

)

) 3 /( 2 ) 3 /( 2 exp ) , , ( ) , , ( ) , ( 1 1 1
1 1 y p x p N p n m SH N ilm N ijn p m n ib p m n a y l x j p                  (5) ) , ( ) , , ( ) , , ( 2 2 yp xp yp xp q m q n q q p m n b p m n a         (6)
If n≠0 and m≠0. Here, Δq xp =Δq x /3 p and Δq yp =Δq y /3 p are the spectral domain grid sizes for sub-harmonic order p. The FFT sub-harmonic phase-screen simulation is no longer periodic and the discontinuities of the phase at the edges of the simulation may produce additional errors. These errors are negligible for beam-wave simulations if the beam is located at the center of the simulation. In Figure 2, phase screens with and without high-frequency components are shown.

In the following example, we explore the impact of the transmitting beam size on the far-field intensity calculated by LAtmoSim. For an atmospheric turbulence condition where the Fried radius is equal to 6 cm, we propagate several collimated beams with an ideal Gaussian shape. The beam diameters are chosen; 1cm < r 0 , 6 cm = r 0 , and 20cm > r 0 . The propagation distance used in this simulation is 2 km. Figure 3 shows the resulting intensity profiles for each beam. Beam with a diameter equal to the Fried parameter seems to be the least affected by turbulence.

C N 2 PREDICTION MODELS

We made simultaneous measurements of C n 2 and other atmospheric parameters that are known causing boundary layer turbulence. From the analysis of C n 2 data gathered from scintillometer measurements and meteorological data, we found that net radiation of the atmosphere, relative humidity, and wind to be the atmospheric parameters that are most correlated with C n 2 . Unfortunately, net radiation data is also not easily available. Therefore, we used Holtslag model [START_REF] Holtslag | A Simple Scheme for Daytime Estimates of the Surface Fluxes from Routine Weather Data[END_REF] which uses the time of the year, latitude, longitude, temperature, cloud cover and albedo as variables to calculate the net radiation. We compared net radiation data from pyrradiometer measurements with the Holtslag model for clear sky conditions and we found quite good agreement between the two (Figure 4).

Next, we used artificial neural networks to have C n 2 as a function of net radiation, relative humidity, and wind speed. Approximately six months of data from measurements of C n 2 , net radiation, relative humidity, and the wind is used as a training set. Later, a user-interface for the C n 2 estimator is developed which allows users to predict C n 2 from date, time, latitude, longitude, temperature, relative humidity, wind speed, cloud cover, albedo and height from the ground level. In Figure 5, a comparison of estimated and measured C n 2 values for the 24-hour duration is shown.

MODELING OF WAVEFRONT ABERRATIONS

The performance of an AO system can be evaluated using several different metrics. The most common performance metric for AO is the wave-front error variance. By minimizing the variance of wave-front error σ φ 2 , according to following Marechal approximation

) exp( 2     S (7)
It is possible to maximize Strehl ratio S. The use of Zernike polynomials for describing the aberrations of atmospheric turbulence is well known [START_REF] Noll | Zernike Polynomials and Atmospheric Turbulence[END_REF] . Beam wander is related to tip and tilt aberrations which are the second and third term of Zernike polynomials. Higher-order terms are the result of scintillation. Far-field phase distribution of a laser beam passing through atmospheric turbulence can be calculated by numerical propagation methods described in Section 3. If we let φ(x, y) be a turbulence-induced phase distribution, we can decompose this function over an orthogonal set of Zernike polynomials Z i (x, y) Where r 0 is the Fried parameter and D is the receiving aperture diameter. At first, LAtmoSim calculates Zernike coefficients a j by using the least square method. The accuracy of the least square method is shown in Figure 6. We also tested the numerical stability of Zernike coefficients against the changing the number of Zernike terms included in the calculations. Although Zernike polynomials form an orthogonal and complete set in the continuous domain, in discrete numerical calculations to have orthogonal and complete set one has to increase the number of sampling. Figure 7 shows that the stability of the Zernike tilt coefficient increases with an increased sampling number. We conclude that we need to increase the sampling size to 512 to have the correct Zernike coefficients. On the other hand, the Gram Schmidt orthogonalization scheme exists as a solution for the stability problem [START_REF] Upton | Gram-Schmidt orthogonalization of the Zernike polynomials on apertures of arbitrary shape[END_REF] .

After finding Zernike coefficients, the total wave-front error in the first J modes is deduced by using Noll's formula

     J j j J a 1 2 2  (10)
This allows us to determine how many modes are needed to be corrected to achieve Strehl ratios close to the diffractionlimited performance. In Figure 8, we see the result of the simulation run for a laser beam with a transmitter aperture size of 12 cm focused at a distance of L=1000 m passing through various atmospheric turbulence conditions. For the simulation, the physical size of the computational grid was set to 0.8 m for both transverse directions and the resolution of the grid was 512 x 512. This corresponds to approximately 1.5 mm pixel size. In the computational domain, pixel size determines the inner scale; Ɩ 0 of atmospheric turbulence. For the case where C n 2 is equal to 1.0e-15 m -2/3 , the Strehl ratio increases from 0.25 to 0.65 with the correction of only the first two modes; tip and tilt. The number of modes corrected is roughly proportional to the number of actuators a deformable mirror used in an AO system has. With the Strehl ratio calculations, we can determine the number of actuators needed for the compensation of atmospheric errors for different engagement scenarios for an ideal AO system. Another important design parameter for DM is stroke; dynamic range of the mirror. The stroke value depends on the PtV of the turbulent phase map. Conventional adaptive optics (AO) systems perform adequately when the phase of a distorted wave-front is continuous across the aperture, which occurs when the turbulence is weak. However, when an optical field is propagated through strong turbulence, it experiences nulls in its amplitude leading to branch points in the phase. It is shown by Fried [START_REF] Fried | Branch point problem in adaptive optics[END_REF] that when branch points are present in the phase of a turbulence-distorted optical field, the ability of an adaptive optics system to sense all of the perturbations is limited. Because in the case of strong turbulence where the value of Rytov number exceeds 0.3, phase values larger than 2π wrapped to the principal-values which are confined to the range [-𝜋, 𝜋], unwrapping of the phase is required to be able to obtain true PtV and RMS wave-front error. We used the software package [START_REF] Zhao | Robust 2D Phase Unwrapping Algorithm Based on the Transport of Intensity Equation[END_REF] which contains a robust phase unwrapping algorithm based on the transport of intensity equation. The aberration coefficients of the unwrapped phase corresponding to different C n 2 values are calculated. The resulting RMS wave front error and PtV values are shown in Table 2.

CONCLUSIONS

In this work, we aimed to determine the design parameters of the AO system that are related to the atmospheric turbulence. We calculated the PtV and RMS wave front errors in the beam for different atmospheric turbulent conditions. Our plan includes making experiments that measure the wave-front errors of the laser beam traveling through the atmosphere and compare measurement results with the results obtained from the simulations.

Figure 1 :

 1 Figure 1: Gaussian beam propagation in LAtmoSim left: only with turbulence, right: turbulence and thermal blooming

Figure 2 :Figure 3 :

 23 Figure 2: Phase screen with only high-frequency component (left), the low frequency component of the phase created by using the sub-harmonic method (middle) and the total phase screen (right)

Figure 4 :

 4 Figure 4: Holtslag model comparison with data

Figure 5 :

 5 Figure 5: Comparison of estimated and measured C n 2 values

  cross-correlations of Zernike coefficients of turbulent phase with Kolmogorov statistics, Noll[START_REF] Noll | Zernike Polynomials and Atmospheric Turbulence[END_REF] deduced that residual wave-front error for piston removed and tip-tilt removed are respectively 3

Figure 6 :

 6 Figure 6: Reconstruction of the phase by Zernike polynomials with 37 modes

Figure 7 :

 7 Figure 7: The change observed in the Z1 coefficient as the number of Zernike terms added in the calculations at different sampling values; 128(red), 256(black) and 512(yellow)

Figure 8 :

 8 Figure 8: Strehl ratio calculated after correction of each mode

Table 2 :

 2 RMS and peak-to-valley values of wave-front error calculated for different levels of atmospheric turbulence

	C n 2 (m -2/3 )	RMS(µm)	PtV(µm)
	1.0e-15	0.02	1.1
	1.0e-14	0.24	2.1
	1.0e-13	1.19	6.2
	2.0e-13	2.09	9.8


text:
Luca Rinaldi

Vincent Michau

Mathieu Boutillier

Coupling of atmospheric perturbed optical beams with guided modes of propagation

Keywords: Adaptive optics, laguerre-gaussian modes, free-space optical communications, atmospheric turbulence

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Optical communications in free space are of great interest because of the throughput they allow, especially for satellite-to-ground links. In addition to the high data rate, optical wavelengths allow slightly divergent beams, thus reducing the SWAP of the devices associated with this technology. Conversely, optical beams are sensitive to atmospheric turbulence. The optical turbulence causes amplitude and phase distortions that lead to fading of the signal coupled into a single-mode optical fibre (SMF). The coupling of such an optical beam into an SMF has already been described, [START_REF] Ruilier | Coupling of large telescopes and single-mode waveguides: application to stellar interferometry[END_REF] including in the presence of adaptive optics (AO) compensation.

Our objective here is to characterize the coupling of waves perturbed by atmospheric turbulence in a Multi Mode fiber (MMF), or in a spatial multiplexer such as a multi-plane light conversion (MPLC) [START_REF] Bade | Fabrication and Characterization of a Mode-selective 45-Mode Spatial Multiplexer based on Multi-Plane Light Conversion[END_REF] or a photonic lantern. [START_REF] Leon-Saval | Modeselective photonic lanterns for space-division multiplexing[END_REF] In all these cases, the coupling may be described by the decomposition of the wave into a set of propagation modes such as Hermite-Gaussian (HG) or Laguerre-Gaussian (LG) modes. In the following, we have chosen the LG modes which does not restrict the generality of our study since any mode of the first family may be expressed as a linear combination of the elements of the other family. [START_REF] O'neil | Mode transformations in terms of the constituent Hermite-Gaussian or Laguerre-Gaussian modes and the variable-phase mode converter[END_REF][START_REF] Cox | The Resilience of Hermiteand Laguerre-Gaussian Modes in Turbulence[END_REF] Used in various applications focused on their orbital angular momentum, [START_REF] Li | Power loss mitigation of orbital-angular-momentum-multiplexed free-space optical links using nonzero radial index Laguerre-Gaussian beams[END_REF][START_REF] Doster | Laguerre-Gauss and Bessel-Gauss beams propagation through turbulence: analysis of channel efficiency[END_REF] such as quantum mechanics or classical telecommunications, [START_REF] Arikawa | Mode diversity coherent receiver with few-mode fiber-coupling for highspeed free-space optical communication under atmospheric turbulence[END_REF] Laguerre-Gaussian modes are also the most popular choice for modal decomposition.

We consider a set of N [START_REF] Tatarski | Wave propagation in a turbulent medium[END_REF] LG modes characterized by their radial and azimuthal degrees, p and l such that 0 ≤ p ≤ N -1 and 0 ≤ l ≤ N -1. [START_REF] Plick | Physical meaning of the radial index of Laguerre-Gauss beams[END_REF] Their complex amplitudes can be written in the paraxial assumption, in circular-cylindrical coordinates (r, φ, z) as:

LG p,l (r, φ, z) ∝ 1 w(z) r √ 2 w(z) |l| L pl 2r 2 w(z) 2 × exp - r 2 w(z) 2 exp ilφ (1)
where L p,l represents the generalized Laguerre polynomial of order p and degree l. Fig. 1 presents the transverse spatial profile of the 16 first Laguerre-Gaussian modes, at waist (z = 0), for four different values of radial and azimuthal degrees. In the following, the modes are denoted as LG i , where i is equal to i = l + N p.

e-mail: luca.rinaldi@onera.fr Finally, the complex amplitude in the pupil plane may be written as :

P 0 Ψ turb = N 2 -1 i=0 a i LG i + Ψ res , (2)
where Ψ turb is the complex amplitude of the incident beam, and P 0 the characteristics function associated to the telescope pupil.

SATELLITE-GROUND OPTICAL LINK: COUPLING EFFICIENCY

Coupling efficiency is defined as the ratio of the average power coupled into a fiber to the average power in the receiver aperture plane. The coupling between the turbulent field and a mode LG i propagating in the spatial multiplexer is given by the normalized overlap integral:

3 S = < P 0 Ψ turb |LG i > (< P 0 Ψ turb |P 0 Ψ turb >< LG i |LG i >) 1/2 (3)
where < | > denotes the inner product. The fields are assumed normalized in the pupil plane:

S =< P 0 Ψ turb |LG i >= a i .
The coupling efficiency, η, between a turbulent field and a MMF, characterized by N [START_REF] Tatarski | Wave propagation in a turbulent medium[END_REF] LG modes, writes:

η = N 2 -1 i=0 |a i | 2 (4)
Section 2.1 is dedicated to the coupling efficiency in absence of atmospheric turbulence for different sets of Laguerre-Gaussian modes is presented. The effects of turbulence on the perturbed optical beam are illustrated in section 2.2. Conclusion are drawn in section 3.

Coupling in case of a plane wave

The coupling efficiency η has been estimated by numerical simulation in case of a plane wave and for sets of growing numbers of LG i modes such that 1 ≤ N 2 ≤ 49. Fig. 2 presents η as function of the ratio D/2w, D being the pupil diameter and w the waist size of the considered set of modes. As expected (Fig. 2 (a)), the coupling efficiency η is equal to 0.82 for N 2 = 1 (Gaussian profile). Increasing the number of LG modes increases the coupling efficiency: the maximum value of 0.95, reached with the set of 49 modes, is obtained for D/2w ≈ 3, 6. Furthermore, in Fig. 2 (b) η is presented for a set of N 2 = 36 Laguerre-Gaussian modes, function of a pupil with a central occultation here noticed by T , from 0 to 50% of the total pupil diameter D. As shown in the figure above, with a set of 36 LG mode the effect of central occultation on the coupling efficiency is not significant as long as the occultation rate is less than 40%.

Coupling in case of atmospheric turbulence

Assuming the atmospheric turbulence described by Kolmogorov's law, its effects on wave propagation can be characterized by the Fried parameter, r o . We have simulated the wave propagation between a satellite and the ground with a wave optics code developed in our laboratory. [START_REF] Védrenne | Turbulence effects on bi-directional ground-to-satellite laser communication systems[END_REF] The propagation is representative of a geostationary (GEO) -ground link at 30 • of elevation. The main parameters characterizing the propagation are: r 0 = 6.9 cm and σ 2 χ = 0.086 at λ = 1550 nm. The diameter of the receiver pupil is equal to D = 25 cm. The complexe amplitudes in receiver pupil plane are decomposed on N 2 = 49 Laguerre-Gaussian modes. We have estimated the covariance matrix of coefficients a i from 40 complex amplitudes obtained with our propagation code. Fig. 3 presents the cumulative term of the eigenvalues of the covariance matrix: it represents the coupling efficiency as function of the number of eigenmodes used for the coupling. A coupling rate of 54% is obtained with a dozen degree of freedom.

CONCLUSION

The coupling efficiency η between an incoming disturbed beam and a set of Laguerre-Gaussian modes was characterized here function of the pupil diameter D normalized by the waist size, also in presence of a central occultation T . A value of η equal to 0.95 is obtained with a combination of 49 LG modes and a circular pupil, while the effect of a central occultation is not significant increasing the number of LG modes. Furthermore, we evaluate the covariance matrix < a i a j >. The average of the coupling efficiency was evaluated function of the number of eigenmodes used for the coupling, obtaining η = 0.54 with a set of 49 LG modes in which the incoming beam is decomposed.

Figure 1 .

 1 Figure 1. Spatial profile (intensity: left, phase: right) of the 16 first Laguerre-Gaussian modes.

Figure 2 .

 2 Figure 2. Coupling efficiency as function of the pupil diameter normalized by the waist, considering (a) a finite set of LG modes of growing number, (b) a set of 36 LG modes with a central occultation T from 0 to 50% of the pupil diameter D.

Figure 3 .

 3 Figure 3. Average of the coupling rate function of the eigenvalues k with a combination of N 2 = 49 LG modes


text:
Nicolas Boehrer

Robert Nieuwenhuizen

Judith Dijk
email: judith.dijk@tno.nl

Using event cameras for imaging through atmospheric turbulence

Keywords: Neuromorphic sensors, event cameras, anisoplanatic turbulence, image reconstruction

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Light traveling through the atmosphere encounters turbulent regions which modify the optical path length [START_REF] Roggemann | Imaging through turbulence[END_REF]. As the light propagates, the effect of turbulent regions accumulate leading to a random phase distortion of the wave front which causes time-varying blurs and shifts in the image recorded by a camera. Atmospheric turbulence therefore limits the effective resolution of optical imaging in many long range observation applications like surveillance or astronomy.

The effect of turbulence can be mitigated by either using hardware capable of measuring and correcting for the wave front distortion while recording (adaptive optics) or by using software: image processing techniques [START_REF] Eekeren | Turbulence compensation: an overview[END_REF]. Most image processing techniques rely on a combination of motion compensation, sharp region (often called lucky region) identification, multiframe data fusion and deconvolution [START_REF] Vorontsov | Anisoplanatic imaging through turbulent media: image recovery by local information fusion from a set of short-exposure images[END_REF][START_REF] Fishbain | Real-time stabilization of long range observation system turbulent video[END_REF][START_REF] Aubailly | Automated video enhancement from a stream of atmospherically-distorted images: the lucky-region fusion approach[END_REF][START_REF] Eekeren | Patch-based local turbulence compensation in anisoplanatic conditions[END_REF][START_REF] Huebner | Turbulence mitigation of short exposure image data using motion detection and background segmentation[END_REF][START_REF] Anantrasirichai | Atmospheric turbulence mitigation using complex wavelet-based fusion[END_REF][START_REF] Oreifej | Simultaneous video stabilization and moving object detection in turbulence[END_REF][START_REF] Zhu | Removing atmospheric turbulence via space-invariant deconvolution[END_REF][START_REF] Anantrasirichai | Atmospheric turbulence mitigation for sequences with moving objects using recursive image fusion[END_REF][START_REF] Nieuwenhuizen | Dynamic turbulence mitigation for long-range imaging in the presence of large moving objects[END_REF]. Those algorithms are applied on classical cameras in which a total frame of pixels is recorded. However, the frame rate of classical cameras is typical too low to capture all the dynamics of the turbulence induced changes in the images. As a result, local regions in single frames combine instances where they are temporarily sharp with instances where they are less sharp, and tip/tilt aberrations are averaged to further blur these regions.

Neuromorphic cameras do not record an entire frame but instead output an asynchronous stream of changes, so-called events. Contrary to classical frame recordings that record a lot of redundant data, i.e. image regions that stay constant, this new recording technique allows to record only the local changes between frames so that the bandwidth and the recording resources are best used to record the local dynamics of the scene. With its low latency and high temporal sampling, the neuromorphic camera is therefore expected to be well suited to record the temporal variations of the atmospheric turbulence otherwise unseen by conventional camera. This additional information may improve the quality of the restored turbulence free image.

Here we investigate how the event stream may be used to for this turbulence mitigation. In section 2 and 3 we present the event camera used and describe our experimental setup. In section 4 we investigate if the events themselves can be related to properties of the local shifts and blurs in the scene caused by turbulence. Subsequently in section 5 we derive an algorithm to reconstruct high resolution frames by combining the event stream with conventional intensity frames. In section 6 we show the results of this initial algorithm. To our knowledge, there is currently no published approach for turbulence correction using a neuromorphic sensor in the image plane of a camera.

EVENT CAMERA

Unlike conventional cameras that record entire frames synchronously, event camera only record logarithmic intensity changes. The event camera encodes the changes as an asynchronous series of spikes called events. Each event 𝑒 = [𝑡 , 𝑢 , 𝑝 ] carries the time stamp 𝑡 , the pixel location 𝑢 = [𝑥 , 𝑦 ] and the polarity (on or off) of the log intensity change. As shown in Figure 1, an event is produced when the difference between the memorized log intensity and the current log intensity exceeds a certain threshold (controlled by the user). When the threshold is passed, the time stamp, the pixel position and the polarity is emitted and a the current log intensity value is memorized for the monitoring. With a temporal sampling that adapts natively to the actual dynamic of the scene, the event camera offers multiple advantages: high temporal resolution (a few microseconds) with a low latency, a very large dynamic range (120dB) and a low power and bandwidth consumption.

Event cameras belong to the family of gradient cameras, which react to intensity changes and thus a need scene with enough contrast (mostly edges) and movement from the camera or the scene. Even though neuromorphic imaging is relatively recent in the field, it gained a lot of attention and event cameras are already available as commercial product. Many applications and research directions have been investigated. See [START_REF] Gallego | Event-based Vision: A Survey[END_REF] for a recent survey. Among all research directions that have been investigated, the ones with strong camera motion or object motion are the most promising. Atmospheric turbulence produces time varying intensity changes in the form of apparent scene motion. Continuous motion recording therefore implies the potential offered by an event camera to outperform conventional frame approaches.

EXPERIMENTAL SETUP

For the research reported here we used the DAVIS346 dynamic vision sensor from iniVation [START_REF] Inivation | Inivation launches next generation event based dynamic vision sensor[END_REF]. Error! Reference source not found. shows the specification of this camera. The DAVIS346 has the particular benefit of recording both events and intensity frames. In order to collect test material, we recorded video sequences through man-made, indoor turbulence. We produced turbulent air flows with a hot plate placed close to the camera and placed a flat test chart with enough contrasted content at a distance of about 4m. We used the camera with the Avenir lens 16-160mm f/2.0 at 100 and 160mm. The lens was focused on a siemens star placed on the target before recording. After a few attempts with different heaters, we decided to use a hot plate without fan, as a fan reduces the turbulence. The resulting flow pattern seems like a reasonable approximation of operationally relevant turbulence. The magnitude of shifts and blurs as well as length scales in the images appeared consistent with previously recorded data in a NATO-SET trial in Quebec [START_REF] Velluet | Data collection and preliminary results on turbulence characterisation and mitigation techniques[END_REF] (NATO Quebec).

EVENT STATISTICS

In section 1 we proposed that event cameras may be beneficial for capturing the dynamics of the turbulence induced changes in images in two respects: identifying lucky patches and measuring the turbulence induced shifts in the images.

In this section we will explore these hypotheses by analysing the spatiotemporal statistics of the event counts.

We first investigated how the event stream relates to the instantaneous sharpness (i.e. luckiness) of the images. To assess this we focused on the variations of a small region of interest (ROI) over 300 frames. The effect of turbulence is a combination of a local spatial shift and a local blur. Our hypothesis is that sharp zones that change only due to shifts produce higher event counts than zones that change due to a variation in blur. The sharp zones can be considered as lucky image patches that once stabilized could deliver sub-pixel sampling to the corrected image estimate. However, when comparing the event counts to the local contrast in an ROI, which is mostly determined by the local sharpness, no evidence of a correlation between these quantities is found as shown in Figure 3. To further analyse the relation between sharpness and event counts we characterized the influence of the lens aperture on the event stream. We measured the event count for different f-numbers and noticed that the event count increased with the f-number. Imaging with a smaller aperture visibly resulted in a video disturbed more by local motion and less by local blur variations. This also resulted in lower event rates as shown in Figure 4. Next we considered if the event camera allows for more accurate measurement of local image shifts due to turbulence. The event stream has the potential to contain high frequent information about local intensity changes. To evaluate if the event stream contains more information than intensity frames about the local shifts induced by turbulence, we measured the frequency at which we could locally estimate the motion. Linearly moving edges typically produce events which correspond to 3d points that reside on a plane in x-y-t space. The slope of the plane corresponds to the motion velocity. Therefore local motion estimation with the event stream consists of fitting a plane on the recent points around the current event location. Figure 5 shows how often a weak plane fit can be done relative to the rate of intensity frames, assuming that we need at least 3 points in a 5x5 neighborhood.

Figure 5: Average number of possible fit per frame using a 5x5 neighborhood Some regions show rates up to 2 motion estimation (plane fit) per frame, but rates generally stay below the 1 estimation per frame rate. Even though the event stream may contain high frequent changes, the turbulence phenomenon does not seem to generate those high frequent variations and the intensity images already carry most of the local shift information.

IMAGE RECONSTRUCTION ALGORITHM

The preceding section did not directly show a relation between event counts and the sharpness of image patches. However, for image reconstruction it is highly desirable to identify image regions (in space and time) with relatively high sharpness and build the final corrected output out of the subset of "lucky" regions. Therefore we implemented an image reconstruction algorithm in which we more directly measure if regions are "lucky" by measuring their instantaneous sharpness. Figure 6 gives an overview of our proposed algorithm. We derive our algorithm by modeling the event camera as a high speed camera that records only changes. The algorithm first extracts a virtual "high speed" frame at a time t by reconstructing the intensity from the previous intensity frame and the sequence of events up to the time tag t.

Given a set of events {k} from 𝑒 = [𝑡 , 𝑢 , 𝑝 ] that occurred at time 𝑡 ∈ [𝑡 , 𝑡] with 𝑡 < 𝑡 , location 𝑢 = [𝑥 , 𝑦 ] with a polarity 𝑝 ∈ [+1, -1] and the corresponding constant relative brightness change C (here identical for positive and negative polarity), one can construct the difference image relative to the last intensity frame:

𝛿𝐼(𝑥) = 𝐶 𝑝 𝛿(𝑥 -𝑢 )
As illustrated in Figure 6, the difference image is used to iteratively update the raw log intensity image and create a virtual frame 𝐼 at any point between intensity images 𝐼 :

𝐼 = 𝐼 + 𝛿𝐼
The contribution of each ON or OFF event is calibrated with the first intensity frames so that updating 𝐼 with the contributions of all events between frame 𝐼 and 𝐼 best matches the intensity measured at frame 𝐼 .

Once the virtual raw frame is available, the algorithm removes the local shifts induced by the turbulence. Before processing the sequence, a reference image is first computed by averaging all intensity frames. The local shifts between each intensity frame and the reference image are estimated using Farnebäck dense optical flow [START_REF] Farnebäck | Two-Frame Motion Estimation Based on Polynomial Expansion[END_REF].

The virtual raw frame is warped to the reference space based on the optical flow interpolated between two intensity frame times.

𝐼 = 𝐼 (𝑊(𝑢, 𝑡))

Here the continuous warp is the linear interpolation between the two optical flows from two subsequent frames:

𝑊(𝑢, 𝑡) = 𝑊(𝑢, 𝑡 ) * 𝛼 + 𝑊(𝑢, 𝑡 ) * (1 -𝛼) 𝛼 = 𝑡 -𝑡 𝑡 -𝑡
With 𝑊(𝑢, 𝑡 ), {𝑖 ∈ ℕ | 0 < 𝑖 < 𝑁}, N the number of intensity frames, the warp computed from the optical flow between each intensity frame and the average reference frame.

In the last step, the algorithm updates the final turbulence corrected image based on a weighting scheme that emphasizes lucky regions. We use the Sobel response as a sharpness (and therefore luckiness) metric.

𝐼 = 𝐼 𝐻 𝐼 + 𝐼 𝐻(𝐼 ) 𝐻 𝐼 + 𝐻(𝐼 )

With

𝐻(𝐼) = 𝐻 (𝐼) + 𝐻 (𝐼) 𝐻 (𝐼) = +1 0 -1 +2 0 -2 +1 0 -1 * 𝐼 𝐻 (𝐼) = +1 +2 +1 0 0 0 -1 -2 -1 * 𝐼

IMAGE RECONSTRUCTION RESULTS

We are interested in applying the lucky imaging technique described in the previous section onto the event stream to:  Take into account the local turbulence motion generating the events. We want to verify if we are able to capture high frequency motion that would make the edge position estimate converge towards the true edge position.

 Benefit from the temporal resolution of the camera. The virtual frame may contain more lucky regions than the ones visible in the intensity frames

We ran the pipeline on multiple targets recorded with the test setup described above. We compared our results to the state of the art turbulence mitigation method based on Hardie's super-resolution algorithm [START_REF] Hardie | High-resolution image reconstruction from a sequence of rotated and translated frames and its application to an infrared imaging system[END_REF][START_REF] Eekeren | Patch-based local turbulence compensation in anisoplanatic conditions[END_REF]. When compared to a simple stabilization, the turbulence correction using events allows to recover some fine details (see Figure 7, windows, wake on the river) but the gain of event based lucky imaging isn't significant when compared to more classical intensity frame based approach e.g. the Hardie approach.

We also looked at the temporal stability of our estimation pipeline. When looking at videos, we observe that a temporally varying correction quality could be more disturbing than the actual turbulence effects it tries to correct. Therefore we inspected the temporal constancy of the output sequences of our pipeline. This confirmed the stability of the iterative correction: compared to the output of a simple geometric stabilization obtained by only warping raw frames, the turbulence correction reduces the effect of local shift and blur induced by the turbulence, leading to a more stable output. Figure 7 shows an example of a region of interest in one of our result videos, which shows that the output frames of our pipeline (the corrected frames) indeed exhibit a more constant image quality than the stabilized raw frames.

CONCLUSIONS & DISCUSSION

In this paper we report on experiments with a neuromorphic camera for turbulence mitigation. First we investigated the properties of the event data stream in relation to the sharpness of the reference intensity frame. Our analysis did not show a correlation between event counts and the sharpness variations of the intensity frames. Nevertheless, we were able to construct a lucky imaging pipeline that makes use of events and therefore benefits from the high temporal resolution of the camera to better record the continuous variations of the atmospheric turbulence. We applied thus lucky imaging pipeline on video's with "artificial" turbulence. The improvements are relatively modest when compared to intensity frame-based correction e.g. Hardie. There may be different reasons for that:

 The turbulence in our experimental setup is generated with a hot plate close to the sensor. This may not be a good approximation to typical atmospheric turbulence found in practical scenarios, which may exhibit faster temporal dynamics due to crosswinds for example. For such high frequency variations the neuromorphic camera would have an advantage over a classic camera.

 The current turbulence motion compensation relies on a dense optical flow computed between the intensity frames and the reference frame [START_REF] Farnebäck | Two-Frame Motion Estimation Based on Polynomial Expansion[END_REF]. The motion to compensate the event position is then interpolated between the frames. This linear interpolation may be inaccurate.

 The sharpness measure used for updating the corrected images relies on a Sobel response. The Sobel operator is designed to measure continuous gradients rather than the spatially discrete event features.

In future research we will record experimental data with natural turbulence in an outdoor imaging setup to further verify our current conclusions. On the image reconstruction side we will explore different possible improvements: spatial localization of the event updates, 2) sharpness measures and 3) weighting schemes for the event-based and intensity framebased updates in the reconstruction algorithm.

Figure 1 :

 1 Figure 1: DVS Pixel operation

Figure 2

 2 Figure 2 Illustration of the experimental setup.

Figure 3

 3 Figure 3 Statistics of event counts in relation to local image sharpness.

Figure 4

 4 Figure 4 Visualization of local event counts between two intensity frames for different aperture settings. Positive events are shown in green, negative events in red.

Figure 6

 6 Figure 6 Schematic overview of the image reconstruction algorithm.

Figure 7

 7 Figure 7 Comparison of image quality of the event based image reconstruction (corrected frames) with the Hardie super-resolution based reconstruction and ground truth.

Figure 8

 8 Figure 8 Comparison of the stability over multiple frames between locally stabilized raw frames and the event based image reconstruction (corrected frames).

Table 1 :

 1 Specifications of the DAVIS346 dynamic vision sensor of iniVation.

	Feature	Value
	Resolution	346 × 260
	Time resolution event stream	1 µs
	Frames	Gray scale
	Dimensions	40×60×25 mm 3 (H×W×D)
	Weight	100 g (without lens)

ACKNOWLEDGEMENTS

The work was funded by the Office of Naval Research Global.