Maxime Pinson

Pablo Caron

Philippe Laurent

Ion Cojocari

Development of a plastic scintillator-based active shield for the ICARE-NG radiation monitor

Keywords: Monte-Carlo simulation, active shield, space physics, geometric factor, decontamination, veto counters, space environment measurements

An active shield using a scintillator and silicon photo-multipliers (SiPMs) has been developed to operate with space environment particle detectors sensitive to both protons and electrons, such as ICARE-NG (Influence sur les Composants Avancés des Radiations de l'Espace-Nouvelle Génération). The method shows a reduction in electron contamination through the sides of the detector, thus increasing energy resolution. Two geometries are studied, one working in coincidence mode with the main detector, the other in anti-coincidence. Performances of both geometries are simulated in proton and electron environments in energy ranges typical of a space environment. Experimental measurements then aimed to validate the Monte-Carlo simulation framework for scintillating materials. Lastly, a case study of electron decontamination is carried out, as well as an error rate estimation in flux reconstruction.

I. INTRODUCTION

M ANY instruments have been flown in the past to measure the proton and electron population in the Earth's radiation belts, as well as galactic cosmic rays. To this end, several detectors such as the Relativistic Electron Proton Telescope (REPT) flown on the Radiation Belt Storm Probes (RBSP) in 2012, have been devised. The REPT uses a stack of silicon diodes in an aluminum casing [START_REF] Baker | The Relativistic Electron-Proton Telescope (REPT) Investigation: Design, Operational Properties, and Science Highlights[END_REF].

Also, the RADMON detector [START_REF] Oleynik | Calibration of RADMON Radiation Monitor Onboard Aalto-1 CubeSat[END_REF] launched in 2017 uses a combination of a silicon detector and a CsI(Tl) scintillator for the detection of trapped protons and electrons.

Similarly, the ICARE-NG detector is composed of five silicon diode detectors arranged in three different casings labelled PE2S, PE1 and PE2. Fig. 1 shows a cutaway of a schematic version of one of the three heads of the ICARE-NG detector, head PE1. The diodes vary in thickness from 300 µm to 700 µm allowing for different energy detection ranges. The instrument has flown on several space missions, including SAC-C (2000), JASON-2 (2008) [START_REF] Boscher | In Flight Measurements of Radiation Environment on Board the French Satellite JASON-2[END_REF], SAC-D (2011) [START_REF] Boscher | In-Flight Measurements of Radiation Environment on Board the Argentinean Satellite SAC-D[END_REF] and JASON-3 (2016).

Bigger instruments such as the Alpha Magnetic Spectrometer (AMS) installed on the ISS (International Space Station) in 2011, uses silicon strip detectors to reconstruct particle trajectories in its Silicon Tracker System [START_REF] Burger | The Alpha Magnetic Spectrometer Silicon Tracker[END_REF]. The AMS also M. Pinson is a PhD student with ONERA-DPHY Toulouse, France (e-mail: maxime.pinson@onera.fr), P. Caron is a research engineer at ONERA DPHY Toulouse, France (e-mail: pablo.caron@onera.fr), P. Laurent is a research engineer with CEA-IRFU Gif-sur-Yvette France (e-mail: philippe.laurent@cea.fr), and I. Cojocari is a PhD student with APC CEA Paris, France (e-mail: ion.cojocari@cea.fr) has an anti-coincidence counter (ACC) [START_REF] Doetinchem | The AMS-02 Anticoincidence Counter[END_REF], designed to vetoout particles using plastic scintillators, wave-length shifting (WLS) fibers and photo-multiplier tubes (PMTs).

The solution studied here was to not use WLS fibers, and the small size of the detector led us to choose Silicon Photomultipliers (SiPMs) over PMTs due to their smaller size and growing commercial availability. Also, recent studies seem to highlight the higher photon detection efficiency of SiPMs over PMTs [START_REF] Hahn | Results of the Parallel Operation of Large-Size SiPM Detector Modules and PMTs in IACTs[END_REF] making them a good candidate for our active shield.

Instruments such as REPT, RADMON, NGRM [START_REF] Desorgher | ESA Next Generation Radiation Monitor[END_REF] and ICARE-NG use multi-channel analysers to sample energy deposits in the sensitive volumes (silicon diodes for the most part). The energy responses of these instruments depend greatly on the type of environment it is placed in, and is a superposition of all the types of particles encountered, mainly electrons and protons. A major hurdle in treating instrument data is the ability to separate the electron and proton contributions and convert the counts in the detector to proton and electron fluxes. One of the objectives of this paper is to try and overcome the difficulties associated with the counts to flux conversion in proton and electron environments. The above cited instruments can all be subject to contamination by misidentified particles as well as high energy particles entering through the lateral aluminum shielding. The following work focuses on one geometry, that of ICARE-NG, to showcase the capabilities of an active shield and calculate the expected performance gain of the detector's measurements. This was made possible by coding end-toend Monte Carlo simulations of the detector and scintillators.

Additionally, experimental studies carried out at ONERA were able to provide validation of the Monte-Carlo code used in this paper.

The most recent version of ICARE-NG (flown on Eutelsat E7C) is sensitive to protons and electrons in the typical kinetic energy ranges : 250 keV to 3 MeV for electrons and 30 to 350 MeV for protons. Specifically, the PE1 head of the instrument is composed of a single silicon diode with a thickness of 700 µm and a radius of 10 mm, encased in a 5 mm thick cylindrical aluminum case. This case has a small 3 mm diameter hole on the front. The following results focus on Monte-Carlo simulations of head PE1 alone because it is the only head containing a single diode, and has a small opening on the front. The other heads contain two diodes. This means that it can be subject to both the electron and proton environments simultaneously. Adding a scintillator and SiPMs can prove to be more cost efficient than adding a second diode for coincidence or anti-coincidence counting.

The diode detectors used are sensitive to an electrical charge created by the incident particle by ionization in the silicon. This electrical charge is proportional to the energy deposited by the incident particle. The ICARE-NG diode acquisition chain sorts the measured electrical charge into channels numbered from 0 to 127, where 127 corresponds to the highest deposited energy: approximately 3 MeV (or 23.5 keV per channel). The exact value of this energy is calculated during pre-flight calibrations.

One of the limits of the ICARE-NG instrument is the high energy electron contamination that occurs as these particles pass through the lateral aluminum shielding of the detector casing. These particles hit the diodes at grazing angles, whereas those passing through the front of the detector hit the diode almost perpendicularly, making their track lengths in the detector longer. This means that for a given incident energy, a particle entering through the sides is likely to deposit more energy than through the front of the detector. This leads to a higher variability in the detector response and reduces the energy resolution of the instrument.

To illustrate this, one can use the geometric factor (Gef), shown in Fig. 2, which is a quantity that expresses the proportion of particles that hit the diode detector out of all the particles simulated at that kinetic energy. The geometric factor is characteristic of the implemented geometry and of the nature of the incident particle, as it translates the effect of the geometry on particle collection in the diode. The geometric factor, G (in cm 2 sr), can be defined as in [START_REF] Pak | A Numerical Method to Analyze Geometric Factors of a Space Particle Detector Relative to Omnidirectional Proton and Electron Fluxes[END_REF] for a given incident kinetic energy by

G = n N 4π 2 r 2 (1)
where n number of detected particles by the diode N number of simulated particles at a given incident energy r radius of the sphere from which particles are initialized, in cm The shapes of these geometric factors can be explained by the sub-atomic physics used in the GEANT4 simulation code. One notices that the proton geometric factor of Fig. 2 (top) is representative of the fact that the linear energy transfer (LET) is inversely proportional to its kinetic energy. The spread of the geometric factor around the main LET curve is due to protons arriving at non-perpendicular angles in the diode, which have a significantly longer track length in the silicon, depositing proportionally more energy. The geometric factor associated to the electrons, which are relativistic in this energy range, deposit an approximately constant amount of energy in the silicon diode. This minimum ionizing particle (MIP) behaviour which protons and electrons exhibit at high energies does not allow for the measurement of the kinetic energy by using the deposited energy in the diode.

Furthermore, electrons are subject to greater multiple scattering angles compared to protons, which are inversely proportional to the momentum of the particles. This greater scattering angle explains the difficulty in mitigating electron contamination in the detector. A low energy limit can also be seen in the electron Gef in Fig. 2 (bottom) due to the limited range of electrons at this energy and the aluminum casing of the detector head.

During a mission, information is gathered in the form of counts in the different channels (0 to 127). These counts can then be used to reconstruct the flux of the radiation belts using the geometric factor of the instrument. An active shield would be able to modify the geometric factors of ICARE-NG in a way that would reduce the spread of the deposited energy, and thus improve its accuracy. In turn, these counts would give more precise information on the energies of the particles in the radiation belts, through the creation of differential channels for protons. This information would allow for the improvement of numerical models of the Earth's radiation belts which rely (partly of fully) on experimental data such as GREEN (Global Radiation Earth ENvironment) [START_REF] Sicard | GREEN: the new Global Radiation Earth ENvironment model (beta version)[END_REF]. The geometric factors presented in Fig. 2 are shown as they are currently, without an active shield. By reducing the spread of these geometric factors, we can hope to dissociate both of these contributions.

II. SIMULATION SETUP

The simulations carried out here aimed to quantify the added-value of an active shield on head PE1 of the ICARE-NG instrument. These simulations used the GEANT4 simulation toolkit [START_REF] Agostinelli | Geant4-a simulation toolkit[END_REF], as well as the GODDeSS [START_REF] Dietz-Laursonn | GODDeSS: a Geant4 extension for easy modelling of optical detector components[END_REF] and G4SiPM [START_REF] Niggemann | G4SiPM: A novel silicon photomultiplier simulation package for Geant4[END_REF] packages.

GEANT4 is a Monte-Carlo based code which relies on statistical and probabilistic methods to simulate the behaviour of sub-atomic particles in materials. The track of each particle is divided into a large number of steps, points at which random numbers are drawn by the simulation, which dictate the type of interaction that occurs. Within the GEANT4 framework, the GODDeSS package simulates the behavior of optical photon production for several different types of scintillators, and the G4SiPM package simulates the electrical response of an SiPM given optical photon stimulation. Along with the basic GEANT4 classes, these extra packages allow for an easier implementation of scintillators, fibers and SiPMs allowing for an end to end simulation of the active shield.

The simulated elements are the plastic scintillator (a BC-408 from Saint Gobain) and the Hamamatsu S10362-11-100U silicon photo-multiplier multi-pixel photon counter (MPPC). The BC-408 scintillator has a peak emission wavelength at 425 nm when ionized [START_REF]BC-408, BC-412, BC-416, solid plastic and liquid scintillation materials[END_REF], with a yield reported to be around 10 000 photons/MeV [START_REF] Ogawara | Signal pulse emulation for scintillation detectors using Geant4 Monte Carlo with light tracking simulation[END_REF], [START_REF] Riggi | Geant4 simulation of plastic scintillator strips with embedded optical fibers for a prototype of tomographic system[END_REF]. The Hamamatsu MPPC is composed of 100 square pixels with a width of 100 µm, which have a peak sensitivity (photon detection efficiency, PDE) of 65 % at around 440 nm. Both the fibers and the scintillator have relaxation times of a few nanoseconds, making them compatible for use in an anti-coincidence counter.

The electron and proton particle tracks are initialized isotropically from a sphere which encompasses the entire geometry in order to simulate a first approximation of a simple space environment.

III. PROPOSED GEOMETRIES OF THE ACTIVE SHIELDS

A. Scintillator and SiPM geometries

Two different active shield geometries were simulated in both electron and proton environments (Fig. 3). The first consists of a hollow cylindrical scintillator placed on the outside of head PE1 of the ICARE-NG instrument. The thickness of the diode was set to 300 µm. The radius of this cylinder is roughly that of the detector (1.5 cm), with a thickness of 0.5 mm. The scintillator must also be covered with a light-tight wrapping made of reflective material such as aluminum (teflon and Vikuiti can also be used) as to trap the photons within the material. Two SiPMs measuring 1 mm by 1 mm are glued to the surface of the scintillator, inside the wrapping.

The second geometry consists of a smaller, square piece of scintillator placed in front of the hole on head PE1. This piece is 3 mm wide and 0.5 mm thick. Similarly, two SiPMs are glued to the sides of the scintillator. This second geometry differs from the first in the sense that it works in coincidence mode, keeping only the particles that hit both the detector and the scintillator.

Several notable differences were found between these two solutions, both having some advantages and drawbacks, mainly because of their very different form-factors and the different counting modes used on each (coincidence/anticoincidence).

B. Use of wave-length shifting (WLS) fibers

Several questions were raised as to whether or not the implementation of WLS fibers inside the scintillators would increase the overall detection efficiency. The role of the fibers is to guide the scintillation photons from the scintillator to the SiPMs which can thus be placed farther away from the scintillator. The fibers are encrusted into the plastic scintillator and one end of the fiber would face the sensitive surface of the SiPMs. This can prove to be useful in order to shield the SiPMs from radiation exposure. In both the cylindrical and the square geometries, the fibers were protruding from the scintillator material so that the SiPMs were 2 cm away from the scintillator. Results show that the use of WLS fibers have a negative effect on photon collection (Fig. 4), compared to a situation where the SiPMs are glued to the surface of the scintillator, as in Fig. 3. The number of photons collected on both SiPMs is reduced by a factor of about 8 for all incident proton energies. This reduced efficiency can be explained by the very small trapping efficiency of WLS fibers. On a more practical level, encrusting fibers inside plastic is quite costly and requires more equipment to make the setup light-tight.

Only a small proportion of the wavelength-shifted light is emitted at an angle that allows for total reflection inside the fiber. This value is said to be around 4 % depending on the type of fiber, for a numerical aperture of 0.74 [START_REF] Saint-Gobain | Scintillating fibers Saint-Gobain Crystals[END_REF]. Although the efficiency of the fibers is low, they still allow for a sufficient number of photons to be collected on the SiPMs for detection. This is especially interesting if the SiPMs need to be placed far from the instrument for electrical reasons and/or radiation exposure.

IV. SIMULATION RESULTS AND DISCUSSION

The simulations were carried out in proton and electron environments separately in the energy ranges 10 MeV to 350 MeV and 500 keV to 7 MeV respectively. The geometric factors were calculated with and without the implementation of the active shield and serve as comparison.

To simulate the anti-coincidence counter, the optical photons are simulated from their creation in the scintillator all the way to their absorption by the SiPM. The electrical response of the SiPM is then simulated and if the peak of the electrical response is above 250 mV, the particle is considered as detected by the active shield. The 250 mV threshold was chosen as to avoid counting single cell pulses due to afterpulsing, cross-talk or thermal noise of the SiPM. The voltage range in which an SiPM is operational is quite narrow and is prone to drift with temperature and dose. In our simulations, the SiPMs were operated at 73 V, at an over-voltage of 1.5 V.

A. Cylindrical geometry

Final results of the cylindrical geometry show that an active shield can reduce the variability of the detector response and create differential modes for several different incident energy ranges for protons. Differential modes are groups of channels which are centered on a certain incident energy, giving information capable of reconstructing a differential flux. These differential modes start at around 50 MeV and can be seen up to approximately 100 MeV. Above this energy, the contributions from higher energy protons are too close to each other to obtain a clear separation. This can be seen on the geometric factor on Fig. 5 as the curve flattens at high energies. Examples of the obtained differential modes are shown on Fig. 6. These curves are horizontal cuts of the geometric factors for certain groups of channels of width 20. Above 200 MeV, integrated modes can be used to determine the contribution of protons with an energy above a certain threshold.

As for the electron environment with this geometry, results show only a minimal change in the number of electrons detected, Fig. 7, and thus does not solve the problem of electron contamination above channel number 20. The results suggest that many electrons are deflected at much greater angles by multiple scattering than the protons in the front and rear sides of the aluminum cover, as mentioned in I. In those cases, the scintillator fails to veto-out these electrons and the high energy contamination remains. These deflections mean that the angle of attack on the diode is not as significantly reduced as for protons. For electrons, the sharp increase in the geometric factor around 1.5 MeV (Fig. 2) is attributed to the greater penetration capabilities of higher energy electrons.

To eliminate this contribution, we need to be able to remove all the counts that pass through the aluminum next to the hole on the front, as well as on the back of the cylinder.

Placing another piece of scintillator on the back side of the detector head would not work because it would mean that we would need to distinguish particles using their time of flight to know which ones pass through the front (diode first and then scintillator) and the back (scintillator first and then diode), this aspect is not studied here. Another geometry has to be considered in order to eliminate the remaining contributions. The second (square) geometry aims to solve the problem of the electron environment while maintaining good performance for protons. The differential channels observed with the cylindrical geometry are conserved, and the associated Gef is not included due to its great resemblance to the cylindrical one (Fig. 5).

B. Square geometry

This geometry succeeds in reducing the number of electrons that pass through the aluminum covers, both front and back, due to its smaller form factor, Fig. 8. As a result, the increase in the geometric factor with incident energy is very limited and a significantly lower plateau can be observed at high energies. A significant advantage of this scintillator geometry is its ability to limit the electron response to lower channels. The square scintillator allows for decontamination above channel 20, which can be used for differential modes in the proton environment, as mentioned above.

V. EXPERIMENTAL VALIDATION OF THE MONTE-CARLO SIMULATIONS

A. Experimental setup

The experimental measurements carried out aimed to quantify the global performance of the Monte-Carlo simulations, for different types of radiation. The setup was devised to be able to characterize the energy spectrum of a 90 Sr + 90 Y electron source. This source has a 2.2 MeV spectrum endpoint, requiring a long irradiation time to measure, which was not possible in our setup.

The experiments were carried out in a vacuum chamber and the read-out was ensured by a pulse shaper and an oscilloscope. The scintillating material was custom made to a thickness of 500 µm with a photon emission spectrum close to that of the BC-408. The SiPM used was a Hammamatsu 13360-6050CS. The scintillator was wrapped with one layer of Vikuiti reflector with optical glue. Several layers of Teflon are then added to ensure a light-tight wrapping. Only the SiPM wiring protrudes from the wrapping. Along the line of sight of the electron source, holes were made in the Teflon and Vikuti in both the front and back, to allow for the electrons to pass through the scintillating material and impinge the silicon diode as well. Otherwise, the range of the electrons was not sufficient for them to pass through the scintillator completely. Because of these holes, it is important for these experiments to be carried out in the total darkness, since the scintillator is no longer fully light-tight. The distance from the electron source to the scintillator was measured to be approximately 6.6 mm and scintillator to diode length was 4.3 mm and stayed the same through out the different experiments. Fig. 9. Schematic of the experimental set-up with an electron source. The entire experiment (except for the read-out) was placed in a vacuum. The electron source was placed in such a way that it did not have a direct line of sight on the diode: it is entirely shadowed by the scintillator.

The silicon diode used was a Canberra PIPS detector of thickness 300 µm, the diode output was captured by a bias-T circuit which allowed for the polarization of the diode to -90 V all while sending the signal to a high-bandwidth current amplifier which delivers short current pulses (from 5 to 50 ns). This reversed polarization configuration allows for full depletion of the diode, maximizing the pulse height. The response was then amplified in a Canberra charge amplifier and read on an 8 GHz oscilloscope.

As for the SiPM, the over-voltage was set to 53 V and the signal was passed through an amplifier and then through a pulse shaper, which returns a short square signal of 100 ns. The activity of the electron source, which can be approximated by a point source, was 2.86 kBq. Being 6.6 mm away from the scintillator, the flux of the electron source was not high enough to cause any sort of pile-up in the SiPM nor diode.

The oscilloscope used has a sampling rate of 20 GS/s and is able to trigger based on pre-determined rules involving two different channels, allowing for a direct implementation of the coincidence counter. To implement a coincidence counter the 90 Sr+ 90 Y electron source was one of the only sources which could be used at ONERA since it is the only source that has a penetration length capable of passing through the entire depth of the scintillator. The other radioactive sources, all of which were alpha emitters, did not have sufficient ranges to impinge both the scintillator and the diode, rendering coincidence counting impossible. The experimental setup was implemented in GEANT4 using the same tool kits (GODDeSS and G4SiPM) as for the active shields. The energy spectrum gathered by the SiPM shows reasonable agreement with the predictions made by GEANT4 on Fig. 10. The experimental curves cut off at low energy (around 75 keV) due to the trigger height set by the oscilloscope.

In the same experiment, it was possible to implement a coincidence count of electrons passing through the scintillator and plotting the energy spectrum from the silicon diode. The simulated counts also show good agreement with the experimental data. To show the efficiency of the square geometry in decontaminating the high-energy electron region, experimental data (Fig. 12) from the RBSP-A MAGEIS and RPS instruments is convoluted with the previously simulated geometric factors (both corrected and uncorrected) for protons and electrons, in Figs. 5,8. This data measures proton flux at 105.56 MeV and electron flux between 1.371 MeV and 1.895 MeV. The result of the convolution is the number of counts per second detected by the simulated geometries had they been subjected to these fluxes of protons and electrons.

In the uncorrected case, the proton and electron fluxes cross around channel 95, which means that the channels below 95

B. Error reduction with differential channels in a proton environment

The differential channels created by the cylindrical geometry are more precise than those of the instrument without an active shield. This can be seen in Fig. 6 where the channels 40 to 60 are skewed towards the high energies in the uncorrected case. The cylindrical geometry reduces this problem since the Gef is centered more tightly and symmetrically around 95 MeV.

The method used here after aims to reconstruct a flux using the geometric factors of the corrected and uncorrected cases separately through a closure test. Using a chosen differential flux, we calculate the count rates in each channel by using the corresponding geometric factors. Then, using Gaussian approximations of the differential channels such as in Fig. 6, the flux can be partially reconstructed wherever these Gaussian approximations can be made. Finally, theses fluxes are compared to the original input flux.

This method uses a simplistic mathematical approach, which does not yield accurate results for the current geometries. Instead, more sophisticated methods are considered, such as singular value decomposition (SVD) [START_REF] Maget | Unfolding JASON-2/ICARE-NG High-Energy Particles Measurements Using a Singular Value Decomposition Approach[END_REF], minimization [START_REF] Matzke | Unfolding procedures[END_REF] and artificial intelligence methods [START_REF] Aminalragia-Giamini | Artificial intelligence unfolding for space radiation monitor data[END_REF]. The fact that our geometric factors are cleaner than those of today's detectors allows us to justify the use of this simplistic inversion method.

For proton fluxes alone, the uncorrected case can be subjected to a higher error rate than the cylindrical geometry when calculating differential counts. If the proton flux were to be relatively flat in the 50 to 350 MeV range, then the channels 40 to 60 would be contaminated by the higher energy protons, due to the tail of geometric factor distribution in the uncorrected case in Fig. 6.

To show this, the proton flux f p in M eV -1 cm -2 sr -1 s -1 was modeled by an exponential power law such that:

f p (E) ∝ E α e -E τ (2)
with parameters α (unitless) and τ in M eV . This form was chosen because the (α, τ ) parameter space covers a large variety of flux shapes.

The count rates C can be obtained by convoluting the flux f p to the geometric factor matrices gef p of the uncorrected and corrected cases (Figs. 2 and5) by:

C(α, τ ) = ∞ 0 f p (E)gef p dE (3)
as done in [START_REF] Maget | Unfolding JASON-2/ICARE-NG High-Energy Particles Measurements Using a Singular Value Decomposition Approach[END_REF] and [START_REF] Sandberg | Unfolding and Validation of SREM Fluxes[END_REF]. From these counts, one can isolate certain groups of channels in order to obtain differential fluxes centered on different incident energies, as demonstrated in Fig. 6. The reconstructed flux is then compared to the initial flux f p . Furthermore, the uncorrected geometry tends to consistently overestimate the input flux, this is because the differential channels associated to this geometry are skewed towards higher energies (Fig. 6) leading to higher energy particles being misinterpreted by the reconstruction as lower energy particles.

Given the variability of the proton flux shape in space, Fig. 15 represents the ratio of the mean relative error of the reconstructed flux compared to the original flux associated to a parameter pair (α, τ ), which scans a wide range of flux shapes.

VII. CONCLUSIONS

To conclude, this study has shown the type of performance improvements that we can expect from an active shield composed of a scintillator for the ICARE-NG space instrument. Such a shield on ICARE-NG would be able to reduce the angle of the particles hitting the diode, and can significantly improve the performance of the detector, as shown by the creation of differential modes for protons, and the electron decontamination above channel 20. Both of the geometries studied (cylindrical and square) manage to greatly modify and accentuate the geometric factors of the proton and electron environments. The main differences between the two geometries is the square geometry's ability to filter out high energy electrons while the cylindrical one does not. On the other hand, the square geometry reduces the total number of counts due to its much smaller form factor. Also, the experimental measurements have managed to validate the Monte-Carlo code that was used to simulate the generation and propagation of photons in scintillating materials.

Fig. 1 .

 1 Fig. 1. Cutaway of the schematic geometry of the ICARE-NG detector head PE1 alone, on which the scintillator active shield is studied. A 3 mm radius hole is seen on the front of the aluminum casing. The rear side, which does not contain a hole, has the same thickness as the front : 5 mm.

Fig. 2 .

 2 Fig. 2. Geometric factor (Gef) of head PE1 of ICARE-NG for protons (top) and electrons (bottom) in the energy range 10 MeV to 350 MeV and 0.5 MeV to 7 MeV respectively, without an active shield.

Fig. 3 .

 3 Fig. 3. The square (left) and cylindrical (right) geometries studied here. In blue, head PE1 of the ICARE-NG detector, and in red the scintillator material. The SiPMs can be seen in yellow glued to the scintillator material.

Fig. 4 .

 4 Fig. 4. Number of photons collected per event on the surface of the SiPM as a function of the number of photons created in the scintillator for 350 MeV protons. A factor of 8 between the collection efficiencies is noted.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Corrected geometric factor for the cylindrical scintillator geometry in a simulated proton environment with incident energies from 25 to 350 MeV

Fig. 7 .

 7 Fig.7. Geometric factors as a function of incident energy for all channels from 9 to 127, before and after implementation of the square and cylindrical geometries. The rise in the geometric factor with electron energy is reduced by a factor of around 100 across the entire spectrum with the square geometry whereas the cylindrical geometry fails to reduce any high energy contamination.

Fig. 8 .

 8 Fig. 8. Corrected geometric factor of the square scintillator geometry in an electron environment

B. 90

 90 Sr + 90 Y electron source measurements

Fig. 10 .

 10 Fig. 10. Comparison of SiPM response of the simulated spectrum for the 90 Sr + 90 Y electron source and experimental data

Fig. 11 .

 11 Fig. 11. Coincidence counts detected on diode with the 90 Sr + 90 Y electron source only

Fig. 12 .

 12 Fig. 12. Electron fluxes (top) and proton fluxes (bottom) in cm -2 s -1 sr -1 M eV -1 for different L* values in 2015 over a 6 month period (February 2015 -July 2015) as measured by the RBSP-A MAGEIS and RPS instruments.

Fig. 13 .

 13 Fig. 13. Proton and electron counts in the different channels for corrected and uncorrected cases, for fluxes measured on 19/03/2015 by the RBSP-A MAGEIS and RPS instruments respectively. The channel to energy conversion is the same as for the previous geometries : 23.5 keV/channel.

Fig. 14 .

 14 Fig. 14. Example of 5 different input fluxes with parameters τ = 125 M eV and different values of α. The reconstructed fluxes for differential channels between 100 and 250 MeV show good agreement with the corrected geometry whereas the uncorrected geometry significantly overshoots the real flux.

Fig. 14

 14 Fig.14shows an example of five different input fluxes and their reconstructions after using the corrected and uncorrected versions of the geometric factors of the cylindrical geometry. Each cross and plus sign represents a group of channels associated to a differential channel. The corrected geometry shows much better agreement with the input fluxes, which shows that differential flux reconstruction is easier with the implementation of an active shield.Furthermore, the uncorrected geometry tends to consistently overestimate the input flux, this is because the differential channels associated to this geometry are skewed towards higher energies (Fig.6) leading to higher energy particles being misinterpreted by the reconstruction as lower energy particles.Given the variability of the proton flux shape in space, Fig.15represents the ratio of the mean relative error of the reconstructed flux compared to the original flux associated to a parameter pair (α, τ ), which scans a wide range of flux shapes.

Fig. 15 .

 15 Fig.[START_REF] Ogawara | Signal pulse emulation for scintillation detectors using Geant4 Monte Carlo with light tracking simulation[END_REF]. Ratio of the mean relative error of five differential channels between the uncorrected and corrected cases for fluxes in the (α,τ ) parameter space. In all fluxes shown here, the corrected cases error rate is lower than the uncorrected case.


Alexandre Ern
email: ern@cermics.enpc.fr

Martin Vohral Ík
email: vohralik@ann.jussieu.fr

A POSTERIORI ERROR ESTIMATION BASED ON POTENTIAL AND FLUX RECONSTRUCTION FOR THE HEAT EQUATION *

Keywords: heat equation, unified framework, a posteriori estimate, discontinuous Galerkin, finite volumes, mixed finite elements, conforming finite elements, nonconforming finite elements AMS subject classifications. 65N15, 65N30, 76S05

We derive a posteriori error estimates for the discretization of the heat equation in a unified and fully discrete setting comprising the discontinuous Galerkin, finite volume, mixed finite element, and conforming and nonconforming finite element methods in space and the backward Euler scheme in time. Our estimates are based on a H 1 -conforming reconstruction of the potential, continuous and piecewise affine in time, and a locally conservative H(div)-conforming reconstruction of the flux, piecewise constant in time. They yield a guaranteed and fully computable upper bound on the error measured in the energy norm augmented by a dual norm of the time derivative. Localin-time lower bounds are also derived; for nonconforming methods on time-varying meshes, the lower bounds require a mild parabolic-type constraint on the meshsize.

1. Introduction. We consider the heat equation

∂ t u -∆u = f
a.e. in Q := Ω × (0, T ), (1.1a) u = 0 a.e. on ∂Ω × (0, T ), (

u(•, 0) = u 0 a.e. in Ω, (

where Ω ⊂ R d , d ≥ 2, is a polyhedral domain, T the finite simulation time, f the source term, and u 0 the initial datum. We assume that f ∈ L 2 (Q) and u 0 ∈ L 2 (Ω).

In the sequel, u is called the potential and -∇u the flux. The purpose of this work is to derive guaranteed (that is, without undetermined constants) and fully computable a posteriori error estimates for the discretization of (1.1a)-(1.1c) by locally conservative methods in space. We consider full discretizations obtained using a backward Euler scheme in time and allow for time-varying meshes. The main focus is on nonconforming methods in space, such as discontinuous Galerkin, cell-centered and face-centered finite volumes, and mixed finite elements. Our framework also covers conforming, locally conservative methods such as vertexcentered finite volumes, and under mild modifications, conforming and nonconforming finite elements.

Following the approach proposed by Verfürth for conforming finite elements [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF], the error is measured in a (broken) energy norm augmented by a dual norm of the time derivative. This yields error upper bounds that are global in space and in time, and error lower bounds that are local in time and global in space. The estimators themselves are local in space and in time and can be used in a space-time adaptive time-marching algorithm.

The present estimates are based on defining a H 1 -conforming reconstruction of the potential, continuous and piecewise affine in time, and a locally conservative H(div)conforming reconstruction of the flux, piecewise constant in time. A salient feature of this approach is that it allows for a unified setting: the a posteriori error analysis is performed under two simple conditions on the potential and flux reconstructions without any specific reference to the underlying discretization scheme in space. Given a certain scheme, it suffices to verify these two conditions to apply the present analysis. One condition exploits the local conservativity of the scheme through the flux reconstruction, while the other links locally the mean values of the potential reconstruction to those of the discrete solution. The potential reconstruction is not needed for conforming methods (e.g., vertex-centered finite volumes), while the flux reconstruction is not needed for flux-conforming methods (e.g., cell-centered finite volumes and mixed finite elements).

The present parabolic potential and flux reconstructions are inspired from those derived in the context of a posteriori error estimates for elliptic problems in [START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF][START_REF] Ern | Guaranteed and robust discontinuous Galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF][START_REF] Vohralík | A posteriori error estimates for lowest-order mixed finite element discretizations of convection-diffusion-reaction equations[END_REF][START_REF]Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients[END_REF][START_REF]Residual flux-based a posteriori error estimates for finite volume and related locally conservative methods[END_REF][START_REF]Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods[END_REF]. The idea to derive parabolic a posteriori error estimates from elliptic estimates on each time level is rather natural. In fact, the residual-based a posteriori error estimates for conforming finite elements derived in [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF][START_REF] Bergam | A posteriori analysis of the finite element discretization of some parabolic equations[END_REF] take this form. We also mention [START_REF] Makridakis | Elliptic reconstruction and a posteriori error estimates for parabolic problems[END_REF][START_REF] Lakkis | Elliptic reconstruction and a posteriori error estimates for fully discrete linear parabolic problems[END_REF][START_REF] De Frutos | A posteriori error estimates for fully discrete nonlinear parabolic problems[END_REF] for the so-called elliptic reconstruction technique allowing for optimal error estimates in higher order norms for conforming finite elements. An important conceptual difference is that we reconstruct the (vector-valued) flux and that this quantity is discrete, is constructed locally by postprocessing, and is directly used to evaluate the estimator. In [START_REF] Babuška | A posteriori error estimation for the semidiscrete finite element method of parabolic differential equations[END_REF][START_REF] Babuška | On one approach to a posteriori error estimates for evolution problems solved by the method of lines[END_REF][START_REF] Strouboulis | Guaranteed a posteriori error estimation for fully discrete solutions of parabolic problems[END_REF], various estimators for elliptic problems are extended to the heat equation in a conforming setting to bound the error measured in the L 2 -norm in the space-time cylinder plus a time-weighted energy-norm; only error upper bounds are considered. Some other results on a posteriori error estimates for parabolic problems in a conforming setting can also be found in [START_REF] Repin | Estimates of deviations from exact solutions of initial-boundary value problem for the heat equation[END_REF] using the so-called functional approach where a flux reconstruction is also considered, but without enforcing any local condition; furthermore, only error upper bounds are derived. Finally, we observe that contrary to conforming finite elements, a posteriori energy-norm error estimates for the heat equation discretized by nonconforming methods are less explored; we mention, in particular, [START_REF] Cascón | Space-time adaptive algorithm for the mixed parabolic problem[END_REF] for mixed finite elements, [START_REF] Nicaise | A posteriori error estimates for a nonconforming finite element discretization of the heat equation[END_REF] for nonconforming finite elements, [START_REF] Georgoulis | A posteriori error control for discontinuous Galerkin methods for parabolic problems[END_REF] for discontinuous Galerkin methods, and [START_REF] Amara | A posteriori error estimator for finite volume schemes[END_REF] for finite volume schemes.

This paper is organized as follows. Section 2 presents the continuous and discrete settings. Section 3 collects the main results, namely the error upper and lower bounds. A space-time adaptive time-marching algorithm is also briefly outlined. Section 4 shows how to apply the present framework to various discretization schemes in space, namely discontinuous Galerkin, mixed finite elements, and various finite volume schemes. Sections 5 and 6 are devoted to the proofs of the error upper and lower bounds, respectively. Finally, Appendix A extends the theory to conforming and nonconforming finite elements.

2. The setting. This section briefly describes the continuous and discrete settings.

The continuous setting

. Since f ∈ L 2 (Q) and u 0 ∈ L 2 (Ω), the exact solution is such that u ∈ X := L 2 (0, T ; H 1 0 (Ω)) with ∂ t u ∈ X ′ = L 2 (0, T, H -1 (Ω)
). For a.e. t ∈ (0, T ) and for all v ∈ H 1 0 (Ω), there holds

∂ t u, v (t) + (∇u, ∇v)(t) = (f, v)(t), (2.1)
x where •, • denotes the duality pairing between H 1 0 (Ω) and H -1 (Ω) and (•, •) the L 2 (Ω)-inner product with associated norm denoted by • . In the sequel, the space H 1 0 (Ω) is equipped with the H 1 -seminorm, and for a region R ⊂ Ω, • R denotes the L 2 (R)-norm with appropriate Lebesgue measure.

Ω t t 0 t 1 t 2 t 3 τ 1 τ 2 τ 3 T 0 T 1 T 2 T 3 u 0 hτ u 1
For y ∈ X, we introduce the space-time energy norm

y 2 X := T 0 ∇y 2 (t) dt. (2.2)
Furthermore, for y ∈ Y := {y ∈ X; ∂ t y ∈ X ′ }, we also introduce the space-time norm

y Y := y X + ∂ t y X ′ , ∂ t y X ′ := T 0 ∂ t y 2 H -1 (t) dt 1/2 .
Observe that ∂ t y X ′ = sup ϕ∈X; ϕ X =1 T 0 ∂ t y, ϕ (t) dt. 2.2. The discrete setting. This section collects some useful notation concerning the discrete setting.

2.2.1. Time steps and time-varying meshes. We consider a strictly increasing sequence of discrete times {t n } 0≤n≤N such that t 0 = 0 and t N = T , together with a set of meshes {T n } 0≤n≤N . The discrete times and the meshes are constructed by a space-time adaptive time-marching algorithm, e.g., that outlined in §3.3 below.

For all 1 ≤ n ≤ N , we define I n := [t n-1 , t n ] and τ n := t nt n-1 . For all 0 ≤ n ≤ N , we assume that T n covers exactly the polyhedral domain Ω. For all T ∈ T n , h T denotes the diameter of T , and we let h n := max T ∈T n h T denote the maximum meshsize of T n . For simplicity, we also assume that the meshes are simplicial and matching. Extensions to general polygonal and nonmatching meshes are possible, but technical; see, e.g., [START_REF] Ern | Flux reconstruction and a posteriori error estimation for discontinuous Galerkin methods on general nonmatching grids[END_REF] for an example of flux reconstruction on such meshes. Furthermore, the initial mesh T 0 is used to approximate the initial condition, while for all 1 ≤ n ≤ N , T n corresponds to the mesh used to march in time from t n-1 to t n , see Figure 2.1. The meshes can be refined or coarsened as time evolves; precise assumptions on the meshes are stated in §3 below. Typically, T n is obtained from T n-1 by refining some elements and coarsening some other ones. For all 1 ≤ n ≤ N , we denote by T n-1,n a common refinement of T n-1 and T n . Let W be a vector space of functions defined on Ω. Then, P 1 τ (W ) denotes the vector space of functions defined on Q such that v(•, t) takes values in W and is continuous and piecewise affine in time. Functions in P 1 τ (W ) are uniquely defined by the (N + 1) functions {v n := v(•, t n )} 0≤n≤N in W . Similarly, P 0 τ (W ) denotes the vector space of functions defined on Q such that v(•, t) takes values in W and is piecewise constant in time; for 1 ≤ n ≤ N , we then set v n := v(•, t)| In . Functions in P 0 τ (W ) are uniquely defined by the N functions {v n } 1≤n≤N in W . Furthermore, we observe that if v ∈ P 1 τ (W ), then

∂ t v ∈ P 0 τ (W ) is such that for all 1 ≤ n ≤ N , ∂ t v n := ∂ t v| In = 1 τ n (v n -v n-1 ).
2.2.2. The discrete solution. For all 0 ≤ n ≤ N , the approximate solution at time t n , say u n hτ , is such that

u n hτ ∈ V n h where V n h := V h (T n
) is a discrete space built on the mesh T n . The spaces V n h consist of piecewise polynomial functions whose degree is uniformly bounded in n. In the sequel, Π n 0 denotes the L 2 -orthogonal projection onto piecewise constant functions on T n , while Π V n h denotes the L 2 -orthogonal projection onto V n h . We introduce the space-time function u hτ : Q → R which is continuous and piecewise affine in time and such that for all 1 ≤ n ≤ N and t ∈ I n ,

u hτ (•, t) := (1 -̺)u n-1 hτ + ̺u n hτ , ̺ = 1 τ n (t -t n-1 ).
More generally, for any function v :

Q → R that is continuous in time, we set v n := v(•, t n ) : Ω → R for all 0 ≤ n ≤ N .
2.2.3. Broken gradients and broken X-norm. Since we allow for nonconforming methods in space, it is convenient to introduce for all 0 ≤ n ≤ N the broken gradient operator ∇ n such that for a function v that is smooth within each mesh element in T n , ∇ n v ∈ [L 2 (Ω)] d is defined as (∇ n v)| T := ∇(v| T ) for all T ∈ T n . The broken gradient operator ∇ n-1,n on the mesh T n-1,n is defined similarly. Because of possible nonconformities in space approximation, the discrete solution u hτ may not be in the energy space X. Thus, we extend the definition (2.2) of the X-norm by setting for all y ∈ X, yu hτ 2

X := N n=1 In ∇ n-1,n (y -u hτ ) 2 (t) dt = N n=1 In T ∈T n-1,n ∇(y -u hτ ) 2 T (t) dt.
Since ∂ t (yu hτ ) X ′ is always well-defined, the quantity yu hτ Y is now welldefined for all y ∈ Y . In the sequel, y is either the exact solution or the potential reconstruction.

2.2.4. Mesh faces. For all 0 ≤ n ≤ N , the mesh faces in T n are collected into the set F n . More specifically, F ∈ F n if F has positive (d -1)-dimensional measure and if either there are distinct mesh elements T ± ∈ T n (arbitrary but fixed once and for all) such that F = ∂T -∩ ∂T + (F is then called an interior face and we write F ∈ F i,n ) or if there is a mesh element T ∈ T n such that F = ∂T ∩ ∂Ω (F is then called a boundary face and we write F ∈ F b,n ). For T ∈ T n , it is convenient to introduce the sets

F n T := {F ∈ F n ; F ⊂ ∂T }, F n T := {F ∈ F n ; F ∩ ∂T = ∅},
that is, F n T collects the faces of T whereas F n T collects the faces having a non-empty intersection with ∂T . We will also use the set F i,n T := F n T ∩ F i,n which collects the interior faces of T .

For F ∈ F i,n and a smooth-enough function v that is possibly two-valued on F , we define its average and jump at F as

{ {v} } := 1 2 (v -+ v + ), [[v]] := v --v + ,
where v ± = v| T ± , and we define n F as the unit normal to F pointing from T -towards T + . For F ∈ F b,n , the above definitions are extended by setting

{ {v} } = [[v]] = v| F ,
while n F coincides with the unit outward normal to Ω. For a subset F ⊂ F n , we define the jump seminorms

|[[v]]| ± 1 2 ,F := F ∈F h ±1 F [[v]] 2 F 1/2
, where h F denotes the diameter of F . Finally, for F ∈ F i,n , we define the jump in the normal derivative of v as n

•[[∇ n v]] := n F •(∇v| T --∇v| T + ).
3. Main results. This section collects the main results of this paper concerning the error upper and lower bounds. We also outline a space-time adaptive timemarching algorithm to be used in conjunction with the present estimates.

Error upper bound.

The approximation error uu hτ will be measured in the Y -norm, while the error upper bound will be formulated in terms of a potential reconstruction s and a flux reconstruction θ.

3.1.1. Assumptions on the potential and flux reconstructions. We assume that

s ∈ P 1 τ (H 1 0 (Ω)), θ ∈ P 0 τ (H(div, Ω)). (3.1)
The potential reconstruction is determined by the (N + 1) functions s n ∈ H 1 0 (Ω) associated with the discrete times {t n } 0≤n≤N , while the flux reconstruction is determined by the N functions θ n ∈ H(div, Ω) associated with the time intervals {I n } 1≤n≤N .

The potential and flux reconstructions must satisfy two important assumptions. Firstly, the mean values on mesh elements of the potential reconstruction s are related to those of the discrete solution u hτ . Specifically, we assume that for all 0 ≤ n ≤ N , (s n , 1)

T ′ = (u n hτ , 1) T ′ , ∀T ′ ∈ T n,n+1 , (3.2)
with the convention that T N,N +1 := T N . An important consequence of (3.2) is the following Lemma 3.1 (Mean values of time derivatives). For all 1 ≤ n ≤ N , there holds

(∂ t s n , 1) T = (∂ t u n hτ , 1) T , ∀T ∈ T n . (3.3)
Proof. Observe first that by definition, for all 1 ≤ n ≤ N ,

∂ t (s -u hτ ) n = 1 τ n [(s n -u n hτ ) -(s n-1 -u n-1 hτ )].
Furthermore, by construction, s n and u n hτ have the same mean values on all the elements of T n and of T n+1 . Similarly, s n-1 and u n-1 hτ have the same mean values on all the elements of T n-1 and of T n . Hence, (s ns n-1 ) and (u n hτu n-1 hτ ) have the same mean values on all the elements of T n , yielding (3.3).

Secondly, the flux reconstruction θ must satisfy the following local conservation property: For all 1 ≤ n ≤ N ,

( f n -∂ t u n hτ -∇•θ n , 1) T = 0, ∀T ∈ T n . (3.4)
Here, we have set

f n := 1 τ n In f (•, t) dt (it is also possible to take f n := f n if f ∈ C 0 (0, T ; L 2 (Ω))). We define accordingly f ∈ P 0 τ (L 2 (Ω)) such that f | In := f n for all 1 ≤ n ≤ N .
The actual design of the flux reconstruction θ exploits the local conservation properties of the numerical scheme in space; examples are given in §4 below.

The error estimators.

For all 1 ≤ n ≤ N and T ∈ T n , we define the residual estimator and the diffusive flux estimator respectively as

η n R,T := C P h T f n -∂ t s n -∇•θ n T , (3.5)
η n DF,T (t) := ∇s(t) + θ n T , t ∈ I n , (3.6)
where C P := 1 π . Observe that only the quantity η n DF,T is time-dependent. Furthermore, still for all 1 ≤ n ≤ N , we define the following nonconformity estimators: For all T ∈ T n ,

η n NC1,T (t) := ∇ n-1,n (s -u hτ )(t) T , t ∈ I n , (3.7)
η n NC2,T := C P h T ∂ t (s -u hτ ) n T . (3.8)
Observe that only the quantity η n NC1,T is time-dependent. We also observe that the four above estimators are local in space and in time. Finally, we define the initial condition estimator as

η IC := 2 1/2 s 0 -u 0 .
(3.9)

3.1.3. Guaranteed and fully computable upper bound. We are now in a position to state our main result concerning the error upper bound.

Theorem 3.2 (Error upper bound). Assume (3.1), (3.2), and

(3.4). Let η n R,T , η n DF,T , η n NC1,T , η n NC2,T
, and η IC be defined by (3.5)-(3.9). Then,

u -u hτ Y ≤ 3 N n=1 In T ∈T n (η n R,T + η n DF,T (t)) 2 dt 1/2 + η IC + 3 f -f X ′ + N n=1 In T ∈T n (η n NC1,T ) 2 (t) dt 1/2 + N n=1 τ n T ∈T n (η n NC2,T ) 2 1/2 .
Remark 3.3 (Conforming methods). For conforming methods, u n hτ ∈ H 1 0 (Ω), and we put s n := u n hτ . Hence, the nonconformity estimators η n NC1,T and η n NC2,T vanish. Remark 3.4 (Flux-conforming methods). The cell-centered finite volume method of §4.2 below, as well as the lowest-order mixed finite element method of §4.3 below, are flux-conforming methods that directly produce u n hτ such that -∇ n u n hτ =: θ n ∈ H(div, Ω) for all 1 ≤ n ≤ N . Consequently, we can drop in such a case the flux reconstruction θ n and only work with u n hτ . Remark 3.5 (Time oscillation of the source term). The quantity ff X ′ can be viewed as an error estimator related to the time-oscillation of the source term.

It is convenient to state a slightly less sharp error estimate by separating the residual and diffusive flux contributions. To shorten the notation, we also introduce the local-in-time and global-in-space version of the four above estimators. We define for all 1

≤ n ≤ N , (η n R ) 2 := 2τ n T ∈T n (η n R,T ) 2 = 2τ n T ∈T n C 2 P h 2 T f n -∂ t s n -∇•θ n 2 T , (3.10)
(η n DF ) 2 := 2 In T ∈T n (η n DF,T (t)) 2 dt = 2 In T ∈T n ∇s(t) + θ n 2 T dt. (3.11) Moreover, let (η n NC,1 ) 2 := In T ∈T n (η n NC1,T (t)) 2 dt = In T ∈T n ∇ n-1,n (s -u hτ )(t) 2 T dt, (3.12) (η n NC,2 ) 2 := τ n T ∈T n (η n NC2,T ) 2 = τ n T ∈T n C 2 P h 2 T ∂ t (s -u hτ ) n 2 T . (3.13)
A straightforward consequence of Theorem 3.2 is the following Corollary 3.6 (Simplified upper bound). Under the assumptions of Theorem 3.2, there holds

u -u hτ Y ≤ 3 N n=1 (η n R ) 2 + (η n DF ) 2 1/2 + N n=1 (η n NC,1 ) 2 1/2 + N n=1 (η n NC,2 ) 2 1/2 + η IC + 3 f -f X ′ .
Remark 3.7 (Spatial and temporal estimators). The diffusive flux estimator η n DF can be further separated into a time and a space contribution. Namely, the triangle inequality yields

(η n DF ) 2 ≤ (η n DF,1 ) 2 + (η n DF,2 ) 2 ,
where

(η n DF,1 ) 2 := 4 In T ∈T n ∇s(t) -∇s n 2 T dt = 4 3 τ n T ∈T n ∇(s n -s n-1 ) 2 T , (η n DF,2 ) 2 := 4τ n T ∈T n ∇s n + θ n 2 T ,
where we have used the fact that s is piecewise affine in time to simplify the expression for η n DF,1 . Then, in an adaptive time-marching algorithm, η n DF,1 can be used as a time error estimator, while η n DF,2 , together with η n R , η n NC,1 , and η n NC,2 can be used as space error estimators. Furthermore, since s-u hτ is piecewise affine in time, the evaluation of η n NC,1 can be simplified into a sum involving the quantities ∇ m (s mu m hτ ) for m ∈ {n -1, n}; see Lemma 6.1 below.

Error lower bound.

The goal of this section is to derive local-in-time upper bounds on the error estimators appearing in Corollary 3.6 in terms of the error uu hτ , the data, and possibly some jump seminorms of the discrete solution u hτ . These bounds are derived for a specific choice of the potential reconstruction s based on a modified Oswald interpolate and a specific approximation property for the flux reconstruction θ.

3.2.1. Assumptions on the space-time meshes. We assume that (M1) Shape regularity: the meshes {T n } 0≤n≤N are shape regular uniformly in n. (M2) Transition condition: for all 1 ≤ n ≤ N , the commonly refined meshes T n-1,n are also shape regular uniformly in n and such that

Ξ := sup 1≤n≤N sup T ∈T n-1 ∪T n sup T ′ ∈T n-1,n ;T ′ ⊂T h T h T ′ < +∞. (3.14) (M3) For all 1 ≤ n ≤ N , (h n ) 2 ≤ Υτ n .
In practice, all the meshes {T n } 0≤n≤N are refinements of a base simplicial mesh T * and their shape regularity parameter is bounded by that of T * . Moreover, for two (open) elements T ∈ T n and T ′ ∈ T m , either T ∩ T ′ = ∅, or T ⊂ T ′ , or T ′ ⊂ T . For instance, the newest vertex bisection procedure [START_REF] Schmidt | ALBERT-software for scientific computations and applications[END_REF][START_REF] Morin | Revised reprint of "Data oscillation and convergence of adaptive FEM[END_REF][START_REF] Stevenson | The completion of locally refined simplicial partitions created by bisection[END_REF] maintains these properties.

Remark 3.8 (Transition condition). The transition condition (3.14) on the meshes is classical in the context of a posteriori error analysis with time-varying meshes and means that meshes cannot be refined or coarsened too quickly; more specifically, T ∈ T n-1 in the second supremum restricts the refinement, while T ∈ T n restricts the coarsening.

Henceforth, C denotes a generic constant whose value can change at each occurrence, and which can depend on the regularity of the meshes, the transition constant Ξ in (3.14), the constant Υ in (M3), and the maximum polynomial degree used to build the approximation spaces V n h , but is independent of the size of the meshes and of the time steps. The inequality A ≤ CB is often abbreviated as A B.

Design of the potential reconstruction. Let 0 ≤ n ≤ N . The Oswald interpolation operator I n

Os on the mesh T n is classically constructed as the H 1 -conforming Lagrange interpolate by prescribing at interpolation nodes averaged values of the piecewise discontinuous function to interpolate; see [START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems[END_REF] for the h-analysis and [START_REF] Burman | Continuous interior penalty hp-finite element methods for advection and advection-diffusion equations[END_REF] for the hp-analysis. Observing that for all 0 ≤ n ≤ N , I n Os (u n hτ ) ∈ H 1 0 (Ω), we set

s n := I n Os (u n hτ ) + T ′ ∈T n,n+1 α n T ′ b T ′ , (3.15)
where for all T ′ ∈ T n,n+1 , b T ′ denotes the standard (time-independent) bubble function supported on T ′ , defined as the product of the barycentric coordinates of T ′ , and scaled so that its maximal value is 1. Moreover, for all 0 ≤ n ≤ N and for all T ′ ∈ T n,n+1 , we set

α n T ′ := 1 (b T ′ , 1) T ′ (u n hτ -I n Os (u n hτ ), 1) T ′ . (3.16)
This important choice guarantees that (3.2) holds.

Approximation property of the flux reconstruction.

We assume that for all 1 ≤ n ≤ N and for all T ∈ T n ,

∇u n hτ + θ n T |n•[[∇ n u n hτ ]]| + 1 2 ,F i,n T + |[[u n hτ ]]| -1 2 ,F n T .
(3.17)

3.2.4. Optimally efficient lower bound. For all 1 ≤ n ≤ N , we introduce the following jump semi-norm (local-in-time and global-in-space)

J n (u hτ ) 2 := 1 2 τ n T ∈T n-1 |[[u n-1 hτ ]]| 2 -1 2 ,F n-1 T + T ∈T n |[[u n hτ ]]| 2 -1 2 ,F n T . (3.18)
We also localize in time the X-and Y -norms as follows:

y 2 X(In) := In ∇ n-1,n y 2 (t) dt, y Y (In) := y X(In) + ∂ t y X ′ (In) , with z 2 X ′ (In) := In z 2 H -1 (t) dt.
Finally, we define for all 1 ≤ n ≤ N the spacetime data oscillation term

(E n f ) 2 := f -f 2 X ′ (In) + τ n T ∈T n h 2 T f n -Π V n h f n 2 T . ( 3
η n R + η n DF + η n NC,1 + η n NC,2 u -u hτ Y (In) + J n (u hτ ) + E n f . (3.20)
Remark 3.10 (Bound on the jumps). Further handling of the jump seminorm J n (u hτ ) depends on the numerical scheme. For conforming methods, this jump actually vanishes. Furthermore, the jump seminorm is bounded by the energy error if the jumps have zero mean on each face of the meshes [1, Theorem 10]; this is the case for mixed finite elements, face-centered finite volumes, and nonconforming finite elements. Finally, for discontinuous Galerkin methods and cell-centered finite volumes, one way to proceed is to add the jump seminorm to the error measure in order to have the same error measure in the upper and lower bounds. Alternatively, one may expect to bound the jump seminorm by the energy error, as done in [START_REF] Ainsworth | A posteriori error estimation for discontinuous Galerkin finite element approximation[END_REF] in the elliptic case.

Remark 3.11 (Condition (M3)). Condition (M3) is only needed for nonconforming methods. For such methods on fixed meshes, the condition can still be avoided provided an additional jump seminorm

τ n T ∈T n |[[∂ t u n hτ ]]| 2 + 1 2 ,F n T 1/2
is added to the right-hand side of (3.20).

A space-time adaptive time-marching algorithm.

To briefly outline a space-time adaptive time-marching procedure to be used in conjunction with the above analysis, let Sol(IC, TS, Mesh) denote the discrete solution produced by the numerical scheme starting from an initial condition IC over a single time step TS and using a mesh Mesh. The initial condition need not be defined on the mesh Mesh; the way this is handled depends on the numerical scheme.

1. Initialization (a) choose an initial mesh T 0 with corresponding space V 0 h := V (T 0 ) and an initial approximation u 0 hτ ∈ V 0 h such that the initial condition estimator η IC is below a prescribed tolerance; (b) select an initial time step τ 0 ; 2. Loop in time: While i τ i < T , (a) set T n * := T n-1 with corresponding space V n * h and set τ n * := τ n-1 ; (b) solve u n * hτ := Sol(u n-1 hτ , τ n * , T n * ); (c) estimate error using the upper bound of Theorem 3.2 at the current discrete time, possibly modify T n * and τ n * and return to step (2b); (d) when the error estimate is below a prescribed tolerance, save approximate solution, mesh, and time step as u n hτ , T n , and τ n . Remark 3.12 (Time-marching evaluation of s). When evaluating the error at the current discrete time t n , the mesh T n+1 is not known yet. However, it suffices to adjust the mean values of s n only at the elements of the current mesh T n (putting temporarily T n+1 := T n ). Then, at the next time step t n+1 , if the mesh is refined, additional bubble functions need to be added to s n , as prescribed by (3.15)-(3.16), before evaluating the error at t n+1 . It is also possible to prescribe (3.2) for all T ′ ∈ T n-1,n ; then the time-marching evaluation of s is straightforward, but the error estimates have to be slightly modified.

4.

Applications. We apply in this section the preceding results to different space discretization schemes. For the upper bound of Theorem 3.2, this is done by specifying the reconstructed flux θ and checking that the conservation property (3.4) holds true. For the lower bound of Theorem 3.9, we need to check the approximation property (3.17). The reconstructed potential s in nonconforming methods is always given by (3.15)-(3.16). In conforming methods, where u n hτ ∈ H 1 0 (Ω) for all 0 ≤ n ≤ N , we put s n := u n hτ , so that (3.2) immediately holds true. The reconstructed flux θ n on all time levels 1 ≤ n ≤ N will belong to the Raviart-Thomas-Nédélec (RTN) spaces of vector functions on the mesh T n ,

RTN l (T n ) := v h ∈ H(div, Ω) ; v h | T ∈ RTN l (T ) ∀T ∈ T n , where RTN l (T ) := P d l (T ) + xP l (T ), l ≥ 0. In particular, v h ∈ RTN l (T n ) is such that ∇•v h ∈ P l (T ) for all T ∈ T n , v h •n F ∈ P l (F ) for all F ∈ F n
T and all T ∈ T n , and such that its normal trace is continuous, cf. [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. In certain cases, a submesh of T n will be used instead. We use the notation P k (T ) for the vector space spanned by polynomials of total degree ≤ k on T and P k (T n ) for the vector space spanned by discontinuous piecewise polynomials of total degree ≤ k on T n . 4.1. Discontinuous Galerkin. The discontinuous Galerkin method for the discretization of (1.1a)-(1.1c) on the time interval I n , 1 ≤ n ≤ N , and the corresponding mesh T n reads: Find

u n hτ ∈ V n h := P k (T n ), k ≥ 1, such that (∂ t u n hτ , v h ) - F ∈F n {(n F •{ {∇ n u n hτ } }, [[v h ]]) F + θ(n F •{ {∇ n v h } }, [[u n hτ ]]) F } + (∇ n u n hτ , ∇ n v h ) + F ∈F n (α F h -1 F [[u n hτ ]], [[v h ]]) F = ( f n , v h ) ∀v h ∈ V n h . (4.1)
Here θ ∈ {-1, 0, 1} and α F are positive parameters. Following [START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF][START_REF] Kim | A posteriori error estimators for locally conservative methods of nonlinear elliptic problems[END_REF], we use the following definition of θ n ∈ RTN l (T n ), l ∈ {k -1, k}: For all T ∈ T n , all F ∈ F n T , and all q h ∈ P l (F ),

(θ n •n F , q h ) F = -n F •{ {∇ n u n hτ } } + α F h -1 F [[u n hτ ]], q h F , (4.2)
and for all r h ∈ P d l-1 (T ),

(θ n , r h ) T = -(∇ n u n hτ , r h ) T + θ F ∈F n T ω F (n F •r h , [[u n hτ ]]) F , (4.3)
where ω F := 1 2 for F ∈ F i,n and ω F := 1 for F ∈ F b,n . Taking v h ∈ P 0 (T n ) in (4.1), yields (3.4) as in [START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF]Theorem 3.1]. Finally, the approximation property (3.17) can be proven from (4.2)-(4.3) as in [START_REF] Ern | Guaranteed and robust discontinuous Galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF].

Remark 4.1 (Residual estimator). A consequence of (4.1) and (4.2)-( 4.3) is that, for l = k, there holds

∂ t u n hτ + ∇•θ n = Π V n h f n . (4.4)
Thus, the residual estimator η n R is in this case the sum of a superconvergent estimator plus the nonconformity estimator associated with ∂ t (su hτ ).

Cell-centered finite volumes.

A general cell-centered finite volume method for the discretization of (1.1a)-(1.1c) on the time interval I n , 1 ≤ n ≤ N , and the corresponding mesh T n reads: Find ūn

hτ ∈ V n h := P 0 (T n ) such that 1 τ n (ū n hτ -u n-1 hτ , 1) T + F ∈F n T S n T,F = ( f n , 1) T ∀T ∈ T n . (4.5)
Here S n T,F are the diffusive fluxes through the faces F of an element T . For our a posteriori error estimates we do not need the specific form of S n T,F , as long as the conservation property S n T -,F = -S n T + ,F holds for all F = ∂T -∩ ∂T + ∈ F i,n . A simple example is the so-called "two-point" scheme, see [START_REF] Eymard | Finite volume methods[END_REF], which requires the "orthogonality condition" of the mesh. Finally, u n-1 hτ , n ≥ 2, is constructed from ūn-1 hτ and the fluxes S n-1 T,F as specified below. The function u n-1 hτ is needed only when the two meshes T n-1 and T n are different; we can replace it by ūn-1 hτ otherwise. We use u n-1 hτ instead of ūn-1 hτ for two reasons. Firstly, under the form (4.5), the scheme enters exactly the present framework; a minor modification would be necessary otherwise. Secondly, since u n-1 hτ is a kind of regularization of the piecewise constant function ūn-1 hτ , we find that it is better suited for being evaluated on the different mesh T n . Following [START_REF]Residual flux-based a posteriori error estimates for finite volume and related locally conservative methods[END_REF], we define θ n ∈ RTN 0 (T n ) by (θ n •n, 1) F := S n T,F for all T ∈ T n and all F ∈ F n T , 1 ≤ n ≤ N . Next, we define u n hτ ∈ V n h := P 1,2 (T n ), where P 1,2 (T n ) is the space P 1 (T n ) enriched elementwise by the parabolas Then, (4.5) immediately yields (3.4) using the Green theorem. By construction, ∇u n hτ + θ n = 0, whence (3.17) is trivial. 4.3. Mixed finite elements. The mixed finite element method for the discretization of (1.1a)-(1.1c) on the time interval I n , 1 ≤ n ≤ N , and the corresponding mesh T n reads: Find

σ n hτ ∈ W n h and ūn hτ ∈ V n h such that (σ n hτ , w h ) -(ū n hτ , ∇•w h ) = 0 ∀ w h ∈ W n h , (4.7a) (∇•σ n hτ , v h ) + 1 τ n (ū n hτ -u n-1 hτ , v h ) = ( f n , v h ) ∀v h ∈ V n h . (4.7b)
The couple W n h × V n h can be any of the usual mixed finite element spaces, cf. [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. In particular, in the RTN method, W n h = RTN l (T n ) and V n h = P l (T n ), l ≥ 0. As for the cell-centered finite volume method, ūn hτ is postprocessed. Following [START_REF] Arbogast | On the implementation of mixed methods as nonconforming methods for second-order elliptic problems[END_REF][START_REF]Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods[END_REF], we define where Π W n h and Π V n h are respectively the L 2 -orthogonal projections onto W n h and V n h . Equation (4.8) thus represents a higher-order polynomial equivalent of (4.6). The spaces V n h are typically P l+1 (T n ) enriched by bubbles in the RTN case and are specified in [START_REF] Arbogast | On the implementation of mixed methods as nonconforming methods for second-order elliptic problems[END_REF]. As in the previous section, u n-1 hτ in (4.7b) can be replaced by ūn-1 hτ whenever the two meshes T n-1 and T n coincide.

u n hτ ∈ V n h such that Π W n h (-∇ n u n hτ ) = σ n hτ , Π V n h (u n hτ ) = ūn hτ , ( 4
In order to use our a posteriori error estimates, we put θ n := σ n hτ . Taking v h ∈ P 0 (T n ) in (4.7b) and considering (4.8) yields (3.4). One does no longer have ∇ n u n hτ + θ n = 0 except for the lowest-order case, see [START_REF]Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods[END_REF], but owing to (4.8), this quantity is expected to be negligible and act as a numerical quadrature.

Remark 4.2 (Residual estimator). Similarly to the relation (4.4) for discontinuous Galerkin methods, there holds for mixed finite elements

∂ t ūn hτ + ∇•θ n = Π V n h f n when T n-1 = T n .
4.4. Vertex-centered finite volumes. The vertex-centered finite volume method for the discretization of (1.1a)-(1.1c) on the time interval I n , 1 ≤ n ≤ N , and the corresponding mesh T n reads: Find

u n hτ ∈ V n h := P 1 (T n ) ∩ H 1 0 (Ω) such that (∂ t u n hτ , 1) D -(∇u n hτ •n D , 1) ∂D = ( f n , 1) D ∀D ∈ D i,n . (4.9)
Here D n is the dual mesh around vertices created using face and element barycenters as indicated in the left part of Figure 4.1; D i,n corresponds to the interior vertices, D b,n to the boundary ones, and n D denotes the outward unit normal to D. We define θ n ∈ RTN 0 (S n ), where S n is the fine simplicial submesh of both T n and D n as indicated in the right part of Figure 4.1. We set

θ n •n F | F := -∇u n hτ •n F | F
on all faces F of S n included in ∂D for some D ∈ D i,n . Note that (3.4) on all D ∈ D i,n is then a direct consequence of (4.9). There are two possibilities to define θ n •n F | F on the remaining faces of S n (that is, those faces located inside some D ∈ D n and those located on the boundary ∂Ω). Firstly, one can directly prescribe

θ n •n F | F := -{ {∇u n hτ •n F } }
on all the remaining faces F of S n ; the present a posteriori error estimates can then be used on the dual meshes D n , but not on the original simplicial ones T n . Alternatively, following [START_REF]Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients[END_REF], one can solve local mixed finite element problems on the simplicial submesh S n D of D. On D ∈ D i,n , a Neumann boundary condition given by -∇u n hτ •n D is prescribed on ∂D. On D ∈ D b,n , zero Dirichlet boundary conditions on ∂D ∩ ∂Ω are supplemented. The right hand-side is given by f n -∂ t u n hτ . This leads to a local linear system solution on all D ∈ D n . The key advantage is that solving these local linear systems to reconstruct the flux more accurately ensures that (3.4) is satisfied

T n D n Fig. 4.2.
Simplicial mesh T n and the associated face-centered dual mesh D n on all T ∈ T n ; see [START_REF]Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients[END_REF] for the details and further alternatives. Finally, (3.17) can be shown using the techniques of [START_REF]Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients[END_REF].

4.5. Face-centered finite volumes. The face-centered finite volume method for the discretization of (1.1a)-(1.1c) on the time interval I n , 1 ≤ n ≤ N , and the corresponding mesh T n reads: Find

u n hτ ∈ V n h such that (∂ t u n hτ , 1) D -(∇ n u n hτ •n D , 1) ∂D = ( f n , 1) D ∀D ∈ D i,n . (4.10)
Here V n h is the Crouzeix-Raviart space of piecewise linear polynomials such that the face jumps are orthogonal to constants, and D n is the dual mesh around faces created using vertices and element barycenters as indicated in Figure 4.2; D i,n corresponds to the interior faces and D b,n to the boundary ones. We define S n D for all D ∈ D i,n as the union of the two simplices sharing the associated face; S n D := {D} for D ∈ D b,n . We again look for θ n ∈ RTN 0 (S n ). We set

θ n •n F | F := -∇ n u n hτ •n F | F
on all faces F of the dual mesh D n . As in the previous section, we then have several ways to define the remaining fluxes, and we easily get (3.4) using the Green theorem from (4.10). Finally, (3.17) can be shown using the techniques of [START_REF]Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients[END_REF].

5. Proof of the error upper bound. The goal of this section is to prove Theorem 3.2. The proof is decomposed into several steps.

Abstract

• Y -norm error estimate. We begin with an abstract • Ynorm error estimate that can be formulated in a rather general setting for the potential and flux reconstructions. In particular, assumptions (3.1), (3.2), and (3.4) are not yet required.

Let s ∈ X ∩ H 1 (0, T ; L 2 (Ω)), so that ∂ t s ∈ L 2 (0, T ; L 2 (Ω)) and hence s ∈ Y . Let θ ∈ L 2 (0, T ; H(div, Ω)). We define the residual R(s, θ) ∈ X ′ such that for all ϕ ∈ X,

R(s, θ), ϕ X ′ ,X := T 0 {(f -∂ t s -∇•θ, ϕ)(t) -(∇s + θ, ∇ϕ)(t)} dt.
(5.1)

At this stage, the contribution of θ can be eliminated since (∇•θ, ϕ) + (θ, ∇ϕ) = 0 owing to the Green theorem; however, the two contributions composing R(s, θ) will be treated separately.

Lemma 5.1 (Abstract • Y -norm error estimate). Let s ∈ X ∩ H 1 (0, T ; L 2 (Ω)). Let θ ∈ L 2 (0, T ; H(div, Ω)). Then, u -u hτ Y ≤ s -u hτ Y + 3 R(s, θ) X ′ + 2 1/2 s 0 -u 0 .
(5.2)

Proof. We first bound us Y . Since v := us is in Y , there holds (see, e.g., [15, Theorem 5.9.3])

1 2 u -s 2 (T ) = 1 2 u 0 -s 0 2 + T 0 ∂ t (u -s), u -s (t) dt.
As a result,

u -s 2 X ≤ 1 2 u -s 2 (T ) + u -s 2 X = 1 2 u 0 -s 0 2 + T 0 { ∂ t (u -s), u -s + (∇(u -s), ∇(u -s))}(t) dt = 1 2 u 0 -s 0 2 + T 0 {(f -∂ t s, u -s) -(∇s, ∇(u -s))}(t) dt,
where we have used (2.1). Adding and subtracting (θ, ∇(us)) in the integrand for a.e. t ∈ (0, T ), using the Green theorem and the definition (5.1) of R(s, θ) yields

u -s 2 X ≤ R(s, θ) X ′ u -s X + 1 2 s 0 -u 0 2 . Since x 2 ≤ ax + b 2 implies x ≤ a + b, it is inferred that u -s X ≤ R(s, θ) X ′ + 2 -1/2 s 0 -u 0 .
Let now ϕ ∈ X with ϕ X = 1 and observe that

∂ t (u -s), ϕ X ′ ,X = T 0 {(f -∂ t s -∇•θ, ϕ) -(∇s + θ, ∇ϕ) -(∇(u -s), ∇ϕ)}(t) dt, whence ∂ t (u -s) X ′ ≤ R(s, θ) X ′ + u -s X , so that u -s Y ≤ 3 R(s, θ) X ′ + 2
R(s, θ) X ′ ≤ N n=1 In T ∈T n (η n R,T + η n DF,T (t)) 2 dt 1/2 + f -f X ′ . Proof. Let ϕ ∈ X with ϕ X = 1. Then, R(s, θ), ϕ X ′ ,X = T 0 {(f -f , ϕ) + ( f -∂ t s -∇•θ, ϕ) -(∇s + θ, ∇ϕ)}(t) dt := T 1 + T 2 + T 3 . Clearly, |T 1 | ≤ f -f X ′ ϕ X = f -f X ′ .
Moreover, owing to (3.1), there holds s ∈ P 1 τ (H 1 0 (Ω)) and θ ∈ P 0 τ (H(div, Ω)), so that

T 2 = N n=1 In ( f n -∂ t s n -∇•θ n , ϕ(t)) dt.
For all 1 ≤ n ≤ N , owing to (3.3) and (3.4),

( f n -∂ t s n -∇•θ n , 1) T = 0, ∀T ∈ T n .
Hence, for a.e.

t ∈ I n , ( f n -∂ t s n -∇•θ n , ϕ(t)) = ( f n -∂ t s n -∇•θ n , ϕ(t) -Π n 0 ϕ(t)) ≤ T ∈T n C P h T f n -∂ t s n -∇•θ n T ∇ϕ T (t),
where we have used the Poincaré inequality on each T ∈ T n stating that ϕ-Π n 0 ϕ T ≤ C P h T ∇ϕ T . Moreover,

T 3 ≤ N n=1 In T ∈T n ∇s(t) + θ n T ∇ϕ T (t) dt.
Collecting the above estimates yields using Cauchy-Schwarz inequalities

|T 2 + T 3 | ≤ N n=1 In T ∈T n (η n R,T + η n DF,T (t)) 2 dt 1/2
.

The conclusion is now straightforward.

s -u hτ Y ≤ N n=1 In T ∈T n (η n NC1,T ) 2 (t) dt 1/2 + N n=1 τ n T ∈T n (η n NC2,T ) 2 1/2 . Proof. It is clear that s -u hτ X = N n=1 In T ∈T n (η n NC1,T ) 2 (t) dt 1/2
.

Let now ϕ ∈ X with ϕ X = 1. Observe that since both s and u hτ are piecewise affine and continuous in time,

∂ t (s -u hτ ), ϕ X ′ ,X = N n=1 In (∂ t (s -u hτ ) n , ϕ(t)) dt.
For all 1 ≤ n ≤ N , it is inferred from (3.3) that the quantity ∂ t (su hτ ) n has zero mean on each element T ∈ T n . Hence, using the Poincaré inequality yields

∂ t (s -u hτ ), ϕ X ′ ,X = N n=1 In (∂ t (s -u hτ ) n , ϕ(t) -Π n 0 ϕ(t)) dt ≤ N n=1 In T ∈T n η n NC2,T ∇ϕ T (t) dt,
whence the desired estimate is inferred using the Cauchy-Schwarz inequality.

5.4. Proof of Theorem 3.2. We observe that using the definition (3.9) for η IC , Lemma 5.1 yields uu hτ Y ≤ su hτ Y + 3 R(s, θ) X ′ + η IC , and we use the bounds on R(s, θ) X ′ and su hτ Y derived, respectively, in Lemmas 5.2 and 5.3 to conclude.

6. Proof of the error lower bound. The goal of this section is to prove Theorem 3.9. The proof is decomposed into several steps. Since the error lower bound is local in time, we keep the integer 1 ≤ n ≤ N fixed in this section.

Equivalent expression for η n

NC,1 . It is convenient to simplify the expression for η n NC,1 . Lemma 6.1 (Equivalent expression for η n NC,1 ). Letting v := su hτ , there holds

τ n 1 6 ( ∇ n-1 v n-1 2 + ∇ n v n 2 ) ≤ (η n NC,1 ) 2 ≤ τ n 1 2 ( ∇ n-1 v n-1 2 + ∇ n v n 2 ).
Proof. Let T ∈ T n-1,n . Since v is piecewise affine in time and smooth in T , it is inferred that

In ∇v 2 T (t) dt = τ n 1 0 ∇v n-1 + τ ∇(v n -v n-1 ) 2 T dτ = τ n 1 3 ( ∇v n-1 2 T + ∇v n 2 T + (∇v n-1 , ∇v n ) T ).
Hence,

τ n 1 6 ( ∇v n-1 2 T + ∇v n 2 T ) ≤ In ∇v 2 T (t) dt ≤ τ n 1 2 ( ∇v n-1 2 T + ∇v n 2 T ).
Summing over T ∈ T n-1,n yields the conclusion.

6.2. Bounds on η n R and η n DF . We introduce the following quantities

(E n R ) 2 := τ n T ∈T n h 2 T Π V n h f n -∂ t u n hτ + ∆u n hτ 2 T , (6.1)
(E n DF ) 2 := τ n T ∈T n |n•[[∇ n u n hτ ]]| 2 + 1 2 ,F i,n T , (6.2)
and observe that E n R and E n DF take the form of usual residual-based a posteriori error estimators for the heat equation. We also recall that the quantity E n f defined by (3.19). Finally, we define the following jump seminorm

J n * (u hτ ) 2 := τ n T ∈T n |[[u n hτ ]]| 2 -1 2 ,F n T . (6.3)
Lemma 6.2 (Bounds on η n R and η n DF ). Assume (M1) and the approximation property (3.17) for the flux reconstruction θ. Then,

η n R E n R + E n DF + E n f + η n NC,2 + J n * (u hτ ), (6.4)
η n DF u -u hτ Y (In) + E n R + E n DF + E n f + η n NC,1 + η n NC,2 + J n * (u hτ ). (6.5) Proof. (i) Bound on η n R . Observe that f n -∂ t s n -∇•θ n = ( f n -Π V n h f n ) -(∂ t s n -∂ t u n hτ ) + (Π V n h f n -∂ t u n hτ + ∆u n hτ ) -(∆u n hτ + ∇•θ n ).
Hence, using the triangle inequality and the above definitions leads to

(η n R ) 2 (E n R ) 2 + (E n f ) 2 + (η n NC,2 ) 2 + τ n T ∈T n h 2 T ∆u n hτ + ∇•θ n 2 T .
Using an inverse inequality to bound the last term above and using (3.17) yields As a result, using the bound (6.4) for η n R leads to

τ n T ∈T n h 2 T ∆u n hτ + ∇•θ n 2 T τ n T ∈T n ∇u n hτ + θ n 2 T (E n DF ) 2 + J n * (u hτ ) 2 , whence (6.
In A n (t) dt u -u hτ 2 Y (In) + (E n R ) 2 + (E n DF ) 2 + (E n f ) 2 + (η n NC,1 ) 2 + (η n NC,2 ) 2 + J n * (u hτ ) 2 + τ n B n . Thus, (η n DF ) 2 u-u hτ 2
Y (In) +(E n R ) 2 +(E n DF ) 2 +(E n f ) 2 +(η n NC,1 ) 2 +(η n NC,2 ) 2 +J n * (u hτ ) 2 +τ n B n ,
and it remains to bound the last term. To this purpose, we use the triangle inequality, the lower bound in Lemma 6.1, and the approximation property (3.17) of the flux reconstruction θ to infer

τ n B n τ n { ∇ n (s n -u n hτ ) 2 + ∇ n u n hτ + θ n 2 } (η n NC,1 ) 2 + (E n DF ) 2 + J n * (u hτ ) 2 .
The conclusion is straightforward. To estimate the bubble contribution, the transition condition (3.14) must be used. For all T ′ ∈ T n,n+1 such that T ′ ⊂ T , using (3.16) and the Cauchy-Schwarz inequality,

∇β n hτ T ′ = |α n T ′ | ∇b T ′ T ′ ≤ h T ′ ∇b T ′ T ′ |T ′ | 1/2 (b T ′ , 1) T ′ h -1 T ′ ι n hτ -u n hτ T ′ .
Owing to the shape regularity of T n,n+1 , the factor between parentheses is bounded uniformly, so that

∇β n hτ T ′ h -1 T ′ ι n hτ -u n hτ T ′ .
We now use the transition condition (3.14), the shape regularity of the meshes, and the approximation properties of the Oswald interpolation operator to infer

∇β n hτ 2 T = T ′ ∈T n,n+1 ;T ′ ⊂T ∇β n hτ 2 T ′ T ′ ∈T n,n+1 ;T ′ ⊂T h -2 T ′ ι n hτ -u n hτ 2 T ′ ≤ Ξ 2 h -2 T T ′ ∈T n,n+1 ;T ′ ⊂T ι n hτ -u n hτ 2 T ′ = Ξ 2 h -2 T ι n hτ -u n hτ 2 T |[[u n hτ ]]| 2 -1 2 ,F n T .
The proof is complete. Lemma 6.4 (Bound on η n NC,2 ). Assume (M1)-(M3) and that the potential reconstruction s is defined by (3.15)- (3.16). Then,

η n NC,2
J n (u hτ ).

Proof. Let v := su hτ . Using the transition condition (3.14) and the triangle inequality yields

T ∈T n h 2 T ∂ t v n 2 T T ∈T n-1 h 2 T (τ n ) -2 v n-1 2 T + T ∈T n h 2 T (τ n ) -2 v n 2 T .
Consider the second term in the right-hand side. Proceeding as in the proof of Lemma 6.3 leads to the bounds

v n T |[[u n hτ ]]| + 1 2 ,F n
T , so that, owing to the shape regularity of the meshes and the condition (M3),

T ∈T n h 2 T (τ n ) -2 v n 2 T T ∈T n h 4 T (τ n ) -2 |[[u n hτ ]]| 2 -1 2 ,F n T T ∈T n |[[u n hτ ]]| 2 -1 2 ,F n T .
The term at t n-1 is treated similarly using again the transition condition in combination with (M3).

Bounds on E n

R and E n DF . The last step consists in bounding the usual residual-based error estimators E n R and E n DF . Lemma 6.5 (Bounds on E n R and E n DF ). Under the assumptions of Lemmas 6.2 and 6.4, there holds

E n R + E n DF u -u hτ Y (In) + E n f + J n (u hτ ).
Proof. Using the technique of element and edge bubble functions introduced by Verfürth [START_REF] Verfürth | A review of a posteriori error estimation and adaptive mesh-refinement techniques[END_REF] and proceeding as in [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF], it can be shown (details are skipped for brevity) that

(E n R ) 2 + (E n DF ) 2 ≤ C( u -u hτ 2
Y (In) + (E n f ) 2 ) + ǫ 2 In ∇ n-1,n (u hτ (t) -u n hτ ) 2 dt,
where ǫ can be chosen as small as needed. The triangle inequality yields ∇ n-1,n (u hτ (t)u n hτ ) ≤ ∇ n-1,n (su hτ )(t) + ∇s(t) + θ n + ∇ n u n hτ + θ n .

Hence, taking square roots

E n R + E n DF ≤ C( u -u hτ Y (In) + E n f ) + ǫ(η n DF + η n NC,1 + E n DF + J n * (u hτ )).
A.2. Nonconforming finite elements. The nonconforming finite element method for the discretization of (1.1a)-(1.1c) on the time interval I n , 1 ≤ n ≤ N , and the corresponding mesh T n reads: Find u n hτ ∈ V n h such that

(∂ t u n hτ , v h ) + (∇ n u n hτ , ∇ n v h ) = ( f n , v h ) ∀v h ∈ V n h . (A.3)
Here V n h is defined as in §4.5; the nonconforming finite element method given by (A.3) is, in fact, very close to the face-centered finite volume method of §4.5. Our estimates can be adapted to the present setting by proceeding exactly as in §A.1.

hτ u 2 hτ u 3 hτFig. 2 . 1 .

 2321 Fig. 2.1. Time-dependent meshes and discrete solutions

d i=1 x 2 i

 2 , such that for all T ∈ T n , -∇u n hτ = θ n , (u n hτ , 1) T = (ū n hτ , 1) T . (4.6)

Fig. 4 . 1 .

 41 Fig. 4.1. Simplicial mesh T n and the associated vertex-centered dual mesh D n (left) and the fine simplicial mesh S n D of D ∈ D n (right)

5. 3 .

 3 Computable upper bound on su hτ Y . The next step is to derive a computable upper bound on su hτ Y since this quantity also involves a dual norm. Lemma 5.3 (Computable upper bound on s-u hτ Y ). Assume (3.2). Let η n NC1,T and η n NC2,T be defined by (3.7)-(3.8). Then,

4 )

 4 readily follows by taking square roots. (ii) Bound on η n DF . Let t ∈ I n . The triangle inequality yields (η n DF ) 2

Proof. 1 2

 1 Owing to the upper bound in Lemma 6.1, it suffices to prove that for m ∈ {n -1, n} and for all T ∈ T m ,∇(s mu m hτ ) T |[[u m hτ ]]| -1 2 ,F mT . The proof is presented for m = n; it is similar for m = n -1. Lettingι n hτ := I n Os (u n hτ ), β n hτ := T ′ ∈T n,n+1 α n T ′ b T ′ ,it is clear that (3.15) yields the decomposition s nu n hτ = (ι n hτu n hτ ) + β n hτ , and we bound, on each T ∈ T n , the two terms in the right-hand side separately. Classical approximation properties of the Oswald interpolation operator yield∇(ι n hτu n hτ ) T |[[u n hτ ]]| -

  .19) Theorem 3.9 (Error lower bound). Assume that the meshes satisfy (M1)-(M3), that the potential reconstruction s is defined by (3.15)-(3.16), and that the flux reconstruction θ satisfies (3.17). Let 1 ≤ n ≤ N . Let E n f be defined by (3.19) and let the jump seminorm J n (u hτ ) be defined by (3.18). Finally, let η n

R , η n DF , η n NC,1 , and η n NC,2 be defined by (3.10)-(3.13). Then,

  1/2 s 0u 0 . The triangle inequality concludes the proof.5.2.Computable upper bound on R(s, θ) X ′ . The dual norm R(s, θ) X ′ in the abstract error estimate (5.2) is not easily computable. We are now going to infer a computable upper bound for this quantity by making use of assumptions(3.

	1),
	(3.2), and (3.4).
	Lemma 5.2 (Computable upper bound on R(s, θ) X ′ ). Assume (3.1), (3.2), and (3.4). Let η n R,T and η n DF,T be defined by (3.5)-(3.6). Then,

  6.3. Bounds on η n NC,1 and η n NC,2 . The next step is to bound η n NC,1 and η n NC,2 . Lemma 6.3 (Bound on η n NC,1 ). Assume (M1)-(M2) and that the potential reconstruction s is defined by (3.15)-(3.16). Then,

	η n NC,1	J

n (u hτ ).

This work was partly supported by the Groupement MoMaS (PACEN/

In ∇s(t) -∇s n 2 dt + τ n ∇s n + θ n 2 =:

Let v be the space-time function on I n such that v(t) := s(t)s n and observe that v ∈ X(I n ). Elementary algebra yields A n (t) = (∇s(t) -∇s n , ∇v(t)) = -{(f, v)(t) -(∂ t s n , v(t)) -(∇s(t), ∇v(t))} + (f (t)f n , v(t))

Owing to (2.1),

, the Cauchy-Schwarz inequality yields

. Collecting these bounds and since

Using the triangle inequality for the first term in the right-hand side and bounding su hτ Y (In) as in Lemma 5.3 yields

Owing to Lemmas 6.2, 6.3, and 6.4,

where we have used the fact that

The conclusion is now straightforward by choosing ǫ sufficiently small to eliminate the terms E n R and E n DF from the right-hand side.

6.5. Proof of Theorem 3.9. The proof is a direct consequence of Lemmas 6.2, 6.3, 6.4, and 6.5.

Appendix A. Conforming and nonconforming finite elements. This appendix extends the above theory to the conforming and nonconforming finite element methods.

A.1. Conforming finite elements. The conforming finite element method for the discretization of (1.1a)-(1.1c) on the time interval I n , 1 ≤ n ≤ N , and the corresponding mesh T n reads: Find

The conforming finite element method given by (A.1) is very close to the vertexcentered finite volume method of §4.4. In fact, as in [START_REF]Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients[END_REF], we find that (A.1) is equivalent to look for u n hτ ∈ V n h such that

Here, V is always the vertex associated with the dual volume D, ψ V is the associated basis "hat" function, and T V is the patch of elements of T n sharing V . Let

on all faces F of S n included in ∂D for some D ∈ D i,n . Then, contrarily to §4.4, (3.4) is not satisfied on D ∈ D i,n ; indeed, the right-hand side of (A.2), clearly equivalent to a numerical quadrature, has to be taken into account. In the proof of Lemma 5.2, the bound on the term T 2 has to be enriched by

where φ is a piecewise linear, vertex-based interpolation of ϕ and φ is a piecewise constant one on the mesh T n . Consequently, two estimators have to be added in Theorem 3.

Ω) by φ and φ. To work on the mesh T n , the flux reconstruction is determined by solving local linear systems, similarly to §4.4. The approximation property (3.17) can be shown using the techniques of [START_REF]Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients[END_REF]. To prove a lower bound for the new estimators C I,T h T ∂ t u n hτ -Π n 0 (∂ t u n hτ ) T , condition (M3) together with the Poincaré inequality is necessary to proceed as in Lemma 6.2, whereas the estimators corresponding to C I,T h T f n -Π n 0 ( f n ) T represent a data oscillation term.