Fabrice Bartolomei
email: <fabrice.bartolomei@ap-hm.fr>

Stanislas Lagarde

Isabelle Lambert

Agnès Trébuchon

Samuel Médina Villalon

Aileen Mcgonigal

Christian Bénar

Agnes Trébuchon

Christian G Benar

Brain connectivity changes during ictal aggression (a

Keywords: SEEG, functional connectivity, aggression, focal seizure *

Minor episodes of aggression are relatively common in patients with epilepsy but they are probably no more frequent than in populations with social disadvantages or brain damage [START_REF] Rantakallio | Association of perinatal events, epilepsy, and central nervous system trauma with juvenile delinquency[END_REF][START_REF] Grafman | Frontal lobe injuries, violence, and aggression: a report of the Vietnam Head Injury Study[END_REF][START_REF] Treiman | Violence and the epilepsy defense[END_REF][START_REF] Van Elst | Aggression and Violence in Patients with Epilepsy[END_REF][START_REF] Fazel | Risk of violent crime in individuals with epilepsy and traumatic brain injury: a 35-year Swedish population study[END_REF][START_REF] Piazzini | Aggressive behavior and epilepsy: a multicenter study[END_REF]. Three different forms of aggressive behaviours have been described in epileptic patients on the basis of their temporal relationship with seizures: interictal, post-ictal, and ictal (Tassinari et al., 2005a;[START_REF] Bronsard | Rhythms, rhythmicity and aggression[END_REF]. Of these, postictal aggression, while the patient is in a confused state, is likely to be the most prevalent. However, reported cases with good quality video-EEG recording are rare (Tassinari et al., 2005a) and as a result, it may be difficult to distinguish between true semiological seizure phenomena and postictal confusional behaviour. In seizures affecting brain regions involved in emotion processing, ictal or postictal aggressive behaviour may occur. Defensive-like behaviour can be F. Bartolomei, et al. observed, particularly during the immediate postictal period, in which the patient is often confused and in an altered mood state [START_REF] Delgado-Escueta | Special report. The nature of aggression during epileptic seizures[END_REF]. In this context, aggressive behaviour is often related to minimally unpleasant stimulation or attempts to restrain the patient (Tassinari et al., 2005a). This phenomenon is rare and has been generally classified as a defensive/aggressive behaviour, related to the interaction between the patient and the caregivers (i.e. nurse during video-EEG monitoring) [START_REF] Bronsard | Rhythms, rhythmicity and aggression[END_REF]. It is generally thought that the aggressive behaviour is not usually truly dangerous. More serious offences can occur during postictal psychosis, but the clinical picture is close to premeditated/predatory aggression [START_REF] Kanemoto | Violence and postictal psychosis: a comparison of postictal psychosis, interictal psychosis, and postictal confusion[END_REF]. The current view of the neurobiology of aggressive behaviour is based on a model of an imbalance between regulatory regions (essentially prefrontal) and regions engaged in emotional and behavioural responses (amygdala, hypothalamus, and brain stem) [START_REF] Siever | Neurobiology of aggression and violence[END_REF][START_REF] Rosell | The neurobiology of aggression and violence[END_REF]. Such a scheme has been proposed for seizures associated with the emergence of stereotyped emotional abnormal behaviour, such as rage or aggression [START_REF] Bartolomei | Acute alteration of emotional behaviour in epileptic seizures is related to transient desynchrony in emotionregulation networks[END_REF], Tassinari et al., 2005a, 2005b). Ictal semiology does not depend on the involvement of a unique region but appears to be related to the involvement of a set of regions within specific anatomical systems [START_REF] Bartolomei | Abnormal binding and disruption in large scale networks involved in human partial seizures[END_REF]. This is particularly true for complex ictal behaviours that necessarily depend upon altered dynamics within widely distributed, rather than focal or regional, neural networks. In this way, considering large-scale network involvement during seizure evolution is a way to investigate the emergence of clinical patterns during seizures. Thus, the study of functional connectivity during seizures has been proposed in the past to be used to quantify the temporo-spatial changes in network properties [START_REF] Bartolomei | Abnormal binding and disruption in large scale networks involved in human partial seizures[END_REF].

In the present study, we report a case in which, during one recorded seizures, we observed violent behaviour, clearly occurring during the ictal phase, culminating in an aggressive act with a strangulation attempt. We describe both SEEG aspects and intracerebral EEG connectivity changes. During the aggressive behaviour part of the seizure, major changes in functional connectivity (Fc) occurred between frontal and temporal regions.

Case study

A 33-year-old male patient underwent presurgical evaluation for his drug-resistant epilepsy in our epilepsy unit. He had had an encephalitic illness at the age 21, presenting with fever and recurrent seizures in the acute phase. He was treated at this time with acyclovir IV and amoxicillin, and outcome was finally favourable with no neurological or significant cognitive deficit. Seizures developed six months later and were resistant to antiepileptic drugs. Non-invasive presurgical evaluation led to the recording of prominent left side fronto-temporal seizures. Psychiatric evaluation did not reveal major personality disorder or an anxio-depressive state. Magnetic resonance imaging (MRI) did not show any structural abnormality. Intracerebral EEG (stereoencephalography [SEEG]) was subsequently performed for presurgical evaluation. Ten electrodes were placed in brain regions according to hypotheses regarding seizure onset and propagation (mainly the temporo-frontal left side), including the right and left temporo-frontal cortices (figure 1A).

Connectivity analysis

Interdependencies between SEEG signals were estimated with Anywave software [START_REF] Colombet | AnyWave: a cross-platform and modular software for visualizing and processing electrophysiological signals[END_REF] which computes a pairwise non-linear regression analysis based on the h 2 coefficient. In summary, a piecewise linear regression is performed between each pair of signals, testing all the shifts of one signal relative to the other within a maximum lag. The h 2 is the coefficient of determination that measures the goodness of fit of the non-linear regression. The h 2 is bounded between 0 (no correlation) and 1 (maximal correlation) [START_REF] Wendling | Modeling EEG signals and interpreting measures of relationship during temporal-lobe seizures: an approach to the study of epileptogenic networks[END_REF]. We used a sliding window of four seconds with an overlap of two seconds, and a maximum delay between signals of 100 ms. Signals were filtered using the Anywave software [START_REF] Colombet | AnyWave: a cross-platform and modular software for visualizing and processing electrophysiological signals[END_REF] (available at http://meg.univamu.fr/wiki/AnyWave), with Butterworth filters of order 4. We used only a 0.5-Hz high-pass filter for all the analyses.

Graph measures

For all the selected channels (each channel being a bipolar derivation), we computed all the pairwise h 2 values. Between two channels, we selected the higher h2 value between chan1-> chan2 and chan2-> chan1. We thus obtained connectivity graphs, with each channel representing a node of the graph, and h2 values the strength of the link between two nodes. We summarized the connectivity graphs with two related graph measures. The first measure, node degrees, first involves thresholding the graphs, and then counting the number of significant links between a given node and the rest of the graph. This has the advantage of removing the links with low connectivity. The threshold was set empirically to 0.25. The second measure, node strength, consists of the mean h2 between one node and all the other nodes. It has the advantage of not relying on a threshold that can be difficult to tune.

Statistical analysis

Node degrees and connectivity strength were compared between a period preceding seizure onset (background and, in particular, the seizure onset) and the period in which the aggression occurred (end of seizure). A comparison was performed using a Wilcoxon non-parametric paired test and a Bonferroni correction was applied. aggressive, showing combative behaviour with groaning, swearing, and threatening the nurse with his fist. Despite verbal reassurance, the patient then grabbed the nurse and attempted to strangle her (figure 1B).

Video-SEEG and intracerebral EEG connectivity

After several seconds and thanks to the intervention of another nurse, the patient stopped and finally returned to his bed. After the seizure, the patient was completely amnesic for the episode. Electro-anatomo-clinical correlations for this seizure are shown in figure 1. The seizure started from the left temporal region (lateral and medial parts), becoming rapidly bilateral and spreading to the prefrontal cortex. Aggressive behaviour coincided with marked slowing of bilateral frontal region activity occurring in the last part of the seizure and was associated with the persistence of ictal activity under the form of burst activities in the right medial temporal lobe (in particular, the medial contacts of electrodes B which were used to investigate the anterior hippocampus/amygdala region).

A connectivity study (figure 2B) was performed between 18 bipolar derivations, recording different brain areas, including the two temporal regions and the left prefrontal cortex. This analysis revealed that seizure onset was associated with connectivity changes (an increase) limited to the left temporal regions. The global connectivity was not significantly changed (figure 2A) at seizure onset (p=0.8), but largely increased at seizure end (p<0.0001). During the aggressive phase (end part of the seizure), a massive change, characterized by increased connectivity, was indeed observed bilaterally, affecting most of the recorded connections between frontal and temporal regions. This hypersynchrony was observed in most of the studied regions, both in terms of strength of connectivity and degrees (figure 2D). We also performed the same analysis during a second temporal seizure for the same patient, starting from the left temporal cortex. In this seizure, despite the same examination by the nurse, no aggression occurred. This seizure started the left mesial and lateral temporal cortex, spreading to the frontal lobe. This seizure was shorter and the significant frontal lobe slowing observed during the seizure with aggression was not evident. Connectivity changes are indicated in figure 3 and did not show the same pattern of fronto-temporal hypersynchrony.

Discussion

We report here seizure-related aggressive behaviour that clearly occurred during the ictal phase of a temporo-frontal seizure, recorded using intracerebral electrodes. The clinical changes observed in our case can be classified as defensive/aggressive ictal behaviour, triggered by the verbal interaction between nurse and patient. However, this was associated with a severe aggressive act (strangulation attempt) which is fairly unusual in this context. In agreement with previous reported cases, our patient had seizures involving the temporo-frontal cortices (Tassinari et al., 2005a;[START_REF] Bronsard | Rhythms, rhythmicity and aggression[END_REF]. The failure of "top-down"control systems in the prefrontal cortex to modulate aggressive acts (with a possible triggering stimulus) appears to play an important role in aggressive behaviour [START_REF] Davidson | Dysfunction in the neural circuitry of emotion regulation: a possible prelude to violence[END_REF][START_REF] Siever | Neurobiology of aggression and violence[END_REF]. An imbalance between prefrontal regulatory influences and hyper-responsivity of the amygdala and other limbic regions involved in affective evaluation is a putative mechanism of aggressive human behaviour [START_REF] Siever | Neurobiology of aggression and violence[END_REF]. The observations in the present case could fit well with this model, since the patient's aggressive behaviour was concomitant with persistent epileptic activity in the limbic system (internal temporal region including the amygdala) and massive slowing of cortical prefrontal regions. Connectivity analysis showed that the brain regions investigated, particularly the frontal-temporal connections, were engaged in massive abnormal hypersynchrony. It could be hypothesized that this pathological hypersynchronous state would lead to a breakdown of normal regulatory function of the cortical regions. Similar patterns of long-range hypersynchronous activity changes have been correlated with ictal altered consciousness, in line with the theory that seizure-related hypersynchrony disrupts the "global workspace of consciousness" [START_REF] Bartolomei | The global workspace (GW) theory of consciousness and epilepsy[END_REF]. Changes in connectivity between amygdala/mesial temporal regions and the frontal lobes have been described in seizures, including profound alterations in emotional behaviour such as screaming, agitation, and facial expression of terror [START_REF] Bartolomei | Acute alteration of emotional behaviour in epileptic seizures is related to transient desynchrony in emotionregulation networks[END_REF]. Some of these clinical symptoms have been interpreted as the emergence of "archaic" or innate motor repertoires (Tassinari et al., 2005b), in which impaired higher cortical control would allow "pre-programmed" and essentially subcortically organised motor behaviours to be inappropriately expressed or "released". This idea of hierarchical control of behaviour can indeed be traced back to Hughlings Jackson [START_REF] Jackson | Selected Writing of John Hughling Jackson Vol 1: On epilepsy and Epileptiform Convulsions[END_REF]. In this context, aggressive behaviour, such as biting, is an interesting finding observed in some patients during seizures [START_REF] Bartolomeil | Fear, anger and compulsive behavior during seizure: involvement of large scale frontotemporal neural networks[END_REF]Tassinari et al., 2005a) that may indeed reflect the stimulation of an innate program because of its similarity to animal behaviours (Tassinari et al., 2005a).

In conclusion, our case is a rare report of intracerebrally recorded ictal aggressive behaviour associated with altered network dynamics within fronto-limbic circuits. Clearly, no firm conclusions can be drawn from analysis of a single seizure, in terms of whether the observed dysfunction within these networks might have an eventual causal role or not. However, the rarity of the observation appears to merit discussion. Current methods of signal analysis of intracerebrally recorded seizures allow estimation of network dynamics correlated with evolving clinical signs, with notably high temporal resolution in comparison to other methods, such as functional MRI. The electrical data obtained is entirely dependent on the recording sites selected according to clinical indication for presurgical evaluation. As such, the vast majority of available data in seizure anatomo-electrical clinical correlation are taken from small-volume samples of cortical and/or limbic structures, albeit precisely defined anatomically and chosen in order to logically sample known networks. Understanding of the neural basis of complex ictal behaviours remains extremely elusive and it can be imagined that in order to make progress in this area, current methods would need to be combined with means of studying dynamic change within more widespread networks, including multiple subcortical structures.

Figure 1 .

 1 Figure 1. (A) SEEG implantation showing the electrodes on a 3D MRI mesh of the cortex. TB and TB' electrodes were used to investigate the temporo-basal cortex to within the entorhinal cortex; B and B' electrodes for the anterior hippocampus/amygdala region (medial leads) and the mid part of MTG (lateral leads); C' electrode for the posterior hippocampus and the lateral temporal cortex; H' electrode for the insula (medial leads) and the anterior part of the superior temporal gyrus (STG) (lateral leads); OF' electrode the frontal operculum; OP' electrode for the parietal operculum; CR' electrode for the prefrontal cortex to within the anterior cingulate region; and OR' electrode for the prefrontal cortex (superior leads) and the orbitofrontal cortex (inferior leads). The prime symbol refers to the left side. (B) Captures from the video-SEEG recordings showing the patient jumping from his bed and grabbing the nurse. (C) SEEG traces showing the seizure onset (SO) and a discharge starting from the left temporal neocortex (blue arrow), before spreading (middle seizure [MS]) to the suprasylvian cortices (H', OF') and the frontal cortices (OR') and affecting the right temporal lobe. The last part (end of seizure [ES]) is characterized by a slowing in the left regions and a burst activity in the right temporal electrodes. L: left; R: right; OFC: orbitofrontal cortex; DLPFC: dorsolateral prefrontal; CG32: cingulate gyrus Brodmann's area 32; OPF: frontal operculum; STG: superior temporal gyrus; TBC: basal temporal cortex; Hip: hippocampus; EC: entorhinal cortex; MTG: middle temporal gyrus.

Figure 2 .

 2 Figure 2. Connectivity analysis. (A) Global changes in connectivity during a seizure with aggression (SO: seizure onset; ES: end of seizure). Mean values over periods of 15 seconds were compared. Changes were significant only for the end of the seizure period (p<0.00001 relative to background [BKG]). (B) Changes in connectivity are illustrated on a 3D mesh of the MRI with the position of electrodes and the graphs thresholded (h 2 >0.25). On the left (seizure onset), the main connected regions are the left temporal contacts, while the end part of the seizure (right) is characterized by connecting links between temporofrontal regions affecting both sides. Colour scale indicates h 2 values. (C, D) For each selected channel, the mean values of h 2 (C) and the degrees (D) are indicated between two conditions (background and seizure end period). *significant interactions (p<0.05 after Bonferroni corrections).

Figure 3 .

 3 Figure 3. Changes in connectivity during a left side-onset temporal seizure without aggressive behaviour. These are illustrated on a 3D mesh of the MRI with the position of electrodes and graphs thresholded (h 2 =0.25). Colour scale indicates h 2 values.

Epileptic Disord, Vol. xxx, No. x, xxx 2017 F.Bartolomei, et al.

Acknowledgements and disclosures.

We thank Pr J Regis for stereotactic procedures and Pr P Chauvel for clinical assessment of the patient. The authors have no conflict of interest to disclose.

Supplementary data.

Summary didactic slides are available on the www.epilepticdisorders.com website.

Author Query

No query


MD Fabrice Bartolomei

Prof Stanislas Lagarde

Samuel Médina Villalon

Aileen Mcgonigal

Christian G Benar

The "Proust Phenomenon": odor-evoked autobiographical memories triggered by direct amygdala stimulation in human

Keywords:

odor-evoked memories: they evoke more emotional and evocative recollections than memories triggered by any other cue and are exceptionally rich in contextual information [START_REF] Larsson | Autobiographical odor memory[END_REF]) [START_REF] Arshamian | The functional neuroanatomy of odor evoked autobiographical memories cued by odors and words[END_REF]Saive, Royet, & Plailly, 2014). The intensity of olfactory recollection and its visceral characteristics have been related to the unique and specific connections of the olfactory system with the neural structures involved in emotion and associative learning [START_REF] Arshamian | The functional neuroanatomy of odor evoked autobiographical memories cued by odors and words[END_REF][START_REF] Herz | Neuroimaging evidence for the emotional potency of odor-evoked memory[END_REF]. The olfactory cortex includes the amygdala, which is also involved in emotional memory and is connected to the hippocampus and in contrast with other sensory modalities, projections from the sensory input onto amygdala do not pass via the thalamus [START_REF] Cahill | The amygdala and emotional memory[END_REF][START_REF] Phelps | Emotional memory: what does the amygdala do?[END_REF]. A main characteristic of the Proust phenomenon is its unpredictability, rendering it particularly difficult to reproduce in an experimental setting [START_REF] Jellinek | Proust remembered: has Proust's account of odor-cued autobiographical memory recall really been investigated?[END_REF].

The present case is the first to analyze the induction of Proust phenomenon by focal electrical stimulation of the amygdala, a rare observation in the context of epilepsy presurgical exploration. Our patient was a 27-year-old, right-handed woman who had suffered from drug-resistant left frontal lobe seizures since the age of 5 years. Cerebral magnetic resonance imaging (MRI) disclosed a small lesion in left fronto-polar cortex. Recorded spontaneous seizures showed complex motor semiology without subjective symptoms, notably no olfactory hallucination or déjà vu phenomena. Intracerebral electroencephalographic (EEG) monitoring was undertaken to precisely localize the epileptogenic zone, with 9 intracerebral electrodes (stereoelectroencephalographic (SEEG) method [START_REF] Bartolomei | Cortical stimulation study of the role of rhinal cortex in deja vu and reminiscence of memories[END_REF] implanted within the temporal lobe and the frontal cortex on the left side. Ten seizures were recorded involving the left prefrontal cortex and a left frontal cortectomy was later performed leading to seizure freedom. Pathological examination revealed a type 2 focal cortical dysplasia. Electrical stimulation was performed (50 Hz, 0.5-2 mA, in a bipolar fashion to each contact in the gray matter during a 3 second period) to map functional cortices and trigger habitual seizures [START_REF] Bartolomei | Cortical stimulation study of the role of rhinal cortex in deja vu and reminiscence of memories[END_REF]. Two stimulations of the amygdala triggered olfactory hallucination with memory reminiscence. signals before and after stimulations. We performed a nonlinear regression analysis based on the h² coefficient (Wendling, Bartolomei, Bellanger, & Chauvel, 2001)(see supplemental material). Figure 1A presents the differences in mean degree (number of significant links) between pre-and post-stim periods. Figure 1B shows connectivity graphs superposed on a 3D rendering of the cortex, with each bipolar channel representing a node of the graph. Link strengths are the h2 values between pairs of nodes, for two positive stimulations (at 1.5 mA and 2.5mA) and one negative stimulation (1 mA) of the amygdala. Positive stimulations induced a significant increase of connectivity values between the amygdala, the temporal pole, and the insular cortex at 1.5 mA. A more extended network was involved at 2.5 mA, including the orbitofrontal and prefrontal cortices. At 1 mA (negative stimulation), only one link between amygdala and temporopolar cortex was elicited without significant change in node degree.

Our observations strongly support the role of the amygdala in inducing the "Proust" phenomenon [START_REF] Bray | Forgetting the madeleine: Proust and the neurosciences[END_REF][START_REF] Jellinek | Proust remembered: has Proust's account of odor-cued autobiographical memory recall really been investigated?[END_REF]Toffolo, Smeets, & van den Hout, 2012), in agreement with the physiological role of the amygdala in processing odor. Autobiographic memories evoked by amygdala stimulation in our case belonged to the adolescent period. A review of mnestic phenomena induced by intracerebral stimulation in patients with temporal lobe epilepsy found that amygdalar stimulation induced déjà-vu more often than hippocampal stimulation (Vignal, Maillard, McGonigal, & Chauvel, 2007) although no olfactory illusion was reported in this series. A large body of data obtained in human beings and animals show that the amygdala participates in various aspects of odor processing, especially in relation to emotion and memory. In particular, during the experience of recollecting an odor-evoked autobiographical memory, the amygdala was more activated than with similar odors that did not evoke a memory [START_REF] Herz | Neuroimaging evidence for the emotional potency of odor-evoked memory[END_REF]. In patients with epilepsy, depth electrode recordings of the amygdala have revealed that odorant stimulation may induce evoked potentials [START_REF] Hudry | Odorants elicit evoked potentials in the human amygdala[END_REF]. Our connectivity study shows that the amygdala stimulation leads to the activation of a network including ipsilateral orbitofrontal cortex and insular cortex, representing "secondary olfactory cortices". In particular, the left anterior insula is likely involved in the evaluation of odor properties [START_REF] Plailly | Involvement of the left anterior insula and frontopolar gyrus in odor discrimination[END_REF] [START_REF] Royet | fMRI of emotional responses to odors: influence of hedonic valence and judgment, handedness, and gender[END_REF]. fMRI studies have found that the amygdala could be activated bilaterally by both positive and negative stimuli (Sergerie, Chochol, & Armony, 2008). A study using direct intracerebral electrical stimulation reported that positive and negative emotional feelings could be triggered by left amygdalar stimulation while only negative feelings were triggered by right amygdalar stimulation [START_REF] Lanteaume | Emotion induction after direct intracerebral stimulations of human amygdala[END_REF]. In addition, neuroimaging studies have revealed that anterior insular cortex can be activated during intensely positive feelings, such as joy or maternal and romantic love as well as seeing or making a smile [START_REF] Craig | How do you feel--now? The anterior insula and human awareness[END_REF]. Thus together with the activation of

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
olfactory regions, the activation of the insula may be responsible for the pleasant feeling associated with the Proust phenomenon.

Saive, A. L., Royet, J. P., & Plailly, J. (2014). A review on the neural bases of episodic odor memory: from laboratory-based to autobiographical approaches. Front Behav Neurosci, 8, 240. doi: 10.3389/fnbeh.2014.00240 Sergerie, K., Chochol, C., & Armony, J. L. (2008). The role of the amygdala in emotional processing: a quantitative meta-analysis of functional neuroimaging studies. Neurosci Biobehav Rev, 32(4), 811-830. doi: 10.1016/j.neubiorev.2007.12.002 Toffolo, M. B., Smeets, M. A., & van den Hout, M. A. (2012). Proust revisited: odours as triggers of aversive memories. Cogn Emot, 26(1), 83-92. doi: 10.1080/02699931.2011.555475 Vignal, J. P., Maillard, L., McGonigal, A., & Chauvel, P. (2007). The dreamy state: hallucinations of autobiographic memory evoked by temporal lobe stimulations and seizures. Brain, 130(Pt 1), 88-99. doi: 10.1093/brain/awl329 Wendling, F., Bartolomei, F., Bellanger, J. J., & Chauvel, P. (2001). Interpretation of interdependencies in epileptic signals using a macroscopic physiological model of the EEG.

Clin Neurophysiol, 112( 7), 1201-1218.

Figure legend:

Connectivity changes associated with Proust phenomenon induced by amygdala stimulations (contacts A'3-4 stimulated at 1.5 mA and 2.5 mA) and in comparison with negative stimulation (1 mA).

A. Node degree connectivity. Connectivity graphs are summarized by a node degrees which consist in thresholding the graph (set empirically to 0.25) and then counting the number of significant links between a given node and the rest of the graph. Comparison were done using a Wilcoxon non-parametric paired test and a Bonferroni correction by comparing a prestimulation period and the stimulation period. * shows significant results.

2

  , Hôpital de la Timone, Service de Neurophysiologie Clinique, Marseille, F-13005 3. AP-HM, Hôpital de la Timone, Service de Neurochirurgie Fonctionnelle, Marseille, F-13005 4. AP-HM, Hôpital de la Timone, Service de Neurochirurgie Pédiatrique, Marseille, F-13005 Correspondance to: Prof Fabrice Bartolomei, MD, PhD Service de Neurophysiologie Clinique, CHU Timone, 264 Rue Saint-Pierre, 13005, Marseille, France Tel: +33491385833 Fax: +33491385826 Email: fabrice.bartolomei@ap-hm.fr Vivid memories triggered by odors were particularly well described by the French writer Marcel Proust in his novel Swann's Way (Du Côté de Chez Swann). The sensorial input provoked by the madeleine cake's odor, flavor and texture immediately transported him into a vivid and rich past childhood episode. Proust constructed a detailed literary description of psychological characteristics of the reminiscence: its unexpected occurrence, the intense positive feeling resembling ecstatic sensation, and the vividness of the memory. Neuroscientific investigation has clarified the characteristics of

  Contacts A'3-4 located in the basolateral region of the amygdala (Fig 1) at 1.5 mA and 2 .5 mA induced a sudden odor of burnt wood, which was immediately followed by the memory of a scene of a campfire on a riverbank during a holiday when she was 15 years old. This scene was associated with a strong feeling of happiness. No electrical after-discharge was noted after stimulations. No other stimulation (including left temporal pole, anterior insula and orbitofrontal cortex) induced such M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 3 phenomena. Connectivity analysis was performed by measuring interdependencies between SEEG

  Selected regions are indicated by bipolar contacts: A'1-2: Amygdala; A'9-10: Lateral temporal cortex FO'2-3: Orbito-frontal cortex; Fp'1-2: Prefrontal cortex (internal); Fp'7-8: Lateral prefrontal cortex; G'1-2: Cingulate gyrus (area 24); K'1-2: Cingulate gyrus (Area 32); OP'1-2: Anterior insular cortex; P'2-3: Parietal cortex; PM'1-2: Internal premotor cortex; TP'1-2: Temporal pole cortex B. Changes of connectivity are illustrated on a 3D mesh of the MRI with the position of electrodes. Each graph indicates the significant changes in connectivity between selected contacts of the electrodes (the color scale express the changes relative to the pre-stimulation period in term of Z-scores, here positive values). Color scale indicate the Z-scores of the h² values during stimulation period relative to pre-stimulation period C. Reconstruction of electrode A' within the left amygdala. In green are indicated the two contacts for which bipolar stimulation induced Proust phenomenon

Acknowledgments

We thank Prof J Regis for stereotactic exploration of the patient


Florent Renac
email: florent.renac@onera.fr

A robust high-order Lagrange-projection like scheme with large time steps for the isentropic Euler equations

Keywords: Lagrange-projection, discontinuous Galerkin method, explicitimplicit, entropy-satisfying, positivity-preserving, isentropic Euler equations Mathematics Subject Classification (2000) 65M12, 65M60

We present an extension to high-order of a first-order Lagrange-projection like method for the approximation of the Euler equations introduced in Coquel et al. (Math. Comput., 79 (2010), pp. 1493-1533). The method is based on a decomposition between acoustic and transport operators associated to an implicitexplicit time integration, thus relaxing the constraint of acoustic waves on the time step. We propose here to use a discontinuous Galerkin method for the space approximation. Considering the isentropic Euler equations, we derive conditions to keep positivity of the mean value of density and to satisfy a discrete entropy inequality in each element of the mesh at any approximation order in space. These results allow to design limiting procedures to restore these properties at nodal values within elements. Numerical experiments support the conclusions of the analysis and highlight stability and robustness of the present method, while it allows the use of large time steps.

Introduction

In this work, we are interested in the design of a robust and accurate numerical method for the description of flows exhibiting multiple space and time scales that can differ by several orders of magnitude. An example is the cooling system of high-pressure gas turbines in turbomachinery. The flow induced by the blade rows is transonic, while the internal cooling channels contain region of very low speed convection. Another examples are multiphase flows where the acoustic time scales strongly differ between liquid and gas, or the Euler equations in the incompressible limit where the speed of acoustic and transport waves differ notably.

In those applications, one may be interested in the transport phenomena associated to slow waves only. Unfortunately, numerical methods are usually designed for the resolution of all waves and suffer from restriction of the fast waves. In standard explicit shock-capturing methods the time step is limited by the fast waves to ensure stability of the numerical scheme. Moreover, their numerical diffusion is proportional to the speed of the fastest waves. Both aspects impose over-resolution in time and space. As a consequence, spurious pressure oscillations that deteriorate the quality of the solution are usually observed.

The present work is mainly based on a Lagrange-projection like scheme introduced in [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] in the context of a first-order finite volume formulation of the Euler equations. The method uses the Lagrange-projection framework [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] for the splitting of acoustic and transport operators, but no mesh movement is applied. The acoustic operator is solved in Lagrange coordinates, while the projection onto the grid is replaced by the transport operator. The Lagrange step is integrated in time with an implicit backward-Euler scheme in order to relax the time step restriction associated to acoustic waves. An explicit forward Euler method is applied for the transport step in order to accurately describe associated unsteady phenomena. This method was then successfully used to design an asymptotic preserving scheme for the discretization of the Euler equations with source terms [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF].

Besides, the work in [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] circumvent the difficulties in the treatment of nonlinearities associated to the equation of state by using a relaxation approximation [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF][START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF]. The latter method approximates the nonlinear system with a linear or a quasi-linear enlarged system with stiff relaxation source terms. In the limit of instantaneous relaxation, the system is consistent with the original system. Relaxation approximations have been applied to the isentropic Euler equations [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF] and to the full system of gas dynamics [START_REF] Chalons | Navier-Stokes equations with several independent pressure laws and explicit predictor-corrector schemes[END_REF][START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] where only the physical pressure is replaced by a relaxation pressure with its own evolution equation.

One attractive feature of the Lagrange-projection like scheme introduced in [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] is the preservation of convex invariant domains and entropy inequality at the discrete level, while allowing the use of large time steps. The objective of the present work is to extend this method to space discretization of arbitrary order by using a discontinuous Galerkin (DG) method [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF][START_REF] Lesaint | On a finite element method for solving the neutron transport equation[END_REF] which has become very popular for the solution of nonlinear convection dominated flow problems [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws II: general framework[END_REF][START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convectiondominated problems[END_REF][START_REF]ADIGMA -A european initiative on the development of adaptive higher-order variational methods for aerospace applications[END_REF][START_REF] Renac | Fast time implicit-explicit discontinuous Galerkin method for the compressible Navier-Stokes equations[END_REF][START_REF] Renac | Aghora: A high-order DG solver for turbulent flow simulations[END_REF]. These are particular aspects that make the DG method well suited. First, it is possible to use the numerical fluxes derived in [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] and therefore the present method may be viewed as a natural extension to high-order of the related numerical method. Then, the effect of the numerical flux on the quality of the approximation is known to decrease as the polynomial degree p in the DG method increases [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convectiondominated problems[END_REF][START_REF] Qiu | A numerical study for the performance of the Runge-Kutta discontinous Galerkin method based on different numerical fluxes[END_REF][START_REF] Renac | Stationary discrete shock profiles for scalar conservation laws with a discontinuous Galerkin method[END_REF]. This avoids the use of local numerical parameters tuned at each interface of the mesh in order to lower the numerical diffusion induced by the first-order approximation [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF][START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF]. This aspect is essential in our analysis to restore the PDE properties at the discrete level and to derive a priori conditions to preserve invariant domains and satisfy an entropy inequality by the present Lagrange-projection DG (LPDG) scheme. These two latter properties are satisfied for the mean value of the numerical solution in the elements of the mesh and present similarities with the positivity preserving scheme in [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF] and the entropy satisfying scheme in [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF]. Besides, they suggest the application of a posteriori limiters introduced in [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF][START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF] that extend the properties to nodal values whithin elements.

The paper is organized as follows. Section 2 presents the model problem with the system of isentropic Euler equations (section 2.1), its relaxation approximation (section 2.2) and the splitting between acoustic and transport operators (section 2.3). The numerical approach for the high-order space discretization is introduced in section 3, while time discretization is described in section 4. The first-order implicit-explicit time integration is described in section 4.1 and the properties of the numerical scheme are analyzed in section 4.2. High-order time discretization and limiting strategies are discussed in sections 4.3 and 4.4, respectively. These results are assessed by several numerical experiments in section 5. Finally, concluding remarks about this work are given in section 6.

2 One-dimensional model problem

Isentropic Euler equations

The discussion in this paper focuses on the Euler equations for an isentropic gas in one space dimension. Let Ω = R be the space domain and consider the following problem

∂ t u + ∂ x f (u) = 0, in Ω × (0, ∞), (1a)
u(•, 0) = u 0 (•), in Ω. ( 1b
)
The vector

u = ρ ρu (2)
represents the conservative variables with ρ the density and u the velocity. The nonlinear convective fluxes in (1a) are defined by

f (u) = ρu ρu 2 + p . (3)
Equations ( 1) are supplemented with an equation of state for the pressure of the form p = p(τ ), with τ = 1/ρ the specific volume. Assuming that p (τ ) < 0 and p (τ ) > 0 for all τ > 0, the system (1a) is strictly hyperbolic over the set of states

Ω a = {u ∈ R 2 : ρ > 0, u ∈ R},
with eigenvalues u ± c(τ ) associated to nonlinear fields. The sound speed is defined by c 2 (τ ) = τ 2 e (τ ) with e(τ ) the specific internal energy defined by e (τ ) = -p(τ ).

Introducing the specific total energy E = e + u 2 /2, the mapping ρE : Ω a u → ρE(u) ∈ R is a strictly convex function [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]. Physically relevant solutions to (1) must hence satisfy an inequality of the form

∂ t ρE + ∂ x (ρEu + pu) ≤ 0. ( 4
)

Relaxation approximation of the Euler equations

Solutions to the isentropic Euler equations (1a) may be approximated by solutions of the following Suliciu relaxation system [6]

∂ t ρ + ∂ x (ρu) = 0, (5a)
∂ t (ρu) + ∂ x (ρu 2 + Π) = 0, (5b)
∂ t (ρΠ) + ∂ x (ρΠu + a 2 u) = -ρ(Π-p(τ )) ε , (5c)
where ε > 0 represents a characteristic relaxation time. Equation (5c) models the evolution of the pressure relaxation, Π, in flows subject to mechanical disequilibrium. The variable Π may be viewed as a linearization of the pressure p around its equilibrium, Π = p(τ ), while a > 0 is a parameter that approximates the Lagrangian sound speed, ρc, and whose value will be specified later. In the limit of instantaneous relaxation, we get

lim ε→0 Π = p(τ ), (6)
and system (5) converges formally toward (1a). We refer the reader to [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF] for an in-depth discussion of the relaxation approximation of the Euler equations for an isentropic fluid.

It may be checked that the homogeneous quasilinear system (5) is strictly hyperbolic over the set of states

Ω r = {w ∈ R 3 : ρ > 0, u ∈ R, Π ∈ R}, (7)
with eigenvalues µ 1 = u -aτ , µ 2 = u, and µ 3 = u + aτ . All the characteristic fields associated with these eigenvalues are linearly degenerate. We note that the relaxation system ( 5) is a dissipative approximation of the Euler equations providing that the following subcharacteristic condition is satisfied

a > max τ -p (τ ), (8)
for all τ under consideration [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF].

Acoustic-transport operator splitting

We decompose equation (1a) between acoustic and transport operators with a sequential splitting: the acoustic step reads

∂ t u + (∂ x u)u + ∂ x 0 p = 0, (9)
then the transport step reads

∂ t u + ∂ x f t (u) -(∂ x u)u = 0, (10)
with f t (u) = uu. Assuming smooth solutions, the acoustic step (9) may be rewritten into the equivalent form

∂ t τ -τ ∂ x u = 0, ( 11a
)
∂ t u + τ ∂ x p = 0. ( 11b
)
In the following, we shall consider a Suliciu relaxation approximation for the Lagrangian gas dynamics applied to [START_REF] Coquel | Second-order entropy diminishing scheme for the Euler equations[END_REF]. Considering the relaxation system for the Euler equations (5), the relaxation system for the acoustic part reads

∂ t w + τ ∂ x f a (w) = s(w), (12)
with

w =   τ u Π   , f a (w) =   -u Π a 2 u   , s(w) =   0 0 -Π-p(τ ) ε   . ( 13
)
With a slight abuse, the whole vector w will be referred to as the Lagrange variables in the following. It may be shown that the relaxation system ( 12) is a dissipative approximation of the acoustic step [START_REF] Coquel | Second-order entropy diminishing scheme for the Euler equations[END_REF] under the subcharacteristic condition [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws II: general framework[END_REF].

Let m be the mass variable defined by dm = ρ 0 (x)dx. Then, approximating the space derivative operator τ ∂ x by τ (x, 0)∂ x with τ (•, 0) = 1/ρ 0 (•) and imposing ε → ∞, the equations in ( 12) may be written in the following conservation form

∂ t τ -∂ m u = 0, ( 14a
)
∂ t u + ∂ m Π = 0, (14b)
∂ t Π + ∂ m (a 2 u) = 0. ( 14c
)
The following results hold for the homogeneous relaxation system [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF].

Theorem 1 (Hyperbolicity of the relaxation system) The homogeneous relaxation system ( 14) is strictly hyperbolic over the set of states [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] with eigenvalues

µ 1 = -a < µ 2 = 0 < µ 3 = a.
All the characteristic fields associated with these eigenvalues are linearly degenerate. Moreover, the characteristic variables associated with these eigenvalues are w= Π + au, J = Π + a 2 τ , and w= Π -au, respectively:

∂ t w +a∂ m w = 0, ( 15a
)
∂ t J = 0, ( 15b
)
∂ t w -a∂ m w = 0. ( 15c
)
Proof Computing the eigenvalues and eigenvectors of ∂f a /∂w is straightforward and is left to the reader. Likewise, the three operations (14c)+a(14b), (14c)+a 2 (14a), and (14c)-a(14b) give (15a,b,c), respectively.

Theorem 2 (Riemann problem for the relaxation system) Consider the Riemann problem defined by the homogeneous relaxation system (14) associated with the initial condition

w(m, 0) = w L , m < 0, w R , m > 0,
where w L and w R are in Ω r . Then, the unique solution to this problem is the self-similar solution W(•; w L , w R ) defined by

W( m t ; w L , w R ) =        w L , m t < -a, w L , -a < m t < 0, w R , 0 < m t < a, w R , m t > a, (16)
where w L = (τ L , u , Π ) , w R = (τ R , u , Π ) , and

u = u L + u R 2 + Π L -Π R 2a , (17a)
Π = Π L + Π R 2 + a(u L -u R ) 2 , (17b)
τ L = τ L + u -u L a , (17c)
τ R = τ R + u R -u a . ( 17d
)
Proof The Rankine-Hugoniot relations associated with the continuity and momentum equations (14a,b) through the steady 2-wave impose u L = u R = u and Π L = Π R = Π . Now, applying the Rankine-Hugoniot relations across the 1-and 3-waves, one obtains

-a(τ L -τ L ) + (u -u L ) = 0, (18a)
a(τ R -τ R ) + (u R -u ) = 0, ( 18b
) -a(u -u L ) -(Π -Π L ) = 0, (18c)
a(u R -u ) -(Π R -Π ) = 0. ( 18d
)
Adding and subtracting (18c) and (18d), one obtains the expressions for u and Π , respectively. Then, (18a) and (18b) lead to the expressions (17c) and (17d).

We end this section by summing up our strategy to solve problem [START_REF] Berthelin | Kinetic schemes on staggered grids for barotropic Euler models: entropy-stability analysis[END_REF]. We split (1a) into acoustic and transport parts by applying the decomposition introduced in section 2.3: (i) the acoustic step is approximated with the homogeneous relaxation system ( 14); (ii) we apply the transport step [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF]. We stress that relaxation mechanisms in the right-hand-side of ( 12) are taken into account through the initial condition (1b) which is transformed into Lagrange variables with data at equilibrium, i.e., w 0 = 1/ρ 0 , ρu 0 /ρ 0 , p(1/ρ 0 )

pointwise. We now propose to mimic this strategy at the discrete level in the next section, where Theorems 1 and 2 will be used to design the numerical flux for the acoustic step.

Discontinuous Galerkin formulation

The DG method consists in defining a discrete weak formulation of problem [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF], [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF] with initial condition (1b). The domain is discretized with a uniform grid

Ω h = ∪ j∈Z κ j with cells κ j = [x j-1/2 , x j+1/2 ], x j+1/2 = (j + 1
2 )h and h > 0 the space step (see Figure 1).

Numerical solution and Lagrange polynomials

We look for approximate solutions in the function space of discontinuous polynomials where P p (κ j ) denotes the space of polynomials of degree at most p in the element κ j . The approximate solutions to systems ( 10) and ( 14) are sought under the form

V p h = {v h ∈ L 2 (Ω h ) : v h | κ j ∈ P p (κ j ), κ j ∈ Ω h }, (19)
u h (x, t) = p l=0 φ l j (x)U l j (t), ∀x ∈ κ j , κ j ∈ Ω h , t ≥ 0, (20a)
w h (x, t) = p l=0 φ l j (x)W l j (t), ∀x ∈ κ j , κ j ∈ Ω h , t ≥ 0, (20b)
where U l j = (ρ l j , ρU l j ) constitute the degrees of freedom (DOFs) in the element κ j and are associated to conservative variables, while W l j = (τ l j , U l j , Π l j ) are coefficients associated to Lagrange variables. The subset (φ 0 j , . . . , φ p j ) constitutes a basis of V p h restricted onto a given element. In this work we will use the Lagrange interpolation polynomials 0≤k≤p associated to the Gauss-Lobatto nodes over the segment [-1, 1]:

s 0 = -1 < s 1 < • • • < s p = 1: k (s l ) = δ k,l , 0 ≤ k, l ≤ p, (21)
with δ k,l the Kronecker symbol. The basis functions in a given element κ j thus write φ k j (x) = k (σ j (x)) where σ j (x) = 2(x -x j )/h and x j = (x j+1/2 + x j-1/2 )/2 denotes the center of the element.

The DOFs thus correspond to the point values of the solution, e.g. given 0 ≤ k ≤ p, j in Z, and t ≥ 0, we have u h (x k j , t) = U k j (t) for x k j = x j + s k h/2. The left and right traces of the numerical solution at interfaces x j±1/2 of a given element hence read (see figure 1):

u - j+ 1 2 (t) := u h (x - j+ 1 2 , t) = U p j (t), ∀t ≥ 0, (22a)
u + j-1 2 (t) := u h (x + j-1 2 , t) = U 0 j (t), ∀t ≥ 0. (22b)
Moreover, assuming equilibrium ε → 0 for the relaxation system [START_REF] Coquel | An entropy satisfying MUSCL scheme for systems of conservation laws[END_REF], the physical and relaxation pressures satisfy

p τ k j (t) = Π k j (t), 0 ≤ k ≤ p, t > 0. ( 23
)
As a consequence, the conservative and Lagrange variables (20) may be related in a weak sense by the relations

W k j (t) = w(U k j (t)) 0 ≤ k ≤ p, t ≥ 0, (24a)
U k j (t) = u(W k j (t)) 0 ≤ k ≤ p, t ≥ 0, (24b)
where

w : Ω a → Ω r ; u → w(u) = 1 ρ , ρu ρ , p 1 ρ , u : Ω r → Ω a ; w → u(w) = 1 τ , u τ
denote, with a slight abuse, the change from conservative to Lagrange variables and its inverse.

Space discretization

We first consider the space discretization of the relaxation approximation of the acoustic step without source term, i.e., equation [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF]. Substitute (20b) in equation ( 14) with ε → ∞, multiply it with a test function v h in V p h with support in a given element κ j and integrate by parts over κ j to obtain

κ j v h ∂ t w h dx - κ j f a (w h )∂ x (v h τ h )dx + v h τ h h a (w - h , w + h ) x j+ 1 2 x j-1 2 = 0, (25)
after replacing the physical flux at interface by the numerical flux h a : Ω r × Ω r → R 3 defined from the solution of the Riemann problem ( 16):

h a (w - j+ 1 2 , w + j+ 1 2 ) = f a W(0; w - j+ 1 2 , w + j+ 1 2 ) =    -u j+ 1 2 Π j+ 1 2 a 2 u j+ 1 2    , (26)
with data at equilibrium, i.e., Π ±

j+ 1 2 = p(τ ± j+ 1 2
). Note that there is no ambiguity in defining the numerical flux at ξ = 0 since the application ξ → f a W(ξ; w L , w R ) is continuous at the origin because of the Rankine-Hugoniot relations associated with the steady wave. Using a second integration by parts, the semi-discrete equation may be equivalently written as

κ j v h ∂ t w h dx + κ j v h τ h ∂ x f a (w h )dx + v h τ h h a (w - h , w + h ) -f a (w h ) x j+ 1 2 x j-1 2 = 0. ( 27
)
We now introduce the space discretization of the transport problem [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] where we consider the approximate solution (20a). Using a similar approach as for the acoustic step, with one integration by parts, the semi-discrete form of the DG discretization in space reads

κ j v h ∂ t u h dx + κ j (v h u h )∂ x u h dx + v h h t (u - h , u + h ) -h u (u - h , u + h )u h x j+ 1 2 x j-1 2 = 0, (28)
where h t : Ω a × Ω a → R 2 and h u : R × R → R denote Lipschitz-continuous numerical fluxes consistent with the physical fluxes of the transport step: h t (u, u) = uu and h u (u, u) = u. In this work, we use upwind numerical fluxes. Introducing

ûj+ 1 2 = U p j , u j+ 1 2 > 0, U 0 j+1 , u j+ 1 2 ≤ 0, ( 29
)
where the quantity u j+1/2 is defined from the first component of the numerical flux [START_REF] Spiteri | A new class of optimal high-order strong-stability-preserving time discretization methods[END_REF] for the acoustic step, we set

h t (u - j+ 1 2 , u + j+ 1 2 ) = u j+ 1 2 ûj+ 1 2 = (u j+ 1 2 ) + U p j + (u j+ 1 2 ) -U 0 j+1 , (30)
where

(u j+ 1 2 ) + = max(u j+ 1 2 , 0), (u j+ 1 2 ) -= min(u j+ 1 2 , 0)
represent positive and negative parts of u j+ 1

2

. Then, we set

h u (u - j+ 1 2 , u + j+ 1 2 ) = u j+ 1 2 . ( 31
)
The integrals in the semi-discrete equations are approximated by using a numerical quadrature. We use a Gauss-Lobatto quadrature rule with nodes associated with the interpolation points of the numerical solution

κ j f (x)dx h 2 p l=0 ω l f (x l j ),
with ω l > 0, p l=0 ω l = 2, x l j = x j +s l h/2 the weights and nodes of the quadrature rule, and s l defined in [START_REF] Qiu | A numerical study for the performance of the Runge-Kutta discontinous Galerkin method based on different numerical fluxes[END_REF]. For p + 1 integration points, this quadrature is exact for polynomials of degree deg(f ) ≤ 2p -1. It is worth noting that the integration by parts used in equations ( 27) and ( 28) are replaced by summation by parts of the form

f, d x v h p j + d x f, v h p j = f (x p j )v h (x p j ) -f (x 0 j )v h (x 0 j ), ∀v h ∈ V p h , where
f, g p j = h 2 p l=0 ω l f (x l j )g(x l j )
represents the discrete inner product in the element κ j . As noticed in [START_REF] Kopriva | On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods[END_REF], this operation holds true in a weak sense for general functions f by considering the interpolation polynomials of degree p of the function, let say f h , at integration points, x 0≤k≤p j .

Time discretization

In section 4.1, we introduce the fully discrete scheme for a one-step first-order implicit-explicit time discretization and analyze its properties in section 4.2, while we briefly discuss high-order time discretization in section 4.3 and limiting strategy in section 4.4.

First order time discretization

Let t (n) = n∆t with ∆t > 0 the time step, the time integration of the Euler equations with Lagrange-projection over a time step is done with a two-step sequential splitting: (i) homogeneous acoustic step ( 27) over (t (n) , t (n+1 -) ] and (ii) transport step ( 28) over (t (n+1 -) , t (n+1) ]. Again, relaxation mechanisms in [START_REF] Coquel | An entropy satisfying MUSCL scheme for systems of conservation laws[END_REF] are taken into account by imposing data at equilibrium at each time step:

Π (n) h = p(τ (n)
h ). We use an implicit backward-Euler scheme for the time discretization of the acoustic step over (t (n) , t (n+1 -) ]. Approximating the mass variable by ∆m = τ k j (t (n) )h at point x k j over (t (n) , t (n+1 -) ] and setting v h = φ k j into (27), the discrete scheme now reads

W k,n+1 - j = W k,n j -λ k τ k,n j ∂ x f a (w n+1 - h ), φ k j p j + δ k,p h a (W p,n+1 - j , W 0,n+1 - j+1 ) -f a (W p,n+1 - j ) -δ k,0 h a (W p,n+1 - j-1 , W 0,n+1 - j ) -f a (W 0,n+1 - j ) , ( 32
)
for all 0 ≤ k ≤ p and j in Z, where we have used the notations W k,n j = W k j (t (n) ) and λ k = 2λ/ω k with λ = ∆t/h. For the sake of clarity, we now use the short notations for the components of the Riemann solver

h a (W p,n+1 - j , W 0,n+1 - j+1 ) =    -u j+ 1 2 Π j+ 1 2 a 2 u j+ 1 2   
without any possible confusion on the time and trace values.

The subcharacteristic condition ( 8) is imposed at the discrete level by requiring that

a > max j∈Z max 0≤k≤p max θ∈[0,1] -p θτ k,n j + (1 -θ)τ k,n+1 - j , n ∈ N. (33)
Likewise, we use an explicit forward Euler time integration of the semi-discrete transport equations [START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF]. Introducing the definitions of the numerical fluxes (30) and (31), we obtain

U k,n+1 j = U k,n+1 - j -λ k u n+1 - h ∂ x u n+1 - h , φ k j p j + δ k,p u j+ 1 2 ûn+1 - j+ 1 2 -U p,n+1 - j )-δ k,0 u j-1 2 ûn+1 - j-1 2 -U 0,n+1 - j ) . (34)
The discrete problem for the Euler equations now reads: for all time t (n+1) with n ≥ 0 find u h (•, t (n+1) ) in (V p h ) 2 such that equations (32) and (34) are satisfied with

W k,n j = w U k,n j , 0 ≤ k ≤ p, j ∈ Z, (35a)
U k,n+1 - j = u W k,n+1 - j , 0 ≤ k ≤ p, j ∈ Z, (35b)
given by transformations [START_REF] Renac | Aghora: A high-order DG solver for turbulent flow simulations[END_REF].

The time sequence u h (•, t (n) ) is associated with the initial condition

Ω h v h w h (x, 0)dx = Ω h v h w u 0 (x) dx, ∀v h ∈ V p h , (36)
which reduces to W k j (0) = w u 0 (x k j ) for all 0 ≤ k ≤ p and j in Z. The following theorem allows us to give the final form (37) of the LPDG scheme.

Theorem 3 (Consistency and conservation) The discrete schemes (32) and (34) with projections (35) constitute a conservative approximation consistent in time and space with the Euler equations (1a) of the form

U k,n+1 j = U k,n j -λ k -f (u n+1 - h ), d x φ k j p j + δ k,p h n+1 - j+ 1 2 -δ k,0 h n+1 - j-1 2 , (37) with h n+1 - j+ 1 2 =   u j+ 1 2 ρn+1 - j+ 1 2 u j+ 1 2 ρu n+1 - j+ 1 2 + Π j+ 1 2   , (38)
evaluated from DOFs at time t (n+1 -) .

Proof We first observe that the first component of the discrete equation (32) reads

τ k,n+1 - j = L k,n+1 - j τ k,n j with L k,n+1 - j = 1 + λ k ∂ x u n+1 - h , φ k j p j +δ k,p (u j+ 1 2 -U p,n+1 - j ) -δ k,0 (u j-1 2 -U 0,n+1 - j ) (39a) = 1 + λ k -u n+1 - h , d x φ k j p j + δ k,p u j+ 1 2 -δ k,0 u j-1 2 . (39b)
Using the change of variables (24b) and assuming that τ k,n j > 0 and τ k,n+1 - j > 0 (see equation (47) in Lemma 1), the two first components of (32) may thus be rewritten under the form

L k,n+1 - j ρ k,n+1 - j = ρ k,n j , (40a)
L k,n+1 - j ρU k,n+1 - j = ρU k,n j -λ k ∂ x Π n+1 - h , φ k j p j + δ k,p Π j+ 1 2 -Π p,n+1 - j -δ k,0 Π j-1 2 -Π 0,n+1 - j .(40b)
Using the expressions of the numerical fluxes (30) and (31) with f t (U k,n j ) = U k,n j U k,n j , the discrete transport step (34) may be rewritten as

U k,n+1 j = U k,n+1 - j -λ k ∂ x (u n+1 - h u n+1 - h ), φ k j p j -∂ x u n+1 - h u n+1 - h , φ k j p j +δ k,p h t (u - j+ 1 2 , u + j+ 1 2 ) -f t (U p,n+1 - j ) + f t (U p,n+1 - j ) -u j+ 1 2 U p,n+1 - j -δ k,0 h t (u - j-1 2 , u + j-1 2 ) -f t (U 0,n+1 - j ) + f t (U 0,n+1 - j ) -u j-1 2 U 0,n+1 - j = L k,n+1 - j U k,n+1 - j -λ k ∂ x (u n+1 - h u n+1 - h ), φ k j p j +δ k,p h t (u - j+ 1 2 , u + j+ 1 2 ) -f t (U p,n+1 - j ) -δ k,0 h t (u - j-1 2 , u + j-1 2 ) -f t (U 0,n+1 - j ) .
Summing the two above schemes, we get

U k,n+1 j = U k,n j -λ k ∂ x f E (u n+1 - h , w n+1 - h ), φ k j p j +δ k,p h n+1 - j+ 1 2 -f E (U p,n+1 - j , W p,n+1 - j ) -δ k,0 h n+1 - j-1 2 -f E (U 0,n+1 - j , W 0,n+1 - j ) ,
with f E (u, w) = uu+(0, Π) consistent with the Euler fluxes f (u) when Π = p(τ ). This is done by imposing instantaneous relaxation [START_REF] Renac | Fast time implicit-explicit discontinuous Galerkin method for the compressible Navier-Stokes equations[END_REF], p τ k,n+1 j = Π k,n+1 j . Now, we use the equivalence between the discrete DG formulations with either one or two integration by parts when using a collocated Gaus-Lobatto quadrature [START_REF] Kopriva | On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods[END_REF]. Applying integration by parts, the above scheme reduces to (37) and constitutes a conservative implicit-explicit discretization of the Euler equations (1a) consistent in space and time.

Then, the numerical scheme (32) provides the following discrete versions of the conservation equations ( 15) for the characteristic variables. Note that from the definition of the numerical flux ( 26) with [START_REF] Kopriva | On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods[END_REF] we have

Π j+ 1 -Π p,n+1 - j = -a(u j+ 1 2 -U p,n+1 - j ), ( 41a
)
Π j-1 2 -Π 0,n+1 - j = a(u j-1 2 -U 0,n+1 - j ). ( 41b
)
Then, setting

W k,n j = Π k,n j + aU k,n j , J k,n j = Π k,n j + a 2 τ k,n j , W k,n j = Π k,n j -aU k,n j , (42)
and using relations (41), the three operations (32c)+a(32b), (32c)+a 2 (32a) and (32c)-a(32b) give respectively:

W k,n+1 - j = W k,n j -aλ k τ k,n j ∂ x w n+1 - h ,φ k j p j -δ k,0 W p,n+1 - j-1 -W 0,n+1 - j ,(43a) J k,n+1 - j = J k,n j , (43b)
W k,n+1 - j = W k,n j +aλ k τ k,n j ∂ x w n+1 - h ,φ k j p j +δ k,p W 0,n+1 - j+1 -W p,n+1 - j . (43c)
Equation (37) constitutes the LPDG scheme where the state u n+1 - h is evaluated from the linear implicit system (43a,c) that may be easily solved for the discrete characteristic variables (42), τ k,n+1 - j being then given explicitly by (43b). Then, the Lagrange variables are obtained by inverting relations (42) at time t (n+1 -) . This result is essential for the performances of the present method.

Properties of the discrete scheme

In this section, we discuss the properties of the LPDG scheme with first-order time integration and arbitrary order for the space discretization. The main results are given in Theorem 4 and prove positivity and an entropy inequality for the mean value of the numerical solution

u n j := 1 h κ j u h (x, t (n) )dx = p l=0 ω k 2 U l,n j . ( 44
)
The entropy inequality applies to the total energy that we introduce at the discrete level via its interpolant

ρE h (x, t) = p l=0 φ l j (x)ρE l j (t), ∀x ∈ κ j , t 0, ( 45
)
with ρE l j (t) = ρE U l j (t) .

Lemma 1 Assume that ρ 0≤k≤p,n j∈Z > 0, then under the CFL condition

λ max j∈Z max 0≤k≤p 1 ω k u n+1 - h , d x φ k j p j -δ k,p (u j+ 1 2 ) -+ δ k,0 (u j-1 2 ) + < 1 2 , ( 46
)
we have ρ k,n+1 - j > 0 ( 47
)
and

u n+1 j = p k=0 ω k 2 -λ u n+1 - h , d x φ k j p j -δ k,p (u j+ 1 2 ) -+ δ k,0 (u j-1 2 ) + U k,n+1 - j -λ(u j+ 1 2 ) -U 0,n+1 - j+1 + λ(u j-1 2 ) + U p,n+1 - j-1 (48)
is a convex combination of DOFs at time t (n+1 -) .

Proof From (40a), we have

L k,n+1 - j ρ k,n+1 - j = ρ k,n j with L k,n+1 - j > 0 from (39b)
and condition (46), hence ρ k,n+1 - j > 0. Then using (34) to evaluate the mean value at time t (n+1) together with the definition of the upwind flux (29), we get

u n+1 j = p k=0 ω k 2 U k,n+1 j = p k=0 ω k 2 U k,n+1 - j -λ k u n+1 - h ∂ x u n+1 - h , φ k j p j +δ k,p (u j+ 1 2 ) -U 0,n+1 - j+1 -U p,n+1 - j ) -δ k,0 (u j-1 2 ) + U p,n+1 - j-1 -U 0,n+1 - j ) .
Then, developing the expression of the volume integral, the above equation reads

u n+1 j = p k=0 ω k 2 U k,n+1 - j -λ p k=0 ω k h 2 U k,n+1 - j p l=0 d x φ l j (x k j )U l,n+1 - j + δ k,p (u j+ 1 2 ) -(U 0,n+1 - j+1 -U p,n+1 - j ) -δ k,0 (u j-1 2 ) + (U p,n+1 - j-1 -U 0,n+1 - j ) .
Inverting indices in the double sum and rearranging terms, one easily obtains (48) and positivity of coefficients follows from (46). Finally, note that all coefficients in (48) are positive from condition (46) with unit sum from p k=0 d x φ k j = 0. Thus (48) is a convex combination.

Lemma 2 Assume that ρ 0≤k≤p,n j∈Z > 0, then the discrete acoustic step satisfies the following discrete entropy inequality

η k,n+1 - j -η k,n j -2a 2 λ k τ k,n j Π n+1 - h u n+1 - h , d x φ k j p j -δ k,p H n+1 - j+ 1 2 +δ k,0 H n+1 - j-1 2 ≤ 0, (49) with η k,n j = ( W k,n j ) 2 + ( W k,n j ) 2 2 = (Π k,n j ) 2 + a 2 (U k,n j ) 2 , (50a)
H n+1 - j+ 1 2 = ( W p,n+1 - j ) 2 -( W 0,n+1 - j+1 ) 2 4a = Π j+ 1 2 u j+ 1 2 . ( 50b
)
Proof Relations (50) follow directly from the definitions in ( 42) and [START_REF] Kopriva | On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods[END_REF]. Then,

multiplying equation (43a) with W k,n+1 - j gives W k,n+1 - j ( W k,n+1 - j -W k,n j ) + aλ k τ k,n j ω k h 2 W k,n+1 - j ∂ w h ∂x n+1 - x k j -δ k,0 W 0,n+1 - j W p,n+1 - j-1 -W 0,n+1 - j = 0, hence ( W k
,n+1 - j ) 2 2 - ( W k,n j ) 2 2 + aλ k τ k,n j ∂ x ( ( w n+1 - h ) 2 2 ), φ k j p j -δ k,0 ( W p,n+1 - j-1 ) 2 2 - ( W 0,n+1 - j ) 2 2 = - ( W k,n+1 - j -W k,n j ) 2 2 -aλ k τ k,n j δ k,0 ( W p,n+1 - j-1 -W 0,n+1 - j ) 2 2 .
Using integration by parts in the above equation and considering the sign of its right-hand-side, one deduces

( W k,n+1 - j ) 2 2 - ( W k,n j ) 2 2 + aλ k τ k,n j - ( w n+1 - h ) 2 2 , d x φ k j p j + δ k,p ( W p,n+1 - j ) 2 2 -δ k,0 ( W p,n+1 - j-1 ) 2 2 ≤ 0.
Likewise, multiplying equation (43c) with W k,n+1 - j and applying similar manipulations give

( W k,n+1 - j ) 2 2 - ( W k,n j ) 2 2 -aλ k τ k,n j - ( w n+1 - h ) 2 2 , d x φ k j p j + δ k,p ( W 0,n+1 - j+1 ) 2 2 -δ k,0 ( W 0,n+1 - j ) 2 2 ≤ 0.
Summing the two last equations gives the desired inequality (49).

Lemma 3 Assume that ρ 0≤k≤p,n j∈Z > 0, then under the CFL condition (46) and subcharacteristic condition (33), the discrete acoustic step satisfies the following discrete entropy inequality

E k,n+1 - j -E k,n j -λ k τ k,n j Π n+1 - h u n+1 - h , d x φ k j p j -δ k,p H n+1 - j+ 1 2 + δ k,0 H n+1 - j-1 2 ≤ 0, (51) with H n+1 - j± 1 2
given by (50b).

Proof Using (50a) and (45), the specific total energy may be written as

E k,n j = e(τ k,n j ) + (η k,n j -(Π k,n j ) 2 )/2a 2 . Hence, we have E k,n+1 - j -E k,n j = e(τ k,n+1 - j ) -e(τ k,n j ) + η k,n+1 - j -η k,n j 2a 2 - (Π k,n+1 - j -Π k,n j ) 2 2a 2 - Π k,n j (Π k,n+1 - j -Π k,n j ) a 2
.

Using (43b) to substitute Π k,n+1 - j -Π k,n j = -a 2 (τ k,n+1 - j -τ k,n j
) and the fact that data at equilibrium impose Π k,n j = p(τ k,n j ) = -e (τ k,n j ), we have

E k,n+1 - j -E k,n j - η k,n+1 - j -η k,n j 2a 2 = e(τ k,n+1 - j ) -e(τ k,n j ) -e (τ k,n j )(τ k,n+1 - j -τ k,n j ) - a 2 2 (τ k,n+1 - j -τ k,n j ) 2 .
Applying a second-order Taylor development with integral remainder of e(τ k,n+1 - j ) about τ k,n j , we obtain

E k,n+1 - j -E k,n j - η k,n+1 - j -η k,n j 2a 2 = (τ k,n+1 - j -τ k,n j ) 2 × • • • 1 0 e (τ k,n j + ξ(τ k,n+1 - j -τ k,n j )) -a 2 (1 -ξ)dξ ≤ 0,
under the subcharacteristic condition (33). Finally, using (49) gives (51).

Theorem 4 Assume that ρ 0≤k≤p,n j∈Z > 0, then under the CFL condition (46) and subcharacteristic condition (33), the LPDG scheme satisfies positivity for the mean value of the solution:

ρ n+1 j > 0, (52)
and the discrete entropy inequality

ρE(u n+1 j )-ρE n j +λ u j+ 1 2 ( ρE n+1 - j+ 1 2 +Π j+ 1 2 )-u j-1 2 ( ρE n+1 - j-1 2 +Π j-1 2 ) ≤ 0, ( 53
)
where ρE

n+1 - j± 1 2
denote upwind fluxes of the form (29) evaluated at time t (n+1 -) .

Proof From assumptions of Theorem 4, the results of Lemmas 1 and 3 hold. We thus infer positivity in (52) by using the convex combination (48) with (47).

Then, the entropy inequality (53) follows from the following arguments. Multiplying (51) with ρ k,n j = L k,n+1 - j ρ k,n+1 - j and using (39a), we get

ρE k,n+1 - j -ρE k,n j + λ k ρE k,n+1 - j ∂ x u n+1 - h , φ k j p j +δ k,p (u j+ 1 2 -U p,n+1 - j ) -δ k,0 (u j-1 2 -U 0,n+1 - j ) +λ k -Π n+1 - h u n+1 - h , d x φ k j p j + δ k,p u j+ 1 2 Π j+ 1 2 -δ k,0 u j-1 2 Π j-1 2 ≤ 0.
The first volume integral may be rewritten as

ρE k,n+1 - j ∂ x u n+1 - h , φ k j p j = ω k h 2 ρE k,n+1 - j ∂u h ∂x n+1 - x k j = ρE n+1 - h ∂ x u n+1 - h , φ k j p j ,
and using integration by parts one thus obtains

ρE k,n+1 - j -ρE k,n j + λ k -u n+1 - h , ∂ x (ρE n+1 - h φ k j ) + Π n+1 - h d x φ k j p j +δ k,p u j+ 1 2 (ρE p,n+1 - j + Π j+ 1 2 ) -δ k,0 u j-1 2 (ρE 0,n+1 - j + Π j-1 2 ) ≤ 0.
Summing over 0 ≤ k ≤ p with weights ω k /2, one obtains

ρE n+1 - j -ρE n j + λ -u n+1 - h , ∂ x ρE n+1 - h p j +u j+ 1 2 (ρE p,n+1 - j + Π j+ 1 2 ) -u j-1 2 (ρE 0,n+1 - j + Π j-1 2 ) ≤ 0. (54)
Now, by convexity of the mapping ρE(u), the convex combination (48) gives

ρE(u n+1 j ) - p k=0 ω k 2 -λ u n+1 - h , d x φ k j p j -δ k,p (u j+ 1 2 ) - +δ k,0 (u j-1 2 ) + ρE k,n+1 - j + λ(u j+ 1 2 ) -ρE 0,n+1 - j+1 -λ(u j-1 2 ) + ρE p,n+1 - j-1 ≤ 0.
Using integration by parts, we obtain

ρE(u n+1 j ) -ρE n+1 - j + λ p k=0 -∂ x u n+1 - h , φ k j p j -δ k,p ((u j+ 1 2 ) --U p,n+1 - j ) + δ k,0 ((u j-1 2 ) + -U 0,n+1 - j ) ρE k,n+1 - j + λ(u j+ 1 2 ) -ρE 0,n+1 - j+1 -λ(u j-1 2 ) + ρE p,n+1 - j-1 ≤ 0.
Using again integration by parts and the definition of the numerical flux (29) applied to ρE h , one obtains

ρE(u n+1 j ) -ρE n+1 - j + λ u n+1 - h , ∂ x ρE n+1 - h p j +u j+ 1 2 ( ρE n+1 - j+ 1 2 -ρE p,n+1 - j ) -u j-1 2 ( ρE n+1 - j-1 2 -ρE 0,n+1 - j ) ≤ 0. ( 55
)
Summing ( 54) and ( 55) gives (53).

We end this section by noting that the inequality (53) should be written as

ρE(u n+1 j ) -ρE(u n j ) + λ u j+ 1 2 ( ρE n+1 - j+ 1 2 + Π j+ 1 2 ) -u j-1 2 ( ρE n+1 - j-1 2 + Π j-1 2 ) ≤ ρE n j -ρE(u n j ),
and is not strictly speaking an entropy inequality since the quantity ρE n j -ρE(u n j ) is positive from a Jensen's inequality. As a consequence, the Lax-Wendroff theorem cannot be applied to prove convergence to an entropy solution. The control of this quantity is a challenging task (see [START_REF] Berthon | An entropy preserving MOOD scheme for the Euler equations[END_REF] and references therein for an in-depth discussion on discrete entropy inequalities for high-order schemes). In the present work, the limiter (57) is an attempt to circumvent this difficulty.

High-order time discretization

The present method is extended to high-order time integration by using strongstability preserving explicit Runge-Kutta methods [START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes[END_REF][START_REF] Spiteri | A new class of optimal high-order strong-stability-preserving time discretization methods[END_REF]. These methods consist in convex combinations of first-order forward Euler methods and thus will keep positivity of Theorem 4 under a given CFL condition. We note however that the first-order time discretization (37) is not an explicit forward Euler method because the residuals are evaluated at an intermediate time step t (n+1 -) . As a consequence, accuracy in time is not guaranteed when using high-order Runge-Kutta schemes. The design of adapted high-order time integration is beyond the scope of the present study. However, the numerical experiments in section 5 tend to indicate that explicit Runge-Kutta time integration do not alter accuracy of the present method and are thus well adapted in practice at least for the present range of applications.

(γ-1) 2 4γ 1.6 RP2 shock-shock (0.5, 1) (1, 1) (2, 1)
(γ-1) 2 4γ 1.6 RP3 rarefaction-shock (1, -0.5) (0.5, -0.25) (γ-1) 2 4γ 1.6 RP4 rarefaction-rarefaction (1, -5) (2, 1)
the highest discretization order p = 3 as visible in the density distributions of RP1 and RP2. The solution for RP4 is made of two symmetric rarefaction waves with formation of near-vacuum in the intermediate region. The positivity limiter is successful to keep robustness of the computation and increasing p reduces the diffusion at the tail of the waves as expected. However, the entropy limiter slightly alters the solution at the head of the rarefaction waves for p = 3. We attribute this effect to the fact that the solution is not smooth in this region, the limiters keeping accuracy for smooth solutions only. We refer to the work in [START_REF] Zhang | A minimum entropy principle of high order schemes for gas dynamics equations[END_REF] for details on the conditions for preserving accuracy for smooth solutions.

Concluding remarks

The LPDG scheme introduced in this work is based on the Lagrange-projection like scheme from [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] derived in the context of a first-order finite volume formulation of the Euler equations. This method is here extended to high-order by using a DG method of arbitrary order for the space discretization and associated to a first-order implicit-explicit time discretization of acoustic and transport operators, respectively. Considering the isentropic Euler equations, a priori conditions on the time step and on the numerical parameter imposing the subcharacteristic condition are derived in order to guaranty positivity and entropy inequality for the mean value of the numerical solution in each mesh element. A posteriori limiters similar to those introduced in [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF][START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF] are then used to extend these properties to nodal values within elements. Strong-stability preserving Runge-Kutta schemes are applied for the time integration in order to keep positivity at any time discretization order.

Numerical experiments in one space dimension highlight high-order approximation of smooth solutions, while the method proves to be robust in the presence of discontinuities or vacuum. Future investigations will consider the full system of gas dynamics with a general equation of state and the extension to several space dimensions.

Fig. 1

 1 Fig. 1 Mesh with definition of left and right traces at interfaces x j±1/2 .

Fig. 2

 2 Fig. 2 RP1: numerical solution for density (left) and velocity (right) at time t = 0.3 for different polynomial degrees p = 1 to p = 3 from top to bottom, h = 1 200 .

Fig. 3

 3 Fig. 3 RP2: numerical solution for density (left) and velocity (right) at time t = 0.3 for different polynomial degrees p = 1 to p = 3 from top to bottom, h = 1 200 .

Fig. 4 Fig. 5

 45 Fig. 4 RP3: numerical solution for density (left) and velocity (right) at time t = 0.4 for different polynomial degrees p = 1 to p = 3 from top to bottom, h = 1 200 .

Table 1

 1 Manufactured smooth solution: different norms of the error at time t = 5 and associated orders of convergence.

	p h	e h L 1 (Ω)	O 1	e h L 2 (Ω)	O 2	e h L∞(Ω)	O∞
		1/4	0.35511e-01	-	0.45109e-01	-	0.85997e-01	-
		1/8	0.97747e-02 1.86 0.12770e-01 1.82 0.31154e-01 1.46
	1	1/16 0.24729e-02 1.98 0.32942e-02 1.95 0.78816e-02 1.98
		1/32 0.61003e-03 2.01 0.81665e-03 2.01 0.19587e-02 2.00
		1/64 0.14369e-03 2.08 0.19544e-03 2.06 0.45383e-03 2.10
		1/4	0.50087e-02	-	0.76828e-02	-	0.26541e-01	-
		1/8	0.42909e-03 3.54 0.65902e-03 3.54 0.24946e-02 3.41
	2	1/16 0.33870e-04 3.66 0.53700e-04 3.61 0.21400e-03 3.54
		1/32 0.50879e-05 2.73 0.79338e-05 2.75 0.25079e-04 3.09
		1/64 0.50163e-06 3.34 0.75007e-06 3.40 0.23383e-05 3.42
		1/4	0.39559e-03	-	0.55401e-03	-	0.19667e-02	-
		1/8	0.23759e-04 4.05 0.36194e-04 3.93 0.13136e-03 3.90
	3	1/16 0.15068e-05 3.97 0.22986e-05 3.97 0.94855e-05 3.79
		1/32 0.94001e-07 4.00 0.14440e-06 3.99 0.61378e-06 3.94
		1/64 0.58836e-08 3.99 0.90435e-08 3.99 0.38706e-07 3.98

Table 2

 2 Initial conditions and physical parameters of Riemann problems

	test	description	left state u L	right state u R	κ	γ
	RP1 shock-shock

Acknowledgements The author would like to thank Frédéric Coquel from École Polytechnique and Christophe Chalons from Université de Versailles-Saint-Quentin-en-Yvelines for valuable discussions and their constructive comments.

Limiting strategy

The properties in Theorem 4 hold only for the mean value in mesh elements of the numerical solution at time t (n+1) , which is not sufficient for robustness and stability of numerical computations. However, these results may motivate the use of a posteriori limiters introduced in [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF][START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF]. These limiters aim at extending preservation of invariant domains [START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF] or maximum-principle [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF] from mean to nodal values within elements. Our strategy differs slightly from the ones in [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF][START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF], so we describe it in the following.

First, we enforce positivity of nodal values of density by using the linear limiter ρk,n+1

with 0 ≤ θ ρ j ≤ 1 defined by ) and strengthen the entropy inequality (53) by observing that the discrete transport step (34) satisfies a maximum principle for any convex function U : Ω a u → U (u) ∈ R. Indeed, using (48) we obtain

we thus impose a maximum principle at nodal values from Ũk,n+1

with 0 ≤ θ s j ≤ 1 defined by

Note that when Ȗk,n+1 j is not in Ω a , there exists a unique 0 ≤ θ s,k j ≤ 1 such that the above relation holds by convexity of U since u n+1 j is in Ω a . In practice, we use the total energy as entropy. Finally, we replace the DOFs at time t (n+1) by the limited values Ũ0≤k≤p,n+1 j∈Z . For high-order time integration, we apply the limiter after each stage of the Runge-Kutta scheme. We stress that the limiters (56) and (57) preserve conservation and accuracy for smooth solutions [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF][START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF].

We end this section by summing up our strategy at the discrete level with the following algorithm applied at each stage of the Runge-Kutta method:

(i) solve the linear system (43) for the characteristic variables (42) with data at equilibrium (35a); (ii) compute the conservative variables with (42); (iii) compute the discrete residuals of the LPDG scheme (37) with these values; (iv) apply the limiters (56) and (57).

Numerical experiments

In this section we present several numerical experiments to illustrate the performances of the LPDG scheme derived in this work. For all experiments, we consider an isentropic polytropic ideal gas with an equation of state of the form p(τ ) = κτ -γ with κ > 0 and γ > 1.

We use strong-stability preserving Runge-Kutta time integration schemes of order p + 1 when using polynomials of degree p for the space discretization: the two-stage second-order Heun method for p = 1, the three-stage third-order scheme of Shu-Osher [START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes[END_REF] for p = 2, and the five-stage fourth-order scheme of Spiteri and Ruuth [START_REF] Spiteri | A new class of optimal high-order strong-stability-preserving time discretization methods[END_REF] for p = 3, respectively.

Finally, the a priori CFL condition (46) and subcharacteristic condition (33) are imposed at time t (n) as was proposed in [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF][START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF].

Manufactured smooth solution

We first consider the convection of a density wave in a uniform flow with Mach number M ∞ . Let Ω = (0, 1), we solve

with periodicity conditions and initial condition ρ 0 (x) = 1 + sin(2πx), u 0 (x) = 1, ∀x ∈ Ω, with = 0.2. The source term s = s(x, t) is such that the exact solution for this problem reads

The parameters of the equation of state are κ = 1/γM 2 ∞ with M ∞ = 0.1 and γ = 1.4. Table 1 indicates different norms of the numerical error on density e h = ρ h -ρ for different polynomial degrees and grid refinements with associated convergence orders in space. The expected p + 1 order of convergence is recovered with the present method.

Riemann problems

We now consider Riemann problems with initial condition

The set of initial conditions is given in Table 2. Problems RP1 and RP3 are taken from [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF], while problem RP4 is taken from [START_REF] Berthelin | Kinetic schemes on staggered grids for barotropic Euler models: entropy-stability analysis[END_REF]. Figures 2 to 5 compare the numerical solution in symbols with the exact solution in lines.

For RP1, RP2 and RP3, results are qualitatively similar. The shock waves are well captured and the increase in the discretization order has a clear positive effect on the approximation of the rarefaction waves in RP3. We observe some spurious oscillations of low amplitude in the neighborhood of strong shocks with