Valentin Dupif
email: valentin.dupif@onera.fr

Marc Massot
email: marc.massot@centralesupelec.fr

Joël Dupays
email: joel.dupays@onera.fr

Frédérique Laurent
email: frederique.laurent@centralesupelec.fr

Clément Le Touze
email: clement.le_touze@onera.fr

On the influence of the numerical strategy on the predictive character of Euler-Euler model for two-phase flow simulations in solid rocket motor instabilities

Keywords: Solid rocket motor, Vortex shedding, Two-way coupling, Euler-Euler simulation, Velocity polydispersion, Anisotropic Gaussian, Realizable high order schemes

Solid Rocket Motors involve strongly coupled two-phase flow. The presence of a polydisperse spray of particles resulting from the combustion of aluminized propellant has been shown to have a strong impact on stability and can eventually yield thrust oscillations. The ability to conduct predictive simulations of such a harsh environment is highly desirable. Euler-Euler models relying on moment methods for the disperse phase constitute interesting approaches due to their efficiency at coupling both phases and their ability for high performance computing. A multi-fluid model coupled to a new numerical strategy for the disperse phase is introduced in order to cope with the natural high stiffness of the resulting systems of conservation laws. The predictive character of the method is strongly related to the possibility of using accurate methods while preserving stability and robustness in the presence of intrinsic singularities occurring in the disperse phase equations. The purpose of this contribution is to stress the impact of several numerics and modeling on the solution. Relevant test cases for solid propulsion involving hydrodynamic instabilities and acoustic coupling are presented. A strategy is proposed in order to produce reliable predictions.

Introduction

Among the space propulsion technologies, Solid Rocket Motors (SRM) are known for their reliability, high trust and ability to be stored during a long period of time. Included their relatively low cost, such engines are first choice possibilities for the first stages of present and future space launchers [START_REF] Guery | Solid propulsion for space applications: An updated roadmap[END_REF]. In order to improve their efficiency, or more particularly their specific impulse, aluminum powder is added to the propellant grain increasing the temperature of the burnt mixture. From this technical solution results the presence of a massive amount of aluminum oxide particles that strongly interact with the internal flow of the motor. Since these particles have a notable impact on thermoacoustic, combustion and hydrodynamic instabilities or combination of these [START_REF] Kuentzmann | Instabilités de fonctionnement des systèmes propulsifs[END_REF], it is a key issue to be able to prevent such phenomenon in the early development of solid rocket motor and thus avoid deleterious effects on the launcher mission.

Due to the harsh conditions existing in the internal flow, measurement possibilities are limited and simulations are mandatory. In addition to the experimental investigations to understand the complex physic of such flows, the ONERA continually refine the modeling of both carrier and disperse phase while improving associated numerical methods of CFD issues. This contribution exhibit the recent progress in the disperse phase modeling including its two-way interaction with the gaseous carrier phase obtained thanks to a long term collaboration with the EM2C laboratory of CentraleSupelec [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF][START_REF] Sibra | Modélisation et étude de l'évaporation et de la combustion de gouttes dans les moteurs à propergol solide par une approche eulérienne Multi-Fluide[END_REF]. This work is focused on the Eulerian representation of polydisperse [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF][START_REF] Laurent | Two-size moment multi-fluid model: a robust and high-fidelity description of polydisperse moderately dense evaporating sprays. To appear[END_REF] spray giving solution for realistic flow of particles taking into account a large spectrum of size of droplets [START_REF] Emre | Eulerian moment methods for automotive sprays[END_REF] and more recently statistical droplet trajectory crossing [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF][START_REF] Sabat | Fully Eulerian simulation of 3D turbulent particle laden flow based on the Anisotropic Gaussian Closure[END_REF]. The ultimate goal of these works is to produce models that can efficiently represent particles from low to high inertia. While being mathematically well posed, such models are exposed to severe non-linearities and have to rely on dedicated numerical methods. In order to aim at industrial configurations, an active research is conducted on accurate and robust numerical methods for unstructured mesh [START_REF] Le Touze | Multislope MUSCL method for general unstructured meshes[END_REF] constantly improving the reliability of computed solutions. In our recent work, a discriminant test case under a realistic configuration for SRMs is investigated and will be the key point of this paper. Highly sensitive to the numerical method, this fictive SRM, namely C1 [START_REF] Lupoglazoff | Numerical simulation of vortex-shedding phenomena in 2D test case solid rocket motors[END_REF][START_REF] Morfouace | Two-phase flow analysis of instabilities driven by vortex-shedding in solid rocket motors[END_REF][START_REF] Kourta | Computation of vortex shedding in solid rocket motors using time-dependant turbulence model[END_REF], possesses an intrinsic purely hydrodynamic instability that can be either excited, reduced or even suppressed by the single presence of particles under a condensed form. Through this new analysis of the problem, we propose to distinguish the spurious phenomena associated to the numerical method from those attributable to the model.

The contribution is organized as follow. Section 2 is devoted to the modeling of the disperse phase according to the briefly described SRM's internal flow conditions. From the closure of the kinetic description of the particles, two Eulerian models are deduced and will form our governing equations solved thanks to the numerical methods proposed in section 3. Special care is taken of the design of robust schemes conserving the mathematical properties of the models while keeping a high precision. Then Section 4 introduces the two-phase C1 test case that will assess the relia-bility of models and methods. The ability of our configuration to be predictive will be the subject final discussions leading to the conclusion of the paper.

Model

The internal flow of a SRM involves complex phenomena that can interact with each other. This includes for instance combustion, turbulence and radiative heat transfer, but in this paper, acoustics and two-phase flow issues will be our only concern. This choice of simplicity is motivated by the two-phase C1 test case that will be introduced in another Section 4. We then focus on the description of the disperse phase and its interaction with the carrier phase. After arguing on the choice of an Euler-Euler model, two kinetic closures are proposed for the disperse phase and the governing equations obtained thanks to a method of moment.

Eulerian model of the flow

Combustion and multi-component gas issues are omitted to simplify the study. The compressible form of the well known Navier-Stokes equations is retained. Since it has been proven that LES models have a marginal impact in the case of a C1 configuration without particles [START_REF] Kourta | Computation of vortex shedding in solid rocket motors using time-dependant turbulence model[END_REF], the equations will be solved under a DNS context.

Combined with their high material density, the small size of the droplets justify the assumption of point-particle. Through a diameter lying between several hundred of nanometers and few hundred of micrometers, their volume fraction in the Navier-Stokes equations can be neglected. Thus only two kinds of interactions between the carrier and carried phase are considered, namely the drag force and the heat exchange. According to [START_REF] Doisneau | Modeling and simulation of polydisperse moderately dense coalescing spray flows with nanometricto-inertial droplets: application to Solid Rocket Motors[END_REF], the Reynolds number of the droplets remain under 500 even through the nozzle, where it is the highest, and classical Schiller-Naumann and Ranz-Marshall corrections are sufficient to describe the dynamics of an isolated particle. These take the following form:

   F St (ug, c, dp) = (ug -c) τ st u , τ st u = ρ l d 2 p 18µg F = F St (1 + 0.15Re 0.687 p ) = (ug -up) τu (1)
H St (Tp, Tg, dp) = c p,l Tg -Tp τ st T , τ st T = 3 2 P rg c p,l cp,g τ st u H = H St (1 + 0.3Re 1/2 p P r 1/3 g ) = c p,l Tg -Tp τ T (2)
where F and H are given by unit of mass, ug and c are the velocities of the gas and the particle, ρ l is the density of the droplet material, dp the diameter of the droplet, µg the dynamical viscosity of the gas, Rep the Reynolds number of the particle, P rg the Prandtl of the gas and finally cp,g and c p,l the heat capacity at constant pressure of the gas and the droplet material.

The number of particles, that can easily reach 10 5 per cubic centimeter, induces a non negligible collisions frequency [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. Such a problem is out of the scope of the paper and will be neglected. The droplets are then consider only to be able to interact between themselves through their coupling with the gaseous phase.

At this point, two classes of models for the resolution of the disperse phase can be considered. It can be described through a Lagrangian method where droplets are tracked, alone or grouped, in their displacement or thanks to a Eulerian framework by describing droplets through local quantities. For phenomena driven by two-way interactions as in the SRMs, one may prefer the Eulerian framework. As a matter of fact such method can accurately calculate droplet-gas interaction with a precision that does not depend on their number while avoiding statistical convergence issues. This last property is a key advantage for two-way coupled unsteady simulations. Moreover, Eulerian closure can be efficiently parallelized through domain decomposition whereas Lagrangian methods quickly suffer from memory exchange overhead in the same situation.

Then, assuming the statistical convergence through an infinite number of particles, the spray can be described thanks to a Number Density Function (NDF) f that depends on the time t, position x, velocity c and temperature T . This simplified form of the Williams-Boltzmann equation assumes a unique particle size without evaporation nor collision terms and take the following form:

∂tf + ∂x.(cf ) free transport + ∂c.(Fnf ) drag force + ∂T ( Hn c p,l f ) thermal transfert = 0 (3)
where Fn and Hn are the drag force and heat transfer applied to a single droplet. This model is not solved directly but approximated through moment methods.

Kinetic based moment method

Equations are written on moments in velocity and temperature. Despite our choice of considering only a discrete size of particles, we wish to be able to extend the presented work to two size-moment hybrid methods [START_REF] Boileau | Two-size moment Eulerian multi-fluid method describing the statistical trajectory crossing: modeling and numerical scheme[END_REF][START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF][START_REF] Laurent | Two-size moment multi-fluid model: a robust and high-fidelity description of polydisperse moderately dense evaporating sprays. To appear[END_REF]. These methods use a size discretization into sections [START_REF] Kuentzmann | Instabilités de fonctionnement des systèmes propulsifs[END_REF] and two size-moments in each section and have reach a mature level but are out of the scope of this paper (see [START_REF] Emre | Eulerian moment methods for automotive sprays[END_REF][START_REF] Laurent | Two-size moment multi-fluid model: a robust and high-fidelity description of polydisperse moderately dense evaporating sprays. To appear[END_REF][START_REF] Boileau | Two-size moment Eulerian multi-fluid method describing the statistical trajectory crossing: modeling and numerical scheme[END_REF] and references therein for a broader view on the literature in this field). Then, the number density n = f dcdT dS and the mass density m are considered, despite their redundancy, as proof of feasibility and straightforward extension of the presented work to such literature. Therefore, the particle diameter dp will not be considered as a fixed parameter but as deduced from n and m thanks to the formula dp = [(6m)/(ρ l πn)] 1/3 .

We propose then to close problem thanks to assumptions on the NDF according to the velocity and temperature at a time t and position x. Assuming that these two local variables can be closed independently, the NDF can be stated as follow:

f (t, x, c, T ) = n(t, x)fc(t, x, c)fT (t, x, T ) (4)
Assuming this explicit form of f , conservation laws can be obtained thanks to the transport equation ( 3): from the integration of (3), one can obtain sets of equations with usual physical meanings. Thus, the moments of order 1 and 2 in velocity lead to the conservation of momentum and kinetic energy, the moment of order 1 in temperature leads to the conservation of the thermal energy and the zeroth order moment to the conservation of particle number and mass.

For the sake of simplicity, a discrete temperature of the particle Tp at a point is assumed, which leads to fT (t, x, T ) = δ (T -Tp(t, x)). After the integration of (3) multiplied by the temperature T and the heat capacity c p,l , one can obtain the equation of the transport of the thermal energy h = c p,l Tp such as:

∂t (mh) + ∂x • (mhu) = mc p,l Tg -Tp τT ( 5
)
where u is the average velocity obtained for the first moment in velocity u = n -1 cf dcdT dS. The form of the velocity distribution will be a key point of this contribution. It is possible to consider locally a unique velocity such as fc(t, x, c) = δ (c -u(t, x)), which will be called the monokinetic (MK) assumption [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. This hypothesis results in the equations of Pressureless Gas Dynamic (PGD) with a drag source terms on the right hand side. Combined with the conservation of energy [START_REF] Dupays | Contribution à l'étude du rôle de la phase condensée dans la stabilité d'un propulseur à propergol solide pour lanceur spatial[END_REF] and particle number, one can obtain the system that follow:

         ∂tn + ∂x • (nu) = 0 ∂tm + ∂x • (mu) = 0 ∂t (mu) + ∂x • (mu ⊗ u) = m ug -u τu ∂t (mh) + ∂x • (mhu) = mc p,l Tg -Tp τ T (6)
The monokinetic assumption finds its limits for relatively inertial droplets for which Particle Trajectory Crossing (PTC) [START_REF] Doisneau | Modeling and simulation of polydisperse moderately dense coalescing spray flows with nanometricto-inertial droplets: application to Solid Rocket Motors[END_REF] can occur. In other word, the particles belonging to the same system of equations cannot cross each other at a given point. This uniqueness of the local velocity created in case of PTC a Dirac delta function in density called δ-shock. To avoid such drawback, recent works inspired by the theory of rarefied gas [START_REF] Levermore | The Gaussian moment closure for gas dynamics[END_REF] proposed a multivariate Gaussian distribution in velocity [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF][START_REF] Sabat | Eulerian models and realizable numerical schemes for the description of low-to-medium inertia polydisperse sprays in turbulent two-phase flows[END_REF][START_REF] Sabat | Fully Eulerian simulation of 3D turbulent particle laden flow based on the Anisotropic Gaussian Closure[END_REF]

fc(t, x, c) = det(Σ) -1/2 (2π) 1 2 N d exp(-1 2 (u -c) T Σ -1 (u -c))
where Σ is the covariance matrix of the velocity. Using the moment in velocity up to the order 2, one can obtain:

             ∂tn + ∂x • (nu) = 0 ∂tm + ∂x • (mu) = 0 ∂t (mu) + ∂x • (mu ⊗ u + P) = m ug -u τu ∂t (mE) + ∂x • ((mE + P) ∨ u) = m ug ∨u-2E τu ∂t (mh) + ∂x • (mhu) = mc p,l Tg -Tp τ T ( 7
)
where E is the energy matrix and P the granular pressure tensor defined by P = mΣ = 2mE -mu ⊗ u. This second system, called AG for Anisotropic Gaussian, is an extension of the MK system. As a matter of fact, the system (6) based on the MK closure is contained in [START_REF] Gottlieb | On high order strong stability preserving Runge-Kutta and multi step time discretizations[END_REF] in the limit of zero velocity dispersion. A first analysis of the left hand side of the system [START_REF] Gottlieb | On high order strong stability preserving Runge-Kutta and multi step time discretizations[END_REF] shows that the nature of this system of this Partial Differential Equations (PDE) is hyperbolic as the compressible Navier-Stokes chosen for the gas. These equations require specific methods of resolution since they contain intrinsic singularities such as shock waves. In the case of the system (6), the left hand side is a weakly hyperbolic PDE, which means that the Jacobian of the matrix is not invertible. From this specificity, result the strong singularity that is the δ-shock. The remaining part of the equations that is the right hand side of the chosen system and its backward action on the carrier phase act as a simple system of ODEs. Since both models possess their own properties, they have to rely on adapted numerical methods in order to produce accurate solutions.

Numerical methods

For the sake of genericity, the involved equations ( 6) or [START_REF] Gottlieb | On high order strong stability preserving Runge-Kutta and multi step time discretizations[END_REF] and the Navier-Stokes equations are split into three operators. The right hand side of the equations (6) ( or ( 7) ) as well as their action on the carrier phase will be treated as a source term operator S and can be written as an ODE dtU = S(U ) where S is the source term function acting on the conservative vector U . In the same way Tg, Tp will be the transport operators corresponding to the usual Navier-Stokes equations and the left hand side of the equations (6) ( or [START_REF] Gottlieb | On high order strong stability preserving Runge-Kutta and multi step time discretizations[END_REF] ) that can take the generic form ∂tU + ∂x • f (U ) = 0, where f is the flux function. After describing our splitting strategy as well as the treatment of the source terms, two methods for the resolution of the transport operators are proposed.

Operator splitting and source terms

Through a Strang splitting, the integration of the source terms is conducted independently from the resolution of the transport operators. Expecting steep source terms, we choose to use the operator S on half a time-step before solving independently the transport operators Tg and Tp since they are uncorrelated at this point and then conclude by half a time of source terms. Under a mathematical formulation, this gives:

U Strang (t+∆t) = S ∆t 2 [Tg +Tp](∆t)S ∆t 2 U (t) (8)
For the integration of the source terms, the use of a third order Strong Stability Preserving Runge-Kutta (SSPRK) time integration [START_REF] Gottlieb | On high order strong stability preserving Runge-Kutta and multi step time discretizations[END_REF] is known to be a reliable choice. Combined with the Strang scheme, it leads to an efficient way to treat the source terms.

MacCormack method

Among the second order methods to solve hyperbolic systems of equations, the MacCormack method, of Lax-Wendroff class, has been very popular during the 80's and 90's thanks to its easiness of implementation and numerical efficiency. In an one dimensional framework, the scheme takes the compact predictorcorrector formulation that follow:

U n i = u n i + ∆t ∆x (f (U n i ) -f (U n i-1 )) U n+1 i = 1 2 U n i + U n i + ∆t ∆x (f (U n i+1 ) -f (U n i ) (9)
where U n i is the solution at the predictor step that is here left oriented.

To our knowledge, every multidimensional application of the MacCormack method is linked to a structured mesh and a direction of the predictor according to the grid directions such as in [START_REF] Lupoglazoff | Numerical simulation of vortex-shedding phenomena in 2D test case solid rocket motors[END_REF]. The experience shows that this predictor should be oriented upstream which for SRMs means toward the propellant grain in the radial direction or toward the head-end in the axial direction. Such choice reduces the spurious oscillations but does not solve this issue. To ensure the quality of the results, an artificial viscosity of Swanson and Turkel [START_REF] Swanson | On central-difference and upwind schemes[END_REF] is chosen. The sensor is based on the pressure for the carrier phase but since this value is absent of the monokinetic closure and irrelevant for the AG, the density field is chosen instead. Such method fits well the requirement for the gaseous phase since purely acoustic phenomenon are expected to be observed and shocks avoided inside the internal flow. However the possibility of some vacuum area and discontinuity in the field of particles can lead to stability issues. Despite the solution should be smooth in most of the chamber, an high artificial viscosity is required to manage local singularities and vacuum. In consequence, the ability for this scheme to solve the disperse phase is ambiguous.

MUSCL method

Form the idea of solving discontinuities at the cell interfaces in a finite-volume context, the Godunov scheme class is proven to be able to efficiently tackle singularities. The procedure consists in estimating in a first time the value at the cell interface, then solving Riemann problems at the cell interfaces to deduce the flux while ensuring the entropy condition. In a one-dimensional framework the scheme takes the form:

U n+1 i = U n i + ∆t ∆x f (U Ri i+1/2 ) -f (U Ri i-1/2 ) U Ri i+1/2 = U Ri ( x t = 0, U l i+1/2 , U r i+1/2 ) (10)
where U Ri i+1/2 is the solution to the Riemann problem between the left state U l i+1/2 and right state U r i+1/2 solved by the procedure U Ri .

At order one, piecewise constant values are assumed in each cell equal to the average quantity contained in the corresponding volume and the extension to multidimensional frameworks is straightforward. The resolution of the Riemann problem finds an explicit solution for the PGD equations [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF]. However, such a procedure is tedious for more complex systems and we will rely on a HLLC approximate Riemann solver for the gas [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF] and the HLL solver proposed in [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF] for the AG equations. Each of them preserves the realizability as long as the reconstructed values are realizable states which is obvious at first order.

To obtain higher order schemes of this kind, the distribution of the value inside the cell have to be estimated. It can be achieved at the second order, for example, through a linear reconstruction leading to methods called MUSCL. For multidimensional frameworks, several procedures claim to be of MUSCL type. Among them the multislope method recently developed at ONERA [START_REF] Le Touze | Multislope MUSCL method for general unstructured meshes[END_REF] provides an efficient and robust Local Extrema Diminishing (LED) reconstruction [START_REF] Jameson | Positive schemes and shock modelling for compressible flows[END_REF]. When CFL independent, this LED property ensures that the reconstructed value will be a convex combination of its neighbors and will not create any artificial extrema. The key issue to obtain the realizability of such LED scheme, once the Riemann solver is proven realizable, is to select the right primitive value for the reconstruction that ensure the realizability of the deduced state.

LED reconstructions of the density and pressure fields ensure the positivity of these variables. Thus for the carrier phase, the density, the velocity and the pressure are usual choices ensuring the robustness. For the monokinetic closure, the choice of density and the velocity as reconstruction variables ensures the realizability [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF]. Compared to the PGD equations, the AG closure needs, in addition, the positivity of the σij component of the covariance matrix as well as the positivity of its determinant σiiσjj -σ 2 ij > 0 ∀i = j. As proposed in [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF] the diagonal terms σii are use for the reconstruction to ensure their positivity. However, we introduce here the new reconstruction value σij/ √ σiiσjj that have to be kept in [0, 1[ to ensure the previous constraints. The value of σij deduced from this variable is then sure to fulfill the realizability requirement thanks to the LED procedure for its reconstruction and the positivity σii and σjj as already provided. Finally, the knowledge acquired in the field of size polydispersion [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF] suggests that the reconstruction over h and n/m prevent from any spurious border effect on the temperature and the size deduced from n and m.

Thanks to their design, these upwind methods are built-in realizable and robust with a L∞ stability proven under CFL condition. This property is preserved through the second order SSPRK [START_REF] Gottlieb | On high order strong stability preserving Runge-Kutta and multi step time discretizations[END_REF] time integration used to keep the convergence orders in space and time equal. In contrast to the MacCormack scheme where a case dependent manual tuning of the artificial viscosity is needed, these method lead to the exact solutions through mesh refinement without any artificial input. This property will be crucial for the analysis of the numerical solutions.

Two-phase C1 analysis

In industrial applications, propellant grain usually have sharp edges. The 2D planar C1 test case has been designed in order to reproduce a critical situation due to this geometrical singularity and asses the robustness of models and numerical methods. After having accurately described the case studied, the specific case where 20µm particles are injected is investigated.

Test case presentation

Downstream geometrical discontinuities on the propellant grain, hydrodynamic instabilities are usually developing and produce a phenomenon called Angle Vortex Shedding (AVS). In SRMs, such unsteady flow can be coupled with the acoustic of the chamber. The original design of the C1 aimed at exciting the first axial acoustic mode in single phase flow thanks to the AVS [START_REF] Lupoglazoff | Numerical simulation of vortex-shedding phenomena in 2D test case solid rocket motors[END_REF]. To do so, a sharp discontinuity is placed at the center of the combustion chamber. Unexpectedly however, various numerical simulations [START_REF] Morfouace | Two-phase flow analysis of instabilities driven by vortex-shedding in solid rocket motors[END_REF][START_REF] Kourta | Computation of vortex shedding in solid rocket motors using time-dependant turbulence model[END_REF] proved that the second axial mode is excited providing a distinctive monochromatic pressure oscillation at the head-end, independently of the numerical method. In a preliminary investigation, both MacCormack and MUSCL methods are used on the single phase case of [START_REF] Lupoglazoff | Numerical simulation of vortex-shedding phenomena in 2D test case solid rocket motors[END_REF] on the 317 × 30 cells mesh presented in Figure 1. In accordance with [START_REF] Lupoglazoff | Numerical simulation of vortex-shedding phenomena in 2D test case solid rocket motors[END_REF][START_REF] Morfouace | Two-phase flow analysis of instabilities driven by vortex-shedding in solid rocket motors[END_REF][START_REF] Kourta | Computation of vortex shedding in solid rocket motors using time-dependant turbulence model[END_REF] we find a monochromatic instability around 2650Hz. Negligible differences on this value are obtained using the Mac-Cormack and MUSCL methods. Moreover, a mesh convergence analysis until a 2536 × 240 mesh leads to the conclusion that the results using the standard mesh already were of high quality and can be considered as converged. In [START_REF] Morfouace | Two-phase flow analysis of instabilities driven by vortex-shedding in solid rocket motors[END_REF], particles are injected such that the particle to gas mass ratio ṁp ṁg = 0.396 corresponds to a propellant using 15 % of aluminum powder. Such high presence of particles have a notable impact on the instability amplitude and frequency. As an example, the injection of 5µm particles creates a shift in frequency to approximately 2250Hz as well as a reduction by almost the half of the oscillation magnitude. The development of the corresponding gas rotational field is presented in Figure 2. Using the two-phase acoustic theory of [START_REF] Temkin | Attenuation and dispersion of sound by particulate-relaxation processes[END_REF], one can estimate that the injection of 16µm particles will attenuate the AVS with the highest efficiency. Numerical experiments using both a MUSCL type of scheme [START_REF] Morfouace | Two-phase flow analysis of instabilities driven by vortex-shedding in solid rocket motors[END_REF] and MacCormack scheme [START_REF] Dupays | Contribution à l'étude du rôle de la phase condensée dans la stabilité d'un propulseur à propergol solide pour lanceur spatial[END_REF] on the mesh presented in Figure 1 reported that the instability completely vanish for this size. This behavior can be observed on a close range of particle diameter around to 16µm. For the numerical experiment presented hereafter, 20µm particles are injected. For such size, the detection and occurrence of AVS is not ensured creating a discriminant situation.

Solutions using the monokinetic closure

In order to validate the results obtained on the standard mesh of the cited paper, we conduct a mesh refinement for the various numerical strategies. To stay consistent with the previous investigation, the equations of pressureless gas dynamics are used. A first test using a MacCormack scheme for both phases as in [START_REF] Dupays | Contribution à l'étude du rôle de la phase condensée dans la stabilité d'un propulseur à propergol solide pour lanceur spatial[END_REF] indicates a steady engine. The flow obtained exhibits a local high mass concentration near the symmetry axis and close to the nozzle inlet.

Among the possibilities at our disposal, we chose the combinations of upwind schemes (first order and MUSCL) for the droplets and MUSCL or MacCormack schemes for the carrier phase. It is important to note that the computation involving the MacCormack scheme for gas has been conducted with the inhouse code SIERRA whereas the semi-industrial software CE-DRE has been used for the fully upwind configurations. The result on the most refined mesh has been conducted with CEDRE since it allows us to finalize the simulation within a reasonable computational time. The characteristics of the AVS, when detected, are given Table 1.

The results shows that the AVS is only detected using an accurate method for the carrier phase and the disperse phase only. If the first order method is known to be highly diffusive, the Mac-Cormack scheme is impacted by the artificial viscosity, which is needed for stability issues. Its effect is especially important close to high density gradient and vacuum due to the sensor chosen. At the same locations, the multislope MUSCL method accurately solves the flow avoiding spurious and dissipative effects. We conclude that combined with the physical dissipation due to the presence of the particles, the resulting numerical dissipation has a strong impact and can make the AVS disappear. Aside the numerical strategy for the disperse phase, the impact of the carrier phase resolution can be observed when the AVS is detected. Since this MUSCL version is LED and thus dissipates at extrema, the pressure amplitude obtained are lower than for the MacCormack scheme. Thanks to the low dissipation of the setup, a tiny instability can be detect thanks to the coarsest mesh but only correctly resolved on a refined one. Besides, the frequency observed indicates that two distinct phenomena are detected. According to the previous study, frequencies close to 2150Hz matching the second axial mode are expected. However the 1450Hz instability does not correspond to any acoustic mode of the chamber. A closer look to the density field in that case indicates the occurrence of δ-shocks that are consequences of PTC, which can not be represented using a monokinetic model. These singularities take their origin in two distinct phenomena. On the one hand, the nozzle geometry leads to the crossing of the particles coming from the head-end of the engine with those deviated by the convergent in a quasi-stationary crossing phenomenon. On the other hand, the AVS causes the ejection of the particles from the vortices and thus their impact on the wall and symmetry axis as well as the dynamical crossing of their trajectories as it can be seen in Figure 3. Even if further studies are to be conducted in order to obtain firm conclusion, we can anticipate that this last phenomenon is the main cause of the shift in frequency.

Solutions using the anisotropic Gaussian closure

In order to go over the influence and limitations of such singularities, which are only due to model limitations, we propose to replace the MK equations by the AG model. Due to the additional realizability constraints compared to the equations of pressureless gas dynamics, only the MUSCL scheme can be used for the resolution of the disperse phase. Moreover since this new model is not yet implemented in CEDRE, the in-house code SIERRA is the only one which can be used and thus the MacCormack scheme solves the carrier phase in our new configuration. We expect reliable results from this combination.

A first test on the coarsest mesh does not conduct to the detection of AVS. The possibility of the new model to take into account trajectory crossing results in a spatial redistribution of particles, thus avoiding artificial particle high concentrations, but thus leading to an effect, which is close to the previous artificial dissipation. The velocity dispersion is mainly observed at the stationary crossing around the nozzle and in a more marginal way at the symmetry axis.

Thanks to the refined mesh (634 × 60), the simulation are conducted starting from an initial condition involving the MK closure. The action of the velocity dispersion when PTC is encountered can be seen in Figure 4. We observe an immediate transition to a new unsteady regime of approximately 2100Hz at a lower average pressure with √ 2RM S = 8mbar. Since the state of the art in this field of research clearly states that the magnitude of oscillations are hard to predict [START_REF] Sibra | Modélisation et étude de l'évaporation et de la combustion de gouttes dans les moteurs à propergol solide par une approche eulérienne Multi-Fluide[END_REF], we are not for now able to clearly assess the obtained improvement. However, the fields resolved provide less artifact of the model, which supposes an higher reliability of these new results, even if further detailed analysis are to be conducted in order to draw some firmer conclusions on the coupled effects of modeling and numerical strategy.

Conclusion

In this contribution, we have introduced both a new numerical strategy as well as a new Eulerian model, able to cope with PTC while preserving a well-posed mathematical system of PDEs for which an entropy inequality guarantees a proper treatment of natural singularities, which are inherent to fully Eulerian modeling of disperse two-way coupled two-phase flows. A discriminant test-case in the field of solid propulsion has been adequately chosen in order to conduct a study on the influence of both modeling and numerical methods.

We conclude this paper by observing the critical influence of both the modeling and the numerics. The interest of realizable methods has been pointed out since they are the only ensuring reliable solution when singularities occur, even locally. As far as the monokinetic model is concerned, the strong singularities resulting from the inertia of the particles and from PTC leads to strong singularities, δ-shocks, the effect of which can not completely be dealt without using adequate numerical methods. The AG level of modeling provides an interesting alternative in order to avoid such singularities but the gain in predictivity for actual SRMs simulations still requires further investigations.

While the detailed studies we are conducting are of real importance in the C1 configuration in order to draw some firm conclusions of the influence of both modeling and numerical methods for the prediction of instabilities, comparisons and validations relying on experimental measurements are necessary. We then have to switch to at least 2D axi-symmetrical configurations or even 3D axi-symmetrical configurations and this requires some further development in terms of models and realizable numerical methods, which is work in progress. It should provide a clear assessment of the improvement of the predictive character of the simulations using a tailored strategy involving the right level of modeling resolved using robust and accurate parameter-free numerical simulations.

Figure 1 :

 1 Figure 1: 317 × 30 mesh of the C1 with inlet indication

Figure 2 :

 2 Figure 2: Gas rotational field evolution with 5µm particles and MK model at time (up to bottom) t = 0µs, t = 147µs, t = 294µs

Figure 3 :

 3 Figure 3: Dynamic PTC in the convergent

Figure 4 :

 4 Figure 4: σ22 velocity dispersion in the nozzle

Table 1 :

 1 Instability characteristics depending on the numerical strategy using the monokinetic closure

	Mesh	Tp	Tg	Main frequency	√	2RM S (mbar)
	(317 × 30)	First order	MacCormack	-		-
		MacCormack	MacCormack	-		-
		Multislope	MacCormack	2121 ± 16 Hz		0.2625
		Multislope	Multislope	-		-
	(634 × 60)	First order	MacCormack	-		-
		MacCormack	MacCormack	-		-
		Multislope	MacCormack	1463 ± 13 Hz		17.8
		Multislope	Multislope	2173 ± 13 Hz		6.21
	(1268 × 120)	Multislope	Multislope	1430 ± 35 Hz		5.57

May 22nd -27th 2016, Firenze, Italy ICMF-2016 -9th International Conference on Multiphase Flow May 22nd -27th 2016, Firenze, Italy

May 22nd -27th 2016, Firenze, Italy

The present research was conducted thanks to a Ph.D. Grant from ONERA. This work relies on a long-term collaboration work between EM2C laboratory and the Fundamental and Applied Energetics (DEFA) department at ONERA. Support from ONERA is gratefully acknowledged.


domain: chim phys

text:
Amara Chrayteh

Chris Ewels
email: chris.ewels@cnrs-imn.fr

Denis Jacquemin
email: denis.jacquemin@univ-nantes.fr

Dual fluorescence in strap ESIPT systems: a theoretical study †

Alkylamine-strapped chromophores based

Introduction

Excited state intramolecular proton transfer (ESIPT) is a fast photo-tautomerization process occurring in the electronically excited state (ES) of molecules possessing a strong intramolecular hydrogen bond (H-bond). 1,2 ESIPT produces a tautomer presenting an entirely different electronic structure from the initial Franck-Condon form. Consequently, the emission from the ESIPT structure is strongly redshifted as compared to the absorption, and a large Stokes shift is observed. Additionally, dual fluorescence originating from both the canonical and the ESIPT structures can also take place when the two tautomers have comparable ES energies. The emission of long-wavelength photons paves the way to many applications including white light-emitting materials, [3][4][5] polymorph-dependent light-emitting materials, 6 and biological fluorescent tags to study DNA-protein interaction, 7 which explains the growing popularity of ESIPT dyes. Most ESIPT processes take place within molecules containing a five or six-member aromatic ring presenting a strong intramolecular H-bond between the proton donor and acceptor groups. In the majority of the previously studied ESIPT molecules, an amino (NH 2 ) or an alcohol (OH) group acts as the proton donor, and a carbonyl group (-CQO) or an heterocyclic nitrogen atom serves as the proton acceptor. Therefore, the ESIPT reported in most compounds is based on a keto-enol or an amino-imino tautomerization (Fig. 1a). [8][9][10][11] The ground state (GS) of these chromophores typically exhibits the enol form, whereas the keto tautomer is formed through an intramolecular proton transfer along the H-bond between the two p-conjugated moieties. However, a few exceptions to such ESIPT structures have been reported, such as 2-(hydroxyphenyl)imidazo [1,2-a]pyridines (HPIPs), in which the ESIPT process yields a formal zwitterionic structure rather than a keto isomer (Fig. 1b). [12][13][14][15][16] A series of new ESIPT chromophores involving a tertiary amine moiety placed in a strap as the proton acceptor and a p-conjugated pyrrole scaffold as the proton donor was recently designed (Fig. 1c). 17 The first member of this new class of ESIPT dyes was reported by Suzuki et al. in 2014: they prepared a dye presenting a dithienylpyrrole skeleton capped by dimesitylboryl groups, 2 in Fig. 2. 17 This dye may show two emission bands originating from the local excited state (LE) and the zwitterionic (Z) tautomeric forms, the latter being generated by ESIPT. The measured fluorescence wavelengths are strongly dependent on solvent polarity. In cyclohexane, a non-polar solvent, 2 shows a single component emission (l fluo = 500 nm) attributed to the LE form, whereas increasing the polarity of the solvent by using acetone leads to a strong bathochromic shift (l fluo = 602 nm) of the emission, attributed to ESIPT (Z) fluorescence. In order to investigate this phenomenon theoretically, Suzuki et al. performed quantum-chemical calculations for 2 in both gas phase and acetone. 17 They employed a Time-Dependent DFT (TD-DFT) approach combined with both the CAM-B3LYP 18 and PBE0 19 exchange-correlation functionals, and accounted for the environmental effects using the Solvent Model Density (SMD) method. 20 A comparative analysis showed that the calculated emission wavelength for the ESIPT state was in good agreement with the experimental value, whereas, the computed l fluo of the LE isomer was redshifted as compared to the measurement. In 2018, the same group proposed a more extensive study evaluating the impact of the substituents at the terminal positions and the influence of the strap length on the optical properties. 21 This second investigation demonstrated that both the nature of the substituents and the strap length significantly affect the fluorescence spectra and the ESIPT efficiency. In more detail, they tested six compounds, with different substitutions at the terminal position, including compounds 1, 3, and 4, and one compound with a longer strap 5 (see Fig. 2). They found that the Z emission depends on the degree of the electron-withdrawing strength of the side group, as this strength affects the acidity of the pyrrole group in the ES. Likewise the chain length of the strap moiety is critical for the formation of the intramolecular H-bond, a too long strap preventing the ESIPT. A compact theoretical study of solvent effects on the LE and Z forms was also included in ref. 21 using the Polarizable Continuum Model (PCM). 22 It was found that the energies of both isomers decrease in a polar medium and, more importantly, that the Z tautomer becomes more stable in the ES when increasing the solvent polarity. Recently, Zhou and Zhao proposed a theoretical study further discussing the ESIPT reaction of 2 in three different solvents. 23 They analyzed the nature of LE and Z tautomers along with locating the ESIPT transition state using TD-DFT as well as the second-order Algebraic-Diagrammatic Construction approach [ADC(2)]. 23 They found that, in non-polar solvent such as cyclohexane, the energy barrier is high and thus, fluorescence originates mainly from the LE structure. In contrast, an increase in solvent polarity decreases the energy barrier for the ESIPT reaction and the fluorescence should mainly result from the Z form. To the best of our knowledge there are no other theoretical nor experimental studies on these systems.

The wide range of applicability of molecules exhibiting redshifted emission motivated us to use advanced theoretical tools to design novel ESIPT systems based on the strap scaffold. To begin with, we consider all previously reported compounds, 1-5, 17,21 displayed in Fig. 2 to study their fluorescence properties using a coherent theoretical model. We aim to: (i) determine if the selected protocol can reproduce the experimental outcomes; (ii) provide theoretical insight into their working mechanism; (iii) investigate unknown derivatives in an effort to design dyes with new spectroscopic signatures. For that reason, we choose the M06-2X 24 hybrid exchange-correlation functional in the present work, since it has been reported to provide accurate transition energies for a large panel of ES and is also suited for various classes of ESIPT dyes. [25][26][27][28][29][30][31] Thereafter, the effect of substituents is evaluated by modifying the terminal functional group and the heteroatom of the side rings in 1, resulting in the design of twelve new derivatives 6-17 (Fig. 3). In addition the solvatochromic effects are analyzed for these molecules by applying various solvent models including both linear response (LR) and state specific (SS) schemes to offer a higher accuracy in the modeling of the environment. The ESIPT mechanism is explored as well by locating the transition state (TS) connecting the LE and Z structures in the ES.

Computational details

Density Functional Theory (DFT) and Time-Dependent DFT (TD-DFT) calculations are carried out using the Gaussian16 package. [START_REF] Frisch | Gaussian16 Revision A.03[END_REF] The impact of the solvent, cyclohexane (e = 2.02) or acetone (e = 20.7), is accounted for using the integral equation formalism (IEF) version of the PCM approach. For all compounds, the GS and the ES geometries are fully optimized using DFT and TD-DFT respectively, combined with the M06-2X 24 hybrid exchange-correlation functional. The 6-31+G(d,p) atomic basis set is used during the geometry optimisations, whereas a larger atomic basis set, 6-311+G(2d,p), is selected to compute the total and transition energies on these optimal geometries using the same functional. Next, vibrational frequency calculations are performed for both states using the optimized geometries to verify that all structures are true minima of the potential energy surface, as well as to obtain the zero-point vibrational energies (ZPVE) for both states. Consistently, the ESIPT transition state geometries are optimized at the TD-DFT level using M06-2X functional and the 6-31+G(d,p) atomic basis set. These TS present a single imaginary frequency corresponding to the proton transfer process. Since we focus on excited-state calculations with PCM, and in order to obtain the best theoretical estimates of solvent effects, three different models are used in this work for computing the solvent impact on the energies: (i) the Linear Response approach [START_REF] Scalmani | [END_REF]34 which is a method well-suited for local bright ES that are characterized by a small change in the electronic density; such model is however less effective for charge-transfer excitations; (ii) the corrected LR (cLR) model that accounts for SS effects in a single-shot perturbative way; cLR is well suited when the transition involves a large density reorganization; 35 and (iii) the combination of these two formalisms together (LR + cLR) [36][37][38][39] which is more general and adequate for, e.g., states that are both bright and involving significant changes of the electronic density. Indeed, the LR + cLR model recovers all major solvation effects and often yields more accurate total and transition energies than the LR and cLR approaches taken separately. Since geometrical relaxation is a slow process, excited-state geometries have been determined in the equilibrium limit of the LR-PCM approach. Thus, vibrational frequencies, including zero point corrections and 0-0 energies, are considered in this limit as well. In contrast, both the absorption and emission energies are given in the nonequilibrium PCM limit, except for LR-PCM fluorescence energies, that are in equilibrium due to the intrinsic limit of the current implementations of this model. In addition, to ascertain for one case our TD-DFT results, we use the second-order Coupled Cluster model, CC2, with the aug-cc-pVDZ basis set. CC2 calculations are done with the Turbomole package [START_REF]TURBOMOLE V6.6 2014, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH[END_REF] applying the RI approach. These CC2 calculations are performed in gas phase and the solvent effects, as computed with TD-DFT, are next added to estimate the CC2 results in solution, following a previously described protocol. 29 We note that condensed phase CC calculations are also possible, [START_REF] Caricato | [END_REF]42 but such calculations have not been carried out here.

Results and discussion

Impact of the solvent model on fluorescence wavelengths

In this section, we perform a comparative analysis of theoretical and experimental fluorescence wavelengths. We compute the emission energies of molecules 1-5 displayed in Fig. 2 in two different solvents (cyclohexane and acetone). Note that compound 5 presents two conformers (5a and 5b, see below), and, we consider the most stable one in the statistics. To determine which solvent model could be viewed as the most appropriate, three PCM models are evaluated (LR, cLR, and LR + cLR). At this stage, it should be reminded that the vertical fluorescence energies are not strictly equivalent to measured wavelengths of maximal emission. 43 However, given the size of the systems, one has to resort to this vertical approximation which is reasonable in a homologous series of compounds. Three different criteria are taken into account to appraise the solvent models: (i) the difference between experimental and theoretical absolute fluorescence wavelengths computed for the two tautomers (LE and Z), (ii) the difference between fluorescence wavelengths originating from LE and Z in acetone, and (iii) the difference between LE emission wavelengths when going from one solvent to another. Our data are collected in Tables S1-S3 in the ESI. † For the first criterion, it can be seen that globally and, especially for the emission originating from the Z tautomer, LR + cLR gives the smallest deviations between experiment and theory. Table 1 provides a statistical analysis confirming this statement. Indeed, it can be noted that even if cLR gives the smallest MSE À0.056 eV/À0.054 eV and MAE 0.080 eV/0.098 eV for the LE fluorescence in cyclohexane/ acetone, LR + cLR provides tiny MSE (À0.020 eV) and MAE (0.050 eV) for the Z emission as well as the overall smallest MAE (0.120 eV). For the second criterion, we notice from Table S2 in the ESI † that the computed difference between LE and Z emission wavelengths are systematically too small regardless of the selected solvent model. It is nevertheless clear that LR + cLR gives slightly better correlation with experiment as compared to LR and cLR taken separately, e.g. a theoretical difference of À56 nm is found between the LE and Z fluorescence wavelengths for compound 2 in acetone with LR + cLR, whereas the experimental data provides a difference of À82 nm. Larger errors are obtained with both LR and cLR (À41 nm and À46 nm respectively). Similar findings are obtained for the other molecules. Let us now move to the third criterion, the magnitude of the solvatochromism for which we note that cLR yields slightly smaller values than experiment, whereas both LR and LR + cLR tend to overshoot the measured values (Table S3 in the ESI †). Therefore, these analyses indicate quite obviously that LR is not suitable for the investigated systems, whereas both cLR and LR + cLR appear to be more reasonable which is not a surprising outcome for ESIPT dyes involving both a large reorganization of the electron density between the GS and the ES, and two tautomers of different polarities (see below). Given that LR + cLR has globally the edge, we apply systematically this model in the following.

Influence of the solvent polarity

Our theoretical studies of compound 1 in cyclohexane and acetone demonstrate a strong dependency between the relative stabilization of LE and Z tautomers and the solvent polarity (Fig. 4). As expected, two local minima are obtained in the ES, one associated with LE and another with Z. A single component emission from LE is clearly predicted for 1 in cyclohexane as the LE form is more stable than Z by 0.1 eV, so there is no driving force for ESIPT. In contrast, the Z tautomer is significantly stabilized in acetone and consequently becomes more stable than LE. Provided ESIPT is faster than emission, which is likely given the fact that the experimentally measured fluorescence lifetimes are of the order of the ns, this potential energy surface indicates a possible ESIPT and, hence, a strong bathochromic shift of the emission when going from LE cyclohexane to Z in acetone (Dl fluo = 102 nm). To understand this large solvent effect, the ES dipole moment of the two tautomers has been determined in both solvents. The dipole moment of the Z tautomer is much larger than that of its canonical counterpart, and, more importantly, it significantly increases with the increase of solvent polarity: m ES = 10.6/13.3 D and m ES = 15.1/19.1 D for LE/Z in cyclohexane and acetone, respectively (Table S4 in the ESI †). This clearly supports the zwitterionic character of the tautomer after ESIPT and explain its very strong stabilization in acetone. For the records, we have also investigated the possibility to form the anion in the ES of 1, i.e., the 1 LE + acetone -1 À + acetoneÁH + . However, this reaction was found to be strongly endergonic (+47 kcal mol À1 ), and it can be likely be ruled out. This is also consistent with the experimental data of Suzuki and coworkers who reported no significant change of the emission even in water. 17 To further analyze the ESIPT process, let us investigate the nature of the ES in 1. Fig. 5 shows the frontier molecular orbitals together with the Electron Density Difference (EDD) plots computed in acetone. As expected, the first excited state presents a typical pp* topology, the changes of densities being localized on the conjugated moiety of the dye. There is a partial quadrupolar like charge-transfer from the central pyrrole towards the external aldehyde acceptors. Although the pyrrole is losing density after absorption, its nitrogen atom seems to undergo a slight p density gain in the lowest ES. Nevertheless this donating nitrogen atom becomes slightly less negative in the ES (q = À 0.02e), whereas, the accepting nitrogen atom is highly negative in both states (q = À 0.33e), hence it is ''willing'' to attract a positive charge, which is favorable for triggering ESIPT.

The theoretical TS for the proton transfer of 1 in cyclohexane was found to be more stable than the Z form by ca. 0.05 eV after inclusion of the ZPVE correction (there is a small maximum on the total energy potential energy surface). This suggests a barrierless ''back-ESIPT'' and hence the infeasibility of the emission from Z tautomer in this solvent. In sharp contrast, the forward process is barrierless in acetone, facilitating the formation of Z tautomer, and consistently only the emission from this structure can be expected to occur according to our calculations (Fig. 4). These findings are in excellent agreement with the experimental observations that show LE emission in cyclohexane but Z fluorescence in acetone. 21 In addition, due to the limits of TD-DFT, that has an error bar significantly larger than 0.05 eV, 26 CC2 calculations were performed for 1 on the geometries obtained from TD-DFT (see Table S5 in the ESI †). However, the CC2-corrected potential energy surface yields a slightly disfavored ESIPT formation in acetone. Therefore, CC2 corrections do not improve correlation with experiment, 21 and we therefore stick to TD-DFT in the following. Similar analyses can also be performed for molecules 2, 3, and 4 (see also Section 3.4). According to our calculations, the fluorescence of 2 is originating only from the canonical tautomer in cyclohexane. Indeed, as can be seen in Fig. 4, the LE form is more stable by ca. 0.25 eV compared to the ESIPT structure. In acetone, the Z state is significantly stabilized and consequently presents an energy comparable to its LE counterpart, indicating that proton transfer becomes a more accessible process. Nevertheless, we notice that for 2, the Z tautomer remains slightly less stable than the LE one, by 0.05 eV in acetone. This is not totally satisfying as Z emission is observed experimentally. This discrepancy is simply due to the limitations of TD-DFT or to the static picture used in this work. We note that a similar problem is also observed in the theoretical calculations performed by Suzuki et al. 17 Despite this small discrepancy, the global trend nicely fits experiment: 17 the theoretical emission wavelengths for compound 2 in both cyclohexane (l LE Fluo = 484 nm) and acetone (l LE Fluo = 523 nm/l ESIPT Fluo = 579 nm) are in reasonable agreement with the measured values in cyclohexane (l LE Fluo = 500 nm) and acetone (l LE Fluo = 520 nm/l ESIPT Fluo = 602 nm), 17 see Table S6 in the ESI. † Thus, our theoretical LE/Z wavelength difference correlates reasonably with experiment. Similar findings are obtained for compound 4, see Fig. S1 and Table S6 in the ESI. † Compound 3 exhibits experimentally a single emission band in cyclohexane at 417 nm originating from LE, which matches our theoretical calculations (l LE Fluo = 440 nm, see Fig. S1 in the ESI †). In acetone, two emission bands are obtained experimentally 21 at 420 nm and 490 nm, but it is not easy to attribute them to the LE or the Z tautomers. Indeed, according to our calculations, the two isomers should show a very small energy difference (ZPVE-corrected DE = À0.02 eV) within the error bar of theory see also (Fig. S1 in ESI †). In addition, regarding the emission wavelengths, the experimental LE/Z differences for 1, 2, and 4 in acetone are quite large, e.g., (Dl 1 = À78 nm, Dl 2 = À82 nm, and Dl 4 = À123 nm) in reasonable agreement with the TD-DFT results (Dl 1 = À32 nm, Dl 2 = À56 nm, and Dl 4 = À78 nm). In contrast, for 3 we compute a very small emission wavelength difference between these two tautomers (Dl 3 = À15 nm, see Table S6 in the ESI †) much smaller than the experimental value (À70 nm). This motivated us to further look at the experimental spectra (Fig. S2 in ref. 21) and, it appeared that the most intense emission peak appears at ca. 445 nm (apparently not around 417 nm) in cyclohexane which is also close to the acetone ''Z band'' that seems to appear at 460 nm (and apparently not around 490 nm) experimentally. Note that such Dl 3 of À15 nm seems to more reasonably fit our predictions. In general the experimental fluorescence spectrum of 3, as displayed in the original paper, bears little relation with those of the other dyes for which ESIPT is very clear. Consequently, we cannot make a definitive conclusion regarding the origin of the emission of 3, and the actual existence of ESIPT or not in that compound would likely deserve more experimental investigations.

Impact of the strap chain length on ESIPT

The effect of the strap chain length on the ESIPT behavior was studied experimentally by comparing compounds 1 and 5 (see Fig. 2). 21 In our calculations, two initial geometries have been considered for the optimization of 5, i.e., the crystallographic structure 5a and a chemically-intuitive structure 5b in which the two nitrogen atoms are closer (Fig. 6). The conformer corresponding to the X-ray like structure, 5a, is slightly more stable than the more compact 5b, by 3.74 kcal mol À1 and 4.50 kcal mol À1 in cyclohexane and acetone, respectively. Interestingly geometric and energy differences are also observed in the ES. For 5a, the LE state is more stable than the ESIPT form (in both solvents), indicating that the emission in 5a occurs from LE only, which fits the experimental data 21 (Table S6 in ESI †). The large distance between the proton donor and acceptor (4.22 Å) makes the identification of the ESIPT TS of 5a effectively impossible. Nevertheless, the stretching of the N-H bond by 0.40 Å results into a dramatic increase in the ES energy (+1.5 eV), hinting that the TS is even higher, confirming that ESIPT is impossible in such a conformation. The LE structure of 5b remains more stable than the Z tautomer in cyclohexane, whereas the Z structure of 5b becomes the most stable in acetone, similar to the case of compound 1, causing a bathochromic shift of the emission band. Additionally, the ESIPT barrier in 5b is quite small (ca. 0.09 eV) indicating that the transfer process is energetically accessible with such conformation. Therefore, theory shows that there can be ESIPT in this compact form, and that the absence of ESIPT experimentally is only due to a GS conformational effect. One can speculate that using another type of long strap, allowing the simultaneous presence of several GS conformers, could lead to ESIPT even if the two nitrogen atoms are distant in some conformations.

Influence of terminal substituents

Compounds 1, 2, and 4 have identical skeleton and strap chain length differ in their terminal substitutions. In order to rationalize these auxochromic effects on the ESIPT phenomenon, we compare their theoretical LE/Z relative energy difference (DE ES ) with their experimental relative emission peak intensities in Table 2. If the level of theory used would be perfectly accurate one would expect the energy difference DE ES = 0.00 eV to be the limit separating dominant LE and Z emissions. Indeed, a positive DE ES indicates that ESIPT is a process favored by the excited-state energy landscape and a strong Z emission is expected, whereas a significantly negative DE ES should correspond to cases where a dominant LE emission is obtained experimentally. As can be noted in Table 2, in cyclohexane the LE form is always more stable than the Z form by at least 0.1 eV for the three compounds. This is perfectly consistent with the experimental observations (100% LE emission). As explained above, 1 shows a clear change in the relative stabilities of Z and LE when going to acetone, hence, the zwitterion becomes favored (DE ES = 0.05 eV). Experimentally, the Z emission is strongly dominant (a relative peak intensity of 98% of Z) and only a trifling emission (2%) is detected for the LE form of 1. 21 For the two other compounds (2 and 4), we notice a strong stabilization of Z when going from cyclohexane to acetone. This stabilization decreases the energy difference between the two species and makes them very close (DE ES 2 = À0.05 eV and DE ES 4 = À0.03 eV). As already discussed, these two compounds exhibit experimentally a strong peak of Z and a weak LE band. Interestingly, the LE emission of 2 and 4 experimentally show higher relative intensities (7% and 8% respectively) than in 1. This indicates that the theoretical trends are indeed accurate, although the TD-DFT energies seem to slightly overstabilize LE as compared to Z. Consequently, according to the relative energy values obtained above and the experimental observations, it appears that a pure Z emission can be obtained for DE ES Z À0.05 eV.

Furthermore, we measured the distance between the two nitrogens (N-N) in the GS for 1, 2, and 4 (Table S7 in the ESI †). The distance between these two nitrogens atoms becomes shorter when going from cyclohexane to acetone, and more importantly, this distance is found to be significantly shorter in 1 than in 2 and 4. This is consistent with an easier ESIPT in 1.

Our theoretical distances are also very similar to the X-ray values. 17,21 In addition, the partial atomic charges of these two nitrogens in the ES are shown in Fig. S2 in the ESI. † The charges are very similar in cyclohexane and acetone. In contrast, the charge of the nitrogen of the strap moiety of 1 appears to be more negative than the one of 2 and 4 in both solvents, and thus, has a higher tendency to accept the proton. These findings are consistent with the fact that the proton transfer is easier in 1 than in 2 and 4, confirming that simple qualitative criteria can provide first insights into the ESIPT mechanism.

Let us now turn towards the emission wavelengths of the three compounds (Table 2). The experimental ranking of the redshifted emission originating from the Z tautomer is 4 4 2 4 1 in acetone. This experimental ranking is rather consistent with the TD-DFT one: 4 4 1 C 2. This ordering can be explained: the electron withdrawing BMes 2 groups in 4 are located in an extended p-conjugated moiety, which yields an improved delocalization of the LUMO (see Fig. S3 in the ESI †). Consequently, 4 exhibits a lower LUMO than that of 1 and 2 and, hence a more redshifted emission is observed. In contrast, 1 and 2, that contain formyl and BMes 2 groups respectively, show very small emission wavelength difference (Dl = 6 nm), indicating that both substituents exert similar effects according to our calculations, consistent with the topologies of the frontier MOs (Fig. S3, ESI †). In short, the theoretical emission wavelengths have the expected accuracy for such TD-DFT simulations, while the computed relative ES energies are found to be in slight favor of the LE form.

Design of novel derivatives substitution effects

A series of unknown compounds, 6-17, containing both electron-donating groups (EDG) and electron-withdrawing groups (EWG) of different strengths at the terminal positions (R) and various heteroatoms (X) have been designed (Fig. 3). In order to predict the experimental outcomes, we list in Table 3 the energy difference (E + ZPVE) between the two tautomers in both cyclohexane and acetone, along with the computed emission wavelengths for the emissive species. Dye 1 is chosen as the reference compound to understand the impact of substituents on ESIPT as well as to determine what kind of emission will likely appear in the spectrum.

Let us first evaluate the impact of the heteroatom on the relative stabilities of the LE and Z forms by replacing the sulfur atom (X = S) in 1 by X = O, X = Se, and X = NH in 6, 7, and 8 respectively. The two former appear to have very similar DE ES values in both solvents as compared to 1. One can therefore predict a LE emission in cyclohexane and a Z emission in acetone. Interestingly, theory predicts a significant redshift for the emission of the seleno derivative in acetone (611 nm), see Table 3. In order to understand the reason behind the similar ESIPT behavior despite the use of atoms with different sizes, we show the distance between the two nitrogen atoms (N-N) in 1, 6, and 7 and the charges of these centers in Fig. 7. No significant change is found in neither the N-N distance nor their charges, indicating that the atomic radius of the X atom has negligible effect on the stability of the LE and Z, hence going down in the chalcogen column should only tune the position of the bands without significantly affecting their relative intensities. Furthermore, for 8, the LE form is found to be largely stabilized in cyclohexane, and hence, a pure LE emission is expected to be observed. In contrast, in acetone the LE and Z tautomers of 8 show comparable energies (DE ES 8 = À0.03 eV), so that, according to the reference energetic window defined previously, the proton transfer in 8 should remain almost quantitative. Note that 8 displays a strongly blueshifted band as compared to 1.

Let us now investigate the effect of auxochromic groups (a set of both EDG and EWG of different strengths and natures is selected) at R position. In cyclohexane, only a pure LE emission is expected to be observed experimentally in almost all cases (DE ES o À0.05, see Table 3). The TS for the ESIPT process in these dyes have also been computed and the values are listed in Table S8 in the ESI. † The TS is never indicating a large barrier for the forward ESIPT in cyclohexane, so it should limit the speed of the process. However, there is a clear change in the energy difference between LE and Z in both 13 and 14 (DE ES 13 = À0.04 and DE ES 14 = À0.01 eV respectively) together with a strong redshift of the LE emission of these two compounds (585 nm and 514 nm respectively), hinting that dual emission is likely to occur even in a non-polar solvent for these two dyes. It is worth noting that 13 and 14 show the strongest EWG that allow to stabilize the negative charge appearing in the Z form. In acetone, the insertion of mesomeric EWG (12, 13, and 14) further stabilizes the Z tautomer making the emission from the latter dominant as in 1. Interestingly, a large bathochromic shift of the emission is predicted by theory for the two latter (749 nm and 687 nm for 13 and 14 respectively). The HOMO/ LUMO and EDD plots of 14 can be found in Fig. S4 (ESI †). It can be seen that the LUMO of 14 is very delocalized and that the density on the nitro groups increases upon photon absorption which yields a strong stabilization of the Z tautomer after ESIPT and hence a large redshift. Adding an inductive EWG (15-17) yields a smaller effect, the Z tautomer being energetically less favored (DE ES 15 = À0.02 eV, DE ES 16 = À0.06 eV, and DE ES 17 = À0.08 eV), and a unique LE emission in acetone is even predicted for 16 and 17, contrasting with 1. For 15 the DE ES is found to be within the reference energetic window for Z emission (DE ES 4 À0.05 eV) and although we predict the proton transfer to be almost quantitative, the emission wavelength predicted for 15 appears to be strongly blueshifted (472 nm) as compared to 1. Clearly there is a large amount of possibilities to separately tune the emission wavelengths and the ESIPT process when using EWG.

Let us now turn to EDG. Globally donor groups tend to stabilize the LE form, and consequently tend to yield a sole LE emission. The fact that EDG makes ESIPT more challenging was already reported before in classical enol-keto systems. 8 We underline that the impact of the EDG on the relative LE/Z energies is much larger than that of EWG, and follows the expected NMe 2 4 NH 2 4 Me ranking. When an EDG is inserted (9, 10, and 11), the LE tautomer is largely stabilized in acetone (DE ES 9 = À0.27 eV, DE ES 10 = À0.28 eV, DE ES 11 = À0.18 eV), and only a single emission originating from this tautomer is expected to be observed, confirming the intuitive prediction. Indeed, adding side donor groups should be detrimental for the formation of the negative charge. In Fig. S4 of the ESI, † we S8 in the ESI †).

Conclusions and outlook

A detailed theoretical study of the ESIPT process occurring in alkylamine-strapped chromophores has been performed with a focus on the impact of chemical substitutions. The ESIPT process in this recently proposed family of dyes yields a zwitterionic structure rather than a keto tautomer. We performed our ab initio calculations by combining TD-DFT method with a continuum solvent approach accounting for both linear-response and state-specific solvent effects (LR + cLR). This model was found more appropriate than the commonly used alternatives, although as all continuum models it does not include the structure of the solvent molecules and hence cannot foresee solvent-assisted proton transfer, which is likely reasonable as the experimental solvents are aprotic. A comparison between the experimental and theoretical values for four compounds 1-4, for which fluorescence spectra have been experimentally measured in both cyclohexane and acetone was performed. According to our calculations, the emission is originating solely from the LE tautomer when cyclohexane is used as a medium. Indeed the energy difference between the two tautomers indicates that ESIPT would yield a less stable tautomer in the ES. When going from cyclohexane to acetone a systematic stabilization of the Z tautomer compared to LE is obtained, making the emission from the latter (much) more accessible. These conclusions are in good agreement with the experimental observations of Suzuki et al., although the computed relative energies seem to indicate that LE is slightly ''overfavored'' by TD-DFT. Additionally, theory was able to deliver reasonably accurate fluorescence wavelengths for the four studied compounds. Several new insights emerge from this current study. First, compound 3 seems to be a less clear-cut case than previously thought. Indeed the experimental fluorescence spectrum as well as the theoretical calculations are likely not clear enough to definitively ascertain the origin of the emissions observed in the polar solvent, and hence the actual existence of ESIPT. This was noticed by the facts that: (i) experimentally, the emission attributed to LE in cyclohexane appears very close to the Z band in acetone; and (ii) very small energy and wavelength differences are computed between the two tautomers making it challenging to determined whether the emission is resulting from LE or Z. Second, by reexamining compound 5 that shows a longer strap, we could locate two different conformers, the most compact being less stable but likely to show ESIPT. This is an example in which exerting pressure on the system could help in developing ESIPT. Third, we proposed and designed 12 new substituted compounds through the introduction of various EWG and EDG at the terminal position (R) and different heteroatoms (X) (Fig. 3). Since compound 1 shows a clear ESIPT behavior, it was chosen as reference to evaluate the effect of these chemical modifications on both ESIPT and the emission wavelengths. We conclude that (i) replacing sulfur by larger heteroatoms affects only the position of the bands, but not the ESIPT process and can therefore be an interesting approach to ''simply'' tune the emission spectra; (ii) adding very strong EWG groups, e.g., NO 2 or (CN) 2 , should allow dual emission even in cyclohexane, thanks to a strong stabilization of the negative charge present in the Z tautomer; (iii) using EDG hampers the ESIPT process even in acetone; and (iv) both EDG and EWG strongly modify the emission wavelengths, e.g., we compute an emission at 749 nm for the Z form of the compound presenting (CN) 2 side groups.

In the future it would be very interesting to: (i) investigate how the ESIPT process would evolve in the solid-state, which can now be done by theoretical simulations; 44,45 (ii) explicitly account for the solvent molecules and their dynamics, as this may significantly impact the transfer process, 46 especially if these dyes are to be used in a protic environment in which both enhancement and quenching of the ESIPT process have been reported. [47][48][49]

Fig. 1

 1 Fig. 1 Representation of the ESIPT process in (a) a benzofuran dye based on the keto-enol tautomerism, (b) an HPIP derivatives that involves zwitterion form, and (c) the systems considered in this work.

Fig. 2

 2 Fig. 2 Molecular structure of the emissive strapped ESIPT compounds experimentally studied in ref. 17 and 21.

Fig. 3

 3 Fig. 3 Molecular structure of the novel derivatives studied here.

Fig. 4

 4 Fig. 4 Potential energy diagrams (E + ZPVE) for compound 1 and 2 obtained in (a) cyclohexane and (b) acetone using the LR + cLR solvent model. Theoretical energies and wavelengths have been obtained from TD-DFT calculations using the M06-2X functional and the 6-311+G(2d,p) atomic basis set. For 2, a very large compound, the determination of the TS was technically beyond reach.

Fig. 5

 5 Fig. 5 Representation of the Highest Occupied Molecular Orbitals (HOMO), Lowest Unoccupied Molecular Orbitals (LUMO) and the Electron Density Difference (EDD) plots (isovalue: 0.02 au for MOs 0.0015 au for EDD) of compound 1 in acetone. The red and blue areas in the EDD represent the increase and decrease of electron density upon photon absorption, respectively. We also provide the GS and ES Mulliken charges computed for the two nitrogen atoms.

Fig. 6

 6 Fig.6Representation of the two different optimized ground-state structures of 5 in cyclohexane: the X-ray-like structure (5a) and a new conformer obtained by molecular modeling (5b), both being optimized using the M06-2X functional and the 6-31+G(d,p) atomic basis set. Hydrogen atoms omitted for clarity.

  report the HOMO/LUMO and EDD plots of 10: the LUMO of 10 is less delocalized than the HOMO, and the density of NMe 2 decreases upon absorption which results in an increase of the basicity of the N-H of the pyrrole and tends to inhibit the proton transfer. The computed ESIPT barriers also support this interpretation (Table

Table 1

 1 Statistical data obtained from the comparison of experimental and theoretical emission energies (Table S1 in the ESI): mean signed error (MSE) and mean absolute error (MAE). All values are given in eV

	Solvent	State	Method	MSE	MAE
	Cyclohexane	LE	LR	0.038	0.102
			cLR	À0.056	0.080
			LR + cLR	0.050	0.098
	Acetone	LE	LR	0.140	0.164
			cLR	À0.054	0.098
			LR + cLR	0.198	0.198
		Z	L R	À0.107	0.107
			cLR	À0.262	0.262
			LR + cLR	À0.020	0.050
		All	LR		0.125
			cLR		0.138
			LR + cLR		0.120

Table 2

 2 Energy difference DE ES (E + ZPVE) between LE and Z tautomers for compounds 1, 2, and 4 in cyclohexane and in acetone together with the computed emission wavelengths for the emissive species in comparison with the experimental values. All the theoretical energies values have been obtained from TD-DFT using the M06-2X functional and the 6-311+G(2d,p) atomic basis set. The LR + cLR PCM model is applied. All the relative LE/Z energies are given in eV, whereas emission wavelengths are given in nm

	Cyclohexane

Table 3

 3 Energy difference DE ES (E + ZPVE) in eV between LE and Z tautomers for the series of novel 6-17 dyes in cyclohexane and acetone together with the computed emission wavelengths in nm for the emissive species. Theoretical energies values have been obtained from TD-DFT calculations using the M06-2X functional and the 6-311+G(2d,p) atomic basis set. The LR + cLR PCM model has been applied

		Cyclohexane		Acetone
		DE ES LR+cLR	l LE /l ESIPT	DE ES LR+cLR	l LE /l ESIPT
	1	À0.10	483/-	0.05	-/568
	6	À0.13	449/-	0.05	-/522
	7	À0.12	501/-	0.04	-/611
	8	À0.16	428/-	À0.03	-/492
	9	a	444/-	À0.27	476/-
	10	a	456/-	À0.28	488/-
	11	a	427/-	À0.18	455/-
	12	À0.10	456/-	0.03	-/525
	13	À0.04	585/697	0.12	-/749
	14	À0.01	514/630	0.14	-/687
	15	À0.14	427/-	À0.02	-/472
	16	À0.18	428/-	À0.06	459/-
	17	À0.19	432/-	À0.08	464/-

a Optimization of Z yields directly to the LE form.

Acknowledgements

The authors thank the RFI LUMOMAT for support in the framework of the Fluo-34 grant. This research used resources of the CCIPL (Centre de Calcul Intensif des Pays de Loire). Fig. 7 Distance between the two nitrogen atoms and partial Mulliken atomic charge of these atoms in Å for compounds 1 (X = S), 6 (X = O) and 7 (X = Se) in their LE forms in acetone. All values obtained for the ES.

Conflicts of interest

There are no conflicts to declare.