Alexandre Ern
email: alexandre.ern@enpc.fr

Iain Smears
email: iain.smears@inria.fr

Martin Vohralík
email: martin.vohralik@inria.fr

Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems *

Keywords: Parabolic partial differential equations, a posteriori error estimates, local spacetime efficiency, polynomial-degree robustness, high-order methods AMS subject classifications: 65M15, 65M60

    Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems

Introduction

We consider the heat equation ∂ t u -∆u = f in Ω × (0, T ), u = 0 on ∂Ω × (0, T ), u(0) = u 0 in Ω, (1.1) where Ω ⊂ R d , 1 ≤ d ≤ 3, is a bounded, connected, polyhedral open set with Lipschitz boundary, and T > 0 is the final time. We assume that f ∈ L 2 (0, T ; L 2 (Ω)), and that u 0 ∈ L 2 (Ω). We are interested here in developing a posteriori error estimates for a class of high-order discretizations of (1.1). In particular, we consider a conforming finite element method (FEM) in space on unstructured shape-regular meshes, and a discontinuous Galerkin discretization in time, where one is free to vary the approximation orders p in space and q in time, as well as the mesh size h and time-step size τ , leading to what we call a hp-τ q method. These methods are highly attractive from the point of view of flexibility, accuracy, and computational efficiency, since it is known from a priori analysis that judicious local adaptation of the discretization parameters can lead to exponential convergence rates with respect to the number of degrees of freedom, even for solutions with singularities near domain corners, edges, and at initial times [START_REF] Schötzau | Time discretization of parabolic problems by the hp-version of the discontinuous Galerkin finite element method[END_REF][START_REF] Schwab | p-and hp-finite element methods[END_REF][START_REF] Werder | hp-discontinuous Galerkin time stepping for parabolic problems[END_REF]. In practice, it is desirable to determine the adaptation algorithmically, which requires rigorous and high-quality a posteriori error control in order to exploit the potential for high accuracy and efficiency of hp-τ q discretizations. We recall that a posteriori error estimates should ideally give guaranteed upper bounds on the error, i.e. without unknown constants, should be locally efficient, meaning that the local estimators should be bounded from above by the error measured in a local neighbourhood, and, moreover, should be robust, with all constants in the bounds being independent of the discretization parameters; we refer the reader to [START_REF]A posteriori error estimation techniques for finite element methods[END_REF] for an introduction to these concepts.

In the context of parabolic problems, the a posteriori error analysis for low-and fixed-order methods has received significant attention over the past decade, with efforts mostly concentrated on fixed-order FEM in space coupled with an implicit Euler or Crank-Nicolson timestepping scheme, leading to estimates for a wide range of norms. These include estimates for the L 2 (H 1 )-norm of the error considered independently by Picasso and Verfürth [START_REF] Picasso | Adaptive finite elements for a linear parabolic problem[END_REF][START_REF] Verfürth | A posteriori error estimates for nonlinear problems: L r (0, T ; W 1,ρ (Ω))-error estimates for finite element discretizations of parabolic equations[END_REF], with efficiency bounds typically requiring restrictions on the relation between the sizes of the time-steps and the meshes. Estimates for the L 2 (H 1 ) ∩ H 1 (H -1 )-norm estimates were first considered by Verfürth in [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF], who crucially proved local-in-time yet global-in-space efficiency of estimators without restrictions between time-step and mesh sizes, see also Bergam, Bernardi, and Mghazli [START_REF] Bergam | A posteriori analysis of the finite element discretization of some parabolic equations[END_REF]. Guaranteed upper bounds were later obtained by Ern and Vohralík in [START_REF] Ern | A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF], with similar efficiency results as in [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF]. There are also upper bounds in L 2 (L 2 ), L ∞ (L 2 ) and L ∞ (L ∞ ) and higher order norms, based on either duality techniques as in Eriksson and Johnson [START_REF] Eriksson | Adaptive finite element methods for parabolic problems. II. Optimal error estimates in L ∞ L 2 and L ∞ L ∞[END_REF] or the elliptic reconstruction technique originally due to Makridakis and Nochetto [START_REF] Makridakis | Elliptic reconstruction and a posteriori error estimates for parabolic problems[END_REF] and later considered in the fully discrete context by Lakkis and Makridakis [START_REF] Lakkis | Elliptic reconstruction and a posteriori error estimates for fully discrete linear parabolic problems[END_REF], see also [START_REF] Lakkis | A comparison of duality and energy a posteriori estimates for L ∞ (0, T ; L 2 (Ω)) in parabolic problems[END_REF] and the references therein. Repin [START_REF] Repin | Estimates of deviations from exact solutions of initial-boundary value problem for the heat equation[END_REF] studied so-called functional estimates. Finally, a posteriori error estimates developed in the context of the heat equation often serve as a starting point for extensions to diverse applications, including nonlinear problems and spatially non-conforming methods among others [START_REF] Di Pietro | Adaptive regularization, linearization, and discretization and a posteriori error control for the two-phase Stefan problem[END_REF][START_REF] Dolejší | A framework for robust a posteriori error control in unsteady nonlinear advection-diffusion problems[END_REF][START_REF] Georgoulis | A posteriori error control for discontinuous Galerkin methods for parabolic problems[END_REF][START_REF] Kreuzer | Reliable and efficient a posteriori error estimates for finite element approximations of the parabolic p-Laplacian[END_REF][START_REF] Nicaise | A posteriori error estimates for a nonconforming finite element discretization of the heat equation[END_REF]. Adaptive algorithms for parabolic problems are studied in [START_REF] Chen | An adaptive finite element algorithm with reliable and efficient error control for linear parabolic problems[END_REF][START_REF] Gaspoz | A convergent timespace adaptive dG(s) finite element method for parabolic problems motivated by equal error distribution[END_REF][START_REF] Kreuzer | Design and convergence analysis for an adaptive discretization of the heat equation[END_REF].

It is apparent from the literature that, even for low-and fixed-order methods, there are remaining outstanding issues, particularly in terms of the efficiency of the estimators. The efficiency of the estimators is significantly influenced by the choice of norm to be estimated, with the strongest available results being attained by Y -norm estimates, where henceforth Y := L 2 (H 1 0 ) ∩ H 1 (H -1 ). However, as mentioned above, even in this norm, the full space-time local efficiency of the estimators is not known. It is helpful to examine here more closely the issue of spatial locality of estimates in order to motivate the approach adopted in this work. For example, let us momentarily consider an implicit Euler discretization in time and a conforming FEM in space, recalling that the implicit Euler method corresponds to the lowest-order discontinuous Galerkin time-stepping method, which uses piecewise constant approximations with respect to time. Since the resulting numerical solution u hτ is discontinuous with respect to time, and thus u hτ / ∈ Y , it is usual to consider a reconstruction, denoted by Iu hτ ∈ Y , obtained by piecewise linear interpolation at the time-step nodes, and it is seemingly natural to seek a posteriori error estimates for u -Iu hτ Y , where • Y is defined in (2.1) below, and where u is the solution of (1.1); for example, this corresponds to the approach adopted in [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF].

The main issue encountered in studying the efficiency of estimators with respect to the error measured by u -Iu hτ Y is that Iu hτ fails to satisfy the Galerkin orthogonality property on the discrete level: instead, Iu hτ satisfies

In (f, v hτ ) -(∂ t Iu hτ , v hτ ) -(∇Iu hτ , ∇v hτ )dt = In (∇(u hτ -Iu hτ ), ∇v hτ )dt, (1.2)
for all discrete test functions v hτ that are constant in time over the given time interval I n and belong to the associated finite element space V n h (see section 3 for complete definitions). It is seen from the right-hand side of (1.2) that a discrete residual arises from the difference between the numerical solution u hτ and its reconstruction Iu hτ . Since it does not appear possible to show that the discrete residual in (1.2) is controlled locally in space and time by the corresponding local space-time components of the Y -norm of the error, one cannot obtain local efficiency of the estimators with respect to this norm; we note that this issue remains essentially independent of the specific construction of the a posteriori error estimators, whether they are residual-type estimators as in [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF] or equilibrated flux estimators as we consider here. Nevertheless, Verfürth [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF] showed in the lowest-order case that the global-in-space local-in-time norm of the discrete residual can be bounded by the corresponding global-in-space local-in-time Y -norm of u-Iu hτ , which leads to time-local yet space-global efficiency of the estimators. In order to overcome the issue of the loss of spatial locality, we observe that the discrete residual in (1.2) is equivalent to the temporal jumps in the numerical solution u hτ , which is an error component in itself since it measures the lack of conformity in Y of the numerical solution u hτ . It is therefore natural to consider a composite norm u -u hτ E Y that includes both u -Iu hτ Y and the norm of jumps of the numerical solution u hτ . As we explain below, we then recover the fully space-time local efficiency of our estimators with respect to u -u hτ E Y : see (1.4) below, and see Theorem 5.2 of section 5.

For hp-FEM discretizations, one of the key issues concerns the robustness of the estimators with respect to the polynomial degree; this issue appears already in the context of elliptic problems, where Melenk and Wohlmuth [START_REF] Melenk | A posteriori error estimation and adaptive computational methods[END_REF] and Melenk [START_REF] Melenk | hp-interpolation of nonsmooth functions and an application to hp-a posteriori error estimation[END_REF] showed that the well-known residual estimators fail to be polynomial-degree robust. In a breakthrough work, Braess, Pillwein, and Schöberl [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF] established the polynomial-degree robustness of estimators based on equilibrated fluxes, in the context of elliptic diffusion problems. These estimators are based on a globally H(div)-conforming flux constructed from mixed finite element approximations of local Neumann problems over vertex-centred patches of the mesh. The polynomial degree robustness of these estimators was then recently generalized to nonconforming and mixed methods for elliptic problems in [START_REF]Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF], to which we refer the reader for further references on the literature of equilibrated flux estimators for elliptic problems. In the context of parabolic problems, we must also address the additional question of robustness of the estimators with respect to the temporal polynomial degrees. In comparison to low-and fixed-order methods, there are comparatively few works on a posteriori error estimates for high-order discretizations of parabolic problems. Building on the earlier work of Makridakis and Nochetto [START_REF]A posteriori error analysis for higher order dissipative methods for evolution problems[END_REF], Schötzau and Wihler [START_REF] Schötzau | A posteriori error estimation for hp-version timestepping methods for parabolic partial differential equations[END_REF] studied the effect of the temporal approximation order of a posteriori estimates for a composite norm of L ∞ (L 2 ) ∩ L 2 (H 1 )-type, in the context of high-order temporal semi-discretizations of abstract evolution equations by the discontinuous and continuous Galerkin time-stepping methods. Otherwise, it appears that a posteriori error estimates for hp-τ q discretizations of parabolic problems remain essentially untouched.

In this work, we present guaranteed, locally space-time efficient, and polynomial degree robust a posteriori error estimators for hp-τ q discretizations of parabolic problems. This is by no means simple, as it requires the treatment of the challenges that have been outlined above. Our main results are the following.

Let the spaces Y := L 2 (H 1 0 )∩H 1 (H -1 ) and X := L 2 (H 1 0 ) be respectively equipped with their standard norms • Y and • X defined in (2.1) below. Let Y + V hτ be the sum of the continuous and approximate solution spaces, recalling that V hτ ⊂ X and that u hτ / ∈ Y due to the temporally discontinuous approximation. Let I : Y + V hτ → Y be the reconstruction operator defined in section 3.5 below, where we note that

Iv = v if and only if v ∈ Y . Let the norm • E Y be defined by v 2 E Y := Iv 2 Y + v -Iv 2 X for all v ∈ Y + V hτ .
Guaranteed upper bounds In Theorem 5.2 of section 5, we show a posteriori estimates in the norm • E Y . In particular, we have Iu = u, so u-u hτ 2

E Y = u-Iu hτ 2 Y + u hτ -Iu hτ 2
X , where we note that u hτ -Iu hτ X is a measure of the temporal jumps of the numerical solution u hτ . In the absence of data oscillation, our bound takes the simple form

u -u hτ 2 E Y ≤ N n=1 K∈T n In σ hτ + ∇Iu hτ 2 K + ∇(u hτ -Iu hτ ) 2 K dt , (1.3)
where σ hτ is the H(div)-conforming reconstruction; see sections 3 and 4 for full definitions of the notation and construction of the estimators.

Polynomial-degree robustness and local space-time efficiency We establish local spacetime efficiency of our estimators with polynomial-degree robust constants, expressed by the lower bound

In σ hτ + ∇Iu hτ 2 K + ∇(u hτ -Iu hτ ) 2 K dt a∈V K |u -u hτ | 2 E a,n Y + oscillation, (1.4)
where K is an element of the mesh T n for time-step I n , where V K denotes the set of vertices of K, and where |u -

u hτ | E a,n Y
is the local component of u -u hτ E Y on the patch associated with the vertex a. Here, and in the following, the notation a b means that a ≤ Cb, with a constant C that depends possibly on the shape-regularity of the spatial meshes, but is otherwise independent of the mesh-size, time-step size, as well as the spatial and temporal polynomial degrees. We stress that this efficiency bound does not require any relation between the sizes of the time-step and the mesh. The full bound is stated in Theorem 5.2 below.

In addition to the above results, the estimators proposed here are advantageous in terms of flexibility, since they do not require restrictions on coarsening or refinement between time-steps that appeared in earlier works, such as the transition condition used in [37, p. 196, 201]. The main tool to avoid this condition is Lemma 8.1 below.

Relation between u -u hτ E Y and u -Iu hτ Y In Theorem 5.1, we simplify and generalize to higher-order temporal discretizations a key result of Verfürth [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF], namely that the jumps in the numerical solution can be controlled locally-in-time and globally-in-space by the Y -norm of u -Iu hτ . Specifically, for arbitrary approximation orders and for any time-step interval I n , we show that

In ∇(u hτ -Iu hτ ) 2 dt In ∂ t (u -Iu hτ ) 2 H -1 (Ω) + ∇(u -Iu hτ ) 2 dt + oscillation,
where the constant, which is in fact known explicitly, is independent of all other quantities, including the temporal polynomial degree. The associated oscillation term involves the minimum of the source term data oscillation and the coarsening error, and thus can be controlled in practice.

In the absence of this oscillation, we have

u -Iu hτ Y ≤ u -u hτ E Y ≤ 3 u -Iu hτ Y , (1.5)
in addition to the upper and lower bounds (1.3) and (1.4). The key implication is that u-u hτ E Y and u -Iu hτ Y are globally equivalent, although their local distributions may differ. We also stress that the equivalence is independent of the polynomial degrees. We offer some more refined equivalence results in section 6 in order to treat the case of possibly non-vanishing oscillation terms. This paper is organized as follows. First, in section 2 we introduce a functional setting for the a posteriori error analysis. We find it worthwhile to provide a complete derivation of the inf-sup analysis of the problem, as we give here quantitatively sharp results that are advantageous for the efficiency of the estimators in practice. Section 3 defines the setting in terms of notation, finite element approximation spaces, and the numerical scheme. Then, in section 4, we define the equilibrated flux reconstruction used in the a posteriori error estimates. In section 5 we gather our main results underlying (1.3), (1.4) and (1.5). The proofs of the main results are treated in the subsequent sections: section 6 establishes the relation between u-u hτ E Y and u-Iu hτ Y ; the proof of the guaranteed upper bound is given in section 7; and the efficiency of the estimators is the subject of section 8.

Inf-sup theory

Recall that Ω ⊂ R d , 1 ≤ d ≤ 3 is a bounded, connected, polyhedral open set with Lipschitz boundary. For an arbitrary open subset ω ⊂ Ω, we use (•, •) ω to denote the L 2 -inner product for scalar-or vector-valued functions on ω, with associated norm • ω . In the special case where ω = Ω, we drop the subscript notation, i.e. • := • Ω . We consider the function spaces

X := L 2 (0, T ; H 1 0 (Ω)) and Y := L 2 (0, T ; H 1 0 (Ω)) ∩ H 1 (0, T ; H -1 (Ω)
), which we equip with the following norms:

ϕ 2 Y := T 0 ∂ t ϕ 2 H -1 (Ω) + ∇ϕ 2 dt + ϕ(T ) 2 ∀ ϕ ∈ Y, v 2 X := T 0 ∇v 2 dt ∀ v ∈ X.
(2.1)

Define the bilinear form B

Y : Y × X → R by B Y (ϕ, v) := T 0 ∂ t ϕ, v + (∇ϕ, ∇v) dt, (2.2)
where ϕ ∈ Y and v ∈ X are arbitrary functions, and •, • denotes here the duality pairing between H -1 (Ω) and H 1 0 (Ω). Then, the problem (1.1) admits the following weak formulation: find u ∈ Y such that u(0) = u 0 and such that

B Y (u, v) = T 0 (f, v) dt ∀ v ∈ X.
(2.

3)

The well-posedness of (2.3) is well-known and can be shown by Galerkin's method [START_REF] Evans | Partial differential equations[END_REF][START_REF] Wloka | Partial differential equations[END_REF].

The following result states an inf-sup stability result for the bilinear form B Y for the above spaces equipped with their respective norms. The inf-sup stability result presented here has the interesting and important property of taking the form of an identity, which is advantageous for the sharpness of a posteriori error analysis.

Theorem 2.1 (Inf-sup identity). For every ϕ ∈ Y , we have

ϕ 2 Y = sup v∈X\{0} B Y (ϕ, v) v X 2 + ϕ(0) 2 . (2.4)
Proof. For a fixed ϕ ∈ Y , let w * ∈ X be defined by (∇w * , ∇v) = ∂ t ϕ, v for all v ∈ H 1 0 (Ω), a.e. in (0, T ), which implies the identity ∇w *

2 = ∂ t ϕ 2 H -1 (Ω) a.e. in (0, T ). Furthermore, we have B Y (ϕ, v) = T 0 (∇(w * + ϕ), ∇v) dt, thus implying that sup v∈X\{0} B Y (ϕ, v)/ v X = w * + ϕ X .
Thus, we obtain the desired identity (2.4) by expanding the square

sup v∈X\{0} B Y (ϕ, v) v X 2 = T 0 ∇(w * + ϕ) 2 dt = T 0 ∇w * 2 + 2(∇w * , ∇ϕ) + ∇ϕ 2 dt = T 0 ∂ t ϕ 2 H -1 (Ω) + 2 ∂ t ϕ, ϕ + ∇ϕ 2 dt = ϕ 2 Y -ϕ(0) 2 , (2.5)
where we note that we have used the identity

T 0 2 ∂ t ϕ, ϕ dt = ϕ(T ) 2 -ϕ(0) 2 .
In order to estimate the error between the solution u of (1.1) and its approximation, we define the residual functional R

Y : Y → X by R Y (ϕ), v := B Y (u -ϕ, v) = T 0 (f, v) -∂ t ϕ, v -(∇ϕ, ∇v) dt, (2.6)
where v ∈ X and ϕ ∈ Y , and •, • denotes here the duality pairing between the dual space X and X. The dual norm of the residuals is naturally defined by R Y (ϕ)

X := sup v∈X\{0} R Y (ϕ),v v X
. Theorem 2.1 implies the following equivalence between the error and dual norm of the residual: for all ϕ ∈ Y , we have

u -ϕ 2 Y = R Y (ϕ) 2 X + u 0 -ϕ(0) 2 .
(2.7)

Finite element approximation

Consider a partition of the interval (0, T ) into time-step intervals

I n := (t n-1 , t n ), with 1 ≤ n ≤ N , where it is assumed that [0, T ] = N n=1
I n , and that {t n } N n=0 is strictly increasing with t 0 = 0 and t N = T . For each interval I n , we let τ n := t n -t n-1 denote the local time-step size. We will not need any special assumptions about the relative sizes of the time-steps to each other. We associate a temporal polynomial degree q n ≥ 0 to each time-step I n , and we gather all the polynomial degrees in the vector q = (q n ) N n=1 . For a general vector space V , we shall write Q qn (I n ; V ) to denote the space of V -valued univariate polynomials of degree at most q n over the time-step interval I n .

Meshes

We associate a matching simplicial mesh T n of the domain Ω for each 0 ≤ n ≤ N , where we assume shape-regularity of the meshes uniformly over all time-steps. This allows us to treat many applications where the meshes are obtained by refinement or coarsening between time-steps. We consider here only matching simplicial meshes for simplicity, although we indicate that mixed simplicial-parallelepipedal meshes, possibly containing hanging nodes, can be also be treated: see [START_REF]hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF] for instance. The mesh T 0 will be used to approximate the initial datum u 0 . For each element K ∈ T n , let h K := diam K denote the diameter of K. We associate a local spatial polynomial degree p K ≥ 1 to each K ∈ T n , and we gather all spatial polynomial degrees in the vector p n = (p K ) K∈T n . In order to keep our notation sufficiently simple, the dependence of the local spatial polynomial degrees p K on the time-step is kept implicit, although we bear in mind that the polynomial degrees may change between time-steps.

Approximation spaces

For a general matching simplicial mesh T with associated vector of polynomial degrees p = (p K ) K∈T , p K ≥ 1 for all K ∈ T , the

H 1 0 (Ω)-conforming hp-finite element space V h (T , p) is defined by V h (T , p) := v h ∈ H 1 0 (Ω), v h | K ∈ P p K (K) ∀ K ∈ T , (3.1)
where P p K (K) denotes the space of polynomials of total degree at most p K on K. For shorthand, we denote

V n h := V h (T n , p n ) for each 0 ≤ n ≤ N . Let Π h u 0 ∈ V 0
h denote an approximation to the initial datum u 0 , a typical choice being the L 2 -orthogonal projection onto V 0 h . Given the collection of timesteps {I n } N n=1 , the vector q of temporal polynomial degrees, and the hp-finite element spaces {V n h } N n=1 , the spatio-temporal finite element space V hτ is defined by

V hτ := v hτ | (0,T ) ∈ X, v hτ | In ∈ Q qn (I n ; V n h ) ∀ n = 1, . . . , N, v hτ (0) ∈ V 0 h . (3.2)
Functions in V hτ are generally discontinuous with respect to the time-variable at the partition points, although we take them to be left-continuous: for all 1 ≤ n ≤ N , we define v hτ (t n ) as the trace at t n of the restriction v hτ | In . Functions in V hτ are thus left-continuous; moreover they also have a well-defined value at t 0 = 0. For all 0 ≤ n < N , we denote the right-limit of

v hτ ∈ V hτ at t n by v hτ (t + n ). Then, the temporal jump operators • n , 0 ≤ n ≤ N -1, are defined on V hτ by v hτ n := v hτ (t n ) -v hτ (t + n ), 0 ≤ n ≤ N -1. (3.3)

Refinement and coarsening

Similary to other works, e.g., [37, p. 196], we assume that we have at our disposal a common refinement mesh T n of T n-1 and T n for each 1 ≤ n ≤ N , as well as associated polynomial degrees

p n = (p K ) K∈ T n , such that V n-1 h + V n h ⊂ V n h := V h ( T n , p n ). For a function v hτ ∈ V hτ , we observe that v hτ n-1 ∈ V n h for each 1 ≤ n ≤ N since v hτ (t n-1 ) ∈ V n-1 h , v hτ (t + n-1 ) ∈ V n h , and V n-1 h + V n h ⊂ V n h .
It is assumed that T n has the same shape-regularity as T n-1 and T n , and that every element K ∈ T n is wholly contained in a single element K ∈ T n-1 and a single element K ∈ T n . We emphasize that we do not require any assumptions on the relative coarsening or refinement between successive spaces V n-1 h and V n h . We note that in the present context, refinement and coarsening can be obtained by modification of the meshes as well as change in the polynomial degrees. Concerning the polynomial degrees, we may choose for example p K = max(p K , p K ). In the case where V n h is obtained from V n-1 h by refinement without coarsening, then we may choose T n := T n and p n := p n so that V n h = V n h . However, we do not need the transition condition assumption from [37, p. 196, 201], which requires a uniform bound on the ratio of element sizes between T n and T n .

Numerical scheme

The numerical scheme for approximating the solution of the parabolic problem (1.1) consists of finding u hτ ∈ V hτ such that u hτ (0) = Π h u 0 , and, for each time-step interval I n ,

In (∂ t u hτ , v hτ ) + (∇u hτ , ∇v hτ ) dt -u hτ n-1 , v hτ (t + n-1 ) = In (f, v hτ ) dt ∀ v hτ ∈ Q qn (I n ; V n h ). (3.4)
Here the time derivative ∂ t u hτ is understood as the piecewise time-derivative on each timestep interval I n . The numerical solution u hτ ∈ V hτ can thus be obtained by solving the fully discrete problem (3.4) on each successive time-step. At each time-step, this requires solving a linear system that is symmetric only in the lowest-order case; this can be performed efficiently in practice for arbitrary orders, see [START_REF] Smears | Robust and efficient preconditioners for the discontinuous Galerkin time-stepping method[END_REF] and the references therein.

Reconstruction operator

For each time-step interval I n and each nonnegative integer q, let L n q denote the polynomial on I n obtained by mapping the standard q-th Legendre polynomial under an affine transformation of (-1, 1) to I n . It follows that L n q (t n ) = 1 for all q ≥ 0, and L n q (t n-1 ) = (-1) q , and that the mapped Legendre polynomials {L n q } q≥0 are L 2 -orthogonal on I n , and satisfy In |L n q | 2 dt = τn 2q+1 for all q ≥ 0. We introduce the Radau reconstruction operator I defined on V hτ by

(Iv hτ )| In := v hτ | In + (-1) qn 2 L n qn -L n qn+1 v hτ n-1 ∀ v hτ ∈ V hτ . (3.5)
It is clear that I is a linear operator on V hτ . It follows from the properties of the Legendre polynomials that

Iv hτ | In (t n ) = v hτ (t n ), and that Iv hτ | In (t + n-1 ) = v hτ (t n-1
) for all 1 ≤ n ≤ N . Therefore, Iv hτ is continuous with respect to the temporal variable at the interval partition points {t n } N -1 n=0 , and thus we have

Iv hτ ∈ H 1 (0, T ; H 1 0 (Ω)) ⊂ Y, Iv hτ | In ∈ Q qn+1 I n ; V n h ∀ v hτ ∈ V hτ , (3.6)
where we recall that

V n-1 h + V n h ⊂ V n h .
We easily deduce the following property of the reconstruction operator I from integration-by-parts and the orthogonality of the polynomials L n qn and L n qn+1 to all polynomials of degree strictly less than q n on the time-step interval I n :

In ∂ t Iv hτ φ dt = In ∂ t v hτ φ dt -v hτ n-1 φ(t + n-1 ) ∀ φ ∈ Q qn (I n ; R), (3.7)
where equality holds in the above equation in the sense of functions in V n h . We may therefore use (3.7) to rewrite the numerical scheme (3.4) as

In (∂ t Iu hτ , v hτ ) + (∇u hτ , ∇v hτ ) dt = In (f, v hτ ) dt ∀ v hτ ∈ Q qn (I n ; V n h ). (3.8)
Note also that Iu hτ (0) = Π h u 0 . Remark 3.1 (Alternative equivalent definitions). The operator I is the Radau reconstruction operator commonly used in the a posteriori error analysis of the DG time-stepping method [START_REF]A posteriori error analysis for higher order dissipative methods for evolution problems[END_REF] and in the a priori error analysis of time-dependent first-order PDEs [START_REF] Ern | Discontinuous Galerkin method in time combined with a stabilized finite element method in space for linear first-order PDEs[END_REF]. Several equivalent definitions of I have appeared in the literature, although it will be particularly advantageous for our purposes to use the definition (3.5) of I due to [START_REF] Smears | Robust and efficient preconditioners for the discontinuous Galerkin time-stepping method[END_REF], to which we refer the reader for further discussion on the equivalence of the various definitions.

Construction of the equilibrated flux

The a posteriori error estimates presented in this paper are based on a discrete and locally computable H(div)-conforming flux σ hτ that satisfies the key equilibration property

∂ t Iu hτ + ∇•σ hτ = f hτ in Ω × (0, T ), (4.1)
where Iu hτ is defined in section 3.5, and f hτ ≈ f is a data approximation defined in (4.4) below. We call σ hτ an equilibrated flux. We consider here the natural extension of existing flux reconstructions for elliptic problems [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF][START_REF] Braess | Equilibrated residual error estimator for edge elements[END_REF][START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF][START_REF]Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF] to the parabolic setting; see also [START_REF] Dolejší | A posteriori error estimates for nonstationary problems[END_REF]. In particular, for each time-step, σ hτ is obtained as a sum of fluxes computed by solving local mixed finite element problems over the vertex-based patches of the current mesh, see Definition 4.1 of section 4.3 below.

Local mixed finite element spaces

We now define the mixed finite element spaces that are required for the construction of the equilibrated flux. For each 1 ≤ n ≤ N , let V n denote the set of vertices of the mesh T n , where we distinguish the set of interior vertices V n int and the set of boundary vertices V n ext . For each a ∈ V n , let ψ a denote the hat function associated with a, and let ω a denote the interior of the support of ψ a , with associated diameter h ωa . Furthermore, let T a,n denote the restriction of the mesh T n to ω a . Recalling that the common refinement spaces V n h were obtained with a vector of polynomial degrees p n = (p K ) K∈ T n , we associate to each a ∈ V n the fixed polynomial degree

p a := max K∈ T a,n (p K + 1). (4.2)
Observe that ψ a ∂ t Iu hτ | K×In is a polynomial function with degree at most q n in time and at most p a in space for each K ∈ T a,n , 1 ≤ n ≤ N . For a polynomial degree p ≥ 0, let the local spaces P p ( T a,n ) and RTN p ( T a,n ) be defined by

P p ( T a,n ) := {q h ∈ L 2 (ω a ), q h | K ∈ P p ( T a,n ) ∀ K ∈ T a,n } RTN p ( T a,n ) := {v h ∈ L 2 (ω a ; R d ), v h | K ∈ RTN p ( K) ∀ K ∈ T a,n },
where RTN p ( K) := P p ( K; R d ) + P p ( K)x denotes the Raviart-Thomas-Nédélec space of order p on K. It is important to notice that whereas the patch ω a is subordinate to the vertices of the mesh T n , the spaces P p ( T a,n ) and RTN p ( T a,n ) are subordinate to the submesh T a,n ; of course, in the absence of coarsening, this distinction vanishes.

We now introduce the local spatial mixed finite element spaces V a,n h and Q a,n h , defined by

V a,n h :=    v h ∈ H(div, ω a ) ∩ RTN pa ( T a,n ), v h • n = 0 on ∂ω a if a ∈ V n int , v h ∈ H(div, ω a ) ∩ RTN pa ( T a,n ), v h • n = 0 on ∂ω a \ ∂Ω if a ∈ V n ext , Q a,n h := q h ∈ P pa ( T a,n ), (q h , 1) ωa = 0 if a ∈ V n int , P pa ( T a,n ) if a ∈ V n ext .
We then define the following space-time mixed finite element spaces

V a,n hτ := Q qn (I n ; V a,n h ), Q a,n hτ := Q qn (I n ; Q a,n h ). (4.3)

Data approximation

Our a posteriori error estimates given in section 5 involve certain approximations of the source term f appearing in (1.1). It is helpful to define these approximations here. First, we define the semi-discrete approximation f τ of f by L 2 -orthogonal projection in time. In particular, the approximation

f τ ∈ Q qn (I n ; L 2 (Ω)) is defined on each interval I n by In (f -f τ , v) dt = 0 for all v ∈ Q qn (I n ; L 2 (Ω)). Next, for each 1 ≤ n ≤ N and for each a ∈ V n , let Π a,n hτ be the L 2 ψa -orthogonal projection from L 2 (I n ; L 2 ψa (ω a )) onto Q qn (I n ; P pa-1 ( T a,n ))
, where L 2 ψa (ω a ) is the space of measurable functions v on ω a such that ωa ψ a |v| 2 dx < ∞. In other words, the projection operator Π a,n hτ is defined by In (ψ a Π a,n hτ v, q hτ ) ωa dt = In (ψ a v, q hτ ) ωa dt for all q hτ ∈ Q qn (I n ; P pa-1 ( T a,n )). We adopt the convention that Π a,n hτ v is extended by zero from ω a × I n to Ω × (0, T ) for all v ∈ L 2 (I n ; L 2 ψa (ω a )). Then, we define f hτ by

f hτ := N n=1 a∈V n ψ a Π a,n hτ f. (4.4)
Remark 4.1 (Definition of f hτ ). The somewhat technical appearance of the definition of f hτ is due to the possible variation in polynomial degrees across the mesh and the particular requirements of the analysis of efficiency, in particular the hypotheses of Lemma 8.1 below. Nevertheless, f hτ has several important approximation properties. First, for any 1 ≤ n ≤ N , any K ∈ T n and any real-valued polynomial φ of degree at most q n , we have

In (f -f hτ , 1) K φ dt = a∈V K In (ψ a (f -Π a,n hτ f ), φ1) K dt = 0, (4.5)
where V K denotes the set of vertices of K, and where we use the fact that the hat functions {ψ a } a∈V n form a partition of unity on Ω. Furthermore, using the orthogonality of the projector Π a,n hτ and the fact that 0

≤ ψ a ≤ 1 in Ω, it is straightforward to show that f -f hτ L 2 (In;L 2 ( K)) ≤ √ d + 1 inf w hτ ∈Qq n (In;Pp K ( K)) f -w hτ L 2 (In;L 2 ( K)) ,
This shows that f hτ defines an approximation of f that is at least of the same order as the one associated with the finite element approximation.

Flux reconstruction

For each 1 ≤ n ≤ N and each a ∈ V n , let the scalar function g a,n hτ ∈ Q qn (I n ; P pa ( T a,n )) and vector field τ a,n hτ ∈ Q qn (I n ; RTN pa ( T a,n )) be defined by τ a,n hτ := -ψ a ∇u hτ | ωa×In , (4.6a)

g a,n hτ := ψ a (Π a,n hτ f -∂ t Iu hτ ) | ωa×In -∇ψ a • ∇u hτ | ωa×In . (4.6b)
We claim that for all a ∈ V n int , (g a,n hτ (t), 1)

ωa = 0 ∀ t ∈ I n , (4.7)
which is equivalent to showing that g a,n hτ ∈ Q a,n hτ for all a ∈ V n . Indeed, we first observe that the construction of the numerical scheme, in particular identity (3.8), implies that, for any univariate real-valued polynomial φ of degree at most q n on I n ,

In (g a,n hτ , φ1) ωa dt = In f, φ ψ a ωa -∂ t Iu hτ , φ ψ a ωa -∇u hτ , ∇(φ ψ a ) ωa dt = 0,
where we have used the orthogonality of the projection Π a,n hτ and the fact that φψ a ∈ Q qn (I n ; V n h ) is a valid test function in (3.8). Since the function g a,n hτ is polynomial in time with degree at most q n , i.e. g a,n hτ ∈ Q qn (I n ; P pa ( T a,n )), we deduce (4.7). Definition 4.1. Let u hτ ∈ V hτ be the numerical solution of (3.4). For each time-step interval I n and for each vertex a ∈ V n , let the space-time mixed finite element spaces V a,n hτ and Q a,n hτ be defined by (4.3). Let g a,n hτ and τ a,n hτ be defined by (4.6). Let σ a,n hτ ∈ V a,n hτ be defined by

σ a,n hτ := argmin v h ∈V a,n hτ ∇•v h =g a,n hτ In v h -τ a,n hτ 2 ωa dt. (4.8)
Then, after extending σ a,n hτ by zero from ω a × I n to Ω × (0, T ) for each a ∈ V n and for each 1 ≤ n ≤ N , we define

σ hτ := N n=1 a∈V n σ a,n
hτ . (4.9)

Note that σ a,n hτ ∈ V a,n hτ is well-defined for all a ∈ V n : in particular, for interior vertices a ∈ V n int , we use (4.7) to guarantee the compatibility of the datum g a,n hτ with the constraint ∇•σ a,n hτ = g a,n hτ . The following key result shows that σ hτ from Definition 4.1 leads to an equilibrated flux. Theorem 4.2 (Equilibration). Let the flux reconstruction σ hτ be defined by (4.9) of Definition 4.1. Then σ hτ ∈ L 2 (0, T ; H(div, Ω)) and we have (4.1), where the discrete approximation f hτ is defined in (4.4).

Proof. After extending each σ a,n hτ by zero from ω a ×I n to Ω×(0, T ), we have σ a,n hτ ∈ L 2 (0, T ; H(div, Ω)) as a consequence of the boundary conditions included in the definition of the space V a,n h . This immediately implies that σ hτ ∈ L 2 (0, T ; H(div, Ω)). To show (4.1), the definition of the flux reconstruction σ hτ in (4.9) implies that for any time-step interval I n and any

K ∈ T n , ∇•σ hτ | K×In = a∈V K ∇•σ a,n hτ | K×In = a∈V K g a,n hτ | K×In = a∈V K ψ a Π a,n hτ f -ψ a ∂ t Iu hτ -∇ψ a • ∇u hτ | K×In = (f hτ -∂ t Iu hτ )| K×In , (4.10)
where V K denotes the set of vertices of K, where we use the fact that the hat functions {ψ a } a∈V n form a partition of unity in order to pass to the last line of (4.10), and where we have used the definition of f hτ in (4.4). This yields (4.1) as required.

For the purposes of practical implementation, it is easily seen that, for each time-step interval I n , the fluxes σ a,n hτ can be computed by solving q n + 1 independent spatial mixed finite element problems, provided only that an orthogonal or orthonormal polynomial basis is used in time over I n . Moreover, the q n + 1 linear systems each share the same matrix, which helps to simplify the implementation and reduce the computational cost.

Lemma 4.3 (Decoupling). Let σ a,n

hτ ∈ V a,n hτ be defined by (4.8). Then σ a,n hτ is equivalently uniquely defined by: let

(σ a,n hτ , r a,n hτ ) ∈ V a,n hτ × Q a,n hτ solve In (σ a,n hτ , v hτ ) ωa -(∇•v hτ , r a,n hτ ) ωa dt = In (τ a,n hτ , v hτ ) ωa dt ∀ v hτ ∈ V a,n hτ , (4.11a
)

In (∇•σ a,n hτ , q hτ ) ωa dt = In (g a,n hτ , q hτ ) ωa dt ∀ q hτ ∈ Q a,n hτ , (4.11b)
Furthermore, for each 1 ≤ n ≤ N , let {φ n j } qn j=0 be an L 2 (I n )-orthonormal basis for the space of univariate real-valued polynomials of degree at most q n . For each a ∈ V n , define the functions {g a,n h,j } qn j=0 and {τ a,n h,j } qn j=0 over the patch ω a by g a,n h,j := In g a,n hτ φ n j dt, τ a,n h,j := In τ a,n hτ φ n j dt. (4.12)

Then, the solution (σ a,n hτ , r a,n hτ ) of (4.11) can be obtained by solving the following spatial problems: for each 0 ≤ j ≤ q n , find σ a,n h,j ∈ V a,n h and r a,n h,j in Q a,n h such that

(σ a,n h,j , v h ) ωa -(∇•v h , r a,n h,j ) ωa = (τ a,n h,j , v h ) ωa ∀ v h ∈ V a,n h , (4.13a)
(∇•σ a,n h,j , q h ) ωa = (g a,n h,j , q h ) ωa ∀ q h ∈ Q a,n h , (4.13b)

and then by defining σ a,n hτ := qn j=0 σ a,n h,j φ n j and r a,n hτ := qn j=0 r a,n h,j φ n j .

Remark 4.2. The analysis in the subsequent sections shows that one particular advantage of the equilibrated flux σ hτ of Definition 4.1 is that it leads to estimators that are robust with respect to coarsening (and refinement) between time-steps. The price to pay is that the size of the linear systems in (4.13) grows with the size of coarsening between two successive time-steps, as (4.13) are defined on the patches ω a partitioned by the common refinement mesh T n for each 1 ≤ n ≤ N . The analysis in [17, Section 6], though, shows that this computational cost can be significantly reduced to the solution of two low-order systems over the patches ω a , followed by local high-order corrections on the sub-patches of T a,n . We refer the reader to [START_REF]Discrete p-robust H(div)-liftings and a posteriori error analysis of elliptic problems with H -1 source terms[END_REF]Section 6] for the full details of this approach.

Main results

In this section, we present the a posteriori error estimate featuring guaranteed upper bounds, local space-time efficiency, and polynomial-degree robustness. Let the norm

• E Y : Y + V hτ → R ≥0 be defined by v 2 E Y := Iv 2 Y + v -Iv 2 X ∀ v ∈ Y + V hτ , (5.1)
where we recall Remark 3.2 on the extension of the linear operator I to Y + V hτ . Since the exact solution u ∈ Y implies that Iu = u, we have the identities

u -u hτ 2 E Y = u -Iu hτ 2 Y + u hτ -Iu hτ 2 X = u -Iu hτ 2 Y + N n=1 τn(qn+1) (2qn+1) (2qn+3) ∇ u hτ n-1 2 , (5.2)
where we have simplified In ∇(u hτ -Iu hτ ) 2 dt = τn(qn+1)

(2qn+1) (2qn+3) ∇ u hτ n-1 2
, which is an identity easily deduced from (3.5) and from In |L n q | 2 dt = τn 2q+1 for all q ≥ 0; see also [START_REF] Schötzau | A posteriori error estimation for hp-version timestepping methods for parabolic partial differential equations[END_REF]. We also introduce the localized seminorms

|•| E a,n Y , for each 1 ≤ n ≤ N and each a ∈ V n , defined by |v| 2 E a,n Y := In ∂ t Iv 2 H -1 (ωa) + ∇Iv 2 ωa + ∇(v -Iv) 2 ωa dt ∀ v ∈ Y + V hτ . (5.3)
Similarly to (5.2), we find that

|u -u hτ | 2 E a,n Y = In ∂ t (u -Iu hτ ) 2 H -1 (ωa) + ∇(u -Iu hτ ) 2 ωa dt + τn(qn+1) (2qn+1) (2qn+3) ∇ u hτ n-1 2 ωa . (5.4)
Although it might not be obvious that u-u hτ E Y is equivalent to the Hilbertian sum of the |u -u hτ | E a,n Y , up to data oscillation, this will come as a consequence of the results shown here and in section 8. We are now ready to state our main results in Theorems 5.1 and 5.2 below. It is helpful to denote the time-localized dual norm of the residual by

R Y (Iu hτ )| In X := sup v∈X, v X =1 In (f, v) -∂ t Iu hτ , v -(∇Iu hτ , ∇v) dt. ( 5
η n C := τn(qn+1) (2qn+1) (2qn+3) ∇ {u hτ (t n-1 ) -P n h [u hτ (t n-1 )]} , (5.6a) [η n osc,τ ] 2 := In f (t) -f τ (t) 2 H -1 (Ω) dt, (5.6b)
where

P n h : H 1 0 (Ω) → V n h denotes the elliptic orthogonal projection onto V n h defined by (∇P n h w, ∇v h ) = (∇w, ∇v h ) for all v h ∈ V n h . Then, we have In ∇(u hτ -Iu hτ ) 2 dt ≤ 8 R Y (Iu hτ )| In 2 X + min [η n C ] 2 , 8[η n osc,τ ] 2 , (5.7)
where R Y (Iu hτ )| In X is defined in (5.5). Furthermore, we have

u -Iu hτ 2 Y ≤ u -u hτ 2 E Y ≤ 9 u -Iu hτ 2 Y + N n=1 min [η n C ] 2 , 8[η n osc,τ ] 2 .
(5.8)

We delay the proof of Theorem 5.1 until section 6 below.

Remark 5.1 (Equivalence). Theorem 5.1 shows that u -u hτ E Y and u -Iu hτ Y are globally equivalent up to the minimum of temporal data oscillation and coarsening errors. In particular, one of our key contributions here is to obtain polynomial-degree independent constants in (5.8).

It is important to note that although u -u hτ E Y and u -Iu hτ Y are essentially globally equivalent, their local distributions may differ.

Remark 5.2 (Relation to [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF]). A similar result to (5.7) was previously obtained in the lowestorder case q n = 0 by Verfürth [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF]; see in particular the bounds of [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF]Section 7] for what is denoted there τn

3 |u n h -u n-1 h | 2
1 , which is equivalent to In ∇(u hτ -Iu hτ ) 2 dt with q n = 0 in our notation. For higher polynomial degrees, we note that Gaspoz, Kreuzer, Siebert and Ziegler [START_REF] Gaspoz | A convergent timespace adaptive dG(s) finite element method for parabolic problems motivated by equal error distribution[END_REF] have obtained independently an inequality of a similar kind as (5.7).

We introduce the following a posteriori error estimators and data oscillation terms:

η n F,K (t) := σ hτ (t) + ∇Iu hτ (t) K , (5.9a) η n J,K := τn(qn+1) (2qn+1)(2qn+3) ∇ u hτ n-1 K , (5.9b) η n osc,h,K (t) := K∈ T n , K⊂K h 2 K π 2 f τ (t) -f hτ (t) 2 K 1 2
, (5.9c)

η osc,τ (t) := f (t) -f τ (t) H -1 (Ω) , (5.9d) η osc,init := u 0 -Π h u 0 , (5.9e)
where t ∈ I n , K ∈ T n , the equilibrated flux σ hτ is defined in Definition 4.1, and where the data approximations f τ and f hτ are respectively defined in section 4.2. The two estimators η n F,K and η n J,K are our principal estimators, where η n F,K measures respectively the lack of H(div)conformity of the gradient of the reconstructed solution Iu hτ , and where η n J,K measures the lack of temporal conformity of the numerical solution u hτ . The term η n osc,h,K represents the data oscillation due to the spatial discretisation, whereas η osc,τ represents the data oscillation due to the temporal discretisation. We define the global a posteriori error estimators as

η 2 Y := N n=1 In K∈T n [η n F,K + η n osc,h,K ] 2 1 2 + η osc,τ 2 dt + [η osc,init ] 2 , (5.10a) η 2 E Y := η 2 Y + N n=1 K∈T n [η n J,K ] 2 . (5.10b)
Notice that in the absence of data oscillation, namely if

f = f τ = f hτ and u 0 = Π h u 0 , then η Y simplifies to η 2 Y = T 0 σ hτ + ∇Iu hτ 2 dt, and η E Y simplifies to η 2 E Y = T 0 σ hτ + ∇Iu hτ 2 + ∇(u hτ -Iu hτ ) 2 dt.
Recall that we write a b for two quantities a and b if a ≤ Cb with a constant C depending only on the shape regularity of T n and T n , but otherwise independent of the mesh-size, time-step size, and polynomial degrees in space and time.

Theorem 5.2 (E Y -norm a posteriori error estimate). Let u ∈ Y be the weak solution of (1.1), let u hτ ∈ V hτ denote the solution of the numerical scheme (3.4), and let Iu hτ denote its temporal reconstruction, where the operator I is defined in (3.5). Let σ hτ denote the equilibrated flux of Definition 4.1. Let • E Y be defined in (5.1), and let the a posteriori error estimators be defined in (5.9), with η E Y defined in (5.10). Then, we have the guaranteed upper bound

u -u hτ E Y ≤ η E Y .
(5.11)

Moreover, for each 1 ≤ n ≤ N and for each K ∈ T n , the indicators satisfy the following local efficiency bound:

In [η n F,K ] 2 dt + [η n J,K ] 2 a∈V K |u -u hτ | 2 E a,n Y + [η a,n osc ] 2 , (5.12)
where

|•| E a,n Y is defined in (5.
3), V K is the set of vertices of the element K, and the local data oscillation term η a,n osc is defined by

[η a,n osc ] 2 := In f -Π a,n hτ f 2 H -1 (ωa) dt. (5.13)
Furthermore, we have the following global efficiency bound for u -u hτ E Y :

N n=1 K∈T n In [η n F,K ] 2 dt + [η n J,K ] 2 u -u hτ 2 E Y + N n=1 a∈V n [η a,n osc ] 2 . ( 5.14)
The proof of Theorem 5.2 is postponed to the following sections: the proof of the upper bound (5.11) is given in section 7, and the proof of the bounds (5.12) and (5.14) is the subject of section 8. Theorem 5.2 shows the local space-time efficiency of the estimators with respect to u -u hτ E Y . As a consequence of the proof of Theorem 5.2, we can also show guaranteed upper bounds and local-in-time and global-in-space efficiency of the estimators with respect to u -Iu hτ Y , thereby generalising the results to [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF] to higher-order approximations, see Corollary 5.3 below. Furthermore, for each 1 ≤ n ≤ N , we have

K∈T n In [η n F,K ] 2 dt + [η n J,K ] 2 In ∂ t (u -Iu hτ ) 2 H -1 (Ω) + ∇(u -Iu hτ ) 2 dt + min [η n C ] 2 , 8[η n osc,τ ] 2 + a∈V n [η a,n osc ] 2 . (5.16)
Remark 5.3 (Temporal data oscillation). The temporal data oscillation term η osc,τ is defined with respect to a negative norm, as usual in the literature [START_REF] Ern | A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF][START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF]. Similarly to [START_REF] Ern | A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF][START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF], this temporal data oscillation term can be of the same order as the error in terms of the time-step size. Since this term already appears in the upper bounds of the residual-based estimates of [37, Eq. (1.5)], it is seen that this issue is not related to the choice of equilibrated flux a posteriori error estimators, but is rather a part of the error estimation in the Y -norm. In practical computations, it is often advisable to determine a minimal temporal resolution for reducing this term to within a prescribed tolerance, in advance of solving the numerical scheme (3.4). Although the negative norm appearing in the definition of η osc,τ is non-computable, there are several possibilities for estimating it. First, we mention that η osc,τ is bounded from above by C Ω f -f τ , with C Ω the constant of the global Poincaré inequality, although this can be pessimistic in practice. If f is a finite tensorial product of spatial and temporal functions, then sharper bounds can be obtained by solving a set of independent coarse and low-order conforming approximations for elliptic problems, followed by equilibrated flux a posteriori error estimates to achieve guaranteed upper bounds. Finally, we also mention that this issue motivates a posteriori error estimators in other norms: in particular, we show in [START_REF] Ern | Equilibrated flux a posteriori error estimates in L 2 (H 1 )-norms for high-order discretizations of parabolic problems[END_REF] that X-norm a posteriori estimates benefit from data oscillation terms that are of higher-order by an additional factor of √ τ + h.

6 Proof of equivalence between u -u hτ E Y and u -Iu hτ Y

In this section, we prove Theorem 5.1, along with some corollary results, which relate u-u hτ E Y with u -Iu hτ Y . Our starting point involves the following two original bounds on the norms of the jumps, which generalize one of the key results of Verfürth [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF] for the lowest-order case q n = 0. In fact, our result sharpens and simplifies the proof of the result of [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF] even in the lowest-order case.

Lemma 6.1. For each 1 ≤ n ≤ N , let P n h : H 1 0 (Ω) → V n h denote the elliptic orthogonal projection to V n h defined by (∇P n h w, ∇v h ) = (∇w, ∇v h ) for all v h ∈ V n h . Then, for each 1 ≤ n ≤ N , the jump u hτ n-1 satisfies τ n 8q n + 4 ∇ u hτ n-1 2 ≤ R Y (Iu hτ )| In 2 X + τ n 8q n + 4 ∇ {u hτ (t n-1 ) -P n h [u hτ (t n-1 )]} 2 , (6.1)
where R Y (Iu hτ )| In X is defined in (5.5). Furthermore, we also have the alternative bound

τ n 8q n + 12 ∇ u hτ n-1 2 ≤ 2 R Y (Iu hτ )| In 2 X + [η n osc,τ ] 2 . (6.2)
Proof. First, note that (u hτ -Iu hτ )| In = (-1) qn 2 (L n qn+1 -L n qn ) u hτ n-1 belongs to the space Q qn+1 (I n ; V n h ). We define the test function v hτ := -(-1) qn

2

L n qn P n h u hτ n-1 , which belongs to Q qn (I n ; V n h ), and we use it in equation (3.8) for the numerical scheme, which yields, by orthogonality of the Legendre polynomials and by the definition of the orthogonal projector P n h , the identity

In ∇v hτ 2 dt = τ n 8q n + 4 ∇P n h u hτ n-1 2 = In (∇(u hτ -Iu hτ ), ∇v hτ ) dt = In (f -∂ t Iu hτ , v hτ ) -(∇Iu hτ , ∇v hτ ) dt. (6.3)
Therefore, we have In ∇v hτ

2 dt ≤ R Y (Iu hτ )| In 2
X . This bound yields the desired result (6.1) once it is combined with (6.3) and the orthogonality relation

∇P n h u hτ n-1 2 = ∇ u hτ n-1 2 -∇ { u hτ n-1 -P n h u hτ n-1 } 2 = ∇ u hτ n-1 2 -∇ {u hτ (t n-1 ) -P n h [u hτ (t n-1 )]} 2 ,
where the last equality above follows from the facts that u hτ n-1 = u hτ (t n-1 ) -u hτ (t + n-1 ) and that u hτ (t + n-1 ) ∈ V n h . This completes the proof of the first bound (6.1). We now turn to the proof of (6.2); the main difference in the proofs of (6.1) and (6.2) is that above we appealed to the numerical scheme using a discrete test function, whereas to establish (6.2), we shall now consider a higher-order polynomial function that is not in the discrete test space. We define v on I n by v| In := (-1) qn 2 L n qn+1 u hτ n-1 , and then we extend v by zero outside of I n , so that v ∈ X. Then, by orthogonality of the Legendre polynomial L n qn+1 to all polynomials of degree at most q n on I n , we have the identities In (f τ , v) dt = 0, In (∂ t Iu hτ , v) dt = 0 and In (∇u hτ , ∇v) dt = 0. Therefore, we obtain

In ∇v 2 dt = τ n 8q n + 12 ∇ u hτ n-1 2 = In (∇(u hτ -Iu hτ ), ∇v) dt = In (f, v) -(∂ t Iu hτ , v) -(∇Iu hτ , ∇v) + (f τ -f, v) dt.
The desired result (6.2) then follows straightforwardly from the above identity.

Proof of Theorem 5.1 The first inequality u -Iu hτ (5.2). Recalling the definitions of η n J,K in (5.9b) and η n C in (5.6a), we deduce from (6.1) and (6.2) that

2 Y ≤ u -u hτ 2 E Y is obvious from the definition of u -u hτ E Y in
K∈T n [η n J,K ] 2 ≤ 4(q n + 1) (2q n + 3) R Y (Iu hτ )| In 2 X + [η n C ] 2 ≤ 2 R Y (Iu hτ )| In 2 X + [η n C ] 2 , (6.4)
and that

K∈T n [η n J,K ] 2 ≤ 8(q n + 1) (2q n + 1) R Y (Iu hτ )| In 2 X + [η n osc,τ ] 2 ≤ 8 R Y (Iu hτ )| In 2 X + 8[η n osc,τ ] 2 . (6.5)
Therefore, we obtain (5.7) by taking the minimum of the right-hand sides of the above bounds. Finally, we get (5.8) by summing the above inequality over all time-steps and noting that

N n=1 R Y (Iu hτ )| In 2 X = R Y (Iu hτ ) 2 X ≤ u -Iu hτ 2
Y which follows from (2.7). It is possible to obtain slightly sharper variants of Theorem 5.1 under more specific assumptions. For instance, the following corollary shows that u -u hτ E Y is equivalent to u -Iu hτ Y , without any additional data oscillation, whenever the mesh coarsening error is kept relatively small to the jumps. Corollary 6.2. Using the notation of Lemma 6.1, assume that there exists a constant θ ∈ [0, 1)

such that ∇ [ u hτ n-1 -P n h u hτ n-1 ] 2 ≤ θ ∇ u hτ n-1 2 for each 1 ≤ n ≤ N . Then, we have u -Iu hτ 2 Y ≤ u -u hτ 2 E Y ≤ 3 -θ 1 -θ u -Iu hτ 2 Y . (6.6)
Proof. The result is a consequence of Lemma 6.1 and the fact that

u hτ n-1 -P n h u hτ n-1 = u hτ (t n-1 ) -P n h [u hτ (t n-1 )] , which leads to τn 8qn+4 ∇ u hτ n-1 2 ≤ 1 1-θ R Y (Iu hτ )| In 2 X .
Proceeding as in the proof of Theorem 5.1 then yields (6.6).

Proof of the guaranteed upper bound

We prove here (5.11) and (5.15). First, it is clear from (5.2) that (5.15) immediately implies (5.11). Therefore, it remains to show (5.15). Keeping in mind the equivalence identity (2.7) between norms of the errors and residuals, we turn our attention to bounds for the residual norm R Y (Iu hτ ) X = sup v∈X\{0} B Y (u -Iu hτ , v)/ v X . To this end, consider an arbitrary function v ∈ X such that v X = 1. Then, we obtain

R Y (Iu hτ ), v = T 0 (f -∂ t Iu hτ -∇•σ hτ , v) -(σ hτ + ∇Iu hτ , ∇v) dt,
where we have inserted the flux σ hτ and used integration by parts over Ω. Next, we use (4.1), and we write f -f hτ = f -f τ + f τ -f hτ . For any K ∈ T n , 1 ≤ n ≤ N , we deduce from (4.5) that the function t → (f τ (t) -f hτ (t), 1) K , which is a real-valued polynomial of degree at most q n on I n , vanishes identically on I n . Therefore, letting v K (t) denote the mean value of v(t) over the element K ∈ T n , which is defined for a.e. t ∈ I n , we deduce from the Poincaré inequality that

|(f τ (t) -f hτ (t), v(t)) K | ≤ h K π f τ (t) -f hτ (t) K ∇v(t) K . Therefore, R Y (Iu hτ ), v can be bounded as follows: R Y (Iu hτ ), v = N n=1 In (f -f hτ , v) -(σ hτ + ∇Iu hτ , ∇v) dt ≤ N n=1 In K∈T n η n F,K ∇v K + K∈ T n h K π f τ (t) -f hτ (t) K ∇v(t) K + η osc,τ ∇v dt ≤ N n=1 In K∈T n [η n F,K + η n osc,h,K ] ∇v K + η osc,τ ∇v dt ≤ N n=1 In K∈T n [η n F,K + η n osc,h,K ] 2 1 2
+ η osc,τ ∇v dt.

Applying the Cauchy-Schwarz inequality leads to an upper bound for R Y (Iu hτ ) X , which we then combine with the identity (2.7) relating errors and residuals to obtain u -Iu hτ Y ≤ η Y .

The corresponding upper bound u-u hτ E Y ≤ η E Y then follows immediately, as explained above.

Proof of local space-time efficiency and robustness

We prove here the bounds (5.12), (5.14), and (5.16).

Preliminary result

The following lemma is a generalisation of important results on polynomial-degree robustness of equilibrated flux estimates from [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF]Thm. 7], in two space dimensions, and [15, Thm 2.3] in three space dimensions. In particular, Lemma 8.1 comes from [17, Thm 1.2] on the existence of a discrete polynomial-degree robust H(div)-lifting of data that are piecewise-polynomials with respect to the submesh T a,n . Note that [3, Thm. 7] and [15, Thm 2.3] only hold for the case where the data are piecewise-polynomials on the elements K ∈ T n of the patch ω a . This generalisation is crucial for allowing arbitrary refinement and coarsening between time-steps.

Lemma 8.1 (Polynomial degree-robust stability bound). For each 1 ≤ n ≤ N and each a ∈ V n , recall that T a,n denotes the restriction of T n to ω a and that ψ a ∈ H 1 (ω a ) ∩ P 1 ( T a,n ) denotes the hat function associated with ω a . Let Γ a = {x ∈ ∂ω a , ψ a (x) = 0}. Then, for any f a,n h ∈ P pa-1 ( T a,n ) and any ξ a,n h ∈ RTN pa-1 ( T a,n ), where it is further supposed that (f a h , ψ a ) ωa = (ξ a,n h , ∇ψ a ) ωa if Γ a = ∂ω a , we have

min v h ∈H(div,ωa)∩RTNp a ( T a,n ) ∇•v h =ψaf a,n h -∇ψa•ξ a,n h v h •n=0 on Γa v h + ψ a ξ a,n h ωa sup ϕ∈H 1 0 (ωa)\{0} (f a,n h , ϕ) ωa -(ξ a,n h , ∇ϕ) ωa ∇ϕ ωa
Proof. The result is directly obtained by applying [START_REF]Discrete p-robust H(div)-liftings and a posteriori error analysis of elliptic problems with H -1 source terms[END_REF]Thm 1.2], where Ω there stands for ω a here, where T there stands for T a,n here, and where ψ † there stands for ψ a here. In applying [17, Thm 1.2], we use the fact that Γ a is the union of the faces of the mesh T a,n on which ψ a vanishes, and we have simplified the constant appearing there by using the fact that h ωa ∇ψ a ∞ ψ a ∞ = 1 by shape-regularity.

Stability of the space-time flux equilibration

For each 1 ≤ n ≤ N and each a ∈ V n , we introduce the patch residual functional R a,n hτ , with R a,n hτ :

L 2 (I n , H 1 0 (ω a )) → R defined by R a,n hτ , v = In Π a,n hτ f -∂ t (Iu hτ ), v ωa -∇u hτ , ∇v ωa dt, (8.1)
for all v ∈ L 2 (I n , H 1 0 (ω a )). We are now ready to state the essential result that forms the starting point for our analysis of the efficiency of the error estimators. where Q qn (I n ; H 1 0 (ω a )) denotes the space of H 1 0 (ω a )-valued univariate polynomials of degree at most q n on I n .

Proof. The definition of σ a,n hτ ∈ V a,n hτ in (4.8) implies that it is enough to show that there exists a v h ∈ V a,n hτ such that ∇•v h = g a,n hτ and such that In v h -τ a,n hτ 2 ωa dt is bounded by the right-hand side of (8.2). Let {φ n j } qn j=0 be an L 2 -orthonormal basis of polynomials on I n , and let the functions {f a,n h,j } qn j=0 and {ξ a,n h,j } qn j=0 be defined by

f a,n h,j := In (Π a,n hτ f -∂ t Iu hτ )φ n j dt, ξ a,n h,j := In ∇u h φ n j dt. (8.3) It will be useful to keep in mind that g a,n hτ = qn j=0 [ψ a f a,n h,j -∇ψ a • ξ a,n h,j ]φ n j and that τ a,n hτ = - qn j=0 ψ a ξ a,n h,j φ n j . Let Γ a := {x ∈ ∂ω a , ψ a (x) = 0}; note that if a ∈ V n int , then Γ a = ∂ω a , whereas Γ a is a strict subset of ∂ω a if a ∈ V n
ext . We will now use Lemma 8.1 to show that, for each 0 ≤ j ≤ q n , there exists v a,n h,j ∈ H(div, ω a ) ∩ RTN pa ( T a,n ) such that

∇•v a,n h,j = ψ a f a,n h,j -∇ψ a • ξ a,n h,j in ω a , v a,n h,j • n = 0 on Γ a , (8.4a) v a,n h,j + ψ a ξ a,n h,j ωa R a,n h,j H -1 (ωa) , (8.4b)
where R a,n h,j ∈ H -1 (ω a ) is defined by R a,n h,j , v := (f a,n h,j , v) ωa -(ξ a,n h,j , ∇v) ωa for all v ∈ H 1 0 (ω a ). To check the hypotheses of Lemma 8.1, we start by observing that the choice of p a in (4.2) implies that f a,n h,j ∈ P pa-1 ( T a,n ) and that ξ a,n h,j ∈ RTN pa-1 ( T a,n ) for all 0 ≤ j ≤ q n . For any interior vertex a ∈ V n int , it is seen from (4.7) that (f a,n h,j , ψ a ) ωa = (ξ a,n h,j , ∇ψ a ) ωa for all 0 ≤ j ≤ q n . Therefore, the hypotheses of Lemma 8.1 are satisfied, and there exists v a,n h,j ∈ H(div, ω a ) ∩ RTN pa ( T a,n ) satisfying (8.4).

Next, we claim that v a,n h,j ∈ V a,n h for all 0 ≤ j ≤ q n . Indeed, the definition of Γ a implies that Γ a = ∂ω a for all a ∈ V n int and that ∂ω a \ ∂Ω ⊂ Γ a for all a ∈ V n ext . Therefore, we have v a,n h,j ∈ V a,n h for all 0 ≤ j ≤ q n . It then follows that the function v a,n hτ := qn j=0 v a,n h,j φ n j ∈ V a,n hτ and that this function satisfies

∇•v a,n hτ = qn j=0 [ψ a f a,n h,j -∇ψ a • ξ a,n h,j ]φ n j = g a,n hτ , (8.5a) In v a,n hτ -τ a,n hτ 2 ωa dt = qn j=0 v a,n h,j + ψ a ξ a,n h,j 2 ωa qn j=0 R a,n h,j 2
H -1 (ωa) , (8.5b)
where the equality in (8.5b) results from the orthonormality of {φ n j } qn j=0 . We now claim that

   qn j=0 R a,n h,j 2 H -1 (ωa)    1 2 ≤ sup v∈Qq n (In;H 1 0 (ωa))\{0} R a,n hτ , v In ∇v 2 ωa dt 1 2 . (8.6)
For each j = 0, . . . , q n , we define

z j ∈ H 1 0 (ω a ) by (∇z j , ∇v) ωa = R a,n h,j , v for all v ∈ H 1 0 (ω a ). It is then straightforward to show that ∇z j 2 ωa = R a,n h,j , z j = R a,n h,j 2
H -1 (ωa) for each j = 0, . . . , q n . Then, we define z * ∈ Q qn (I n ; H 1 0 (ω a )) by z * := qn j=0 z j φ n j . It follows from the orthonormality of the temporal basis that In ∇z * Hence, the above identities immediately imply (8.6). Therefore, we combine (8.5) and (8.6) to deduce that v a,n hτ ∈ V a,n hτ satisfies ∇•v a,n hτ = g a,n hτ and In v h -τ a,n hτ 2 ωa dt is bounded by the right-hand side of (8.2). This implies (8.2) as explained above.

Local efficiency

We can now prove the local efficiency bound (5.12).

Proof of the local efficiency bound (5.12) Consider a time-step I n and an element K ∈ T n . First, note that [η n J,K ] which in combination with (8.7), yields the desired result (5.12).

Global efficiency

We finally prove the global efficiency bounds (5.14) and (5.16).

Proof of (5. The proof is essentially a counting argument after local Riesz mappings are introduced to evaluate the negative norms ∂ t (u -Iu hτ ) 2 H -1 (ωa) for all a ∈ V n , see also [START_REF] Blechta | Localization of the W -1,q norm for local a posteriori efficiency[END_REF]. Summing (5.12) over K ∈ T n and using (8.10)

Conclusion and outlook

We have studied a posteriori error estimates for hp-τ q discretizations of parabolic problems based on arbitrarily high-order conforming Galerkin spatial discretizations and discontinuous Galerkin temporal discretizations. The equilibrated flux reconstructions lead to guaranteed upper bounds for the norm u -u hτ E Y . Furthermore, the estimators have the key property of being unconditionally locally space-time efficient with respect to the local errors |u -u hτ | E a,n Y , with constants that are fully robust with respect to both the spatial and temporal approximation orders. The estimators are flexible in the sense that they do not require restrictive transition conditions on the refinement and coarsening between time-steps. We also showed that the composite norm of the error u-u hτ E Y is globally equivalent to u-Iu hτ Y up to the minimum of coarsening error and data oscillation, with polynomial degree-robust constants in the equivalence. Finally, the analysis given here can be extended in various directions: in [START_REF] Ern | Equilibrated flux a posteriori error estimates in L 2 (H 1 )-norms for high-order discretizations of parabolic problems[END_REF], we show that the equilibrated flux reconstruction employed here can also be used for obtaining a posteriori estimates for the X-norm of the error, with guaranteed upper bounds, and local space-time efficiency under the natural parabolic condition that h 2 τ . Furthermore, the adaptation of Lemma 8.1 to the case of residual-based estimators is currently under investigation.

. 5 )Theorem 5 . 1 (

 551 Note that R Y (Iu hτ )| In X can always be bounded from above by the restriction of the Y -norm of the error u -Iu hτ to the time-step interval I n . Equivalence of norms). Let the norm • E Y be defined by (5.1), and, for each 1 ≤ n ≤ N , let the temporal data oscillation η n osc,τ and the coarsening error indicator η n C be defined by

Corollary 5 . 3

 53 (Y -norm a posteriori error estimate). Let the estimator η Y be defined by (5.10a). Then, we have u -Iu hτ Y ≤ η Y .(5.15)

Lemma 8 . 2 (

 82 Space-time stability bound). Let σ a,n hτ denote the patch-wise flux reconstructions of Definition 4.1, and let R a,n hτ denote the local patch residual defined by (8.1). Then, we have In σ a,n hτ + ψ a ∇u hτ

2 HR a,n h,j 2 H

 22 -1 (ωa) . Fubini's theorem and(8.3) imply thatR a,n hτ , z * = qn j=0 In (Π a,n hτ f -∂ t Iu hτ , φ n j z j ) ωa -(∇u h , φ n j ∇z j ) ωa dt = qn j=0 In (Π a,n hτ f -∂ t Iu hτ φ n j dt, z j ωa -In ∇u h φ n j dt, ∇z j ωa = qn j=0 (f a,n h,j , z j ) ωa -(ξ a,n h,j , ∇z j ) ωa = qn j=0R a,n h,j , z j = qn j=0 -1 (ωa) .

  [START_REF]Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]) and (5.16) Recalling the definition (5.3) of the localized seminorms|•| E a,n u -Iu hτ ) 2 H -1 (Ω) + ∇(u -Iu hτ ) 2 dt + In ∇(u hτ -Iu hτ ) 2 dt.(8.10)

  Remark 3.2 (Extensions of I to Y +V hτ ). In what follows, it will be helpful to extend I to a linear operator over Y + V hτ . Note that the definition of the jump operators (3.3) can be naturally extended to Y + V hτ , and therefore the definition (3.5) also extends naturally to Y + V hτ . In particular, I : Y + V hτ → Y , and we have Iϕ = ϕ if and only if ϕ ∈ Y , since the jumps of any ϕ ∈ Y vanish identically.

  2 ≤ a∈V K |u -u hτ | 2 E a,nY trivially, where we recall that V K denotes the set of vertices of K. Hence, it remains only to bound In [η n F,K ] 2 dt. To this end, observe that σ hτ | K×In = a∈V K σ a,n hτ | K×In , and thatIn a∈V K (σ a,n hτ + ψ a ∇Iu hτ ) 2 K dt ≤ (|V K | + 1)| is the number of vertices of the element K, which equals d + 1 for simplices and where we have used that • K ≤ • ωa and the definition of η n J,K in the last line. Keeping in mind Lemma 8.2, we therefore turn our attention to bounding the dual norm of the patchwise residuals R a,n hτ for each a∈ V K . Consider an arbitrary v ∈ Q qn (I n , H 1 0 (ω a )) such that In ∇v 2 ωa dt = 1; then (2.3) implies that R a,n hτ , v = In Π a,n hτ f, v ωa -∂ t (Iu hτ ), v ωa -∇u hτ , ∇v ωa dt =The definition of the • H -1 (ωa) -norm and the Cauchy-Schwarz inequality then yield|E 1 + E 2 + E 3 | |u -u hτ | E a,nY. Finally, we find that |E 4 | ≤ η a,n osc where η a,n osc is defined in (5.13).Recalling (8.2) of Lemma 8.2, we deduce that, for each a ∈ V K ,

	[η n F,K ] 2 dt =
	In
		In a∈V K	σ a,n hτ + ψ a ∇u hτ	2 K + ∇(u hτ -Iu hτ ) 2 K dt
	≤ (|V K | + 1)	In a∈V K	σ a,n hτ + ψ a ∇u hτ	2 ωa dt + [η n J,K ] 2 ,	(8.7)
	where |V K Therefore,
	we find that
	sup v∈Qq n (In,H 1 0 (ωa))\{0}	R a,n hτ , v In ∇v 2 ωa dt	1 2	|u -u hτ | 2 E a,n Y	+ [η a,n osc ] 2	1 2 .	(8.8)
						E a,n Y	+ [η a,n osc ] 2 ,	(8.9)

In ∂ t (u -Iu hτ ), v + (∇(u -Iu hτ ), ∇v) ωa dt + In (∇(Iu hτ -u hτ ), ∇v) ωa -(f -Π a,n hτ f, v) ωa dt =: E 1 + E 2 + E 3 + E 4 . In σ a,n hτ + ψ a ∇u hτ 2 ωa dt |u -u hτ | 2

  then leads to In ∂ t (u -Iu hτ ) 2 H -1 (Ω) + ∇(u -Iu hτ ) 2 dt + In ∇(u hτ -Iu hτ ) 2 dt + Summing the bound (8.11) for all 1 ≤ n ≤ N immediately yields (5.14), whereas (5.16) results from (8.11) after invoking (5.7) and observing that R Y (Iu hτ )| In 2 X is bounded from above by In ∂ t (u -Iu hτ ) 2 H -1 (Ω) + ∇(u -Iu hτ ) 2 dt.

		[η n F,K ] 2 dt + [η n J,K ] 2
	K∈T n	In
		[η a,n osc ] 2 . (8.11)
		a∈V n

* This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 647134 GATIPOR).


Measurement of the intrinsic damping constant in individual nanodisks of Y3Fe5O12 and Y3Fe5O12|Pt

 1)Service de Physique de l' État Condensé (CNRS URA 2464), CEA Saclay, 91191 Gif-sur-Yvette, France 2) Unité Mixte de Physique CNRS/Thales and Université Paris Sud 11, 1 av. Fresnel, 91767 Palaiseau, France 3) Institute of Physics, Kazan Federal University, Kazan 420008, Russian Federation 4) Instituto de Sistemas Optoelectrónicos y Microtecnología (UPM), Madrid 28040, Spain 5) Instituto de Microelectrónica de Madrid (CNM, CSIC), Madrid 28760, Spain (Dated: 18 February 2014)

We report on an experimental study on the spin-waves relaxation rate in two series of nanodisks of diameter φ =300, 500 and 700 nm, patterned out of two systems: a 20 nm thick yttrium iron garnet (YIG) film grown by pulsed laser deposition either bare or covered by 13 nm of Pt. Using a magnetic resonance force microscope, we measure precisely the ferromagnetic resonance linewidth of each individual YIG and YIG|Pt nanodisks. We find that the linewidth in the nanostructure is sensibly smaller than the one measured in the extended film. Analysis of the frequency dependence of the spectral linewidth indicates that the improvement is principally due to the suppression of the inhomogeneous part of the broadening due to geometrical confinement, suggesting that only the homogeneous broadening contributes to the linewidth of the nanostructure. For the bare YIG nano-disks, the broadening is associated to a damping constant α = 4 • 10 -4 . A 3 fold increase of the linewidth is observed for the series with Pt cap layer, attributed to the spin pumping effect. The measured enhancement allows to extract the spin mixing conductance found to be G ↑↓ = 1.55 • 10 14 Ω -1 m -2 for our YIG(20nm)|Pt interface, thus opening large opportunities for the design of YIG based nanostructures with optimized magnetic losses.

Yttrium iron garnet (Y 3 Fe 5 O 12 ), commonly referred as YIG, is the champion material for magneto-optical applications as it holds the highest figure of merit in terms of low propagation loss. It is widely used in high-end microwave and optical-communication devices such as filters, tunable oscillators, or non-reciprocal devices. It is also the material of choice for magnonics 1 , which aims at using spin-waves (SW) (or their quanta magnons) to carry and process information. The development of this emerging field is presently limited by the damping constant of SW. Recently it was proposed that spin-current transfer generated by spin Hall effect from an adjacent layer can partially or even fully compensate the intrinsic losses of the traveling SW beyond the natural decay time. The achievement of damping compensation by pure spin current in 1.3µm thick YIG covered by Pt was reported by Kajiwara et al. 2 , although attempts to reproduce the results have so far failed 3,4 . Given that the spin-orbit torque is purely interfacial in such hybrid system, it is crucial to work with nanometer-thick films of epitaxial YIG. Primarly, because the interface spin-current transfer scales inversely with the YIG thickness 5 . Secondly, because it permits nano-patterning of the YIG and thus the engineering of the spin-wave (SW) spectra through spatial confinement [6][7][8] . To the best of our knowledge, however, there is no report yet 9 on the measurement of the dynamical properties on submicron size nanostructure patterned out of YIG ultrathin films.

Benefiting from our recent progress in the growth of very high quality YIG films by pulsed laser deposition (PLD) 4 , here we will demonstrate that these ultra-thin YIG films (thicknesses ranging from 20 to 4 nm) can be reliably nano-patterned into sub-micron size nano-disks. In the following, we will perform a comparative study of the linewidth on these nanostructures. We will analyze the different contributions to the damping by separating the homogeneous from the inhomogeneous broadening and quantify the spin-pumping contribution when an adjacent metallic layer is added. We will also evaluate the consequences of the damages produced by the lithographic process on the dynamical response of these devices. It will be shown that the different alterations produced by chemical solvent, heat treatment, amorphization and redeposition of foreign elements during the lithography, edge roughness etc... do not produce any increase to the linewidth, offering great hope for incorporation of these YIG films in magnonics.

The 20 nm thick YIG has been grown by PLD on a (111) Gd 3 Ga 5 O 12 substrate following the preparation described in Ref. 4 . A 13 nm thick Pt layer was then deposited by sputtering on half of the YIG surface. One 10 resonance values for the SW modes. These predictions are uncorrected by the additional stray field of the MRFM probe. The last row display the SWs eigen-values obtained by adjusting both the radius of the disks (respectively φ = 700, 520 and 380 nm) and the probe sample separation (respectively h = 3.0, 2.0 and 1.5µm). The SWs are labeled by their azimuthal and radial number ( , n) slab of 1 × 5 mm was cut to perform standard magnetic hysteresis and ferromagnetic resonance (FMR) measurements. The measured magnetic parameters of the bare YIG film are summarized in Table I. On the remaining piece, a series of YIG nanodisks has been subsquently patterned using standard electron lithography and dry ion etching. After the YIG lithography, an insulating layer of 50 nm SiO 2 was deposited on the whole surface and a 150 nm thick and 5 µm wide microwave Au-antenna was deposited on top. Fig. 1 shows an optical image of the sample and the antenna pattern. The series of decreasing diameters bare YIG nanodisks is placed on the left. The spacing between the disks is 3 µm. The series of disks on the right side mirrors the first one and has a 13 nm Pt layer on top. Here we concentrate on the disks enclosed in the rectangular area, with nominal diameter φ =700, 500 and 300 nm.

To measure the FMR-spectra of the nanodisks buried under the microwave antenna, we use a ferro-magnetic resonance force microscope (f-MRFM) 11 . It is based on measuring the deflection of a MFM-cantilever with a magnetic Fe particle of about 800 nm diameter affixed to the tip. The tip magnetic dipole moment senses the stray field produced by the perpendicular component M z of the magnetization of the magnetic nanodisks, which is modulated by the exciting microwave power at the mechanical frequency of the cantilever. In Fig. 2 we show f-MRFM spectra recorded for different diameters (by row) on both the YIG and YIG|Pt nanodisks (by column). The spectra correspond to the spinwave eigenmodes of the disks biased by a perpendicular magnetic field H 0 = 4.99 kOe. The largest peak at lowest frequency stems from the lowest energy FMR-mode, the so-called uniform mode. The smaller peaks at higher energy correspond to higher order modes. The splitting corresponds to the quantization of the SW wavenumber ∝ nπ/φ (n being an integer) in the radial direction 11 . One can thus infer from the peak separation the lateral size of the disk. Using the literature [START_REF] Sparks | Ferromagnetic relaxation theory[END_REF][START_REF] Damon | [END_REF] value for the YIG exchange length Λ ex = 15 nm, a fit of the peak separation leads to an effective confinement of respectively 700, 520, 380 nm for our 3 disks, assuming total pinning at the disk edge (see Table II). This is in very good agreement with nominal sizes targeted by the patterning. Differences with the nominal value are due both to imperfection of the lithographic process and the dipolar pinning condition 14 , which scales as the aspect ratio. Confining a spin wave in a smaller volume leads also to an overall increase of the exchange and self-dipolar energy, and thus a shift of the fundamental mode with increasing energy. We have checked that the position of the peak is compatible with the magnetic parameters shown is Table I assuming that the center of the probe is approached from 3.0 to 1.5 µm when moved from the largest to the smallest disk.

From the comparative measurements presented in Fig. 2, we see that the peak position is not altered by the addition of the Pt layer ontop of the YIG disk whereas its linewidth is clearly increased. To emphasize this result, Fig. 3 shows the linear f-MRFM spectra of the 700 nm bare YIG (b, blue dots) and of the YIG|Pt (c, red dots) disks. For comparison the spectrum measured on the extended YIG thin film is shown in (a, black). It was recorded with a transmission line Au-antenna of 500 µm width and an in-plane external field oriented perpendicularly to the exciting microwave field, similarly to the configuration used in Ref. 3 .

A first striking result of Fig. 3 is obtained by comparing the spectra measured in extended film (a) and the nanostructure (b) for the bare YIG: the nanopatterning improves the linewidth 7 . One can observe that, while the lineshape of the resonance has a Lorentzian shape in the nanostruture (continuous line), the peak shape is asymetric in the YIG film. We attribute it to inhomogeneous broadening. This is confirmed by performing experiments on different slabs of the reference film, which actually yield different lineshapes (not shown). The usual method of separating the homogeneous contribution from the inhomogeneous one is to study the frequency dependence of the half-linewidth ∆f /2. The slope gives α, the Gilbert damping constant, while the zero frequency intercept gives the inhomogeneous contribution. In Fig. 3(d), we have thus plotted the full-linewidth ∆f of the extended reference film at different frequencies using the same color as in Fig. 3(a). The linewidth varies almost linearly with the frequency. We extract from the slope the damping 4• 10 -4 , while the intercept at zero frequency indicates the amount of inhomogeneity of the resonance: ∆f =2.5 MHz (or ∆H=1 Oe).

On the same Fig. 3(d), we show using red dots the frequency dependence of the linewidth of the YIG|Pt nanostructure. The fit of the slope yields α YIG|Pt = 13 • 10 -4 . The striking feature is that a linear fit now intercepts with the origin of coordinates. It means that the linewidth measured in the nanostructures directly yields the homogeneous contribution. For the bare YIG nanodisk it was only possible to reliably extract the linear linewidth at one point at 8.2 GHz (blue dot). Interestingly enough the slope of the straight line from this point to the origin is exactly that of the line fitted to the extended film data. We speculate by analogy to the Pt/YIG case that the linewidth in the YIG nanostructure is purely homogeneous in nature, while the intrinsic part of the damping has been unaffected by the lithographic process.

In Fig. 3, we also see that the linewidth of the 700 nm YIG|Pt disk is 22 MHz (c) i.e. about three times wider than that of the 700 nm YIG disk, which is 7 MHz (b). The influence of an adjacent Pt layer is shown to increase the damping threefold through spin-pumping effect 15 . The characteristic parameter for the efficiency of spin transfer across the interface and the accompanying increase of damping in YIG is the spin mixing conductance G ↑↓ . One can directly evaluate the increase of damping induced by the presence of the Pt layer. As we find α YIG = 4•10 -4 and α YIG|Pt = 13•10 -4 , we deduce a large spin-pumping contribution α sp = 9 • 10 -4 that adds to the intrinsic damping : α YIG|Pt = α YIG + α sp . From this, one can calculate the spin mixing conductance according to 16 :

G ↑↓ = α sp 4πM S t gµ B G 0 , (1)
where t = 20 nm is the YIG thickness, g the electron Landé factor, µ B the Bohr magneton, and G 0 = 2e 2 /h the quantum of conductance. This corresponds to: G ↑↓ = 1.55 • 10 14 Ω -1 m -2 for our YIG|Pt interface.

In summary, we have conducted a study of spin wave spectra of individual submicron YIG and YIG|Pt disks. We find that the litography process does not broaden the linewidth and on the contrary, the linewitdh decreases compared to the extended film. The influence of an adjacent Pt layer on the YIG through the spin pumping effect is investigated and quantified to increase the damping 3 fold. As a non zero spin mixing conductance is determined, these experiments pave the way for observation of inverse spin Hall effects in a YIG|Pt nanodisk accessing its individual spin wave modes. Most importantly, thanks to the small volume and purely intrinsic damping of the YIG sample, we will address in future studies the influence of direct spin Hall effect on the linewidth of these YIG nano-disks as it was demonstrated in allmetallic NiFe/Pt dots 17 .

FIG. 1 .

 1 FIG. 1. Optical image of the YIG disks, with (red) and without (blue) Pt on top, placed below a microwave antenna. The direction of the microwave magnetic field is indicated by a blue arrow. The external bias magnetic field is oriented perpendicularly to the surface.

FIG. 2 .

 2 FIG. 2. Mechanical-FMR spectra at H0=4.99 kOe of the 700, 500 and 300 nm diameter YIG disks arranged in rows by decreasing lateral size. The spectra of the pure YIG disks are shown in the left column (blue), while the ones covered with 13 nm of Pt are shown in the right column (red).

FIG. 3 .

 3 FIG. 3. (a) Low-power spectrum of the reference YIG film recorded at f0 = 8.2 GHz. Lineshape of the uniform mode measured at H0 = 4.99 kOe on the 700 nm disk bare (b, blue circles) and with (c, red circles) Pt. (d) Dependence of the linewidth on the resonance frequency.

TABLE I .

 I Magnetic parameters of the 20 nm thick YIG film.

	4πMs (G)	α	Λex (nm)	γ (rad.s -1 .G -1 )
	2.1 • 10 3	4 • 10 -4	15	1.79 • 10 7

TABLE II .

 II Comparative table of the measured and predicted (analytical and SpinFlow simulations)

  10 .

	φ (nm)	mode	fexp (GHz)	simu.	analyt.	fit
	700	(00)	8.74	8.54	8.69	8.76
	700	(01)	9.14	9.04	9.10	9.16
	700	(02)	9.65	9.66	9.68	9.72
	500	(00)	9.10	8.73	8.90	9.05
	500	(01)	9.72	9.51	9.63	9.69
	500	(02)	10.52	10.51	10.72	10.65
	300	(00)	9.58	9.19	9.47	9.51
	300	(01)	10.57	10.66	11.25	10.63

ACKNOWLEDGMENTS

This research was supported by the French Grants Trinidad (ASTRID 2012 program) and by the RTRA Triangle de la Physique grant Spinoscopy. We also acknowledge usefull contributions from C. Deranlot, A.H. Molpeceres and R. Lebourgeois.