Nicolas Chauchat
email: nicolas.chauchat@univ-pau.fr

Roland Becker
email: roland.becker@univ-pau.fr

Eric Schall
email: eric.schall@univ-pau.fr

Performance of DG methods based on different variables for low Mach number flows

Keywords: Discontinuous Galerkin finite elements method (DGFEM), Euler equations, Entropy variables, Primitive variables, Incompressible limit, Low Mach number

Based on its good theoretical properties, the use of entropy variables is an excellent choice for computing compressible flows at low Mach number. In this paper, we discuss the use of entropy variables in a discontinuous Galerkin discretization of the compressible Euler equations and generalize the numerical flux proposed by Barth to physical and conservative variables. Next, we compare the DG0 discretization based on the entropy variables with several other DG0 discretizations, and also with a standard finite volume method. Comparisons of DG1 discretization with the different sets of variables give hope in an all-Mach number solver.

Introduction

The discretization of the Euler equations [1] is generally based on conservative variables, since they allow the equations to be written as a hyperbolic system. However, several formulations have been proposed for various reasons.

In the context of finite element methods, Barth [START_REF] Barth | Numerical methods for gas dynamic systems on unstructured meshes[END_REF] used entropy variables [START_REF] Bova | An entropy variable formulation and applications for the two-dimensional shallow water equations[END_REF] in order to construct an entropy-stable discretization. The reason for this was that the discrete equations can be tested with these variables, whereas the nonlinearity of the change of variables leads to additional interpolation errors, unless the approximation is piecewise constant, which is the case for cell-centered finite volumes and discontinuous Galerkin (DG) at the lowest order.

It is also well known that the system written in conservative variables degenerates in the incompressible limit and the change of variables becomes singular, see Hauke & Hughes [START_REF] Hauke | A comparative study of different sets of variables for solving compressible and incompressible flows[END_REF] and Pesch & van der Vegt [START_REF] Pesch | A discontinuous Galerkin finite element discretization of the Euler equations for compressible and incompressible fluids[END_REF]. In view of this difficulty, the latter articles developed stable methods for continuous and discontinuous finite element approximations, both based on a combination of entropy variables and physical variables. We should, however, mention that the discretization of the diffusion in the compressible Navier-Stokes equations is considerably easier in physical variables than in entropy variables.

The purpose of this article is to present numerical experiments illustrating different aspects of the choice of variables in the context of the discontinuous Galerkin method. Our method is based on the entropy-stable flux [START_REF] Winters | Affordable, entropy conserving and entropy stable flux functions for the ideal mhd equations[END_REF] of in combination with the change of variables developed in [START_REF] Pesch | A discontinuous Galerkin finite element discretization of the Euler equations for compressible and incompressible fluids[END_REF]. This leads to identical discrete solutions in the lowest-order case DG0 and allows for fair comparisons for higher-order discretizations DG1. Since this numerical flux is different from standard approaches, our first concern is the comparison of the proposed method with the HLLC flux [START_REF] Kim | Robust HLLC Riemann solver with weighted average flux scheme for strong shock[END_REF] and the Roe-based [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF] finite volume solver in FLUENT [START_REF] Fluent | Ansys fluent[END_REF]. Our second topic is the comparison of the different sets of variables: conservative, entropy, and physical variables. Although the discrete solutions coincide in the lower case (i.e. DG0), the nonlinearities are different and have an influence on the practical behavior of the methods. Finally, our third topic is the behavior of the different methods at low Mach numbers [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF], [START_REF] Guillard | On the behavior of upwind schemes in the low mach number limit: II. Godunov type schemes[END_REF]. Since the theory of numerical methods for the Euler equations in the incompressible limit [START_REF] Lions | Incompressible limit for a viscous compressible fluid[END_REF] is not complete, computational experiments are important.

The article is organized as follows. We close this section with some notation. Sections 2 and 3 describe the proposed DG discretization with different sets of variables. The Jacobians required for the numerical fluxes are given in Appendix. Section 4 presents comparisons between the proposed numerical flux and standard approaches based on the Roe-FVM and HLLC-DGFEM discretizations. For a fair comparison, we use the lowest order DG approximation DG0, which is very close to cell-centered finite volumes. Then in Section 5 we compare the use of different variables in the DG method using well-known test problems for low Mach number flows. We also consider the higher-order approximations DG1, but restrict ourselves to cases without shocks, in order to focus on the role of the choice of variables. Higher-order DG discretizations for shock problems would require additional slope-limiters [START_REF] Hoteit | New two-dimensional slope limiters for discontinuous Galerkin methods on arbitrary meshes[END_REF], obscuring the comparison.

Finally, we draw some conclusions on the numerical experiments presented in this article in section 6.

Notation

We denote by v • w the Euclidian scalar product for v, w ∈ R n (which are column vectors), and by |•| the associated norm. We frequently use the notation v T w = v • w. We denote the transpose of a matrix A ∈ R n×m by A T .

For a differentiable function f : A ⊂ R n → R m we denote by Df its Jacobian matrix, (Df ) ij = ∂fi ∂uj . In the special case m = 1 we let ∇f = Df T . In particular, we have for g

(u) = f (T (u)) that ∇g = (Df DT ) T = DT T ∇f . Let φ : R n → R, f : R n → R n and φ(u) := u•f (u). Then Dφ = f (u)+u T Df , since ∂φ ∂u i = n k=1 ∂(u k f k (u)) ∂u i = n k=1 δ ik f k (u) + u k ∂f k (u) ∂u i = f i (u) + (u T Df ) i .

The Euler equations as a hyperbolic system

We consider the classical Euler equations for an ideal gas p/ρ = RT [START_REF] Schochet | The compressible Euler equations in a bounded domain: existence of solutions and the incompressible limit[END_REF].

Let ρ, V , p, T and E be the density, velocity, pressure, temperature and total energy respectively. The Euler equations expressed as a system of mass and momentum and energy conservation equations read:

∂ρ ∂t + div (ρV ) = 0, ∂ρV ∂t + div (ρV ⊗ V ) + ∇p = 0, ∂ρE ∂t + div (ρV E) + div (pV ) = 0. (1)
We write (1) as a system of hyperbolic equations in the conservative variables

u =      ρ ρV ρE      (2)
with the help of smooth flux functions

f i : A → R m , (A ⊆ R m ) 1 ≤ i ≤ d ∂u ∂t + d i=1 ∂f i (u) ∂x i = 0 (3)
combined with the initial condition

u(0) = u 0 (4)
and appropriate boundary conditions. The quasi-linear form of (3) uses the flux

Jacobians A i (u) = Df i (u) : ∂u ∂t + d i=1 A i (u) ∂u ∂x i = 0 (5)
Equation ( 3) is understood in the weak sense, since the solutions are in general not differentiable and even discontinuous (formation of shocks). From a mathematical point of view, a weak formulation is therefore considered, where all the derivatives are put on smooth test functions via integration by parts.

For the Cauchy problem (Ω = R n ), this reads for example:

Find u ∈ L 1 ([0; +∞[×R n )
such that for all smooth compactly supported test functions φ

∞ 0 R n u ∂φ ∂t + d i=1 f i (u) • ∂φ ∂x i = R n u 0 φ(0)
The mathematical theory of hyperbolic systems in several space dimensions is still incomplete. The situation is much more satisfactory for scalar equations.

In order to obtain uniqueness of the weak solution, they must satisfy an additional entropy condition.

If u is a piecewise smooth solution with a surface of discontinuity Γ, the weak formulation implies the Rankine-Hugoniot condition, which determines the shock speed s:

s[u] = [f n (u)] (6)
where [.] is the jump.

The system is hyperbolic if the Jacobian

A n = Df n (u) = d i=1 A i (u)n i admits a spectral decomposition A n (u) = R n (u)Λ n (u)R n (u) -1 (7)
with a real diagonal matrix of eigenvalues Λ(u) and a regular matrix R n (u) composed of the eigenvectors of A n for all vectors n ∈ R d .

Since hyperbolicity is conserved under a smooth change of variables u = Φ(w), the quasilinear system can be transformed into

DΦ(w) ∂w ∂t + d i=1 A i (u(w))DΦ(w) ∂w ∂x i = 0. ( 8
)
It is called symmetric if DΦ(w) is symmetric positive definite and matrices A i (u(w))DΦ(w) = g i (w), g i (w) := f i (u(w)) are symmetric.

The system admits a mathematical entropy if there is a pair of smooth functions

(U, F ) U : A → R, F i : A → R (1 ≤ i ≤ d, A ⊂ R m ) such that U is strictly convex and for all i DU (u)Df i (u) = DF i (u). ( 9
)
The existence of an entropy for the Euler equations is crucial and follows from thermodynamics. It implies hyperbolicity and the possibility to symmetrize the system.

We impose that the entropy is a decreasing function of time (the physical en-tropy, which is more or less its negative, is increasing):

∂U (u) ∂t + d i=1 ∂F i (u) ∂x i ≤ 0 (10)
Since solutions develop shocks, [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF] has again to be understood in a weak sense.

Similarly as before we get the Rankine-Hugoniot condition [START_REF] Kennel | MHD Intermediate Shock Discontinuities. Part 1. Rankine-Hugoniot Conditions[END_REF] for entropy:

s[U (u)] ≤ [F n (u)] (11)
For a smooth solution we do have an equality in [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF] thanks to [START_REF] Fluent | Ansys fluent[END_REF]:

∂U (u) ∂t + d i=1 ∂F i (u) ∂x i =DU (u) ∂u ∂t + d i=1 DF i (u) ∂u ∂x i =DU (u) T • ∂u ∂t + d i=1 DU (u)Df i (u) ∂u ∂x i =DU (u) T • ∂u ∂t + d i=1 ∂f i (u) ∂x i = 0.

Entropy variables

The Legendre-Fenchel transform [START_REF] Corrias | Fast Legendre-Fenchel transform and applications to Hamilton-Jacobi equations and conservation laws[END_REF] of U (u), U * (v), is defined by

U * (v) = sup u (u • v -U (u)) = u(v)•v -U (u(v)) with DU (u(v)) = v
T and is called entropy potential. Then we have

DU * (v) = u(v) T + v T ∂u ∂v -DU ∂u ∂v = u(v) T . (12)
Since the conjugate of a smooth convex function is a smooth convex function,

it follows that M := M (v) = ∂u ∂v = D 2 U * (v) is symmetric positive-definite and A n is similar to M -1/2 A n M 1/2 = M -1/2 (A n M ) M -1/2 . ( 13
)
which is symmetric, since A n M is symmetric as we will see below, and the hyperbolicity of the system follows.

The mapping v → u(v) is a special change of variables and we call v = DU T the entropy variables in contrast to the conservative variables u.

Clearly V (u) = U (u(v)) is convex since D 2 V = I. Let g i (v) := f i (u(v)), and G i (v) := F i (u(v)).
Then

DG i (v) = DF i (u(v))M (v) = DU (u(v))Df i (u(v))M (v) = DU (v)Dg i (v). ( 14
)
This change of variables transforms the quasilinear form (5) into

M (v) ∂v ∂t + d i=1 C i (v) ∂v ∂x i = 0 at Ω, (15)
with M (v) symmetric positive definite and

C i (v) := Dg i (v) = Df i (u(v))M (v)
symmetric, as will be seen below. Therefore, the entropy variables symmetrize the system. The equivalence between the existence of a generalized entropy and symmetrizability of the equations has been shown in Mock [START_REF] Mock | Systems of conservation laws of mixed type[END_REF] and Godunov [START_REF] Godunov | An interesting class of quasilinear systems[END_REF], [START_REF] Godunov | The problem of a generalized solution in the theory of quasilinear equations and in gas dynamics[END_REF].

The symmetry of the Jacobian in entropy variables remains to be seen. The

entropy potential F * i (v) is defined such that F i (u) + F * i (v) = f i (u) • v resp. G i (v) + F * i (v) = g i (v) • v (16)
It follows from eq. ( 14) that

DF * i (v) = g i (v) T + v T C i (v) -DG i (v) = g i (v) T (17)
From this identity follows the symmetry of

C i (v) = Dg i (v) = D 2 F * i (v).
In addition to symmetrizing the system, the use of entropy variables avoids problems of Jacobian definition when the Mach number tends to 0. Indeed, Hauke and Hugues [START_REF] Hauke | A comparative study of different sets of variables for solving compressible and incompressible flows[END_REF] showed that the Jacobian limit is not defined in conservative variables when the Mach number tends to 0. For example, in 1D:

lim M 0 A 1 (u) =      0 1 0 +∞ ∓ 0 0 ± 0 0 ±∞ -0 0 ± 0 0      , (18)
while we do not encounter this problem in entropy variables.

In the following, the entropy vector will correspond to the following vector:

v =      μ-1 2 |V | 2 T V T -1 T      , (19)
with μ = h -T s the chemical potential of the fluid, V the velocity, T the temperature, h the enthalpy and s the entropy. For a perfect gas, the entropy variables can be written as:

v =      cpT (1-ln(T ))+RT (ln(p)-ln(R))-1 2 |V | 2 T V T -1 T      (20)
with p the pressure, c P the specific heat at constant pressure and R the perfect gas constant.

Primitive variables

Next to the entropy variables, the primitive variables [START_REF] Fletcher | A primitive variable finite element formulation for inviscid, compressible flow[END_REF] also have a welldefined Jacobian when the Mach number tends to 0, but the primitive variables do not symmetrize the system. The primitive variables are defined as:

y =      p V T      . ( 21
)

Approximation with discontinuous finite elements

In this section the written formalism is the one that the mathematical community used. Symbols and operators are described in the appendix.

Let H be a family of admissible regular meshes in the usual sense. A mesh h consists of cells K h . We denote by S int h the set of interior faces and we put

S h = S int h ∪ S ∂ h . The measure of K (or S ∈ S h ) is denoted by |K| (|S| respec- tively).
For an interior face S, n S is a fixed unit vector normal to S. If the face S lies on ∂Ω, we set n S = n Ω , the outward pointing unit normal vector.

We define the spaces of piecewise polynomial functions of degree k ∈ N on elements:

D k h := v h ∈ L 2 (Ω) : v h | K ∈ P k (K) ∀K ∈ K h . Let u ∈ D k h .
We define for a given interior face S ∈ S int h and x ∈ S u in S (x) := lim

ε 0 u(x -εn S ), u ex S (x) = lim ε 0 u(x + εn S ).
Next, we define the jump and mean for x ∈ S by

[u](x) := u in S (x) -u ex S (x), {u}(x) := 1 2 u in S (x) + u ex S (x) ,
such that n S is oriented from K in towards K ex , the two triangles containing the side S.

For the sake of brevity, we use the following notation for piecewise integration:

K h = K∈K h K and S h = S∈S h S .
We will first consider discretization of the Euler equations in conservative variables based on the approximation space (D k h ) m , i.e., we are looking for a piecewise polynomial discontinuous vector of unknowns u h , the components of which are approximations to ρ, ρV and ρE. Note that in the lowest-order case k = 0 this yields approximations ρ h , V h , E h (and p h ) to the physical quantities. For k > 0 we generally do not have polynomial approximations to v h and E h .

First, we define a form a

h : (D k h ) m × (D k h ) m → R by a h (u h )(φ h ) := - K h d i=1 f i (u h ) • ∂φ h ∂x i + S int h f Σ n S (u in h , u ex h ) • [φ h ] + S ∂ h f ∂ n S (u in h , u D )φ h , (22)
with an interior flux f Σ n and a boundary flux f ∂ n . Since the sign of the jump term (and the definition of • in and • ex ) depends on the orientation of n S , we require

f Σ -n (u ex h , u in h ) = f Σ n (u in h , u ex h ), (23)
in order for a h to be independent of the choice of the normal.

In addition, the internal flux is supposed to verify the consistency condition

f Σ n S (u, u) = f n (u). ( 24
)
We now consider the semi-discrete problem

∂u h ∂t , φ h + a h (u h )(φ h ) = 0 ∀φ h ∈ D k h . ( 25
)
Complete discretization is achieved if we use an approximation of the time derivative by one of the classical methods. In general, we will use the implicit Euler method which amounts to

u n h -u n-1 h ∆t , φ h + a h (u n h )(φ h ) = 0 ∀φ h ∈ D k h , (26)
for the time steps n = 1, 2, . . . and an approximation u 0 h of the initial data. The definition of the internal fluxes is crucial for the entropy-stability of the method (E-flux, see Tadmor [START_REF] Tadmor | Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems[END_REF]).

We consider upwind fluxes of the form

f Σ n (u in , u ex ) = {f n (u)} + D Σ (u in , u ex )[u], (27)
which trivially fulfill the conditions ( 23) and [START_REF] Barth | Numerical methods for gasdynamic systems on unstructured meshes[END_REF]. The dissipation matrix D Σ can be defined in many ways; the simplest choice is the diagonal matrix scaled by the maximum absolute eigenvalues (Lax-Friedrich). Another classical example is given by the Roe flux [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF]. In the lowest-order case k = 0, the DG method is closely related to cell-centered finite volumes, and any of the fluxes developed in this context, such as approximate Riemann solvers (HLLC), can also be used.

In our implementation, we use an integral of the absolute value matrix function, see below.

For the boundary fluxes, we follow a similar approach with one of the states defined by the boundary conditions.

Change of variables

Next we consider the discretization for a change of variables. Let w be another set of variables, either the entropy variables v = DU (u) T or the primitive variables y = (p, V, T ). In order to simplify notation, the change of variables to the conservative variables is simply written as u(w).

The semi-discretization looks for

w h ∈ D k h such that ∂u(w h ) ∂t , φ h + a h (w h )(φ h ) = 0 ∀φ h ∈ D k h , (28)
with

a h (w h )(φ h ) := - K h d i=1 f i (u(w h )) • ∂φ h ∂x i + S int h f Σ n S (w in h , w ex h ) • [φ h ] + S ∂ h f ∂ n S (w in h , w D )φ h (29)
The first observation is that the contribution involving the time derivative is nonlinear:

∂u(w h ) ∂t , φ h = ∂u(w h ) ∂w ∂w h ∂t , φ h . (30)
In the following simulations, an implicit Euler method is used and the solution for each time step will be calculated using a Newton method. The linear system used in the Newton method is expressed in the dimensionless form as will be explained in Section 3.2.

The second observation concerns the internal flux function the dimensionless form of the equations is obtained using a recurrent set of four reference values; we choose the set {ρ ∞ , a ∞ , T ∞ , L}, where a ∞ is the free-stream speed of sound defined as:

f Σ n (w in , w ex ) = {f n (u(w))} + D Σ (w in , w ex )[w], (31)
a ∞ = γp ∞ ρ ∞ ( 32
)
and L is the characteristic length scale of the problem.

Table 1 gives the relation between the dimensional and dimensionless quantities (denoted with a tilde), for example Using these formulas, the dimensionless equations can be obtained by replacing the dimensional quantities by their expression in terms of the reference values and the dimensionless quantities. Equation (33) is the dimensionless equations for the 1D Euler equation with the conservative variables.

p = ρ ∞ (a ∞ ) 2 p. Quantity Symbol Magnitude Formula mass m [M ] ρ ∞ L 3 length x [L] L time t [T ] L/a ∞ absolute temperature 1 T [θ] T ∞ density ρ [M L -3 ] ρ ∞ velocity V [LT -1 ] a ∞ total internal energy E [L 2 T 2 ] (a ∞ ) 2 pressure p [M L -1 T -2 ] ρ ∞ (a ∞ ) 2 specific heats c v , c p [L 2 T -2 θ -1 ] (a ∞ ) 2 /T ∞
     ρ∞a∞ L 0 0 0 ρ∞(a∞) 2 L 0 0 0 ρ∞(a∞) 3 L           ∂ ∂ t      ρ ρ Ṽ ρ Ẽ      + ∂ ∂ x      ρ Ṽ (ρ Ṽ 2 + p) Ṽ (ρ Ẽ + p)           = 0 ( 33
)
1 Do not confuse the symbol T for temperature with the magnitude time [T ]

For inviscid aerodynamics applications, the following six parameters define the flow problem [START_REF] Kok | An industrially applicable solver for compressible, turbulent flows[END_REF]:

• the characteristic length scale L,

• the free-stream density ρ ∞ , velocity V ∞ and temperature T ∞ ,

• the gas constant R and the specific heat at constant pressure c p .

Since the equations contain four fundamental magnitudes, two dimensionless Pi groups can be formed according to Buckingham's Pi theorem. There are:

• the free-stream Mach number M ∞ = V ∞ /a ∞ ,
• the ratio of specific heats γ = c p /c v , with the set of four reference values, these are the inputs for our computations.

The values of the other parameters required in the simulation are calculated using:

V ∞ = M ∞ a ∞ (34) p ∞ = ρ ∞ (a ∞ ) 2 γ (35) R = p ∞ ρ ∞ T ∞ (36) c p = γ γ -1 R (37)

Numerical flux

In order to obtain entropy-stability, we wish to choose (for general variables) the test function φ h = v(w h ). This form is particularly interesting for the choice of entropy variables, since φ h ∈ D k h . For other variables this is only possible for k = 0. Next, choosing w = v, we immediately obtain the numerical dissipation for the entropy, and a theoretical sound choice, see Barth [START_REF] Barth | Numerical methods for gasdynamic systems on unstructured meshes[END_REF][START_REF] Barth | Simplified discontinuous Galerkin methods for systems of conservation laws with convex extension[END_REF], is

D Σ (v in , v ex ) = 1 0 R n (v s )|Λ n (v s )|R n (v s ) T ds, (38)
where v s := (1 -s)v in + sv ex . In practice, the integral is approximated by numerical integration.

In the case of primitive variables (and also conservative variables), we obtain an entropy flux using the transformation to entropy variables

D Σ (w in , w ex ) = D Σ (v in , v ex ) 1 0 ∂v(w s ) ∂w ds. (39)
The numerical integration is performed with the Simpson method which uses three integration points to calculate (38) and with a two-point integration method to calculate With entropy variables and the Simpson method for each integral, the solution obtained in fig. 1 shows the presence of this phenomenon. On the other hand, using the proposed integration method, the solution in fig. 2 does not have any Carbuncle.

Comparison of finite elements method DG0 and finite volume method (FLUENT)

In order to compare the discontinuous Galerkin finite elements method and the numerical flux presented in Section 3. Here, we study the same test case with a finite volume method and discontinuous Galerkin finite elements method.

All computations are done on a mesh with 7124 triangles. The dimensionless inflow velocity and density are set equal to one. The inflow pressure is chosen such that the Mach number is 10 -1 and 10 -3 , respectively. A comparison of the different solutions is shown for the contour of the normalized pressure p norm = (p -p min )/(p max -p min ). The pressure extrema for each simulation are shown in Table 3.

Figures 6 and10 show the residuals for each calculation. In each case we converge as much as possible in order to reduce the residual to values below 10 -12 .

For both numerical fluxes, the discontinuous Galerkin method converges for the Mach number 10 -1 and the simulations converge in 50 iterations. The finite volume method, however, converges only after 5000 iterations. For the Mach number 10 -3 , convergence is better with the entropy variables-based flux than with the HLLC flux. The finite volume method requires 10 times more iterations to reach a converged solution. High order DG implicit methods are well know to be efficient in computing iteration [START_REF] Schall | Implicit method and slope limiter in AHMR procedure for high order discontinuous Galerkin methods for compressible flows[END_REF].

The solutions presented in Figs. 5 and9 are with conservative variables and DG0 and the flux presented. In both cases, the pressure contours are smooth and correspond to the incompressible solution [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF].

Similar to the results of Guillard and Viozat, the finite volume method (see Figs. 3 and7) and the discontinuous Galerkin finite element method with the HLLC flux (see Figs. The nozzle studied here (see fig. 11) is 1m long and 0.1m in diameter at the inlet and outlet. The diameter at the neck is 0.09m. The mesh used for the computation is a structured mesh composed of 2352 triangular cells, symmetrical with respect to the horizontal axis. The residuals for these simulations are plotted in figs. 14 and 17. Like for the Naca test case, we try to reach a normalized residual lower than 10 -12 . With the finite elements methods, convergence is achieved in 20 iterations and the finite volume method needs 100 times more iterations. However, the cost of an iteration is not the same with the finite volumes and finite elements methods.

A more complete study should be conducted with the computational time to know the quickest method.

Comparison of different sets of variables at low Mach number

In order to evaluate the contribution of the use of primitive and entropy variables compared with conservative variables, the new formulation is benchmarked on previous test cases with higher-order DG discretization (i.e. DG1).

The test cases cover a range of low Mach number, from 10 -4 to 10 -1 . All these test cases are non-viscous flows, without thermal diffusion and stationary. There is a description in table 4.

In order to compare the methods and see the differences between them, the simulations are performed on a coarse mesh and without mesh convergence. A relative L 2 norm of the errors with respect to the analytical solution is presented in tables 5 and 6. Primitive 1.0357276e -3 (-0.0001%) 1.85658e -2 (0.002%) Entropy 1.0357214e -3 (-0.0007%) 1.85669e -2 (0.007%)

Mach number Pressure

(% with respect to conservative variables) Conservative 1.0361418e -3 9.90094082908e -1

Primitive 1.0361418e -3 (0%) 9.90094082908e -1 (0%) Entropy 1.0361436e -3 (0.0002%) 9.90094082908e -1 (0%)

Conclusion

In this article we have generalized the numerical flux for the discontinuous Whatever the set of variables with high order DG1 discretization, results are the same at low Mach number for the compressible Euler equations. This leads us to believe that we are on track to reach an all-Mach number solver.

The Jacobian with conservative variables.

The Jacobian of the change of variables:

A 0 (u) = ∂u ∂u =         1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1         (40)
The Jacobian of the Euler flux:

A 1 (u) = ∂f 1 (u) ∂u =         0 1 0 0 1 2 R c V |V | 2 -V 2 1 2 -R c V V 1 -R c V V 2 R c V -V 1 V 2 V 2 V 1 0 V 1 1 2 R c V |V | 2 -H H -R c V V 2 1 -R c V V 1 V 2 cp c V V 1         (41)
A 2 (u) = ∂f 2 (u) ∂u =         0 0 1 0 -V 1 V 2 V 2 V 1 0 1 2 R c V |V | 2 -V 2 2 -R c V V 1 2 -R c V V 2 R c V V 2 1 2 R c V |V | 2 -H -R c V V 1 V 2 H -R c V V 2 2 cp c V V 2         (42)
The Jacobian with primitive variables.

The Jacobian of the change of variables:

A 0 (y) = ∂u(y) ∂y = 1 RT         1 0 0 -p T V 1 p 0 -pV1 T V 2 0 p -pV2 T E pV 1 pV 2 -p|V | 2 2T         (43)
The Jacobian of the Euler flux:

A 1 (y) = ∂f 1 (y) ∂y = 1 RT         V 1 p 0 -pV1 T V 2 1 + RT 2pV 1 0 -pV 2 1 T V 1 V 2 pV 2 pV 1 -pV1V2 T V 1 H p(H + V 2 1 ) pV 1 V 2 -pV1|V | 2 2T         (44)
A 2 (y) = ∂f 2 (y)

∂y = 1 RT         V 2 0 p -pV2 T V 1 V 2 pV 2 pV 1 -pV1V2 T V 2 2 + RT 0 2pV 2 -pV 2 2 T V 2 H pV 1 V 2 p(H + V 2 2 ) -pV2|V | 2 2T         (45)
The Jacobian with entropy variables.

Calculating the entropy Jacobian is easier after calculating that of the change in variables from primitive ones to entropy.

We start with its inverse:

∂v(y) ∂y =         R p -V1 T -V2 T -H+|V | 2 T 2 0 1 T 0 -V1 T 2 0 0 1 T -V2 T 2 0 0 0 1 T 2         (46) ⇒ ∂y(v) ∂v =         p R pV1 R pV2 R pH R 0 T 0 V 1 T 0 0 T V 2 T 0 0 0 T 2         (47)
The Jacobian of the change of variables:

A 0 (v) = ∂u(v) ∂v = p R 2 T         1 V 1 V 2 E V 2 1 + RT V 1 V 2 V 1 H V 2 2 + RT V 2 H sym HE + 1 2 |V | 2 RT         (48)
The Jacobian of the Euler flux:

A 1 (v) = ∂f 1 (v) ∂v = p R 2 T         V 1 V 2 1 + RT V 1 V 2 V 1 H V 1 (V 2 1 + 3RT ) V 2 (V 2 1 + RT ) V 2 1 (RT + H) + RT H V 1 (V 2 2 + RT ) V 1 V 2 (RT + H) sym V 1 H(RT + H)         (49) A 2 (v) = ∂f 2 (v) ∂v = p R 2 T         V 2 V 1 V 2 V 2 2 + RT V 2 H V 2 (V 2 1 + RT ) V 1 (V 2 2 + RT ) V 1 V 2 (RT + H) V 2 (V 2 2 + 3RT ) V 2 2 (RT + H) + RT H sym V 2 H(RT + H)         (50)
Description of symbols and operators

3 3. 2 .

 32 Different choices are possible for the dissipation matrix D Σ . The choices made in this work are presented in Section 3.Non-dimensionalization of the Euler equation The magnitudes of the physical quantities involved in the Euler equations can be expressed in terms of the four fundamental magnitudes: mass [M ], length [L], time [T ] and temperature [θ]. According to Buckingham's Pi theorem [22],

  If we compute the flow around a circular cylinder at the Mach number 3 on a quadrilateral mesh with 1680 elements and k = 0 and do not use sufficiently accurate quadrature rules, as indicated in the previous paragraph then an unphysical Carbuncle solution is obtained.

Figure 1 :

 1 Figure 1: Flow around a circular cylinder at Mach number 3 with DG0 and entropy variables and Simpson method for each integral of the numerical flux.

Figure 2 :

 2 Figure 2: Flow around a circular cylinder at Mach number 3 with DG0 and entropy variables and the proposed integration method for the numerical flux.

4. 1 .

 1 Naca0012 airflowIn[START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF], Guillard and Viozat studied the flow around a NACA0012 airflow at several low Mach numbers. They showed by a first-order accurate finite volume method that with a Roe approximate Riemann flux, the pressure field shows fluctuation at low Mach numbers. They proposed a preconditioned dissipation matrix to obtain acceptable results.

4 and 8 )

 8 , the pressure contours are not accurate for very low Mach numbers and the solution at Mach number 10 -1 is the closest to the incompressible solution.

Figure 3 :Figure 4 :

 34 Figure 3: Pressure field NACA0012 airfoil at Mach 10 -1 , Fluent, p min = 71.2294P a, pmax = 71.7836P a

Figure 5 :Figure 6 :

 56 Figure 5: Pressure field NACA0012 airfoil at Mach 10 -1 , DG0 CONS, p min = 71.4286P a, pmax = 71.5148P a

Figure 7 : 20 Figure 8 :Figure 9 : 21 Figure 10 : 3 22 4 . 2 .

 720892110342 Figure 7: Pressure field NACA0012 airfoil at Mach 10 -3 , Fluent, p min = 714285P a, pmax = 714286P a

Figure 11 :

 11 Figure 11: Computational domain for ILSN nozzle.

Figure 12 :Figure 13 :Figure 14 :Figure 15 : 0036 Figure 16 :Figure 17 :

 1213141500361617 Figure 12: Evolution of the Mach on the central axis of the nozzle. Inlet Mach number : 0.036

5. 1 .

 1 Internal low speed nozzle (ILSN) at Mach 0.036 and 0.0036The internal low-speed nozzle studied here is the one in section 4.2. The results obtained for both infinite upstream Mach number 0.036 and 0.0036 are all in good agreement with the analytical solution and superimposed. In each case, we try to reduce the normalized residual to values below 10 -12 . They are plotted in figs.18 and 19.

Figure 18 :Figure 19 :

 1819 Figure 18: Normalized residuals for computation of flow in ISLN nozzle at Mach 0.036

5. 2 .

 2 Naca0012 airflow Like in section 4.1, we studied a flow around a NACA0012 airflow at several low Mach numbers. The dimensionless inflow velocity and density are set equal to one. The inflow pressure is chosen such that the Mach number is between 10 -1 and 10 -4 . Figures 20 to 23 show the residual of the simulations. In every case we converge as much as possible in order to reduce the normalized residual to values below 10 -12 . Here, the comparison of the solution is made on the entropy. For this incompressible test case, the entropy must be constant in the whole field. The relative L 2 norms of the errors with respect to the infinite upstream value are shown in table 7. Whatever the set of variables and whatever the Mach number, results are the same. However, the lower the Mach number the smaller the difference.

Figure 20 : 1 30Figure 21 :

 20121 Figure 20: Normalized residuals for computation of flow around NACA0012 airfoil at Mach 10 -1

Figure 22 : 3 31Figure 23 :

 22323 Figure 22: Normalized residuals for computation of flow around NACA0012 airfoil at Mach 10 -3

  Galerkin method proposed byBarth to physical and conservative variables in the case of the low Mach number limit. A new numerical flux consistent with all the possible set of variables used in CFD has been performed. Comparisons of different variables are carried out with entropy variables in particular and results are satisfactory. The presented High DG method combined with a dedicated flux shows accurate results whatever the set of variables. Comparison with the standard numerical HLLC flux in conservative variables shows satisfactory behavior. Furthermore, we have shown that our approach outperforms the standard finite volume methods in FLUENT at very low Mach numbers.

  the mean v • w = v T w the Euclidian scalar product for u, v ∈ R n (which are column vectors) |.| the norm of a column vector A T the transpose of a matrix A ∈ R n×m Df the Jacobian matrix of a differentiable function f H a family of admissible regular meshes in the usual sense K h cells of a mesh h S int h the set of interior faces of a mesh h S ∂ h the set of boundary faces of a mesh h |K| the measure of a cell K |S| the measure of a face S n S a fixed unit vector normal to a face S f the change of variable from w to u

Table 1 :

 1 Physical quantities in the Euler equations, their symbol, magnitude and formula in terms of the recurrent set {ρ∞, a∞, T∞, L}

Table 2 :

 2 Summary of test cases.

	3 with the existing method, two sta-

(i.e. DG0), which is very close to cell-centered finite volumes. These test cases consider 2D inviscid flows without thermal diffusion. Because of the low Mach number considered in these test cases, which ranges from 10 -3 to 10 -1 , the flow can be considered nearly incompressible. It is worth noting that the set of equations used to make the computations is that of the compressible model. This demonstrates the robustness and the accuracy of the presented method at low mach number. A description of the test cases can be found in table 2.

These test cases allow a comparison between the finite volume method with a

Roe flux (without any low Mach number preconditioner) and two discontinuous Galerkin finite element formulations, one with the HLLC flux and the other with the new flux presented in Section 3.3. The finite volume method with a Roe flux used here is implemented in the ANSYS-Fluent code.

Table 3 :

 3 Pressure extrema for each simulation.

Table 4 :

 4 Summary of test cases.

	Test case	Mach number	Compressibility State	Analytical solution
	Naca0012 airflow	10 -4 to 10 -1	Incompressible	Steady No
	Internal low spedd nozzle 0.0036 and 0.036 Incompressible	Steady Yes

Table 5 :

 5 Relative L 2 norms of the errors with respect to the analytical solution for ILSN.

	Simulation of a flow in a nozzle at Mach number 0.036 with DG1 finite elements method and
	different sets of variables.

Table 6 :

 6 Relative L 2 norms of the errors with respect to the analytical solution for ILSN.

	Simulation of a flow in a nozzle at Mach number 0.036 with DG1 finite elements method and
	different sets of variables.

Table 7 :

 7 Relative L 2 norms of the errors with respect to the infinite upstream value for

	Naca0012 airflow.

Acknowledgments

The authors would like to thank professor Jaap van der Vegt from the Department of Applied Mathematics of the University of Twente for his collaboration on this work.

Appendix

Here we present the Jacobians of the Euler equations expressed in the different sets of variables. The expressions are given for the 2D case with ideal gas law.

The following notations are used for physical quantity : ρ, V = (V 1 , V 2 ) T , T and p for respectively : the density, the velocity, the temperature and the pressure.

The total energy is noted E = e + 1 2 |V | 2 with e the internal energy and enthalpy H. c p , c V and R are respectively specific heat at constant pressure, specific heat at constant volume and the ideal gas constant.


Bérangère Godde

Abdelaziz Jouaiti

Matteo Mauro

Roberto Marquardt

Alain Chaumont

Vincent Robert

The Motion of an Azobenzene Light-Controlled Switch: A Joint Theoretical and Experimental Approach

Keywords: density functional calculations, diazo compounds molecular dynamics, molecular machines, photoswitches

To gain further insight into the internal motion of molecular objects, we have synthesized a molecular turnstile AzoT composed of a rotor based on flexible tetraethyleneglycol (TEG) chains grafted on aromatic moieties and a stator containing a photoswitchable azobenzene (Azo) fragment. The control of the reversible light-induced E-AzoT.Z-AzoT isomerization is supported by both NMR spectroscopy and photophysical investigation, which show that the system exhibits a fatigueless isomerization switching process. Furthermore, 2D NMR spectroscopy points to the fact that the free internal motion is triggered by the E-AzoT.Z-AzoT isomerization. Using molecular dynamics simulations and DFT calculations we have investigated the nature of the internal motions. An internal rotation characterized by an energy barrier of 23 kJ/mol is found for the Z-AzoT isomer. In contrast, this barrier reaches 151 kJ/mol for the E-AzoT isomer, excluding any "classical" rotation at room temperature. This rotational movement could in principle occur via tunneling. A simple model calculation, however, excludes tunnelling to take place before 20 ms.

The ability to make use of energy sources, either electric, chemical or radiative, to generate motion is at the fundamentalessence of biological machines and life. [1][2][3] Even though molecular and biochemical systems that satisfy the conditions to establish controlled and directional movements have been synthesized, [4][5][6] the preparation and characterization of systems that collectively operate as molecular machines remains a challenging quest. One prerequisite is to build up architectures that exhibit multi-stability characterized by distinguishable molecular structures -local minima in the potential energy surface (PES) of the electronic ground state, referred to as stations. Stations are ideally populated at will by means of external inputs. Thus, both the relative energies, i.e. thethermodynamics, and the energy barriers that separate the stations, i.e. the kinetics, must be taken into account for the design of molecular machines. In addition, one would set up the conditions to drag the system away from its thermodynamic equilibrium. Even more refined attempts would monitor the molecular dynamics and the time-dependent molecular structure via the coherent control of amplitudes and phases of the quantum mechanical state describing this structure. [7][8][9] One example of this is the laser induced inhibition [10][11] or enhancement [11][12] of tunnelling.

Following the pioneering works of various groups, [13][14][15][16][START_REF] Balzani | Proc. Natl. Acad. Sci[END_REF][START_REF] Feringa | [END_REF] light-activated molecular machines are the object of extensive investigation. [START_REF]Molecular devices and machines: concepts and perspectives for the nanoworld[END_REF][START_REF] Kassem | [END_REF][21] The use of light input introduces several advantages as it can be remotely applied with high spatiotemporal and energetic resolutions and a variable intensity of the perturbation. Photochemical reactions provide a straightforward approach for displacing a chemical system out of its thermodynamic equilibrium.

Prompted by the recent demonstration of unidirectional transit arising from the E.Z isomerization, [22] we wish to status on the accepted picture of rotating pieces as well as to explore whether light can be used to control the kinetics of such intramolecular motion. Tamaoki and co-workers produced a wealth of systems based on azobenzenophanes where the photosensitive Azo fragment is introduced in the more flexible macrocycle. [23][24] A light-controlled molecular brake was reported, supporting a possible suppression of internal rotation. Still, we felt that theoretical inspection might be desirable to support the mechanistic views and to rationalize the actual molecular rotation. In the meantime, Hosseini and co-workers reported on molecular turnstiles, constructed on a rotor (typically a tetraethylene glycol (TEG) chain grafted onto a porphyrin, [25][26][27] an organometallic [28] or an aromatic [29][30] moiety) and a stator. With the possibility to introduce a photo-sensitive group almost at will on the stator, we designed a system where the macrocycle structure would be maintained upon irradiation. Not only does it offer another system to test the concept of motion slowdown but the inherent flexibility of the turnstile would be maintained identical in both isomers. Thus, the microscopic inspection would not suffer from deep structural changes whereas the discriminating factor should be attributed to the stator modification. We believe that such an engineering that maintains the structure of the macrocyclic loop might produce a molecular slowdown that complements previous molecular brakes. How much a simplistic steric view can be adapted to properly describe modifications in the kinetics of the geared internal rotation is at the heart of the present study.

We herein report on a synthetic strategy combining a macrocycle with a stator in which an Azo moiety can be converted from its E (E-AzoT) to Z (Z-AzoT) form by light as a wireless trigger between stations that possess different steric demand. As a proof of principle, the length of the TEG chains has been selected in such a way that the resulting system can be reversibly switched between an ideally blocked (hindered "classical" rotation) and an unlocked (allowed rotation) station. The potential of the E-and Z-AzoT isomers to yield the targeted light-controlled kinetics is jointly supported by experimental and theoretical investigations. The key element is the identification of the structural deformation and the heights of the energy barriers that must be overcome to travel along the PES. Our inspections not only offer a microscopic description of the pathway for the process, but also afford structural and energetic evidence, questioning the traditional views on free rotation. The design of future molecular machineries based on the rotation of two portions of the molecule relative to each other should benefit from the combined experimental and theoretical investigation.

The detailed procedure employed for the synthesis of AzoT is described in Scheme S1 of the Supporting Information (SI). Two multistep strategies were used, affording global yields typically larger than 38%, which differ in the order of functionalization of the stator.

The design of the AzoT macrocycle (Scheme 1) stems from the molecular turnstiles previously reported, [25][26][27][28][29][30] to which a pyridyl moiety was introduced to better follow its dynamical behavior by 1 H NMR spectroscopy. On the other side the light responsive Azo moiety can be isomerized between E (thermodynamically stable) and Z (metastable) using different excitation wavelengths, [22,[31][32] with concomitant large geometrical change of the compound.

The dynamical behavior of the AzoT macrocycle was then investigated by both 1D and 2D 1 H NMR spectroscopies and the experimental findings are supported by theoretical investigations (vide infra). The AzoT appears to be a mixture of E-and Z-isomers, the former being the major component from the 1 H NMR spectrum. This finding is most likely due to photo-isomerization occurring under ambient light. The 1 H NMR spectrum is displayed in Figure S1, where attributions of the 1 H resonances were made on the basis of similar Azo derivatives. [29][30] Expectantly, full conversion of the E/Z-mixture into the thermodynamically stable E isomer can be attained by thermal treatment of a CD3CN solution (Figure S2), a procedure that was systematically followed for all samples. The 2D 1 H-1 H NOESY NMR spectrum of E-AzoT in CD3CN at 298 K is displayed in Figure 1 (for atom labeling see Scheme S2). Correlations between the rotor and stator protons are restricted to Hs, Ht (pyridine moiety), Hk (Azo) and the protons of the TEG chains Hu (Hu'), Hv (Hv'). In combination with the absence of correlations between Hl, Hq, Hp, Ho and the protons of the TEG chains, one may first conclude that the Azo moiety lies outside the loop. Besides, the phenyl ring of the Azo is distant from to the macrocycle chains. Hence the E-AzoT isomer may not be adapted to accomplish the desired internal rotation, a feature which shall be confirmed by molecular simulations. To investigate the optical and photo-switching properties, electronic absorption spectra were recorded at a concentration of 1.0×10 -5 M in CH3CN and 298 K (Figure 2a). The UV-visible absorption spectra of E-AzoT feature three absorption bands. The narrow and intense (ɛ=2.0×10 4 M -1 cm -1 ) band centred at λ abs =283 nm can be ascribed to a 𝜋 -𝜋* transition. The λ abs = 352 nm band (ɛ=2.37×10 4 M -1 cm -1 ) can be ascribed to a 𝜋 -𝜋* transition on the Azo moiety. The much less intense (ɛ ≈200 M - 1 cm -1 ) absorption band at λ abs ca. 450 nm corresponds to the formally symmetry-forbidden n ! p* transition on the Azo fragment. [22,33] Continuous irradiation of E-AzoT at λ exc =365 nm induces a time-dependent change of the electron absorption features with a decrease of the 𝜋 -𝜋* band (λ abs =352 nm) and an increase of the lowest-energy lying n -𝜋* band (λ abs = 450 nm). This indicates the light-promoted E-AzoT:Z-AzoT isomerization until the photostationary state (PSS 365 ) is finally reached (Figure 2a). To determine the PSS 365 composition, a concentrated sample of E-AzoT (1.0×10 3 M) in CD 3 CN was irradiated at 365 nm to reach PSS 365 . Signal integration of the 1 H NMR spectrum points to a 36:64 E-AzoT:Z-AzoT ratio (Figure S3). Subsequent irradiation of the initial sample at λ exc=455 nm yields back the UV-vis spectrum, with a partial recovery of the intensity of the 𝜋 -𝜋* and a decrease of the n -𝜋* band (Figure 2b). The resulting PSS 455 is characterized by a 73:27 ratio (Figure S4-S5). As previously reported, [22] the switching between PSS 365 and PSS 455 can then be repeated with almost no fatigue (Figure 2c), a prerequisite to fabricate a light-controlled rotation system. At this stage, one may question whether the timeframe needed to carry out NMR experiments is suitable for investigating a sample rich in Z-AzoT species as well, since Azo photoswitches are known to undergo thermal Z!E isomerization. Thus, this process was studied by means of spectrophotometry in CH3CN at 298 K. After irradiation to reach PSS 455 , a sample of AzoT in CH3 CN was left in the dark and thermal back isomerization investigated by means of UV-vis spectrophotometry (Figure S6a). Treatment of the time-dependent absorbance at λ abs=340 nm (Equation S1 in SI) allowed us to estimate the back-isomerization kinetic rate constant to be 2.0×10 -6 s -1 (Figure S6b). As a conclusion, the long-lived (half-life=96 hours) Z isomer is compatible with the NMR recording time scale. The resulting E-AzoT:Z-AzoT mixture in PSS 365 was then studied by 1 H-1 H NOESY NMR spectroscopy. Given the sizeable chemical shift differences in the aromatic region between the E and Z isomers, two different sets of correlations can be easily evidenced in the 2D NMR spectrum (Figure 3). First, correlations can be attributed to the presence of the E isomer, by comparison with Figure 1. Then, the predominance of Z-AzoT over the thermodynamically most stable E isomer is apparent on the spectrum. As compared to E-AzoT, a new set of cross correlation peaks between HO , HP and protons of the TEG chains appear for Z-AzoT. This finding indicates that the phenyl ring of the Azo lies in close proximity to the macrocycle chains. Such an observation supports the hypothesis that the Z-AzoT isomer is able to perform a rotational movement through the macrocycle. To further strengthen this picture, extensive computational investigations were carried out.

The large flexibility of the TEG chains makes them particularly interesting from a synthetic point of view, but any theoretical examination calls for particular care. Thus, such architectures were examined using both classical molecular dynamics (cMD) and density functional theory (DFT) calculations to analyse the mechanical picture (see SI for technical details). First, 20 μs cMD simulations of both AzoT isomers were performed at 298 K in the presence of 1000 CH 3 CN molecules. Using Boltzmann statistics, we derived a two-dimensional Gibbs energy surface (Figure 4), as a function of the distance R (defined as the distance between the centre of N=N bond of the azo group and the centre of mass of the rotor) and the dihedral angle θ (see Figure 5 for the geometrical parameters definition). The latter can be seen as the angle between the plane of the rotor and the plane of the stator. The θ=0° and θ=180° values correspond to the stator crossing the TEG macrocycle via its Azo and pyridyl groups, respectively. Three structures emerge from these calculations, named here OUT-(θ~180°), MIN-(θ~130°) and MIN IN -stations (θ ~0°) of the internal rotation and are represented along "idealized" internal rotation path is depicted schematically in Figure 6. The light green region in Figure 4 defines a flat zone of the Gibbs energy surface which contains the equilibrium structure (MIN IN -station). Indeed, the E-and Z-AzoT molecules are characterized by dihedral angle θ values of ca. 100° and R~10 Å. For both Z-and E-isomers an identical energy barrier kJ/mol) is found for the crossing of the pyridyl group of the stator through the rotor. Hence this process is not discriminat-ing between the Eand Z-isomers and the rotational barrier is easily overcome at ambient temperature. This is to be contrasted with a crossing where the Azo group enters the centre of the rotor via the MIN IN-station that is characterized by θ~0° and R~2 Å. This region is sampled during the dynamics of the Z-AzoT isomer, while it is not for the E-AzoT isomer. Figure 7 shows the crossing of the stator through the centre of the rotor at three different times. From these pathways (Figure 4), the free energy barrier for Z-AzoT rotation is found to be 23.0 kJ/mol. In sharp contrast, the coordinates of the transition states of the E-AzoT isomer lie outside the regions that are sampled during the dynamics for this isomer (red dots in Figure 4), pointing out that the internal rotation is virtually forbidden in this case. DFT calculations were then performed using as starting points the geometries of the MIN IN -, OUT-and MIN-stations, which were extracted from cMD. For each isomer, full geometry DFT optimizations lead to energetically competing structures within 10-13 kJ/mol (see Table S1 in the SI). Structurally, the rotation involves dihedral angle modifications along the rotor and stator link (α1 and α2 angles, Figure 5). Starting from a ring shape in the MIN IN -station, a conrotatory motion generates a much more compact OUT form (Figure 5). Ultimately, a significant twisting is observed, featured by a helicoidally organization of the TEG chains in the MIN station. The DFT calculations yield that for the E-AzoT isomer these three structures lie lower in energy than the corresponding ones for the Z-AzoT isomer. Besides, the energy difference between the OUT-and MIN-stations is about 0.07 eV (see Table S1). Thus, the kinetics of the internal rotation must be governed by the energy of the MIN IN -station or that of a nearby lying transition structure. As displayed in Figure 6, the MIN IN -station lies 11.4 kJ/mol above the MIN-station in Z-AzoT, suggesting that thermodynamics directly controls the Z-AzoT rotation. In contrast, the MIN IN -and MIN-stations in E-AzoT are quasi-degenerate (see Table S1). By forcing the entry of the benzyl group in the TEG loop, a barrier of 151 kJ/mol was found from this constrained geometry. Despite its approximate character, this value supports the hypothesis of a kinetically locked internal rotation involving the two quasidegenerate structures in the E-AzoT isomer. Beyond the agreement between molecular dynamics simulations and DFT calculations, one may question the possibility for the E-isomer to overcome the barrier by tunnelling. Thus, a quantum dynamical model calculation was elaborated (see SI). From a two-level model calculation, we can exclude that tunneling takes place before 20 ms which is much longer than the typical duration of the classical trajectories discussed in this paper. In conclusion, the traditional picture of an internal rotation between two stable molecular structures characterized by local minima is revisited from NMR, UV-vis spectroscopy, molecular dynamics simulations and DFT calculations. Z/E photo-switching is evidenced using UV-vis and 2D NMR spectroscopies. A counterintuitive route that connects two local minima (MIN IN and OUT) is suggested with the existence of a third well (MIN) in the potential energy surface. From numerical inspections, the overall kinetics is governed by the passage of the Azo group inside the PEG chains loop, supporting a rotation picture. Classical molecular dynamics simulations suggest an internal rotation (μs timescale) of Z-AzoT governed by an estimated barrier of 23 kJ/mol. They further suggest a blockade for E-AzoT, with a DFT estimate of the barrier for the internal rotation of 151 kJ/mol.

Scheme 1 .

 1 Scheme 1. Chemical structures of the AzoT macrocycle in its E-and Z-forms. Red and black portions of the molecule represent the stator and the rotor, respectively.

Figure 1 .

 1 Figure 1. Aromatic region of the 1 H-1 H NMR NOESY spectrum of E-AzoT recorded in CD3CN at 298 K.

Figure 2 .

 2 Figure 2. a) UV/Vis absorption spectra of E-AzoT at 1.0×10 -5 M in CH3CN before (black trace) and after (violet trace) 30 s UV light irradiation (λ exc=365 nm) at 298 K defining the photostationary state PSS365 . Grey traces indicate intermediate irradiation times. Black arrows indicate the variation of absorbance over irradiation time; b) UV/Vis absorption spectra of AzoT at 1.0×10 5 M in CH3CN at the PSS 365 before (violet trace) and after (blue trace) visible light irradiation (λ exc=455 nm) until PSS 455 is reached. For the sake of comparison, the E-AzoT species is also displayed (black trace). c) Variation of the absorption at λ abs=340 nm upon different cycles of UV (λ exc=365 nm) and visible (λ exc=455 nm) light irradiation showing the reversibility of the photoswitching process.

Figure 3 .

 3 Figure 3. 2D 1 H-1 H NOESY NMR spectrum of AzoT at the PSS 365 recorded in CD3CN at 298 K.

Figure 4 .

 4 Figure 4. Gibbs energy (kJ/mol) surface for E-AzoT (left) and Z-AzoT (right) as a function of the dihedral angle θ (°; see Figure 5 for definition) and the distance R (defined as the distance between the center of the N=N bond of the azo group and the center of mass of the stator (Å)). Pathways for the three different crossings via the MIN IN station in the case of the Z-AzoT are reported in black (crossing after 1.5 μs), blue (crossing after 13.7 μs) and magenta (crossing after 16.2 μs). The shown pathways start 50 ps before and finish 50 ps after the crossing. Snapshots along these trajectories are given in Figure 7.

Figure 5 .

 5 Figure 5. Definitions of the geometrical parameters α1, α2 and θ angles discussed in the text.

Figure 6 .

 6 Figure 6. Schematic profile of the Gibbs energy along a hypothetical path for the internal rotation for compound Z-AzoT together with the corresponding stations and estimated energy barriers. The profile for the E-AzoT is given in the Supporting Information.

Figure 7 .

 7 Figure 7. Z-AzoT isomer: Snapshots along the dynamics showing the crossing of the stator via the azo moiety through the centre of the rotor after 1.5 μs (top), 13.7 μs (middle) and 16.2 μs (bottom).

Acknowledgments

We thank the University of Strasbourg and the Ministère de l'Enseignement Supérieur et de la Recherche for financial support and a scholarship to B.G.


Alessandro Cotronei

On Wider Output Sieves

We describe and compare several novel sieve-like methods. They assign values of several functions (i.e., the prime omega functions ω and Ω and the divisor function d) to each natural number in the considered range of integer numbers. We prove that in some cases the algorithms presented have a relatively small computational complexity. A more detailed output is indeed obtained with respect to the original Eratostene's Sieve.

Introduction

The use of sieve theory is one of the most used way to find primes since the third century B.C.E. Despite faster methods are available [START_REF] Atkins | Prime Sieves using Binary Quadratic Forms[END_REF], this old method is still in wide use today because of the simplicity of its implementation. Howevere there is a relevant downside in particular from its use: the big number of repeated operations. We will start by showing an algorithm that solves this problem and introduce several novel ways of sieving natural numbers with a wider output (hence the title). We will show also several other algorithms that work in a similar way. Sections of this paper are organized with respect to the output for each algorithm.

We refer to the original sieve in the following form, that uses O(n log log n) basic operations (i.e., additions) and O(n) bits of storage [START_REF] Sorenson | An introduction to prime number sieves[END_REF]:

Algorithm 1 Eratostene's Sieve 1: Input: n.

2: Consider an array of natural numbers 1, 2, ..n -1, n. Mark the lowest unmarked number of the list p (2 in the first case) as prime.

6:

for i multiple of p not exceeding n do 7:

Mark i Composite.

Alternatively cancel out i.

8: Output: Characteristic Prime function or equivalently, all prime numbers ≤ n.

We discuss several characteristics for n = p a1 1 p a2 2 ...p a k k with p i distinct primes [START_REF] Hardy | An Introduction to the Theory of Numbers[END_REF]:

• the little omega prime function ω(n) = k. ω(1) = 0.

• the big omega prime function Ω(n) = a 1 + a 2 + ... + a k . Ω(1) = 0.

• the divisor function d(n) = (a 1 + 1)(a 2 + 1)...(a k + 1). d(1) = 0.

We make also use of several other functions as the prime counting functions π(x) and the generalized prime counting function π * k (x), defined as the number of prime numbers ≤ x (resp. positive integers ≤ x that can be written as product of not necessarily distinct of k prime numbers) [START_REF] Nathanson | Elementary Methods in Number Theory[END_REF].

We use some theorems as the Prime Number Theorem in the form:

π(x) = O x log x
, and a generalization

in the form π * k (n) = O x(log log x) k-1 (k -1)! log x [7].
We also make use of Abel's Identity [START_REF] Nathanson | Elementary Methods in Number Theory[END_REF]:

y≤x a(n)f (n) = A(x)f (x) -A(y)f (y) - x y A(t)f (t)dt,
where A(x) = 0 if x < 1, f has continuous derivative over the considered interval and

A(x) = n≤x a(n).
We give several estimates in terms of the big O notation [START_REF] Nathanson | Elementary Methods in Number Theory[END_REF] for computational complexity and bit use (memory consumption) for the algorithms.

We split this paper in three sections, each of them describes algorithms with a different function (resp. ω, Ω, d) as output.

3 Method with output ω

First Method

The first method uses a slight modification in the original sieve to evaluate the funcion ω(n) for each number in an interval of consecutive natural numbers that includes 1. It has the same computational cost with the algorithm from which it derives (if we consider the addition step as expensive as the canceling out step of the original algorithm). At the same time the memory required for the implementation is higher because the algorithm stores integer numbers (corresponding to the ω function) instead of logical values (sometimes zeroes and ones). This will be true for each algorithm we present. We want to remark that this algorithm does not discard the repeated operations of canceling out as in the original Eratostene's Sieve:

Algorithm 2 Modified Eratostene's Sieve 1 (With output ω)

1: Input: n.

2: Consider an array with indices in the natural numbers 2, 3, ..., n -1, n and an array of natural numbers and size n, we call this A.

3: a ← 0. Mark the lowest number of the list p such that a(p) ∈ {0, 1} (2 in the first case) as prime.

6:

for i multiple of p not exceeding n do 7:

a(i) ← a(i) + 1. If we assume that that memory allocation is without unused bits for each stored value, then the number of bits used is:

n k=2 log 2 (ω(k)) .
It is not our interest to give a more precise estimate for this quantity since an actual implementation of the corresponding method would not be practical.

4 Methods with output Ω

Second Method

In this procedure we determine for each number the corresponding Ω function by multiplying numbers with known Ω for prime numbers (so that the Omega is increased by one by this procedure). We mainly work with intervals of natural numbers between consecutive powers of two. We call these sets I n :

Definition 1 (I n ). We define the set

I n as {2 n , 2 n + 1, ..., 2 n+1 -2, 2 n+1 -1} or equivalently {x ∈ N : 2 n ≤ x < 2 n+1 -1}. Remark 2. It is obvious that a ∈ I n ⇐⇒ log 2 a = n.
We have

I 0 = {1}, I 1 = {2, 3}.
The following proposition gives a hint on how multiplied numbers between intervals behave:

Proposition 1. If m ∈ I q and n ∈ I r , then mn ∈ I q+r ∪ I q+r+1 .
Proof. From the Hypotheses, we have 2 q ≤ m < 2 q+1 and 2 r ≤ n < 2 r+1 . If follows that 2 q+r ≤ mn < 2 q+r+2

and the latter is equivalent to the thesis.

We will use the following Lemma to obtain numbers with m prime factors:

Lemma 1 (Main Lemma for Method 2). If m ∈ I r , and m has k prime factors (k ≤ r) and q = Q(r, k), one of the following two statements holds:

• there exists a prime number in p ∈ I s , s ≤ q and a number with k -1 prime factors t ∈ I r-s such that pt = m.

• there exists a prime number in p ∈ I s , s ≤ q and a number with k -1 prime factors t ∈ I r-s-1 such that pt = m.

Proof. Let m ∈ I r , m with k prime factors (2 ≤ k ≤ r) and let q = Q(r, k). We want to prove first that it has at least a prime factor p ∈ I s , s ≤ q. First of all we have that l = max{p prime : p ≤ 2 We can now enunciate the following algorith:

Algorithm 3 Sieve-Like Method (With output Ω)

1: Input: n.

2: Consider a vector B of size 2 n+1 .

3: for i=1 to n do 4:

for j=i-1 to 2 do This step must be skipped for i = 1

5:
r ← Q(i, j).

6:

for k=1 to r do 7:

For each p prime in I k and for each t with j -1 factors in I i-k , store pt as number with j factors.

Result is in I i or I i+1 .

8:

Store all remaining number in the interval I i as primes.

9: Output: B ≡ Ω for all numbers < 2 n and some numbers m:

2 n ≤ m < 2 n+1 .
If we assume that step 7 takes only one operation as set multiplication, the computational complexity is in the same order of magnitude with respect to:

n i=2 2 j=i-1 1 + i+1 j k=1 1 ∼ n i=2 2 j=i-1 1 + i + 1 j ≤ n i=2 (i -2 + (i + 1) log(i)) = O(n 2 log n).
Where we used again Abel's summation formula to evaluate the last sum. This computational cost is very convenient, since the algorithm assigns the function ω to all the numbers between 1 and 2 n . In practical configurations, the implementation could be much more expensive since set multiplication could be not feasable for regular computers. If we consider Prime Number Theorem and its generalization concerning π * k (x) (see introduction) to estimate roughly the number of prime numbers and numbers with distinct prime factors in the proper intervals (in particular we estimate the prime numbers in I k by approximating them by π(2 k+1 ) -π(2 k ) and numbers with j -i prime factors in I i-k in a similar way), we obtain:

n i=2 2 j=i-1 1 + i+1 j k=1 2 i (k -1) (log 2) 2 (k 2 + k)(j -2)! (log log 2 i-k+1 ) j-2 i -k + 1 - (log log 2 i-k ) j-2 i -k .
Being this sum really complex to estimate, we estimate a bigger sum doing several majorizations:

n i=2 2 j=i-1 1 + i+1 j k=1 2 i (k -1) (log 2) 2 (k 2 + k)(j -2)! (log log 2 i-k+1 ) j-2 i -k + 1 ∼ n i=2 i-1 j=2 1 + 2 i (log 2) 2 (j -2)! i+1 j k=1 1 ki (log(i log 2)) j-2 ∼ n i=2 i-1 j=2 1 + 2 i (log 2) 2 (j -2)! i+1 j k=1 1 ki (log i) j-2 * ∼ n i=2   (i -2) + 2 i i i-1 j=2 (log i) j-2 (j -2)! (log(i + 1) -log j)   * * ∼ n i=2   (i -2) + 2 i i i-1 j=2 (log i) j-1 (j -2)!   ∼ n i=2 i -2 + 2 i = O(2 n+1
).

Where we have used several majorations and approximations, as the logarithmic approximation of the harmonic series (*, see [START_REF] Nathanson | Elementary Methods in Number Theory[END_REF]) and Taylor expansion for the function e x (**, see [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]).

We want to remark that if we consider m = 2 n , the computational complexities are respectively O (log 2 m) 2 log log 2 m or equivalently O (log m) 2 log log m and O(2m), those amounts are smaller then the corresponding complexities of Eratostene's Sieve. However this algorithm uses more multiplications, that can be more intrinsecally expensive then additions, as used in the original sieve.

The memory consumption is given by the sum of the memory consumption for storing the values of Ω and conveniently, the number of the interval interval I n for each number. Both of these values can be approximated by 2 n+1 log 2 (n + 1), i.e., the number of numbers to be stored times the maximum number of bits needed to store the corresponding Ω function (resp. the interval), At most n + 1, see proposition 4 (resp. Remark 2).

The overall memory consumption is therefore proportional to O(2 n+1 log(n+1)). Again if we consider m = 2 n , the overall memory consumptions are proportional to O (m log (1 + log 2 m))) and O (m log log m)).

A few properties can make the implementation of the algorithm even more efficient. We suggest some of them, but other ideas could go in the direction of using the additional Ω values in the biggest interval computed by the algorithm to evaluate the function in the whole interval or avoiding redundant calculations by checking the range of the operations before they are performed:

Proposition 2. If m + 1 ≤ 2n then Q(m, n) = 1.
From the hyphothesis we have m + 1 n ≤ 2 and 2

m+1 n ≤ 2 2 = 4. Therefore l = max{p prime : p ≤ 2 m+1 n ≤ 4}
can contain only the primes 2 and 3, both belonging to I 1 . We have reached the thesis.

The following Proposition is a generalization of the previous one and can be proved similarly:

Proposition 3. If m + 1 ≤ n(t + 1) then Q(m, n) ≤ t.
This Proposition can avoid some calculations for the numbers with n not necessarily distinct factors in I n : Proof. By use of the Fundamental Theorem of Arithmetic we notice that 2 n is the smallest number with the required characteristics because the existence of a smaller number with this properties would imply the existence of a prime number smaller than 2. If we now notice that 2 n-1 • 3 is another number with the required number of factors in the correct interval, we must prove that the interval cannot contain numbers with factors bigger than 3 repeated more then once and factors bigger than 5. But this is a consequence of the following inequalities that consider p i , i ∈ {1, ..., n} as (not necessarily distinct) prime numbers:

p 1 • p 2 • ... • p n ≥ 2 n-2 • 3 • 3 = 2 n-2 • 9 > 2 n-2 • 8 = 2 n+1 , p 1 • p 2 • ... • p n ≥ 2 n-1 • 5 > 2 n-2 • 4 = 2 n .
We finally give some tables that show some data on how numbers are repetedly counted as output of the algorithm (as 12 = 2 • 6 = 3 • 4, where both 2 and 3 are indeed in I 1 ):

Third and Fourth Method

In this section we describe some modifications of Eratostene's method that rely on the calculation of the number of prime divisors for each number in the interval 1, ..., n. The output is again the function Ω:

Algorithm 4 Modified Eratostene's Sieve 2 (With output Ω)

1: Input: n.

2: Consider an array of natural numbers and size n, we call this C.

3: a ← 0.

4: p ← 2. 5: while p ≤ √ n do 6:
Mark the lowest number of the list p such that a(p) ∈ {0, 1} (2 in the first case) as prime.

7:

Evaluate the multiplicity of p for each number (which is multiple of p) up to at least n ( n p numbers).

8:

for i multiple of p do 9:

a(i) ← a(i) + multiplicity of p in i.

10: Output: C ≡ Ω for all numbers ≤ n.

By construction this method yields as output the Ω function for all the numbers in the array. We need a convenient method to evaluate the multiplicity v p (n) of p in n and complete the algorithm. One possible way is to use the following properties of v p (n):

Proposition 5. v p (n) = 1 for each p prime, n multiple of p such that n < p 2 .

Proof. This is obvious as a multiple of p has v p (n) ≥ 1 and v p (n) cannot exceed 2 as n < p 2 . Proposition 6. v p (n) = p(n + p l ) for each p prime, l, n ∈ N, n < p l .

Proof. let v p (n) = l, then there exists b coprime with p such that n = p l b. If now we consider n + p k = p l b + p k = p l (b + p k-l ) and since the last multiplicand is coprime with p in any case with n < p l , we have the thesis.

The following algorithm follows from the last Propositions. We can generate indeed the prime multiplicity of p for each multiple of p by copying properly some values. A small number of other operations is needed:

Algorithm 5 Prime Multiplicity Generation 1

1: Input: p,n.

2: Consider an array A of natural numbers and size n p .

3: b ← log p (n) .

4: assign 1 to the first p -1 entries of A.

Unless they exceed the size of A, in that case A ← 1 and exit. Assign the correct value (log p s) to the first unassigned entry s of A.

7:

Copy the first assigned numbers in the interval [1, p i ] p times.

8:

Delete the last number of the array A. change the last entry of the array to l.

12: Output: Multiplicity of p for all multiples of p smaller than n.

Another possibility, that relies on the multiplication of prime powers of every prime numbers in the interval by numbers which are not divisible by the prime number itself and store them with the exponent of the respective powers is the following:

Algorithm 6 Prime Multiplicity Generation 2

1: Input: p,n.

2: b ← log p (n) .

3: Consider an array Ā of natural numbers and size n p .

4: for i from 1 to b do This step could be empty.

5:

for j not multiple of p, p i-1 j ≤ n p do 6:

Ā(p i-1 * j) = i.

7: Output: Multiplicity of p for all multiples of p smaller than n.

The computational complexity of the incomplete algorithm Modified Eratostene's Sieve 2 is again in the order of O(n log log n), as the cost of the original Eratostene's method. The two methods can be overlapped. In the first case (using Prime Multiplicity Generation 1 as step 7) the additional complexity is given by:

p≤n   1 + log p (n) i=2 (2 + p)   = p≤n 1 + (2 + p)( log p (n) -1) = p≤n -1 -p + (2 + p) log p (n) = O n 2 log n ,
which makes the first part complexity of the algorithm negligible. This estimate has been made by using Abel's summation formula and Prime Number Theorem [START_REF] Apostol | Introduction to Analytic Number Theory[END_REF]. In particular the dominant term of the sum is

p log p (n) .
In the second case (Second algorithm as step 7), the computational complexity is:

p≤n   1 + log p (n) i=1 j =kp,p i-1 j≤ n p 1    ∼ p≤n   1 + np p -1 log p (n) i=1 1 p i    * * * ∼ p≤n   1 + np p -1    1 -1 p (1+log p n) 1 -1 p -1       ∼ p≤n 2 + -p + 1 p - 1 n = O n 2 log n .
We have used an identity concerning geometric series (***), several approximations and again Abel's Identity.

The amount of bits used is given by the logarithm of the maximum size of Ω in the considered interval times the dimension of the interval. If we consider Ω(n) = O(log n) [START_REF] Mathworld | DistinctPrimeFactors[END_REF], the final bit consumption is O(n log log n).

These final estimates take into account the fact that the memory used from Modified Eratostene's Sieve 2 is at most in the same order of magnitude.

Methods with output d

We can modify Algorithm 4 so that the output is the function d, we only change step 9 so the final output is multiplied instead of being summed:

Algorithm 7 Modified Eratostene's Sieve 3 (With output d)

1: Input: n.

2: Consider an array of natural numbers and size n, we call this C.

3: a ← 0. Mark the lowest number of the list p such that a(p) ∈ {0, 1} (2 in the first case) as prime.

7:

Evaluate the multiplicity of p for each number (which is multiple of p) up to at least n ( n p numbers).

8:

for i multiple of p do 9:

a(i) ← a(i) • (1+multiplicity of p in i).

10: Output: C ≡ d for all numbers ≤ n.

The computational complexity is comparable with those used from the original method.

The maximum memory consumption is given by Gronwall's approximation [START_REF] Gronwall | Some asymptotic expressions in the theory of numbers[END_REF] d(n) = O(n log log n) and it is in the same order of magnitude with O(n log n).

3 :

 3 p ← 2.

  in particular we have p ≤ 2 r+1 k and if we take p = min{s prime, s > p}, we have from the definition of p, p > 2 r+1 k . If by absurd the smallest prime of m is p, we have m = p 1 • ... • p k ≥ pk > 2 r+1 k •k = 2 r+1 . We have obtained an absurd since m would not belong anymore to the interval of the Hypotheses. To prove the remaining part of the Lemma, we first notice that if we use the Fundamental Theorem of Arithmetic is obvious that if pt = m then t must have exactly k -1 prime factors. The final part of the lemma follows from Proposition 1. It is also useful to make the following definition: Definition 3 (Q function). Let m, n ∈ N, we define the function Q(m, n) as follows: let l = max{p prime : p ≤ 2 m+1 n }, then there exist unique r ∈ N such that l ∈ I r , we define Q(m, n) := r or equivalently (Remark 2) Q(m, n) := log 2 l .

Proposition 4 .

 4 Each set I n contains exactly 2 numbers with n prime factors (Numbers m such that Ω(m) = n), namely 2 n and 2 n-1 • 3.The same set I n contains no numbers with at least n + 1 prime factors.

5 :

 5 for i from 2 to b do 6:

9:

  Copy the first assigned numbers in the interval [1, p b ] until the number n p is reached.10: if n p = p l for some l ∈ N then 11:

4: p ← 2 . 5 :

 25 while p ≤ √ n do 6:

  It is clear that at the end of the algorithm the array a is coinciding with the function ω because every number of the array is increased by one (instead of being marked as true or false or canceled like in the original sieve) for each prime number lower or equal. This method still uses O(n log log n) additions as it does not contain modification with respect to the original algorithm, but a bigger amount of memory has to be assigned. We have to store ω(k) ∀i ∈ 1, 2, ..., n so the required maximum amount is n log 2 max

			1≤k≤n	ω(k)	that can be
	bounded by O n log	log(n) log log(n)	where we have used the fact that ω(n) is at most O	log(n) log log n	[8].
	8: Output: A ≡ ω ∀l ∈ {1, 2, ..., n}.

Table 1 :

 1 Repeated output from Algorithm 2.

	Numbers with repeated output 12 18 24 30 36 42 45 48 50 54 60 ...
	Number of repetitions	2	2	2	2	2	2	2	2	2	2	2 ...
	Corresponding Ω	3	3	4	3	4	3	3	5	3	4	4 ...

Table 3 :

 3 Smallest numbers with amount of repetitions from Algorithm 2

	Smallest number with higher number of repetitions Number of repetitions Corresponding Ω
	2	1	1
	12	2	3
	70	3	3
	1050	4	5
	182490	5	6
	1771770	6	7

Last affiliation: University of Kiel (Germany). ORCID: https://orcid.org/0000-0002-5495-6863.

Acknowledgements

The author of this paper wishes to thank Roberto Amato (University of Messina), Francesco Belardo (University of Naples), Enzo Maria Li Marzi (University of Messina) and Thomas Slawig (University of Kiel) for their valuable mentorship. The author wishes to thank also Daniela Giumbo (University of Messina) for the valuable suggestions that led to the improvement of this manuscript.

Conclusions

The shown methods can evaluate several functions in a way that the primality of any number is obtained, together with other informations. Of course we can have algorithms that yield only prime numbers as output, as in classical Sieve Theory, by considering that prime numbers are characterized by ω(p) = 1, Ω(p) = 1 and

In the end we have achieved additional precision at a bigger cost for memory consumption. Modern machines will be able to front this cost more and more easily. These methods, however, are not published with the aim to be the best possible, but to give new directions in the study of sieve theory.