Eduardo De

Los Santos

Ana Maria Alonso Rodríguez

Francesca Rapetti

Construction of a spanning tree for high-order edge elements

Keywords: High-order edge elements, tree-cotree gauge, small simplices' graph, incidence matrix

The well-known tree-cotree gauging method for low-order edge finite elements is extended to high-order approximations within the first family of Nédélec finite element spaces. The key point in this method is the identification of degrees of freedom for edge and nodal finite element spaces such that the matrix of the gradient operator is the transposed of the all-node incidence matrix of a directed graph. This is straightforward for low-order finite elements and it can be proved that it is still possible in the high-order case using either moments or weights as degrees of freedom. In the case of weights, the geometrical realization of the graph associated with the gradient operator is very natural. We recall in details the definition of weights and present an algorithm for the construction of a spanning tree of this graph. The starting point of the algorithm is a spanning tree of the graph given by vertices and edges of the mesh (the so-called global spanning tree, that is the one used in the low-order case). This global step, interpreted in the high-order sense, is enriched locally, with a loop over the elements of the mesh, with arcs corresponding to edge, face and volume degrees of freedom required in the high-order case.

INTRODUCTION

The aim of this work is to present the construction of a spanning tree of the graph associated with the degrees of freedom (dofs) of high-order nodal and edge finite elements to extend the well-known tree-cotree gauging method. The idea of the tree-cotree gauge is to identify a set of dofs for edge finite elements that spans the kernel of the curl operator. If the domain is simply connected, these dofs are identified on a spanning tree of the graph whose nodes (arcs) are the dofs of nodal (edge) elements, respectively. Using the weights [START_REF] Bossavit | Generating Whitney forms of polynomial degree one and higher[END_REF][START_REF] Rapetti | Geometrical localisation of the degrees of freedom for Whitney elements of higher order[END_REF][START_REF] Rapetti | Whitney forms of higher degree[END_REF] as dofs for high-order edge and nodal finite elements, the graph can be thought on a particular refinement (associated with the principal lattice ) of the elements of the initial mesh. The small tetrahedra generated by this (virtual) refinement are affine contractions of those of the initial mesh and are called small elements. They are not constructed actually (and in this sense they are virtual), however they allow for a geometrical localization of the previously referred graph and thus for a suitable visualization of the dofs for high-order nodal and edge elements for any polynomial degree. The construction of a spanning tree of this graph is not more involved than the construction of a spanning tree of the well-known vertex-edge graph of the original mesh. Indeed, the visualization o the new selected arcs, namely those that do not result from the splitting of edges of the mesh, is the same in all the faces and at the interior of all the tetrahedra of the mesh. If the domain is not simply connected, the kernel of the curl operator is bigger than the image of the gradient operator. The spanning tree has to be enriched of one edge (called, the fastener) for each generator of the first homology group, thus resulting on a "belted" spanning tree [START_REF] Kettunen | Formulation of the eddy current problem in multiply connected regions in terms of h[END_REF][START_REF] Bossavit | Computational Electromagnetism[END_REF] .

Whitney edge elements [START_REF] Bossavit | Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism[END_REF] are widely recognized as a useful tool for computational electromagnetism [START_REF] Bossavit | Computational Electromagnetism[END_REF] . They offer a simple construction of polynomial differential forms on simplicial complexes [START_REF] Nédélec | Mixed finite elements in R 3[END_REF][START_REF] Barton | New vector finite elements for three-dimensional magnetic field computation[END_REF] . Their associated dofs have a very clear physical meaning (circulation of the approximated field on the edges of the element) [START_REF] Bossavit | Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism[END_REF] . Because they model the curl kernel well, the matrix associated with the curl-curl operator is singular, and this singularity is related to the graph whose geometrical realization is the finite element mesh. Specifically when the domain is simply connected, a tree-cotree decomposition of the arcs of this graph is used to classify dofs into dependent and independent ones for the curl operator [START_REF] Albanese | Integral formulation for 3D eddy-current computation using edge elements[END_REF] . In fact, a vector potential with fixed values of its independent dofs (identified by means of the tree) is unique [START_REF] Albanese | Magnetostatic field computations in terms of two-component vector potentials[END_REF][START_REF] Manges | A generalized tree-cotree gauge for magnetic field computation[END_REF] . This was successfully applied in a three-dimensional magnetostatic problem [START_REF] Webb | A single scalar potential method for 3D magnetostatics using edge elements[END_REF] and also in a nonlinear eddy current problem [START_REF] Albanese | On the numerical solution of the nonlinear three-dimensional eddy current problem[END_REF] .

In the development of high-order extensions of Whitney edge elements, the classical approach [START_REF] Nédélec | Mixed finite elements in R 3[END_REF][START_REF] Nédélec | A new family of mixed finite elements in R 3[END_REF][START_REF] Monk | Finite element methods for Maxwell's equations[END_REF] introduces moments to define the needed dofs. Such a construction includes dofs associated with faces and tetrahedra that can not be interpreted as field circulations along edges and consequently it does not allow for an immediate geometrical realization of its associated graph. The approach proposed in [START_REF] Bossavit | Generating Whitney forms of polynomial degree one and higher[END_REF][START_REF] Rapetti | Comments on a High-Order Whitney Complex for Simplices[END_REF][START_REF] Rapetti | High order edge elements on simplicial meshes[END_REF] generalizes to higher orders the methodology of Whitney: if one finds a way to represent a manifold of dimension k by a k-chain (i.e., formal sum of mesh simplices of dimension k) then by duality one knows how to represent a field by a finite weighted sum of basis fields (the Whitney forms).

In this work, we follow this approach that yields new dofs for high-order approximations, the so-called weights [START_REF] Rapetti | Whitney forms of higher degree[END_REF][START_REF] Rapetti | Weights computation for simplicial Whitney forms of degree one[END_REF][START_REF] Rodríguez | Some remarks on spanning families and weights for high order Whitney spaces on simplices[END_REF] . In the case of nodal elements, they are the classical functionals that associate with a function its values at the small vertices, that are the nodes of the principal lattices of the mesh elements. In the case of edge elements, they are functionals that associate with a field its circulations along the small edges, that are the edges of the small tetrahedra [START_REF] Rapetti | Geometrical localisation of the degrees of freedom for Whitney elements of higher order[END_REF][START_REF] Rapetti | High order edge elements on simplicial meshes[END_REF] . These weights are relevant, since they pop up naturally when we generalize Whitney construction to the high-order case. Note that, each small edge can be localized either (i) on an edge of the mesh, or (ii) at the interior of a face of the mesh or (iii) at the interior of an element of the mesh. The small edges for high-order extensions of Whitney 1-forms are always linking adjacent points of the principal lattice of the same order. They have the same orientation and are parallel to edges of the mesh. When dofs are weights, the matrix representative of the gradient is an incidence matrix (between small edges and small vertices) for any polynomial degree [START_REF] Rodríguez | The discrete relations between fields and potentials with high order Whitney forms[END_REF] . On the other hand, the duality property between unisolvent dofs (circulations on small edges) and generators of the finite element space is lost in the high-order case if all the edges of the small tetrahedra are considered. Duality is recovered by choosing an appropriate subset [START_REF] Rodríguez | Minimal sets of unisolvent weights for high order Whitney forms on simplices[END_REF] of the small edges. The resulting graph, with the principal lattices of the mesh elements as nodes and this subset of small edges as arcs, is an oriented connected graph. Let us consider the canonical basis of the nodal and edge finite element spaces for the weights supported on the nodes and arcs of this graph respectively. In this context, we are interested in the matrix that represents the gradient as a linear operator between finite dimensional spaces. This means, the matrix that computes the circulations along the small edges of a field that is the gradient of a nodal finite element function, knowing the values of this function at the nodes of the principal lattices of the mesh elements. It is easy to see that this matrix is the transpose of the all-node incidence matrix of the graph. In the simply connected case, the gradients of the nodal finite elements fill-in the kernel of the curl operator acting on edge finite element functions. Then the selection of a spanning tree of this graph is a useful tool for the identification of high-order potentials [START_REF] Rapetti | High Order Whitney Forms on Simplices and the Question of Potentials[END_REF] or for the tree-cotree gauge within high-order approximations. In this contribution we propose an efficient algorithm for the construction of a spanning tree of this graph.

The outline of this work is as follows. In Section 2 we introduce the fundamental notion of small simplices. In Section 3 we explain how to use this notion to define basis and dofs for high-order Whitney finite elements, and we introduce the high-order graph induced by the supports of these dofs. We then relate this graph with the matrix that represents the gradient operator on the canonical bases. Finally in Section 4 we present the algorithm for the construction of a spanning tree of the high-order graph and some considerations on its particularities.

SMALL SIMPLICES

Let T = (V, E, F, T ) be a tetrahedral mesh of the polyhedral domain Ω of R 3 , where V is the set of vertices, E is the set of edges, F is the set of faces, and T is the set of tetrahedra of T . Let ∆ d (T ) denote the set of d-simplices of the mesh, for For S ∈ ∆ d (T ), we denote by ∆ (S) the set of -simplices contained in S, for ∈ {0, . . . , d}. For each Σ ∈ ∆ (S) with ∈ {0, . . . , d -1}, if S = [v m S (0) , . . . , vm S(d) ] and Σ = [vm Σ(0) , . . . , vm Σ( ) ], then there exists a unique increasing map m S Σ : {0, . . . , } → {0, . . . , d} such that, for m S (0) (x), . . . , λ S m S (d) (x)) the vector of its barycentric coordinates with respect to the vertices {v m S (0) , . . . , vm S(d) } of S. Each function λ S m S ( ) : S → R satisfies the property λ S m S ( ) (v m S (j) ) = δ j, , with , j ∈ {0, . . . , d}. The functions λ S m S ( ) are positive at the interior of S, and non negative on S. The barycentric coordinates also satisfy the properties of representation of points and partition of unity, namely

d ∈ {0, 1, 2, 3}, thus ∆ 0 (T ) = V , ∆ 1 (T ) = E, ∆ 2 (T ) = F ,
each i ∈ {0, . . . , }, m Σ (i) = m S (m S Σ (i)).
x = d j=0 λ S m S (j) (x)v m S (j) and d j=0 λ S m S (j) (x) = 1 . Let t ∈ T , for each i ∈ {1, . . . , n V } we denote λ i : Ω → R λ i : x → λ i (x) =    λ t m t ( ) (x) if x ∈ t and i = mt( ) with ∈ {0, . . . , 3} 0
otherwise.

The set of real valued polynomials defined on S of degree less than or equal to r is denoted by Pr (S) and the subspace of Pr (S) of homogeneous polynomials of degree r is denoted by Pr (S). We have dimPr (S) = r+d . When S is reduced to a single point, then Pr (S) = R for all r ≥ 0. For any S, if r < 0, then Pr (S) = {0}. It is known that the barycentric coordinate functions λ S m S (i) constitute a basis for P 1 (S). Given S ∈ ∆ d (T ) and α ∈ I(d + 1, r), we define the continuous function

λ α S : x → λ α S (x) = Π d j=0 λ S m S (j) (x) α j .
It holds that Pr (S) = Span α∈I(d+1,r) {λ α S }.

Definition 1 (Principal lattice of order r + 1). Given S = [m S (0), . . . , m S (d)] and an integer r ≥ 0, the principal lattice of order r + 1 in the d-simplex S is the set of points

L r+1 (S) = x ∈ S : λ S m S (j) (x) ∈ 0, 1 r + 1 , 2 r + 1 , . . . , r r + 1 , 1 , 0 ≤ j ≤ d .
Equivalently,

L r+1 (S) =    x ∈ S : x = 1 |β| d j=0 β j v m S (j) , ∀ β ∈ I(d + 1, r + 1)    .
In what follows, to give a general definition (in the sense of the ambient dimension q) of the small -simplices (with 0 ≤ ≤ q), we consider a mesh T of a q-dimensional domain Ω, where S ∈ ∆q( T ) is a triangle (resp. tetrahedron) when q = 2 (resp. q = 3).

Definition 2 (Small -simplex). Let S ∈ ∆q( T ) and α ∈ I(q + 1, r), for an integer r ≥ 0. The small q-simplex {α, S} is the q-simplex with barycenter 1 r + 1

q j=0 1 q + 1 + α i v m S (j)
and 1/(r + 1)-homothetic to S. For any Σ ∈ ∆ (S), the small -simplex {α, Σ} is the -simplex which is parallel, equally oriented to Σ and belongs to the small q-simplex {α, S}.

It is worth noting that the small q-simplices do not pave the q-simplex S. In other words, depending on r, geometric shapes not homothetic to S, such as octahedra or reversed tetrahedra (resp. reversed triangles) for q = 3 (resp. q = 2), appear together with the small q-simplices in each S ∈ ∆q( T ) to completely pave S (as detailed in Table 1).

Number and type of shapes associated with

L r+1 (S)
in a tetrahedron S small edges TABLE 1 Number of small -simplices, for > 0, and other shapes in a -simplex S as a function of the polynomial degree r + 1, with r ≥ 0.

In Figure 1 is presented a triangle f = [n 0 , n 1 , n 2 ] together with its principal lattice L r+1 (f ) for r = 2. In grey we have the corresponding family of small simplices. In Figure 2 we have the small simplices associated with the principal lattice

L 3 (t) in a tetrahedron t = [n 0 , n 1 , n 2 , n 3 ].
The resulting decomposition of t includes 10 small tetrahedra, and also 4 octahedra and 1 reversed tetrahedron. On each face of the boundary of t we can see 6 small faces and 3 reversed triangles. Note that each edge of t is decomposed in 3 small edges.

The small simplices are of key importance in the generalization of Whitney finite elements to higher order. Indeed, they allow to make a list of the generators of any polynomial degree of the finite element spaces and constitute the support of new dofs, the weights, for the representation

n0 = {(2, 0, 0), [n0]} n1 = {(0, 2, 0), [n1]} n2 = {(0, 0, 2), [n2]} {(0, 1, 1), [n0]} {(0, 1, 1), [n1]} {(0, 1, 1), [n2]} {(0, 1, 1), [n0, n1, n2]} {(0, 1, 1), [n0, n2]} FIGURE 1 Principal lattice L 3 (f ) in the face f = [n 0 , n 1 , n 2 ],
with the visualization of all the small simplices. The label of the small simplex

{(0, 1, 1), [n 0 , n 1 , n 2 ]
} and of some of its sub-simplices are indicated for illustration.

FIGURE 2 (Taken from reference [25]). Principal lattice

L 3 (t) in the tetrahedron t = [n 0 , n 1 , n 2 , n 3 ],
with the visualization of all the small simplices.

Each face of the boundary of t is decomposed in 6 small faces and 3 reversed triangles (highlighted on the left with red solid lines). The decomposition of t includes 10 small tetrahedra, 1 reversed tetrahedron (highlighted at the center with red dashed lines) and 4 octahedra (one is highlighted on the right in green).

of fields in these spaces (together with a direct geometrical visualization of the tree-cotree structure that we discuss in the last part of the present work).

HIGH-ORDER WHITNEY EDGE ELEMENTS AND THE HIGH-ORDER GRAPH

For a more compact presentation, we rely on the notation of the finite element exterior calculus. It is well known (see [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF][START_REF] Rodríguez | Some remarks on spanning families and weights for high order Whitney spaces on simplices[END_REF] ) that the proxy fields of the space of trimmed polynomial differential forms (usually denoted by P - r+1 Λ k (T )) are the finite element space of Lagrange of degree r + 1 ≥ 1 if k = 0, the first family of Nédélec finite elements of order r + 1 conforming in H(curl; Ω) if k = 1, the first family of Nédélec finite elements of order r + 1 conforming in H(div; Ω) if k = 2, and discontinuous elements of degree ≤ r if k = 3. Let d be the exterior derivative for differential k-forms for k < d.

Definition 3 (Whitney k-forms).

In a tetrahedron t, Whitney k-forms associated with S ∈ ∆ k (t) are the fields

w S =    λ m S (0) , with k = 0, S = [n m S (0) ], k j=0 (-1) j λ m S (j) d w S-[n m S (j) ] , with k > 0, S = [n m S (0) , . . . , n m S (k) ].
In particular, for k = 1 the exterior derivative d corresponds to the gradient operator and Definition 3 yields

w S = λ m S (0) grad λ m S (1) -λ m S (1) grad λ m S (0) .
Whitney k-forms of higher degree in each tetrahedron t can be associated with the principal lattice in t as recalled in the next Definition.

Definition 4 (High-order Whitney k-forms). In a tetrahedron t, Whitney k-forms of higher polynomial degree are the fields

ws = λ α t w S ,
for all the small k-simplices s := {α, S}, with α ∈ I(4, r) for r > 0, and S ∈ ∆ k (t), being w S the Whitney k-form associated with S. The space of high-order Whitney k-forms in t is

P - r+1 Λ k (t) = Span {ws : s = {α, S}, α ∈ I(4, r), S ∈ ∆ k (t)} .
Fields in P - r+1 Λ k (t) are defined as products between a field in P - 1 Λ k (t) and the continuous function λ α t , with α ∈ I(4, r). Therefore, fields in P - r+1 Λ k (t) enjoy the same conformity properties as those in P - 1 Λ k (t) [START_REF] Rapetti | Whitney forms of higher degree[END_REF] .

Definition 5 (Weights on small k-simplices). For a polynomial k-form u ∈ P - r+1 Λ k (t), with 0 ≤ k ≤ 3 and r ≥ 0, the weights are the real numbers

{α,S} u, (1)
for all the small k-simplices {α, S}, with α ∈ I(4, r) and S ∈ ∆ k (t).

The weights of a high-order Whitney k-form u ∈ P - r+1 Λ k (t) along all the small k-simplices {α, S} of a mesh T are unisolvent 25, Proposition 3.14 , namely if u ∈ P - r+1 Λ k (t) is such that {α,S} u = 0 for all the small k-simplices s = {α, S}, then u = 0. The small k-simplices can thus support the dofs for fields u ∈ P - r+1 Λ k (t). Since the result on unisolvence holds true also by replacing t with Σ ∈ ∆ (t) (for k ≤ < 3), then the respective trace of u over any Σ ∈ ∆ (t) (trace that lives in P - r+1 Λ k (Σ)) is uniquely determined by the weights localized on small k-simplices in Σ. We thus can use the weights on the small k-simplices {α, S} as dofs for the fields in the finite element space P - r+1 Λ k (T ), being aware that their number is greater than the dimension of the space.

The weights given in Definition 5 have a meaning as cochains and this relates directly the matrix describing the exterior derivative with the matrix of the boundary operator. The key point is the Stokes theorem C du = ∂C u, where u is a (k -1)-form and C a k-chain. More precisely,

if u ∈ P - r+1 Λ k-1 (T ) then du ∈ P - r+1 Λ k (T ) and {α,S} du = ∂{α,S} u = {β,Σ} B {α,S},{β,Σ} {β,Σ} u, (2)
being B the boundary matrix with as many rows as small k-simplices, and as many columns as small (k -1)-simplices. The small k-simplex {α, S} inherits the orientation of the k-simplex S, so the coefficient B {α,S},{β,Σ} is equal to the coefficient B S,Σ of the boundary of the simplex S if β = α and zero otherwise. Therefore, B S,Σ is 1 or -1 depending of whether the orientation of Σ matches or not with the orientation induced by that of the k-simplex S on its (k -1)-face Σ, and B S,Σ is zero when Σ is not part of the boundary of S. This fact is straightforward if dim(Σ) > 0, but also when dim(Σ) = 0 21 . For k = 1, note that the matrix B {α,S},{β,Σ} is an incidence matrix between small edges and small vertices; moreover, it is an all-node incidence matrix of a (non connected) graph, whose geometrical realization is shown on the right of Figure 3.

Remark: For k = 1, in order to have B an incidence matrix between small edges and small vertices, some points in the principal lattice have different labels. The point

1 |η| 3 j=0 1 4 + η j v m t (j) ∈ L 3 (t)
, with η = (2, 0, 1, 0) has a label as the small vertex {(2, 0, 0, 0), [n 2 ]} (the brown point in the fragmented visualization) and also as the small vertex {(1, 0, 1, 0), [n 0 ]} (the black point in the fragmented visualization). However, the 0-forms associated with the different labels are the same. Indeed, {(2, 0, 0, 0), [n 2 ]} corresponds with the 0-form λ 2 m t (0) wn 2 = λ 2 n 0 λn 2 and {(1, 0, 1, 0), [n 0 ]} with the 0-form λ m t (0) λ m t (2) wn 0 = λn 0 λn 2 λn 0 = λ 2 n 0 λn 2 , and we can see that "both" 0-forms are actually the same.

FIGURE 3 (Taken from reference [21]

). Visualization (on the left) of all the small simplices associated with the principal lattice L 3 (t) in the tetrahe-

dron t = [n 0 , n 1 , n 2 , n 3 ].
The same set of small simplices in a fragmented visualization (on the right), a small tetrahedron in blue, a small face in green, a small edge in red and small vertex in brown and also in black (note the different labels for the same point in the principal lattice L 3 (t)).

We say that a set of unisolvent degrees of freedom for high-order Whitney finite elements is minimal if its cardinality is equal to the dimension of the finite element space. In order to use a set of minimal unisolvent degrees of freedom for high-order Whitney edge finite elements, we have to disregard some small edges. In fact, it can be proved [START_REF] Rodríguez | Minimal sets of unisolvent weights for high order Whitney forms on simplices[END_REF] that to be unisolvent it is enough to consider on each face small edges parallel to two of the three edges of the face, and on each tetrahedron small edges parallel to three of the six edges of the tetrahedra. Proceeding in this way, the corresponding set of weights is unisolvent and minimal.

We use the term active to indicate all the small k-simplices such that λ α t w S belongs to a basis of P - r+1 Λ k (T ). The dimension of the space P - r+1 Λ k (T ) coincides with the number of active small k-simplices in the mesh T [START_REF] Rapetti | High Order Whitney Forms on Simplices and the Question of Potentials[END_REF] . We recall the following result 26, Proposition 6.4 : Proposition 1 (A basis of P - r+1 Λ k (t)). The set λ α t w S : α ∈ I(4, r), S ∈ ∆ k (t) and α j = 0 if j < m t S (0) is a basis of P - r+1 Λ k (t).

A possible rule for selecting the active small vertices {α, [n i ]} and set of active small edges {α, [n me(0) , n me(1) ]}, that induces minimal and unisolvent weights, is the same given in Proposition 1 to select a basis of P - r+1 Λ k (t). We denote by d L and d N the dimension of P - r+1 Λ 0 (T ) and P - r+1 Λ 1 (T ), respectively. We consider the graph G ho with nodes the set {q }

d L =1
of active small vertices (the points of the principal lattices of the mesh elements) and arcs the set {em}

d N m=1 of active small edges. Let {φ h,i } d L i=1
and {ω h,j } d N j=1 be the canonical bases of the nodal and edge finite element spaces for the weights supported on the nodes and arcs of this graph, respectively. This means that φ h,i (q ) = δ i, and em ω h,j = δ j,m , where δ•,• is the Kronecker delta. For each active small edge em, we denote by q m,ini and q m,f in its initial and final small nodes, respectively. It is well known that if the function ψ h belongs to the space of nodal Lagrangian finite elements of degree r + 1, then grad ψ h belongs to the first family of Nédélec finite elements of degree r + 1, and clearly em grad ψ h = ψ h (q m,f in ) -ψ h (q m,ini ). Therefore when using these canonical bases, the matrix that gives the circulations along the active small edges of the gradient of a nodal finite element function (knowing the values of this function at the nodes of the principal lattices in each mesh element) is the transposed of the all-node incidence matrix of the graph G ho .

Its geometrical realization in a tetrahedron is shown (in a fragmented version for a friendly visualization) in Figure 4. It is worthwhile to mention that the set of active weights is proved to be unisolvent in [START_REF] Rodríguez | Minimal sets of unisolvent weights for high order Whitney forms on simplices[END_REF] . The structure shown in the Figure 4 has to be replicated on each tetrahedron t ∈ T .

We thus obtain the high-order graph, G ho , over which we can construct the high-order spanning tree, S ho . FIGURE 4 (Taken from reference [25]). A fragmented visualization of the geometrical realization of the connected graph in a tetrahedron with the active small edges (in red) as arcs and active small vertices as nodes.

ALGORITHM FOR THE CONSTRUCTION OF THE SPANNING TREE (THE HIGH-ORDER CASE)

Thanks to the additional information given by the geometrical localization of the connected graph G ho , formed by the active small vertices and the active small edges (shown in Figure 4), we propose in the Algorithm 2 a suitable set of rules to construct a spanning tree in the referred graph G ho for each tetrahedron; then we use this set of rules to construct a convenient global spanning tree S ho , giving in this way a classification rule in the set of active small edges of the mesh T . The Algorithm 1 starts by reading the arcs of the spanning tree S. We recall that, in a topologically trivial domain Ω, a spanning tree S can be obtained by applying algorithms such as the Breadth-First search or Depth-First search in the connected and directed graph G= (V, E). The arcs of G that are not in the spanning tree compose the cotree S c . For domains with no trivial topology, Algorithm 1 has to start from the arcs of a belted-spanning-tree S b , that can be obtained by applying Algortihm 3.

In Algorithm 1 the information on S is saved by mean of a flag associated with each edge of the mesh (if an edge is part of the spanning tree S, the flag is set to true, otherwise is set to false); when the algorithm is initialized, the flag associated with each face is also set to false (we use this flag to control the previously visited face in the loop over all elements of the mesh).

The information saved in the flags is used within a loop over all the tetrahedra of the mesh, where the spanning tree S ho is constructed; in each cycle of the tetrahedra loop, the previously proposed set of rules is used conveniently. The spanning tree S gives a selection criteria to put some active small edges of the global graph in the spanning tree S ho . More precisely, all the small edges on arcs of the spanning tree S and all but one (in particular the last one) the small edges on arcs of the cotree S c . Then the global spanning tree S ho is further enriched with some active small edges at the interior of each face (only those faces of the tetrahedra that have not been visited previously), and finally with some active small edges

Algorithm 1 Initialization

Require: The mesh T (V, E, F, T ), and the first Betti number g ≥ 0 of Ω.

Ensure: Logical tree-flag, visited-edge-flag and visited-face-flag arrays.

1: if g == 0 then 2:

Construct a spanning tree S * of the graph (V, E) by using a Breadth-First search.

3:

Set S = S * . If Ω is simply connected 4: else 5:
Construct a belted-spanning-tree S b of the graph (V, E) by using Algorithm 3.

If Ω is not simply connected visited-face-flag(f ) ← false.

17: end loop FIGURE 5 The spanning tree S ho (in green) on a tetrahedron t = [n 0 , n 1 , n 2 , n 3 ]. On the left, the subgraph corresponding to active small edges on the edges and in the interior of the face f = [n 0 , n 1 , n 2 ], and the corresponding part of the spanning tree S ho . The red arcs belong to the spanning tree S ho , if and only if, its corresponding edge e belongs to the spanning tree S. On the right, the active small edges in the interior of t belonging to S ho are visualized, (note that for a friendly visualization, some active small edges were omitted in the figure). In this example we are considering

r = 4.
at the interior of each tetrahedron. The proposed algorithm can be executed as part of pre-processing step. The algorithm requires as input the mesh T and the first Betti number g of the domain Ω; it gives as output the previously mentioned flags.

At a first sight, the Algorithm 2 seems to work on a graph G ho with much more nodes and arcs than G (the one in the low-order case). In reality the high-order graph G ho (and its incidence matrix) is not realized. We have used it only for visualization purposes. However, each arc in the high-order graph G ho corresponds with a dof (and hence with an element of the canonical basis). Therefore any action on the set of dofs can be "visualized" as an action on this fictive graph. In Algorithm 2 we use only the low-order graph and make a local choice of dofs, this is done once on a reference element (as in Figure 5) and then applied repeatedly on the elements of the mesh. The local choice has to guarantee that the merge of all local spanning trees is a spanning tree of the (virtual) high-order graph, as illustrated in Figure 6. Algorithm 2 implements the choice illustrated in the figures.

In Figure 5 we can see the structure of the spanning tree S ho (in green) given by the Algorithm 2, considering the case r = 4 in a mesh with a unique tetrahedron t = [n 0 , n 1 , n 2 , n 3 ]; a red arc is in S ho if its respective (big) edge is in the spanning tree S, otherwise is out of S ho .

In Figure 6 we can see the structure of the spanning tree S ho (in green and red) given by the Algorithm 2, considering the case r = 4 in a mesh with two tetrahedra. It is worth noting that the red arcs in this figure correspond to the red arcs in Figure 5 belonging to S ho .

Algorithm 2 A spanning tree of active small edges Require: The mesh T = (V, E, F, T ) and logical tree-flag, visited-edge-flag and visited-face-flag arrays.

Ensure: A (belted-) spanning tree S ho of the graph G ho with nodes: the active small vertices; and arcs: the active small edges. (1) = 0.

All the small edges {α, e} over the edge e. first is always intended in the reversed lexicographical order. Proof. For r = 0, the small edges coincide with the mesh edges and their number in S is (N V -1). Moreover, N V = d L , thus the property holds true for r = 0.

For r > 0, we recall that d L can be rewritten as

d L = N V + N E dim P r-1 (e) + N F dim P r-2 (f ) + N T dim P r-3 (t)
where N V (resp., N E , N F , N T ) is the number of nodes (resp., edges e, faces f , tetrahedra t) in the mesh over Ω. By the proposed algorithm, the number d S ho of active small edges in the spanning tree S ho can be counted as follows

d S ho = (r + 1) #S + r #S c + N F r-1 q=1 (r -q) + N T r-2 q=1 [ r-(q+1) s=1 r -(q + s) ] with cardinalities #S = (N V -1) and #S c = N E -(N V -1). Note that r-1 q=1 (r -q) = r r-1 q=1 1 - r-1 q=1 q = r(r -1) - r(r -1) 2 = r(r -1) 2 = dim Pr-2(f ).

FIGURE 6

The spanning tree S ho (in green and red) on a mesh T with two tetrahedra. The arcs of the spanning tree S in blue. It is worth noting that the red arcs in this figure correspond to the red arcs in Figure 5 belonging to S ho . In this example we are considering r = 4.

Moreover, for the last term in d S ho we have

r-2 q=1 [ r-q-1 s=1 (r -q -s) ] = r-2 q=1 [ (r -q) r-q-1 s=1 1 - r-q-1 s=1 s ] = 1 2 r-2 q=1 [ (r -q)(r -q -1) ] = 1 2 [ r-2 q=1 q 2 -(2r -1) r-2 q=1 q + r(r -1) r-2 q=1 1 ] = 1 2 [ 1 6 (r -2)(r -1)(2r -3) - 1 2 (2r -1)(r -1)(r -2) + r(r -1)(r -2) ] = 1 2 (r -1)(r -2)[ r 3 - 1 2 + 1 2 -r + r] = 1 6 r(r -1)(r -2)
that is indeed dim P r-3 (t). Summing up, we get

d S ho = N V -1 + N E dim P r-1 (e) + N F dim P r-2 (f ) + N T dim P r-3 (t) = d L -1,
thus the property.

Remark: If Ω is not simply connected, let g be the first Betti number of Ω, namely, the number of homologically independent non bounding cycles in Ω. The algorithm can be easily generalized to this case by replacing the global spanning tree, S, with a global belted-spanning-tree [START_REF] Kettunen | Discrete spaces for div and curl-free fields[END_REF][START_REF] Kettunen | Formulation of the eddy current problem in multiply connected regions in terms of h[END_REF] , S b that has g additional arcs. Hence, the number of arc in the high-order belted-spanning-tree

S b ho is d S b ho = d L -1 + g.
For more details on the computation of g and the construction of a belted-spanning-tree, see [START_REF] Hiptmair | Generators of H 1 (Γ h , Z) for triangulated surfaces: construction and classification[END_REF][START_REF] Rapetti | Discrete vector potentials for non-simply connected three-dimensional domains[END_REF][START_REF] Dlotko | Automatic generation of cuts on large-sized meshes for the T -Ω geometric eddy-current formulation[END_REF][START_REF] Rodríguez | Construction of a finite element basis of the first de Rham cohomology group and numerical solution of 3D magnetostatic problems[END_REF][START_REF] Pellikka | Homology and cohomology computation in finite element modeling[END_REF] . Note that Algorithm 2 has to start from the logical arrays tree-flag, visited-edgeflag and visited-face-flag filled in considering the belted-spanning-tree S b . The latter can be constructed as detailed in Algorithm 3. For example,

if Ω is an empty torus, we have g = 2 and a belted-spanning-tree is presented in Fig. 7 (right).

FIGURE 7

The mesh (left), the 2 cycles (center) and the belted-spanning-tree (right) computed from Algorithm 3.

Algorithm 3 Computing a belted-spanning-tree when Ω is not simply connected (g > 0)

Require: The mesh T = (V, E, F, T ).

Ensure: A belted-spanning-tree S b . 1: Compute the g homologically independent not bounding cycles γ 1 , γ 2 , ..., γg in T .

We assume that the g cycles γ 1 , γ 2 , ..., γg are mutually disjoint and without self-intersection. This can be done by relying on the Smith normal form in Z as explained in [START_REF] Rapetti | Discrete vector potentials for non-simply connected three-dimensional domains[END_REF][START_REF] Pellikka | Homology and cohomology computation in finite element modeling[END_REF] select one edge, say e i ∈ γ i , that does not belong to 1≤j≤i-1 γ j 5:

define γi = γ i \ {e i } 6: end for 7:

complete the acyclic sub-graph corresponding to g i=1 γi to get a spanning tree S of the graph (V, E)

8:
define the belted-spanning-tree S b = S g i=1 {e i }.

It has to be said that in the lowest order case, the conditioning of the linear system gets worse once the gauge is applied and some work has been done to try to avoid it [START_REF] Ren | Boundary edge elements and spanning tree technique in three-dimensional electromagnetic field computation[END_REF] . For the high-order case, some tests on three-dimensional eddy current problems are done in [START_REF] Ledger | hp-Finite element simulation of three-dimensional eddy current problems on multiply connected domains[END_REF] where the magnetic vector potential based formulation is gauged by regularization instead of by the imposition of dofs associated with the arcs of a (belted) spanning tree.

Some numerical tests will be done to see how badly the fact of gauging affects the "effective conditioning" of the system matrix in the high-order case, where effective conditioning means the ratio of eigenvalues of the matrix but not considering the (nearly) zero ones.

CONCLUSIONS

The weights, namely integrals over mesh sub-simplices, are natural dofs to adopt when the solution is a cochain on a simplicial mesh. They allow to extend straightforwardly the classical tree-cotree techniques from the low-order to the high-order case, since there is a correspondence of dofs with geometrical entities, the small simplices, that constitute a graph. In this work we have presented an algorithm for the construction of a (belted) spanning tree of this graph. The construction is based on a (belted) spanning tree of the graph given by vertices and edges of the mesh (the socalled global (belted) spanning tree). Then a loop by elements enriches this initial (belted) spanning tree with arcs corresponding to face dofs (only those faces of the tetrahedra that have not been visited previously) and with the arcs corresponding to the volume dofs. The enrichment in each element is the same for all elements. Gauging through such a tree-cotree approach (applied to the dual graph with the small triangles as vertices and the small edges in common as arcs) has been successfully applied in the frame of high-order reconstructions of a Darcy flow in a 2D domain starting from its weights on the small edges [START_REF] Rodríguez | New degrees of freedom for high-order Whitney approximations of Darcy's flows[END_REF] . Further numerical investigations are however necessary to analyze how this gauging influences the conditioning of the final algebraic system associated with an electromagnetic problem.

  and ∆ 3 (T ) = T . Let us fix an orientation on each edge, face, and tetrahedron of T . This can be done by choosing a total ordering of the vertices in V = {v i } N V i=1 and by associating with each d-simplex of the mesh, S ∈ ∆ d (T ), an increasing function m S : {0, 1, . . . , d} → {1, . . . , N V }. The oriented d-simplex S is hence given by S = [v m S (0) , . . . , vm S(d) ].

  For r, d ∈ N, let us set I(d + 1, r) := {η = (η 0 , . . . , η d ) ∈ N d+1 : |η| = r} , being |η| := d i=0 η i . The cardinality of I(d + 1, r) is equal to r+d d = (r+d)! r!d! . For a point x in a d-simplex S ∈ ∆ d (T ) (with d ∈ {1, 2, 3}), we denote by λ S (x) = (λ S

d

  and dim Pr (S) = r+d-1 d

6 :Set S = S b . 7 : end if 8 : 15 :

 67815 loop over the edges e ∈ E loop over the faces f ∈ F 16:

8 : else 9 :

 89 add to S ho the small edges {α, e} for the r -1 first indices α ∈ I(4, r) with α m t

  .

2: for i ← 1 to g do 3 :

 3 consider the i-th disjoint cycle γ i = e∈E c i e e, with c i e ∈ {-1, 0, 1} 4:

1 :

 1 loop over the tetrahedra t ← [na, n b , nc, n d ] ∈ T 2: loop over the faces f ← [n i , n j , n k ] ∈ ∆ 2 (t)

	3:	if visited-face-flag(f ) == false then
	4:	loop over the edges e ← [np, nq] ∈ ∆ 1 (f )
	5:	if visited-edge-flag(e) == false then
	6:	if tree-flag(e) == true then
	7:	add to S ho the small edges {α, e}, ∀ α ∈ I(4, r) with α m t t-e	(0) = 0 and α m t t-e

  S ho the small edges {α,[n m f (0) , n m f (1) ]} for the r -α m t f (2) first indices α ∈ I(4, r) with α m t Adding small edges {α, [n m f (0) , n m f (1)]} of the interior of f to S ho . ho the small edges {α, [n m f (0) , n m f (1) ]} for the r -(α 2 + α 3 ) first indices α ∈ I(4, r). small edges {α, [n m t (0) , n m t (1) ]} of the interior of t to S ho . If Ω is simply connected, the number d S ho of small edges in S ho is (d L -1).

	10:	end if
	11:	visited-edge-flag(e) ← true
	12:	end if
	13:	end loop
	14:	for α m t f (2) ← 1 to r -1 do
	15:	add to t-f	(0) = 0.
	16:	end for
	17:	visited-edge-flag(f ) ← true
	18:	end if
	19:	end loop
	20:	for α 3 ← 1 to r -2 do
	21:	for α 2 ← 1 to r -(α 3 + 1) do
	22: add to S Adding 23: end for
	24:	end for
	25: end loop
	Proposition 2.

Abbreviations: Spanning tree for high-order edge elements

ACKNOWLEDGMENTS Financial disclosure

Financial support from the French MathIT program ANR-15-IDEX-01 and the Italian project PRIN 201752HKH8.

Conflict of interest

The authors declare no potential conflict of interests.


Philippe Lestaevel

Christelle Huet

Veronique Lejeune

Carole Moreno

Carmen Villagrasa

Juliette Feuardent

François Trompier

Véronique Lejeune

Cosmic radiation exposure of airline crews in France over the period 2015-2019

Keywords: aircraft crew exposure, occupational, dosimetry, solar cycle, gender

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

Introduction

Galactic cosmic rays are highly energetic particles

Ionizing radiation originating from various phenomena travels through space and therefore continuously bombards the Earth's atmosphere. Overall, all particles composing the cosmic rays (CR), mainly energetic charged particles (mostly proton and He nuclei, up to 10 21 eV) have sufficient energy to penetrate earth's magnetosphere and to produce secondary particle shower deep in the atmosphere until sea level. There are two main sources of energetic charged particles): (1) the everpresent galactic and extra galactic components that likely originate from explosions of supernovae and galactic nuclei, and (2) transient solar events, known as flares or coronal mass ejections (CME) whose radiative effects in our atmosphere usually last several hours. At 35 000 feet above sea level (about the cruising altitude of subsonic commercial aircrafts), the effective dose rate typically ranges from 2 to 10 µSv/h depending primarily on latitude, altitude, and the level of solar activity (ICRU 84, 2010). Because of the decreasing shielding from the magnetosphere (towards charged particles) with latitude, the dose rate increases with latitude up to 60°. The dose rate, because of atmosphere attenuation, decreases from about 20-25 km (altitude of the Regener-Pfotzer maximum corresponding to the maximal dose rate in atmosphere) to sea level. The impinging charged particle components of the galactic and extragalactic rays are modulated by the sun's magnetic field protecting our solar system according to sun activities. The intensity of solar activity waxes and wanes according to an average eleven-year cycle.

The regulation is transposed for airline flight crews

In 1990, the International Commission on Radiological Protection (ICRP) recommended that flying crews (FCs) be classified as occupationally exposed (ICRP, 1991). This ICRP recommendation was transposed into European regulations by the European Directive No. 96-29 EURATOM of May 13, 1996(European Commission, 1996), which substantially amended the standards governing the European health protection of workers against the effects of ionizing radiation. Concerning the protection of FCs, article 42 states that: 'Each Member State shall make arrangements for undertaking operating aircraft to take account of exposure to cosmic radiation of aircrew who are liable to be subject to exposure over 1 mSv per year'. This directive was transposed into French law by the order no. 2001-270 of March 28, 2001(Ordonnance, 2001). A French decree of December 2003 (J Off, 2003) states that: (1) the individual exposure of any aircrew likely to be exposed to an effective dose >1 mSv per year must be assessed, (2) the effective dose must be assessed for each flight using a calculation method based on flight data provided by the airlines (date, duration, location, etc.) and (3) the calculation method must take into account the influence of the normal solar activity and effects of solar energetic particle (SEP) events classified as ground-level events (GLE). The SIEVERT system was then developed by IRSN to provide a dosimetry service to French airlines and was the first world-wide operational system starting in 2001 [START_REF] Bottollier-Depois | An operational approach for aircraft crew dosimetry: the SIEVERT system[END_REF] and including also GLE dose evaluation [START_REF] Lantos | Semi-empirical model to calculate potential radiation exposure received on board airplane during solar particle events[END_REF]. The French airlines were in charge of the assessment of the aircrew effective dose from the dos e calculated for each flight with SIEVERT system and of their transmission to the French national dose registry. French regulations have evolved through the decree of July 17, 2013, applicable on July 1, 2014, which imposes a significant change in the system. Under this new regulation, effective doses must no longer be assigned to FCs by the airlines but by an approved organization which is also in charge of their transmission to the national dose registry. In this context, IRSN has carried out a complete redesign of the SIEVERT system to a new application, SIEVERTPN, to meet the new regulatory requirements.

The dose monitoring is performed by calculation method

Although monitoring of the occupational exposure of FCs could be possibly performed by measurement means such as personal passive dosimetry or in situ monitoring with fixed instrumentation on board, monitoring based on calculations is mandatory in French regulation. As a matter of fact, since the radiation field at aviation altitude is very complex in terms of particle composition and energy range, dose assessment based on measurement remains a challenging task and presents some drawbacks compared to calculation methods. Individual dosimeters, in addition to the costs, can be forgotten or even lost by FCs, or frequently exposed to X-rays luggage control system inducing either over or under dose evaluation. Dosimeters are calibrated against an imperfect estimator of the effective dose, either the personal or ambient dose equivalent. For all that reasons, in countries applying the recommendation on FC dose monitoring, dose monitoring is performed by calculation. Effective doses are usually calculated monthly using computer codes that take into account flight parameters such as altitude, latitude, and solar activity. Several web-based platforms, such as PANDOCA, AVIDOS, PCAIRE or SIEVERTPN, have been developed in different countries for dosimetry services. These dose calculator systems use a 3D grid of the atmosphere in terms of effective dose rate which allows, knowing the flight route, to estimate the integrated dose on a given flight based. The 3D grid is either generated by simulation using Monte Carlo particle transport codes (Fluka, PHITS, PLANETOCOSMICS), an analytical calculation code such as LUIN or PLOTINUS or using a database of measurements converted into effective dose (FDOScalc and PCAIRE). In those codes, the geomagnetic shielding is evaluated through the effective vertical cut-off rigidity geomagnetic either calculated with codes such as PLANETOCOSMICS or from literature data. The solar activity modulation is assessed with heliocentric potential or ground neutron monitor data. Additional models are used to take into account the solar activity. All these codes have been validated against large set of measurements performed with a tissue-equivalent proportional counter (TEPC) [START_REF] Unscear ; Wissmann | A mobile TEPC-based system to measure the contributions to H*(10) at flight altitudes[END_REF][START_REF] Bottollier-Depois | Exposure of aircraft crew to cosmic radiation: on-board intercomparison of various dosemeters[END_REF] and also through code intercomparison exercises [START_REF] Eurados | Comparison of Codes Assessing Radiation Exposure of Aircraft Crew due to Galactic Cosmic Radiation[END_REF][START_REF] Eurados | Comparison of Codes Assessing Radiation Exposure at Aviation Altitudes in Case of Solar Particle Events report[END_REF]. The relative differences between measured and calculated values were mostly found to be within ± 30% [START_REF] Bottollier-Depois | Exposure of aircraft crew to cosmic radiation: on-board intercomparison of various dosemeters[END_REF], Lillok et al., 2007). The overall agreement between the codes for the galactic and extra galactic components was found to be better than 20 % from the median [START_REF] Eurados | Comparison of Codes Assessing Radiation Exposure of Aircraft Crew due to Galactic Cosmic Radiation[END_REF]. The agreement between the codes is considered quite satisfactory; dose estimates in radiation protection generally have uncertainties from  20 % to  30 %. For GLE dose estimation, larger discrepancies were found among codes, up to 70 %, with median values of the calculated route doses for H*(10) of 161 µSv and 123 µSv [START_REF] Eurados | Comparison of Codes Assessing Radiation Exposure at Aviation Altitudes in Case of Solar Particle Events report[END_REF]. This paper describes the cosmic radiation exposure of aircrew in French airline companies over the period 2015-2021 (i.e., most of the decreasing phase of the solar activity in solar cycle 24) obtained from SIEVERTPN (System of information and evaluation per flight of exposure to cosmic radiation in air transport) calculations and recorded in SISERI (Information system for monitoring occupational exposure to ionizing radiation: French national dose registry). In addition, possible causes of exposure increase during the studied period are analyzed.

Materials and methods

In France, aircrew members are subject to occupational radiation protection surveillance if (i) they are employed by an airline company under the French law, and (ii) they are likely to exceed the effective dose of 1 mSv per year. The obligation to assess the doses received by personnel lies with the employer. The airline must also take measures to keep the doses of the personnel low by appropriate flight planning and crew allocation. The studied population consisted of flight personnel working in France for an airline company with at least one effective dose registered in SISERI between 2015 and 2021. The effective doses were estimated using SIEVERTPN. As introduced above, SIEVERTPN is a calculation tool that allows assessing the doses due to cosmic radiations received during flights according to the routes taken that are made available on a server dedicated to airlines access. In short, the company prepares monthly two files: one containing the information of flights completed and one listing the presence of aircrew members on these flights. Each month, they are sent to the SIEVERTPN application via FTP or by download. The doses are evaluated for each flight according to flight characteristics. The calculation of the flight dose is based on a 3D grid of effective dose rates of cosmic radiation up to an altitude of 80 000 ft computed by EPCARD software [START_REF] Schraube | EPCARD-European program package for the calculation of aviation route doses[END_REF]. In SIEVERTPN, airspace is divided into cells of altitude, longitude, and latitude to form a 3D grid consisting of a mesh of 265 000 cells. Each cell is 1000 feet in altitude, 10° in longitude and 2° in latitude. An effective dose rate value is monthly assigned to each cell. The SIEVERTPN application evaluates the time spent by the aircraft in each cell and deduces the dose received. In the event of a significant solar flare, a specific grid is created using the SiGLE model [START_REF] Lantos | Semi-empirical model to calculate potential radiation exposure received on board airplane during solar particle events[END_REF]. Once the flight file is processed by SIEVERTPN, the monthly effective dose received by each crew member is calculated thanks to the information on their presence on the flights. Then, a file containing the monthly effective dose of each crew member is automatically generated by SIEVERTPN and sent to SISERI. Over the 2015-2021 period, the data from 7 to 13 French airlines per year have been transmitted to SISERI. In the SISERI database, workers are identified by their company, their last name, their first name and by their registration number in the national identification directory of natural persons. Some additional information is available in SIEVERTPN and/or in SISERI, such as the number of flight hours, the categories of FCs, i.e., technical flight crews (TFCs) or commercial flight crews (CFCs), the occupation (airline technical aircrews, aircraft captains, stewards, flight attendants, cabin managers) or the type of aircraft.

Results

The number of monitored FCs increased over the 2015-2019 period and has decreased since 2020, probably due to the COVID-19 crisis (Fig. 1). Likewise, the average annual effective dose steadily increased between 2015 and 2019 (1.98 mSv to 2.19 mSv). In 2020, it decreased by 54% due to the decrease in air traffic linked to COVID-19. In 2021, the average annual effective dose slightly increased compared to 2020 (12%), due to improved health conditions and the re-increased of air traffic. In 2019, the FCs account for nearly 48% of the total occupational collective dose in France whereas "only" 31% in 2020 and 28% in 2021 [START_REF] Lestaevel | La radioprotection des travailleurs -Exposition professionnelle aux rayonnements ionisants en France : bilan[END_REF].

This study focuses on the most FCs exposed and is limited to the 2015-2019 period, to avoid biases linked to the decline in activity in the sector in 2020 and 2021 due to the COVID-19 pandemic.

Between 2015 and 2019, the average annual individual effective dose of the 100 most exposed FCs increased (5.2 mSv in 2019 versus 4.0 mSv in 2015, i.e., 30%); over the same period a slight increase in annual flight time was observed (753 hours in 2019 versus 740 hours in 2015, i.e., 1.7%) (Fig. 2).

In 2015, no FCs received an annual effective dose above 5 mSv. Over the period 2016 to 2019, 169 FCs received an effective dose above 5 mSv by year. Most of these FCs were technical flight crews (TFCs): 109 TFCs versus 60 commercial flight crew (CFCs). Between 2016 and 2019, the number of CFCs tends to increase year over year (Fig. 3).

The ratio of male to female FCs does not change over time (approximately equal to 1 over the 2015-2019 period) (data not shown). In contrast, the proportion of female staff who received a dose greater than 5 mSv changes over time. In 2018 and 2019, this proportion of female staff was 25% compared to 8% in 2017 and 0% in 2016 (Fig. 4).

Of the 109 TFCs who received an annual effective dose above 5 mSv between 2016 and 2019, 73 % are airline technical aircrews and 14% are aircraft commanders (Fig. 5A). In French airlines companies, all airline technical aircrews represent about 56 % of TFCs (data not shown). Of the 60 CFCs who received a dose above 5 mSv, 52% of them are stewards, 41% are flight attendants and 7% are cabin managers (Fig. 5B).

Over the period 2016 to 2019, most FCs with an effective dose above 5 mSv worked mainly on longhaul aircraft with flight profiles (altitude) and destinations from Europe to North America and Asia, resulting in higher exposure. To a lesser extent, they have worked on medium-haul airliners (Fig. 6).

Discussion

The annual exposure of FCs, is, in France, among the highest occupational exposures compared to workers in other activities such as nuclear power plant or industry. Air transport in France has been continuously increasing over the period of study, with a drastic drop in air traffic which will not be compensated by 2021 due to the worldwide health crisis linked to COVID-19 pandemic. Collective dose of French FCs represents nearly 48 % of total occupational collective dose in 2019. This result can be compared to those observed in other countries, where collective dose of FCs represented nearly 58 % of total collective dose in 2018 (data available on ESOREX platform: https://esorexplatform.org/charts/g2). Our results show that the average annual effective dose of French FCs (2.19 mSv in 2019) is between the standard permissible limits set for the public (1 mSv) and the limit set for occupationally exposed persons (20 mSv). The comparison with other countries is not straightforward for different reasons (geographical location and extent of the country, airline practices, commercial and tourist development of the country ...). Nevertheless, these results are in the same range that values from the literature. According to the responses to the UNSCEAR Occupational Exposure Survey from 22 countries, the average annual effective dose of FCs was 2.7 mSv for the period 2010-2014 (UNSCEAR, 2022). In a German cohort study, covering 12 981 aircrew members from 2004 until 2015, Wollschläger et al. showed that the average annual effective dose was 2.25 mSv (range 0.01 -6.39 mSv) [START_REF] Wollschläger | Estimated radiation exposure of German commercial airline cabin crew in the years 1960-2003 modeled using dose registry data for 2004-2015[END_REF]. In Canada, the survey of cosmic radiation exposure of Air Canada pilots in 2009 showed that 93.8% of the pilots received more than the public limit of 1 mSv (versus 83.2 % for French FCs in 2019) and 54% received between 2.5 and 3.5 mSv [START_REF] Bennett | A survey of the cosmic radiation exposure of Air Canada pilots during maximum galactic radiation conditions in 2009[END_REF]. However, due to the dip in the vertical cutoff rigidity over Canada, it places Canada in the highest area of exposure to cosmic radiation in the world.

In this study, different factors that could explain the increase in the effective dose received by the most exposed French FCs during the period 2015-2019 have been investigated, such as the solar cycle, changes in the staff flight time or the aircraft model and thus the route. Some explanations are given below.

First, the total time spent in flight could influence the rise of average annual effective dose observed in the present study. Air transport has been continuously increasing during these last few years in France (Statistiques du trafic aérien | Ministères Écologie Énergie Territoires (ecologie.gouv.fr)) as well as in other countries. In Ireland, during the years 2002-2005, an increase in the total number of flights per year was also observed, resulting in a 75% increase in the number of FCs receiving doses above 1 mSv [START_REF] Colgan | Individual and collective doses from cosmic radiation in Ireland[END_REF]. In the present study, we showed that the average number of flight hours increased slightly over 2015-2019 (Figure 2). Therefore, it can't explain the consequent effective dose increase.

Apart from the total time spent in flight, the solar cycle can influence the FCs individual annual dose. The magnitude of cosmic radiation getting to the Earth depends on the sun's activity at a particular period. The sun constantly ejects particles with an intensity which varies according to a cycle of about eleven years (solar cycle). The cosmic radiation reaching the Earth is less when solar activity is strong and vice versa. In the present study, the period, which was analyzed regarding the most FCs exposed, i.e., 2015-2019, covers the decreasing phase of solar cycle 24 (its maximum was in the year 2014 and its minimum in 2019). Consequently, the ambient dose rate from high altitude radiation increased during the monitoring period. The average annual effective dose of French aircraft crew personnel increased from 2015-2019 by 11%. Previously, for German aircrew, the average annual effective dose during the middle of the decreasing phase of solar cycle 23 ( 2004) was 1.9 mSv and increased during its minimum (2009) to 2.3 mSv [START_REF] Frasch | Radiation exposure of German aircraft crews under the impact of solar cycle 23 and airline business factors[END_REF] thus 21%. Also, high altitudes and latitudes are associated with higher dose values. The comparison of the effective route doses between 2004 and 2009 showed an increase from 25 to 30 µSv (20%) on the trans-equatorial route Frankfurt-Johannesburg, and from 50 to 78 µSv (56%) on the polar route Frankfurt-New York [START_REF] Frasch | Radiation exposure of German aircraft crews under the impact of solar cycle 23 and airline business factors[END_REF]. From SievertPN data, an increase from 48 µSv to 62 µSv (29%) is observed between 2014 and 2019 for the route Paris-New-York whereas it remains stable over the same period for the route Paris-Johannesburg (27 µSv).

It can also be noted that in absolute value, more TFCs received a dose greater than 5 mSv over the 2016 to 2019 period than CFCs. TFCs have a qualification for a given type of aircraft (for example Airbus A380 or Airbus A320, etc.). The TFCs having received a dose greater than 5 mSv are those working mainly on long-haul jumbo jets and, to a lesser extent, workers on medium-haul airliners. In 2018, Los Angeles, San Francisco, New York, Washington, Mexico City, Shanghai, Abidjan, and Johannesburg were the possible destinations on these jumbo jets from Paris. Unlike TCFs who only have a qualification for one type of aircraft (long-haul type, medium-haul type), CFCs fly on several types of aircraft and quite several of them are trained on large aircraft long-haul carriers. Thus, it is easier for the airline company to adapt the CFC flight planning. The differentiation between long-and mediumhaul flights is difficult since CFCs are distributed to all categories of flights. They perform about one or two long-haul flights per month (compared to about five flights per month for TFCs) (data not shown). Previous studies have demonstrated, for example in Finnish airline personnel, that the median cumulative individual estimate dose was 18.5 mSv among cockpit (i.e., TFCs) and 14.5 mSv among cabin crew (i.e., CFCs) from [START_REF] Oksanen | 2001 relative à la transposition de directives communautaires dans le domaine de la protection contre les rayonnements ionisants[END_REF]-2014(Tolibov et al., 2017)). Oksanen et al. obtained annual doses of 2.19 mSv (1.08-2.83 mSv) for pilots and 2,27 mSv (0.72-3.1 mSv) for cabin crew in Finland [START_REF] Oksanen | 2001 relative à la transposition de directives communautaires dans le domaine de la protection contre les rayonnements ionisants[END_REF]. Thus, the difference of dose exposure between TFCs and CFCs is not always clear, and it depends, on the one hand, on the destinations of the companies, and on the other hand, on the practices of the company.

The proportion of female workers was about 52 % for the whole French monitored FCs and the male/female ratio remained relatively unchanged over the period 2016 to 2019 (data not shown).

According to an ILO report, the gender ratio in the civil aviation has, on average, remained relatively also unchanged [START_REF] Seligson | Women and aviation: quality jobs, attraction and retention[END_REF]. Between 2016 and 2019, the whole French male FCs monitored had a higher average annual effective dose (2.20 mSv) than female (2.02 mSv) (data not shown). In Germany, over the period 2004-2009, male cabin crew had also a higher average annual effective dose (2.2 -2.4 mSv) than female cabin crew (1.9 -2.4 mSv) [START_REF] Frasch | Radiation exposure of German aircraft crews under the impact of solar cycle 23 and airline business factors[END_REF]. Wollschläger et al. show that female aircrew members' drop in average effective dose around age 35 may correspond to an increase in social or family responsibilities with potentially more short-haul flights or part-time employment [START_REF] Wollschläger | Estimated radiation exposure of German commercial airline cabin crew in the years 1960-2003 modeled using dose registry data for 2004-2015[END_REF]. In the light of these results, the annual doses of FCs vary significantly with gender, possibly due to shorter flight times for female. Despite this, our study indicated that the number of female FCs having received a dose greater than 5 mSv increased over the period 2016 to 2019. These female FCs are exclusively CFCs. In the United States, between 1992 and 1996, about 84% of flight attendants were female with an annual average effective dose of 0.36 mSv [START_REF] Grajewski | Measurement of dose equivalent distribution on board commercial jet aircraft[END_REF]. Among the 24 female FCs employed by a French airline with a dose above 5 mSv in 2019, nine performed flights, but not only, between France and Japan.

Limitations of the study

As mentioned, aircrew effective doses were estimated by calculation means. Thus, some parameters that may affect the assessment of the effective dose are not considered.

First, the shock wave effect associated with CME leading to a significant decrease in the continuous galactic and extra galactic components impinging Earth is not considered. This effect is known as Forbush decreases (FDs) [START_REF] Forbush | Cosmic-ray intensity variations during two solar cycles[END_REF]. The maximum observed decrease in counting rate on ground neutron monitors typically ranges from a few percent to 20 % and exceptionally up to 30 %, depending on the size, the strength, and the proximity of the CME to the Earth. The decrease occurs rather suddenly, within a few hours or less, but the return to the previous level takes days or even weeks. These types of events are more frequent and longer than solar energy particles events. As an example, during the solar cycle 23 (1997-2007), 123 FDs were detected against 16 GLE events (GLE 55 to GLE 70). The impact of FDs on individual exposure has not been yet evaluated. Regarding GLE, if the maximal dose rate can reach tens of mSv/h at the flight altitudes (as for example for the Feb. 1956 event, the most intense GLE of the 20th century), since 2000, only four events have been considered in SIEVERTPN dose evaluation leading at maximum to double the flight dose on a transatlantic flight. For the worst case GLE, a maximum effective dose was estimated to 1.2 mSv for US aircrew personnel [START_REF] Anderson | Flight attendant radiation dose from solar particle events[END_REF].

In addition, Monte Carlo simulations [START_REF] Ferrari | Evaluation of the influence of aircraft shielding on the aircrew exposure through an aircraft mathematical model[END_REF] and measurements (Kubancak et al., 2014) of the ambient dose equivalent H*(10) at different positions inside the aircraft have shown that it is up to 19% higher at the front (and back) of the cabin (where TFCs are located) than in the central part (where CFCs are located). These differences are due to the aircraft shielding, the luggage load, the amount of fuel, etc. As the location in the plane is not considered in SIEVERTPN calculations, it is therefore likely that the difference in dose between the two categories of personnel (TFC and CFC) is even more pronounced than the one observed in that study.

Expected future trends

Several factors lead us to expect a decrease in the doses received by French aircrew staff in the coming years. First, although air traffic has recovered from the drastic decline due to the pandemic, it still appears to be below 2019 levels (QlikView Printing (ecologie.gouv.fr)). Moreover, the 25 th eleven-year solar cycle began after a minimum of solar activity observed in December 2019. The peak activity of this new cycle is expected in July 2025. This means that the dose rates should decrease until 2025. In addition, the Covid-19 pandemic marked the end of the operation of the A380 by French airlines on which the doses were higher, even if other types of aircraft, such as 777, were able to take over. Finally, the SIEVERTPN system has been updated with the new ICRP 103 weighting coefficients in place of the ICRP 60 weighting coefficients from July 2021 (ICPR, 2007). The modification of the proton and neutron weighting coefficients leads, under identical conditions to cosmic radiation exposure, to a decrease in the effective dose value of the order of 20 to 30% depending on the flight characteristics (Meier and Matthia, 2019).

However, it should be kept in mind that unexpected events (as was the case with the Covid pandemic) can alter these predictions. For instance, the war in Ukraine has had an impact on international air transport and has forced to modify some routes to bypass Russian airspace. For a Paris to Tokyo flight, it is possible to fly south to Turkey or Kazakhstan and return north to the pole. This itinerary has consequences on the flight time (count about 14 hours instead of 11) and on the doses received. If the outbound flight is less dosing, the return is up to twice as much. We can therefore observe a possible effect on the average annual effective dose in 2022.

Conclusion

After an increase in aircrew exposure over the 2015-2019 period, several factors suggest that cosmic radiation exposure for these workers should decrease in the coming years. The evaluation of the occupational exposure for the year 2022, scheduled for the first half of 2023, should already provide us with elements of comparison.

The health effects of these low dose exposure on FCs remains to be delineated. Except for the study by [START_REF] Wollschläger | Estimated radiation exposure of German commercial airline cabin crew in the years 1960-2003 modeled using dose registry data for 2004-2015[END_REF], no large epidemiological study on the long has term has been conducted on this specific exposed population. This would be of interest, especially with the significant although quite narrow range of exposure, as demonstrated in the present study.

Figure 1 :Figure 2 :EffectiveFigure 3 :Figure 4 :

 1234 Figure 1: Distribution of numbers (histogram bars) and average annual effective doses (± SD) (in mSv) of FCs (curve) over the period 2015 to 2021

Figure 5 :Figure 6 :

 56 Figure 5: Number of TFCs (A) and CFCs (B) who received an effective dose greater than 5 mSv by occupation over the period 2016 to 2019 (A)

In any case, another category of workers, for whom no dosimetric monitoring is performed, deserves to be studied: "the frequent flyers", i.e., the persons who, for professional duties, regularly travel by plane on commercial flights. Exposures for these individuals may exceed 1 mSv per year (Tate et al., 2021). A pilot study to objectivize these exposures and possibly move towards a modification in the regulations could be also proposed.

Acknowledgments

We would like to thank all people who helped us in any way during this study.

Funding

This work was supported by Institute of Radioprotection and Nuclear Safety

Informed consent

This article contains no experimental studies involving human subjects.

Ethical approval

Ethical approval was not required.

Conflict of interest

The authors declare they have no conflict of interest.

Authors' contribution

P Lestaevel: study design, analysis, writing original draft. C Huet: analysis, writing and reviewing original draft. V Lejeune: data collection and analysis. C Moreno: investigation, reviewing original draft. Carmen Villagrasa: writing and reviewing original draft. J Feuardent: reviewing original draft. F Trompier: writing and reviewing original draft.


Tam Khac

Minh Truong

Roland Faure

Rumen Andonov

Assembling close strains in metagenome assemblies using discrete optimization

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


Loic Le Guennec
email: loic.leguennec@aphp.fr

Alexander Balcerac

Nicolas Weiss

MD Loïc Le Guennec

A Letter Concerning A Role for Blood-brain Barrier Dysfunction in Delirium following Non-Cardiac Surgery in Older Adults #

Keywords:

To the Editor, Devinney et al. [START_REF] Devinney | A Role for Blood-brain Barrier Dysfunction in Delirium following Non-Cardiac Surgery in Older Adults[END_REF] highlight a link between blood-brain barrier (BBB) dysfunction and post-surgery delirium in older adults. Using the MADCO-PC and INTUIT studies, they analyzed cerebrospinal fluid (CSF) samples from before and after surgery. CSF-to-plasma albumin ratio (CPAR) was used as a surrogate for BBB dysfunction. In their study, 12.6% of subjects experienced postoperative delirium, with a significant rise in CPAR values post-surgery. This change correlated with delirium and longer hospital stays.

It should be specified that CPAR is an imperfect marker for BBB dysfunction for several reasons. First, CPAR evaluates the blood-CSF barrier (BCSFB) rather than the BBB 2, 3 . It's essential to distinguish between the BBB and BCSFB, as they have different diffusion pathways and CSF flow rates. The authors suggest a continuous relationship between CPAR and BBB permeability, referencing a study by Montagne et al. [START_REF] Montagne | Blood-brain barrier breakdown in the aging human hippocampus[END_REF] which shows a correlation between CPAR and a specific MRI sequence called "Ktrans". This sequence reflects the efflux rate of gadolinium from blood plasma into the brain's extravascular-extracellular space. Unfortunately, the generalization of these findings established on a population of normal aging subjects or subjects with various causes of cognitive impairment to patients undergoing surgery lacks evidence. Also, as discussed by Reiber and others [START_REF] Reiber | Proteins in cerebrospinal fluid and blood: barriers, CSF flow rate and source-related dynamics[END_REF][START_REF] Reiber | Flow rate of cerebrospinal fluid (CSF)--a concept common to normal blood-CSF barrier function and to dysfunction in neurological diseases[END_REF] , CPAR doesn't consider CSF flow rate, and factors like reduced CSF production that can alter its values. Since CPAR focuses on albumin (70kDa), lower levels of increased permeability can be missed [START_REF] Dadas | Biomarkers in traumatic brain injury (TBI): a review[END_REF][START_REF] Shlosberg | Blood-brain barrier breakdown as a therapeutic target in traumatic brain injury[END_REF] . Olsson et al. [START_REF] Olsson | Repeated lumbar punctures within 3 days may affect CSF biomarker levels[END_REF] also showed that consecutive lumbar punctures within 24h-48h raised CPAR, possibly due to inflammation from the initial puncture. Furthermore, Devinney's "2-hit" theory for delirium onset lacks data on peri-operative drugs. Indeed, some drugs, especially CYP3A or P-glycoprotein inhibitors [START_REF] Skrobik | Factors predisposing to coma and delirium: fentanyl and midazolam exposure; CYP3A5, ABCB1, and ABCG2 genetic polymorphisms; and inflammatory factors[END_REF] , are tied to delirium. Our team found that in conditions with increased BBB permeability, such drugs often appear in patient CSF [START_REF] Weiss | Cerebrospinal fluid metabolomics highlights dysregulation of energy metabolism in overt hepatic encephalopathy[END_REF] . Finally, the study omits details on the orthopedic surgeries performed, even though 27.1% of participants had them. Differentiating between spinal surgeries, which could be "neurosurgery" or "orthopedic surgery", is vital. Spinal surgeries might cause inflammation, affecting CSF-albumin levels.

In summary, Devinney's study provides really new insights in the pathophysiology of post-operative delirium. Nevertheless, it has gaps and potential biases. Given the complexities surrounding BBB and BCSFB, a singular approach may not capture the full extent of their dynamics and future studies should consider a multimodal assessment approach: First, incorporate a broader range of blood biomarkers such as Protein S-100 Beta, Neuron specific enolase, and GFAP. These markers can provide a more comprehensive understanding of BBB and BCSFB integrity and function. Second, utilize advanced MRI sequences allowing to measure "Ktrans". However, performing an MRI on a patient experiencing delirium is challenging, as their agitation can compromise the quality and safety of the imaging procedure. CT-Perfusion imaging [START_REF] Wintermark | Cerebral perfusion CT: technique and clinical applications[END_REF] is feasible for agitated patients because it is fast, and allows BBB assessment. Third, account for the potential bias of drugs administered during the perioperative period. By integrating these additional measures, future research can offer a more holistic and accurate understanding of BBB and BCSFB dynamics in the context of postoperative delirium. It's essential to interpret the current data with caution, recognizing its limitations, and the potential for various biases. Only with a comprehensive approach can we truly unravel the complexities of BBB and BCSFB in relation to postoperative complications in older adults.