Jérémy Payan
email: jeremy-payan@etud.univ-tln.fr

Claude Jauffret
email: jauffret@univ-tln.fr

Annie-Claude Pérez
email: annie-claude.perez@univ-tln.fr

Maximization of the Fisher Information in PDA

In a cluttered environment, the probabilistic data association (PDA) model allows constructing efficient estimators. In this case, the Fisher information matrix (FIM) is equal to the FIM in the clean environment multiplied by the socalled information reduction factor. This factor depends implicitly on the detection and false alarms probabilities, hence on the threshold prior to the estimation step. The topic of this paper is to seek the optimal threshold, the one that maximizes the information reduction factor. Whereas the Poisson law is used as an approximation in PDA model, here we consider the binomial law for the large false alarms probabilities. An example, coming from signal processing, illustrates our analysis.

I. INTRODUCTION

Depending on the architecture of a surveillance system, target tracking takes place at three different levels:

• "Track-before-detect" (TBD) based on unthresholded measurements (for example in [START_REF] Kirubarajan | Passive ranging of a low observable ballistic missile in a gravitational field[END_REF]).

• "Track-after-detect": after the detection step, the available data are composed of false alarms and correct detections. The tracking can be made by the Probabilistic Data Association (PDA) algorithm. We are facing a cluttered environment.

• "Track-after-detect" and extraction of the line of measurements: the input of the tracking function is a line of supposedly 'good' measurements, which is an optimistic and unrealistic situation. Such tracking is said to be in a clean environment. The topic of this paper takes place in the second scheme. We are concerned with the determination of the 'optimal' threshold in the detection step. Indeed, a poor choice of the threshold can lead to too high a number of false alarms or a very low level of detection probability. In both cases, we intuitively understand that no information may be drawn from these situations. We deduce that the information might be maximal at a particular threshold. Fortunately, we have a tool to quantify the information contained in the available measurements: the Fisher Information Matrix (FIM).

In the PDA model, the FIM in the presence of clutter equals the FIM in a clean environment multiplied by a scalar factor [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF]. The so-called q 2 factor lies between 0 and 1, and explains the information loss due to a certain amount of false alarms and non-unity detection probability [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Bar-Shalom | Tracking and Data Association[END_REF]. This factor is a function of the detection probability and the amount of false alarms in a restricted space. With all these considerations, the existence of an optimal detection threshold that would maximize q 2 is a fundamental question. This paper deals with the search for an optimal threshold that maximizes q 2 . Optimal thresholds have been investigated in [START_REF] Fortmann | Detection thresholds for tracking in clutter-A connection between estimation and signal processing[END_REF], [START_REF] Li | Stability evaluation and track life of the PDAF for tracking in clutter[END_REF], with a specific criterion at each time (target position error in [START_REF] Fortmann | Detection thresholds for tracking in clutter-A connection between estimation and signal processing[END_REF], track lifetime and track loss in [START_REF] Li | Stability evaluation and track life of the PDAF for tracking in clutter[END_REF]).

Our paper is composed of four main sections. Section II presents briefly the PDA with its uses and necessary statistical assumptions, as well as the computation of the FIM in PDA. Section III shows how to optimize q 2 , by changing the false alarms repartition model. After detailing a simple example to solve our problem, we present in Section IV some numerical results. Then, the conclusion, appendix and references follow.

II. THE FISHER INFORMATION MATRIX IN PROBABILISTIC DATA ASSOCIATION

A. Introduction to Probabilistic Data Association

Probabilistic Data Association (PDA) is well described in [START_REF] Bar-Shalom | Tracking and Data Association[END_REF]. PDA is used in realistic tracking contexts. Indeed, it takes into account the presence of false alarms (which are the element of the clutter), as well as a non-unity detection probability. The main aim is to estimate the parameters of interest when true detections are drowned in false alarms, and don't surely appear. Historically, two estimation techniques have been studied:

• A recursive technique : Probabilistic Data Association Filter (PDAF) allowing real-time estimation, based upon a modified Kalman filter. • A batch technique : PDA Likelihood Maximization (PDA-MLE). In our study, the latter is used. Some statistical assumptions are necessary to establish our model, following [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Bar-Shalom | Tracking and Data Association[END_REF] :

• The detections due to the target are corrupted by an additional zero-mean white Gaussian noise. The power of the noise is σ 2 . Subsequently, these detections will be called "true detections".

= (z 1,k , • • • , z m k ,k ) T •
The vectors z k are independent conditionally to X for k = 1 : K. Applying the total probability theorem, we end up with the likelihood 1 of X given the measurement vector z k :

L(X|z k ) = 1 -P d u m k µ f a (m k ) + P d u m k -1 µ f a (m k -1) m k × m k i=1 1 √ 2πσ exp - 1 2 z i,k -s k (X) σ 2 (1)
with s k (X) the state model at time k. Thanks to the independence assumption, we get

L(X|z) = K k=1 L(X|z k ).

B. The Fisher Information Matrix

We recall the well-known expression for the FIM in a clean environment, called J:

J = K k=1 1 σ 2 ∇ X s k (X)∇ T X s k (X) (2)
In the presence of clutter, the computation of the FIM is more complicated, but can be done with the above assumptions. The technique consists in restricting the measurement space to a gate around the true detection [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF]- [START_REF] Fortmann | Detection thresholds for tracking in clutter-A connection between estimation and signal processing[END_REF]. Its size, denoted by v g , depends on the standard deviation of the measurements:

v g = 2gσ (3)
Since the true detections have a Gaussian distribution, the most informative measurements are located in this gate when g = 5.

Indeed, more than 99 % of a Gaussian population are within five standard deviations of the mean. Conversely, measurements outside the gate contain no additional information. As a consequence, the likelihood function can be approximated by the one computed only with measurements in the gate. Doing this, the FIM in a cluttered environment is proven to be equal to the FIM in clean environment multiplied by a factor q 2 :

F = q 2 J (4)
When µ f a is a Poisson law (the usual assumption), q 2 is expressed with respect to parameters λv g and detection probability P d : [START_REF] Jauffret | Trajectographie passive, observabilité et prise en compte des fausses alarmes[END_REF].

q 2 = q 2 (P d , λv g ) [1]-[3],

III. OPTIMIZATION OF THE INFORMATION REDUCTION FACTOR

A. Poisson law versus binomial law

The use of the Poisson law in (1) is justified by the fact that the threshold of the prior test is chosen for a low probability of false alarms P f a [START_REF] Sittler | An optimal data association problem in surveillance theory[END_REF]. Indeed, if M is the number of 1 We implicitly use the notation proposed in [START_REF] Mendel | Lessons in estimation theory for signal processing, communications, and control[END_REF]. measurement cells on which the test is carried out, we have λu M P f a , provided P f a is small enough. Consequently, λv g N g P f a , where N g is the number of measurement cells in the gate. Hence, q 2 is a function of the couple (P d ; P f a ). Finally, using the ROC curves of the detector, we end up with a new expression for q 2 , depending only on the threshold t: q 2 (P d , λv g ) = q 2 (t)

The search for the maximum q 2 is nothing else than the search for the optimal t. But, the optimal t could define a probability of false alarm which could be large. In this case, the use of the Poisson law is no longer justified. This is why we have to introduce the 'natural' false alarm probability into (1): the binomial law with N g and P f a as parameters. Generally speaking, N g depends on an analysis of the methods of measurement, and the way the data are processed and quantified.

B. Information Reduction Factor Derived from a Binomial Law

The reduction information factor q 2 may then be derived from the binomial distribution:

q 2 (t) = Ng m=1 2P d (t)B(m -1) √ 2πg m-1 × g 0 • • • g 0 exp -ξ 2 1 ξ 2 1 a(t) βm(t) + m i=1 exp -1 2 ξ 2 i dξ 1 • • •dξ m (5)
q 2 is then a sum of m-fold integrals, with B(m -1) being the binomial distribution with parameters N g and P f a , t the detection threshold, and

a(t) = 1 -P d (t) (6)
β m (t) = 2gP d (t) [1 -P f a (t)] √ 2πP f a (t)(N g -m + 1) (7)
g = 5 (8)
The derivation of the above is detailed in the Appendix, and based on [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Jauffret | Trajectographie passive, observabilité et prise en compte des fausses alarmes[END_REF].

We note that under the conditions of the convergence of the binomial towards the Poisson distribution, q 2 expressed with both distributions has the same values. Fig. 1 gives us an example of the two q 2 , coming from signal processing and detailed hereafter. So, we have generalized the information reduction factor concept, by extending it to cases with high P f a .

IV. APPLICATION TO SIGNAL PROCESSING

A. Problem Formulation

The context we describe serves as an example for optimal threshold research. The framework is the following: we have a sine wave signal, with unknown amplitude, frequency and phase [START_REF] Rife | Single tone parameter estimation from discrete-time observations[END_REF]. This is our signal of interest. This sine wave is corrupted by an additive zero-mean white Gaussian noise. There are two detection hypotheses:

• H 0 "there is only noise" • H 1 "there are signal of interest and noise" Ref. [START_REF] Kay | Fundamentals of Statistical Signal Processing: Detection Theory[END_REF] shows that the Neyman-Pearson optimal test is equivalent to comparing the periodogram of the signal to a threshold. The probability distribution of the periodogram under each previous hypothesis is:

• H 0 : central χ 2 2 distribution • H 1 : noncentral χ 2 2 distribution
The noncentrality parameter γ depends on the SNR and on the length of the samples as follows [START_REF] Kay | Fundamentals of Statistical Signal Processing: Detection Theory[END_REF]:

γ = N A 2 2σ 2 (9)
with N the number of samples, A the amplitude of the sine wave, and σ 2 the power of the noise. So, if we call the Neyman-Pearson test Λ and the threshold t, we have by definition

P d (t) P (Λ > t| H 1 ) (10)
P f a (t) P (Λ > t| H 0 ) (11)
In this application, the temporal signal is segmented into K non-overlapped scans. The periodogram of each scan is compared to the threshold. The frequency cells of scan k greater than the threshold compose the measurement vector z k . By the knowledge of probability density function (pdf), we deduce easily that

P f a (t) = exp - t 2 (12)
The expression of P d is a numerical series detailed in [START_REF] Kay | Fundamentals of Statistical Signal Processing: Detection Theory[END_REF], in which an algorithm numerically computing P d is presented. So, we are able to express both P d and P f a with respect to the threshold. In this context, N g is the number of frequency cells in the gate.

B. Numerical Results

With the previous knowledge of the pdf under both hypotheses, and the new expression for q 2 , we may henceforth compute q 2 with respect to the detection threshold. We compute q 2 numerically using Monte-Carlo runs [START_REF] Robert | L'analyse statistique bayésienne[END_REF] with 500 000 simulations, and N g set to 5. The search for the maximum is based on a grid-search method. The optimal threshold in the sense of q 2 may then be placed on a ROC curve, to find the optimal operating point, for different values of γ.

Fig. 1 shows the importance of choosing the binomial law instead of the Poisson law to search for the maximum q 2 . The Poisson law may indeed prevent the finding of the maximum in certain cases, as shown by the considerations from previous section. Most of the time, the search for the maximum would be possible under the Poisson distribution if the optimal threshold is not too low. Fig. 2 shows q 2 with respect to the threshold for different values of the noncentrality parameter γ, as well as the optimal threshold in each case. The higher the noncentrality parameter, the higher the optimal threshold. This underscores the differences between the pdf under both hypotheses. Fig. 3 represents the ROC curves for each previous case and the operating point from the optimal threshold, highlighting the powerful performance, even in the case of relatively low γ.

V. CONCLUSION

In this paper, we have posed and solved the problem of choosing the optimal detection threshold in a surveillance system. The criterion is the amount of information relative to the parameter of interest and contained in the retained binary measurements. Based upon the assumptions of the PDA, the optimal threshold is the one for which the information reduction factor is the largest. An example of signal processing illustrates our approach. In the future, a similar problem will be studied when the measurements are no longer binary, but coupled with their own energy as in [START_REF] Kirubarajan | Passive ranging of a low observable ballistic missile in a gravitational field[END_REF].

APPENDIX DERIVATION OF THE INFORMATION REDUCTION FACTOR WITH BINOMIAL LAW

We derive q 2 with a binomial law, using the same principles as explained in [START_REF] Jauffret | Track formation with bearing and frequency measurements in clutter[END_REF], [START_REF] Jauffret | Trajectographie passive, observabilité et prise en compte des fausses alarmes[END_REF]. First, we express the likelihood at time k with (1) restricted to the validation region gate of volume v g . Moreover, µ f a is replaced with the binomial distribution with parameters N g and P f a . For the sake of clarity, the subscript k is dropped:

L(X|z) = 1 -P d (t) v m g B(m) + P d (t) v m-1 g B(m -1) m × m i=1 1 √ 2πσ exp - 1 2 z i -s(X) σ 2 (13)
where v g = 2gσ is the volume of the gate. After factorization and highlighting that

B(m -1) B(m) = m(1 -P f a (t)) P f a (t)(N g -m + 1) (14)
we get

L(X|z) = B(m) v m g 1 -P d (t) + P d (t)v g (1 -P f a (t)) P f a (t)(N g -m + 1) × m i=1 1 √ 2πσ exp - 1 2 z i -s(X) σ 2 . (15)
Then, it is possible to compute FIM in clutter from (15). The latter is then the FIM in a clean environment multiplied by the scalar factor q 2 , which is now

q 2 (t) = Ng m=1 σ 2 1 2πσ 2 P d (t)v g (1 -P f a (t)) P f a (t)(N g -m + 1) 2 1 σ 2 B(m) v m g × vg • • • vg 1 Φ(X|z) m i=1 m j=1 exp - 1 2 z i -s(X) σ 2 - 1 2 z j -s(X) σ 2 z i -s(X) σ z j -s(X) σ dz 1 • • •dz m (16)
with

Φ(X|z) = 1 -P d (t) + P d (t)v g (1 -P f a (t)) P f a (t)(N g -m + 1) × m i=1 1 √ 2πσ exp - 1 2 z i -s(X) σ 2 . (17)
This way, q 2 is a finite sum of m-fold integrals. We define

ξ i = z i -s(X) σ (18)
Reintroducing in (16) the notations of a(t) and β m (t) given in ( 6) and ( 7), we get

q 2 (t) = Ng m=1 β 2 m (t) B(m) v m g σ m × g -g • • • g -g 1 a(t) + β m (t) m i=1 exp -1 2 ξ 2 i × m i=1 m j=1 exp - 1 2 ξ 2 i - 1 2 ξ 2 j ξ i ξ j dξ 1 • • • dξ m (19)
For ξ j fixed, q 2 (t) is an odd function with respect to ξ i . Then, the cross-terms vanish from the integral because of the symmetry of the integration domain: We use the parity of the function with respect to ξ i and note that the integrals are the same for each considered index. After some manipulations, we obtain q 2 (t) = we finally get the information reduction factor expressed with binomial distribution given in [START_REF] Li | Stability evaluation and track life of the PDAF for tracking in clutter[END_REF].

q 2 (t) =

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. q 2 with both Poisson and binomial laws, γ = 2

  True detections appear at most one time at each sampling time, with a detection probability P d .• The available measurement vector at scan k is z k

	• The clutter is distributed according to an uniform law in
	a finite scan space u. The clutter is composed of the false
	alarms.
	• The amount of false alarms at a specific sample k follows

a discrete probability law µ f a . Usually, µ f a = P, the Poisson law of parameter λu (expected number of false alarms).

•


D Passiatore
email: donatella.passiatore@poliba.it

L Sciacovelli

P Cinnella

G Pascazio

Thermochemical non-equilibrium effects in turbulent hypersonic boundary layers

Keywords: compressible boundary layers, hypersonic flow, turbulent reacting flows

A hypersonic, spatially evolving turbulent boundary layer at Mach 12.48 with a cooled wall is analysed by means of direct numerical simulations. At the selected conditions, massive kinetic-to-internal energy conversion triggers thermal and chemical non-equilibrium phenomena. Air is assumed to behave as a five-species reacting mixture, and a two-temperature model is adopted to account for vibrational non-equilibrium. Wall cooling partly counteracts the effects of friction heating, and the temperature rise in the boundary layer excites vibrational energy modes while inducing mild chemical dissociation of oxygen. Vibrational non-equilibrium is mostly driven by molecular nitrogen, characterized by slower relaxation rates than the other molecules in the mixture. The results reveal that thermal non-equilibrium is sustained by turbulent mixing: sweep and ejection events efficiently redistribute the gas, contributing to the generation of a vibrationally under-excited state close to the wall, and an over-excited state in the outer region of the boundary layer. The tight coupling between turbulence and thermal effects is quantified by defining an interaction indicator. A modelling strategy for the vibrational energy turbulent flux is proposed, based on the definition of a vibrational turbulent Prandtl number. The validity of the strong Reynolds analogy under thermal non-equilibrium is also evaluated. Strong compressibility effects promote the translational-vibrational energy exchange, but no preferential correlation was detected between expansions/compressions and vibrational over-/under-excitation, as opposed to what has been observed for unconfined turbulent configurations.

Introduction

High-speed flows represent a challenging topic of interest for manifold configurations, including objects entering a planetary atmosphere or for atmospheric hypersonic flight [START_REF] Gnoffo | Computational aerothermodynamic design issues for hypersonic vehicles[END_REF]). In such flows, conversion of massive amounts of kinetic energy into internal energy causes a sudden rise of the flow temperature. The effects triggered by the high temperatures include chemical reactions and vibrational relaxation phenomena on characteristic time scales comparable to the flow time scales. The complex thermochemical state induced by such conditions may affect quantities of interest for the design of high-speed vehicles significantly [START_REF] Bertin | Critical hypersonic aerothermodynamic phenomena[END_REF][START_REF] Candler | Rate effects in hypersonic flows[END_REF]). In the past, thermochemical effects caused by hypersonic conditions have been investigated in more or less simple configurations. Most of the efforts focused on stagnation-point flows, of interest for bluff bodies re-entering the atmosphere [START_REF] Fay | Theory of stagnation point heat transfer in dissociated air[END_REF][START_REF] Armenise | Nonequilibrium vibrational kinetics in the boundary layer of re-entering bodies[END_REF][START_REF] Bonelli | An MPI-CUDA approach for hypersonic flows with detailed state-to-state air kinetics using a GPU cluster[END_REF][START_REF] Colonna | Impact of fundamental molecular kinetics on macroscopic properties of high-enthalpy flows: the case of hypersonic atmospheric entry[END_REF]. In later stages of re-entry trajectories or for suborbital flight conditions characterized by higher free-stream densities, the flow description is further complicated by the transition from a laminar to a turbulent regime. Such a transition may be triggered by free-stream disturbances, erosion and ablation of the thermal protection systems (TPS) or due to surface defects, leading to sharp rise of the skin friction and heat fluxes at the walls. For this reasons, laminar-to-turbulent transition has been identified as a subject of major concern for the accurate prediction of the aerothermodynamic flow fields around objects flying at hypersonic speeds. A significant number of contributions in the literature have addressed the linear and nonlinear stability of flat-plate boundary layers in chemical non-equilibrium, sometimes up to the initial stages of transition. The pioneering studies of [START_REF] Malik | Real gas effects on hypersonic boundary-layer stability[END_REF], [START_REF] Hudson | Linear stability of hypersonic flow in thermochemical nonequilibrium[END_REF] and [START_REF] Perraud | Studies of laminar-turbulent transition in hypersonic boundary layers at ONERA[END_REF] pointed out that chemical reactions tend to promote the so-called second mode instability, similarly to strongly cooled walls. More recently, [START_REF] Marxen | A method for the direct numerical simulation of hypersonic boundary-layer instability with finite-rate chemistry[END_REF] and [START_REF] Marxen | Direct numerical simulations of hypersonic boundary-layer transition with finite-rate chemistry[END_REF] used direct numerical simulations (DNS) to evaluate the maximum streamwise velocity disturbance caused by large-or small-amplitude waves in presence of chemical reactions. [START_REF] Kline | Multiple boundary layer instability modes with nonequilibrium and wall temperature effects using lastrac[END_REF] and [START_REF] Knisely | Significant supersonic modes and the wall temperature effect in hypersonic boundary layers[END_REF] investigated the effects of thermochemical non-equilibrium (TCNE) on boundary-layer stability; the effect of ablation has also been considered [START_REF] Mortensen | Real-gas and surface-ablation effects on hypersonic boundary-layer instability over a blunt cone[END_REF][START_REF] Miró Miró | Decoupling ablation effects on boundary-layer stability and transition[END_REF]. On the other hand, while wall-bounded turbulence under low-enthalpy conditions (i.e. with air behaving as a calorically perfect gas) has been carefully scrutinized over the years, the interplay between TCNE conditions and turbulence has received much less attention. Compressible turbulent boundary-layer (TBL) configurations have been investigated with different turbulence injection techniques and wall treatments; for a thorough overview of DNS of TBLs at relatively low free-stream Mach numbers (M ∞ ≤ 3), the reader may refer to the work of [START_REF] Wenzel | DNS of compressible turbulent boundary layers and assessment of data/scaling-law quality[END_REF] and references therein. Many efforts have been devoted to extending the validity of scaling laws derived for incompressible configurations [START_REF] Morkovin | Effects of compressibility on turbulent flows[END_REF]) to compressible ones. Inspection of the hypersonic regime by means of high-fidelity simulations was initiated by the studies of [START_REF] Duan | Direct numerical simulation of hypersonic turbulent boundary layers. Part 3. Effect of Mach number[END_REF], who performed DNS of temporally evolving TBLs with nominal free-stream Mach numbers up to 12 to prove the robustness of Morkovin's hypothesis even in high-speed flows. A DNS of a transitional spatially evolving boundary layer at Mach 6 was performed by [START_REF] Franko | Breakdown mechanisms and heat transfer overshoot in hypersonic zero pressure gradient boundary layers[END_REF], whereas [START_REF] Zhang | Direct numerical simulation database for supersonic and hypersonic turbulent boundary layers[END_REF] produced a database of cooled boundary layers up to M ∞ = 14, with flow conditions representative of those encountered in hypervelocity wind tunnels. [START_REF] Hirschfelder | Simulation and modeling of cold-wall hypersonic turbulent boundary layers on flat plate[END_REF] performed DNSs of zero-pressure-gradient hypersonic TBLs and conducted a priori and a posteriori analyses of the performances of different Reynolds-averaged Navier-Stokes (RANS) models. A handful of works focused on high-enthalpy flows, at free-stream conditions for which dissociation and recombination reactions are triggered [START_REF] Martin | Temperature fluctuation scaling in reacting boundary layers[END_REF][START_REF] Duan | Effect of finite-rate chemical reactions on turbulence in hypersonic turbulence boundary layers[END_REF]. A comparative study between low-and high-enthalpy, reactive, temporally evolving boundary layers was performed by Duan & Martín (2011b); the authors found that the two closely resemble each other, since the scaling laws validated under low-enthalpy conditions still hold or could be generalized for high-enthalpy applications. The interaction of finite-rate chemistry with transitional and turbulent spatially developing boundary layers was recently investigated by [START_REF] Passiatore | Finite-rate chemistry effects in turbulent hypersonic boundary layers: a direct numerical simulation study[END_REF] and Di Renzo & Urzay (2021), who isolated the effects of chemical non-equilibrium on pseudo-adiabatic and cooled TBLs, respectively. In both studies, the effect of thermal non-equilibrium was neglected since the selected edge conditions were not expected to promote significant vibrational relaxation phenomena. The influence of these processes on turbulent flow behaviour has been investigated only for unconfined flow configurations so far, such as isotropic turbulence decay [START_REF] Neville | Effect of thermal non-equilibrium on decaying isotropic turbulence[END_REF][START_REF] Khurshid | Decaying compressible turbulence with thermal non-equilibrium[END_REF][START_REF] Zheng | Vibrational relaxation in compressible isotropic turbulence with thermal nonequilibrium[END_REF]) and mixing layers [START_REF] Neville | Thermal non-equilibrium effects in turbulent compressible shear flows[END_REF][START_REF] Fiévet | Numerical investigation of the coupling of vibrational nonequilibrium and turbulent mixing using state-specific description[END_REF].

The present work aims at bridging this knowledge gap with the investigation of a TBL exposed to both thermal and chemical non-equilibrium conditions. For that purpose, we perform a DNS of a flat-plate boundary layer, subjected to post-shock conditions of a 6 • sharp wedge flying at Mach 20. The edge thermodynamic state is such that vibrational relaxation phenomena are significant within the boundary layer, while a mild chemical activity is developed.

The paper is organized as follows. The governing equations, the numerical method, the outline of the simulation and the main parameters are described in § 2. Numerical results are presented in § 3; focusing mostly on the fully turbulent regime, we illustrate the behaviour of the dynamic field and compare different transformations for the streamwise velocity. Afterwards, thermal non-equilibrium effects are inspected by highlighting the tight coupling with turbulent mixing mechanisms. Correlations deriving from the strong Reynolds analogy are also evaluated, as well as classical closures of the new terms arising in the RANS equations. Lastly, the coupling between compressibility effects and thermal non-equilibrium is investigated. Conclusions are then drawn in § 4.

Methodology

Governing equations and thermodynamic models

The fluid under investigation in the current study is air at high temperature, modelled as a five-species mixture of N 2 , O 2 , NO, O and N. When considering a gas under thermal non-equilibrium conditions, the vibrational energetic levels of the molecules in the mixture are partially excited and no longer equilibrated with the roto-translational ones. A direct consequence is that, even assuming that particle populations follow a Boltzmann distribution, the utilization of a single static temperature to represent all the energetic modes is no longer valid [START_REF] Anderson | Hypersonic and High Temperature Gas Dynamics[END_REF]. A classical approach to deal with such conditions, referred to as multi-temperature in the literature, consists of taking into account the vibrational levels separately, by means of additional 'vibrational' temperatures for each molecule. In order to keep a reasonable number of equations to be integrated, the two-temperature model of [START_REF] Park | Two-temperature interpretation of dissociation rate data for N 2 and O 2[END_REF] is adopted for the present calculations. Such a model, widely used in previous works (see, e.g. [START_REF] Hudson | Linear stability of hypersonic flow in thermochemical nonequilibrium[END_REF][START_REF] Franko | Effects of chemistry modeling on hypersonic boundary layer linear stability prediction[END_REF][START_REF] Bitter | Stability of highly cooled hypervelocity boundary layers[END_REF], assumes that the vibrational energy states of each molecule satisfy a Boltzmann distribution characterized by only one vibrational temperature T V , common to all diatomic species in the mixture (that is, N 2 , O 2 and NO).

A single additional conservation equation is, therefore, needed for the total vibrational energy e V . Thus, the behaviour of such flows is governed by the compressible Navier-Stokes equations for multicomponent chemically reacting and thermally relaxing gases, which read ∂ρ ∂t + ∂ρu j ∂x j = 0 (2.1)

∂ρu i ∂t + ∂ ρu i u j + pδ ij ∂x j = ∂τ ij ∂x j (2.2) ∂ρE ∂t + ∂ (ρE + p) u j ∂x j = ∂(u i τ ij ) ∂x j - ∂(q TR j + q V j ) ∂x j - ∂ ∂x j NS n=1 ρ n u D nj h n (2.3) ∂ρ n ∂t + ∂ ρ n u j ∂x j = - ∂ρ n u D nj ∂x j + ωn (n = 1, .., NS -1) (2.4) ∂ρe V ∂t + ∂ρe V u j ∂x j = ∂ ∂x j -q V j - NM m=1 ρ m u D mj e Vm + NM m=1 (Q TVm + ωm e Vm ) .
(2.5)

In the preceding formulation, ρ is the mixture density, t the time coordinate, x j the space coordinate in the jth direction of a Cartesian coordinate system, with u j the velocity vector component in the same direction, p is the pressure, δ ij the Kronecker symbol and τ ij the viscous stress tensor, modelled as

τ ij = μ ∂u i ∂x j + ∂u j ∂x i - 2 3 μ ∂u k ∂x k δ ij , (2.6)
with μ the mixture dynamic viscosity. In (2.3), E = e + 1 2 u i u i is the specific total energy (with e the mixture internal energy), q TR j and q V j the roto-translational and vibrational contributions to the heat flux, respectively; u D nj denotes the diffusion velocity and h n the specific enthalpy for the nth species. In the species conservation equations (2.4), ρ n = ρY n represents the nth species partial density (Y n being the mass fraction) and ωn the rate of production of the nth species. The sum of the partial densities is equal to the mixture density ρ = NS n=1 ρ n , NS being the total number of species. To ensure total mass conservation, the mixture density and NS -1 species conservation equations are solved, while the partial density of the NSth species is computed as ρ NS = ρ -NS-1 n=1 ρ n . In the following, we set such species as molecular nitrogen, since it is the most abundant one throughout the computational domain. In (2.5), e V represents the mixture vibrational energy, given by

e V = NM m=1 Y m e Vm ,
(2.7)

with e Vm the vibrational energy of the mth molecule and NM their total number. In the same equation, Q TV = NM m=1 Q TVm represents the energy exchange between vibrational and translational modes (due to molecular collisions and linked to energy relaxation phenomena) and NM m=1 ωm e Vm the vibrational energy lost or gained due to molecular depletion or production.

Each species is assumed to behave as a thermally perfect gas; Dalton's pressure mixing law leads then to the thermal equation of state

p = ρT NS n=1 RY n M n = T NS n=1 ρ n R n , (2.8)
with R n and M n the gas constant and molecular weight of the nth species, respectively, and R = 8.314 J mol -1 K -1 the universal gas constant. The thermodynamic properties of high-T air species are computed considering the contributions of translational, rotational and vibrational modes; specifically, the internal energy reads

e = NS n=1 Y n h n - p ρ , with h n = h 0 f ,n + T T ref (c T p,n + c R p,n ) dT + e Vn .
(2.9)

Here, h 0 f ,n is the nth species enthalpy of formation at the reference temperature (T ref = 298.15 K), c T p,n and c R p,n the translational and rotational contributions to the isobaric heat capacity of the nth species, computed as

c T p,n = 5 2 R n and c R p,n = R n for diatomic species 0 for monoatomic species (2.10a,b)
and e Vn the vibrational energy of species n, given by

e Vn = ⎧ ⎨ ⎩ θ n R n exp (θ n /T V ) -1
for diatomic species 0 for monoatomic species

(2.11) with θ n the characteristic vibrational temperature of each molecule (3393 K, 2273 K and 2739 K for N 2 , O 2 and NO, respectively). After the numerical integration of the conservation equations, the roto-translational temperature T is computed from the specific internal energy (devoid of the vibrational contribution) directly, whereas an iterative Newton-Raphson method is used to compute T V from (2.7).

The heat fluxes are modelled by means of Fourier's law, q TR j = -λ TR (∂T/∂x j ) and q V j = -λ V (∂T V /∂x j ), λ TR and λ V being the roto-translational and vibrational thermal conductivities, respectively. Pure species' viscosity and thermal conductivities are computed using curve fits by [START_REF] Blottner | Chemically reacting viscous flow program for multi-component gas mixtures[END_REF] and Eucken's relations, respectively [START_REF] Hirschfelder | Simulation and modeling of cold-wall hypersonic turbulent boundary layers on flat plate[END_REF]. The corresponding mixture properties are evaluated by means of Wilke's mixing rules [START_REF] Wilke | A viscosity equation for gas mixtures[END_REF]. Mass diffusion is modelled by means of Fick's law

ρ n u D nj = -ρD n ∂Y n ∂x j + ρ n NS n=1 D n ∂Y n ∂x j , (2.12)
where the first term on the right-hand side represents the effective diffusion velocity and the second one is a mass corrector term that should be taken into account in order to satisfy the continuity equation when dealing with non-constant species diffusion coefficients [START_REF] Poinsot | A theoretical decomposition of mean skin friction generation into physical phenomena across the boundary layer[END_REF]. Specifically, D n is an equivalent diffusion coefficient of species n into the mixture, computed following Hirschfelder's approximation [START_REF] Hirschfelder | Simulation and modeling of cold-wall hypersonic turbulent boundary layers on flat plate[END_REF], starting from the binary diffusion coefficients which are curve fitted in [START_REF] Gupta | A review of reaction rates and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30 000[END_REF]. Note that the molecular weight gradient contribution is neglected in (2.12), which therefore represents a rather simple model (albeit allowing for variable mass diffusion coefficients and non-constant Lewis numbers). The five species interact with each other through a reaction mechanism consisting of five reversible chemical steps [START_REF] Park | Nonequilibrium Hypersonic Aerothermodynamics[END_REF])

R1 : N 2 + M ⇐⇒ 2N + M R2 : O 2 + M ⇐⇒ 2O + M R3 : NO + M ⇐⇒ N + O + M R4 : N 2 + O ⇐⇒ NO + N R5 : NO + O ⇐⇒ N + O 2 , ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (2.13)
with M being the third body (any of the five species considered). Dissociation and recombination processes are described by reactions R1, R2 and R3, whereas the shuffle reactions R4 and R5 represent rearrangement processes. The mass rate of production of the nth species is governed by the law of mass action

ωn = M n NR r=1 ν nr -ν nr × k f ,r NS n=1 ρY n M n ν nr -k b,r NS n=1 ρY n M n ν nr , ( 2.14)
where ν nr and ν nr are the stoichiometric coefficients for reactants and products in the rth reaction for the nth species, respectively, and NR is the total number of reactions. Furthermore, k f ,r and k b,r denote the forward and backward reaction rates of reaction r, modelled by means of Arrhenius' law. The tight coupling between chemical and thermal non-equilibrium, due to their concurrent presence in such flows, is taken into account by means of a suitable modification of the temperature values used for computing the reaction rates. A geometric-averaged temperature is considered for the dissociation reactions R1, R2 and R3 in (2.13), computed as T avg = T q T 1-q V with q = 0.7 [START_REF] Park | Two-temperature interpretation of dissociation rate data for N 2 and O 2[END_REF]). Lastly, the vibrational-translational energy exchange is computed as

Q TV = NM m=1 Q TV,m = NM m=1 ρ m e Vm (T) -e Vm (T V ) t m , (2.15)
where t m is the mth molecular relaxation time evaluated by means of the expression of [START_REF] Millikan | Systematics of vibrational relaxation[END_REF]. Specifically, the relaxation time of the mth molecule with respect to the nth species writes (2.16) where p is the pressure, p atm = 101 325 Pa and a mn and b mn are coefficients reported in [START_REF] Park | Review of chemical-kinetic problems of future NASA missions. I-Earth entries[END_REF]. Since this expression tends to underestimate the experimental data at temperatures above 5000 K, a high-temperature correction was proposed by [START_REF] Park | Assessment of two-temperature kinetic model for ionizing air[END_REF] t mn = t MW mn + t c mn with t c mn = φ mn M m σ .

t MW mn = p p atm exp a mn (T -(1/3) -b mn ) -18.42 ,
(2.17)

Here, φ mn = M m M n /(M m + M n ) and σ = 8R/Tπ(7.5 × 10 -12 NA/T), NA being Avogadro's number. The mean value is then evaluated with a weighted harmonic average

t m = NS n=1 ρ n M n NS n=1 t mn ρ n /M n .
(2.18)

The complete formulation of transport coefficients laws and thermochemical models is provided in Appendix A.

Numerical method

The Navier-Stokes equations are integrated numerically by using a high-order centred finite-difference scheme [START_REF] Sciacovelli | Assessment of a high-order shock-capturing central-difference scheme for hypersonic turbulent flow simulations[END_REF]. The convective fluxes are discretized by means of central tenth-order differences, supplemented with a higher-order adaptive artificial dissipation. This consists in a blend of a nineth-order accurate dissipation term based on tenth-order derivatives of the conservative variables (used to damp grid-to-grid oscillations) along with a low-order shock-capturing term. A highly selective sensor, based on Ducros' extension of Jameson's pressure-based sensor [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF]) is used to turn on shock capturing in the immediate vicinity of flow discontinuities for all equations except the vibrational energy equation. For the latter, a shock sensor based on second-order derivatives of the vibrational temperature was preferred to ensure appropriate damping of spurious oscillations. Standard fourth-order differences are used for the viscous fluxes. Time integration is carried out by means of a third-order total variation diminishing (TVD) Runge-Kutta scheme [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF]. More details about the present numerical technique, as well as a complete assessment for a variety of highly compressible flow problems including chemically reacting hypersonic boundary layers can be found in [START_REF] Sciacovelli | Assessment of a high-order shock-capturing central-difference scheme for hypersonic turbulent flow simulations[END_REF].

Computational set-up

The configuration under investigation is a spatially evolving, zero-pressure-gradient flat-plate boundary layer, sketched in figure 1. The prescribed edge conditions of M e = 12.48, T e = 594.3 K and p e = 4656 Pa are representative of those downstream of a shock wave generated by a 6 • sharp wedge flying at M = 20 at an altitude of approximately 36 km. The stagnation enthalpy at the edge of the boundary layer is H e = h e + u 2 e /2 = 18.66 MJ kg -1 , a value comparable to those of the high-enthalpy cases considered by Duan & Martín (2011b) and Di Renzo & Urzay (2021). Air at such free-stream conditions is supposed to be in thermochemical equilibrium (X N 2 = 0.79, X O 2 = 0.21, X n being the nth species molar fraction). Of note, a similar scenario has been considered by [START_REF] Kline | Multiple boundary layer instability modes with nonequilibrium and wall temperature effects using lastrac[END_REF] for stability studies. The computational domain is a rectangular box enclosed within the shock layer, highlighted in red in figure 1. The extent of the domain is L x × L y × L z = 3000δ in × 120δ in × 30πδ in , with x, y and z denoting the streamwise, wall-normal and spanwise directions, respectively, and δ in the displacement thickness of the boundary layer at the inlet section, computed as δ = δ 0 (1 -ρu/ρ e u e ) dy (δ being the boundary-layer thickness at 99 % of the edge velocity). The computational grid is N x × N y × N z = 9660 × 480 × 512, for a total of approximately 2.4 billions grid points. The grid spacing is uniform in the streamwise and spanwise directions, whereas a constant grid stretching of 0.7 % is applied in the wall-normal direction. The calculation is initiated in the laminar region at a distance x 0 from the plate leading edge. The profiles of ). Unfortunately, no information is available for turbulent flows. Nevertheless, the contribution of the vibrational heat transfer to the total wall heat flux will be shown to be small with respect to the roto-translational counterpart (as later detailed in § 3), justifying a posteriori the choice of a Dirichlet boundary condition as a first approximation. The selected wall temperature value, on the other hand, is representative of realistic flight data [START_REF] Park | Stagnation-point radiation for apollo 4[END_REF]).

Transition to turbulence is induced by means of a suction and blowing strategy, similar to the one adopted in [START_REF] Passiatore | Finite-rate chemistry effects in turbulent hypersonic boundary layers: a direct numerical simulation study[END_REF]. Specifically, the following wall-normal velocity is imposed along a wall strip close to the inflow:

v wall u ∞ = e -0.4g(x) 2 A sin (2πg(x) -ωt) + cos(βz) + 0.05 sin 2πg(x) -ωt + π 4 cos(βz) , (2.19) where g(x) = (x -x forc )/(δ in √ 2σ
), with √ 2σ = 0.85(2π/ω) and x forc is the centre of the Gaussian-like distribution which modulates the forcing function. The Reynolds number at the forcing location is Re δ forc = 10 4 , corresponding to x forc = 7.4 × 10 -2 m. In (2.19), the non-dimensional pulsation ω = ωδ in /c ∞ (c ∞ being the free-stream speed of sound) corresponds to a dimensional frequency of f = ω/2π = 75 kHz and β = 0.4/δ in is the spanwise wavenumber. Lastly, A is the forcing amplitude, which has been set equal to 5 %. The presence of stationary modes at a strong amplitude in the suction-and-blowing injection strip mimics the effects of roughness spots and is used to speed up breakdown to turbulence within the prescribed computational domain.

Data collection and analysis

In the results presented below, flow statistics are computed by averaging in time and in the spanwise homogeneous direction, after that the initial transient is evacuated. Firstand second-order moments of various flow quantities will be presented and discussed. For a given variable f , we denote with f = ff the standard time and spanwise average, with f the corresponding fluctuation, whereas f = ff denotes the density-weighted Favre averaging, with f the Favre fluctuation and f = ρf / ρ. The sampling time interval is constant and corresponds exactly to 300 samples per each period of the forcing harmonic; specifically, Δt + stats = Δt stats (u 2 τ ρ w /μ w ) = 0.16, with u τ = τ w /ρ w the friction velocity based on the averaged wall shear stress τ w at the end of the computational domain. Data are collected for more than three turnover times, corresponding to T stats (u 2 τ /μ w /ρ w ) ≈ 6700 for a total of ≈40 000 temporal snapshots. In addition to flow statistics, instantaneous planes and specific meshlines are extracted with a frequency twenty and forty times higher with respect to the fundamental mode in the forcing function, for a total of 3200 and 6400 samples, respectively. The analysis will mainly focus on five selected streamwise stations; one located in the laminar region (for purpose of comparison), one in the transitional region and the last three in the turbulent portion of the domain. Table 1 reports the values of the Reynolds number based on the distance from the leading edge Re x = ρ e u e x/μ e at the selected stations and some boundary-layer properties. The Reynolds number based on the local momentum thickness Re θ = ρ e u e θ/μ e (with θ = δ 0 ρu/ρ e u e (1 -ρu/ρ e u e ) dy) reaches values close to ≈ 6000 in the fully turbulent region, corresponding to Re inc θ = Re θ μ e /μ w of approximately 3000. The displacement-thickness-based Reynolds number Re δ = ρ e u e δ /μ e and the friction Reynolds number Re τ = ρ w u τ δ/μ w reach values up to ≈16 × 10 4 and ≈1100, respectively. Of note, the grid spacings ensure a DNS-like spatial resolution everywhere. Here, the notation '• + ' denotes normalization with respect to the viscous length scale l v = μ w /(ρ w u τ ). Unless otherwise specified, the wall-normal evolution of statistics is displayed in inner semi-local units y = ρu τ y/ μ, with u τ = √ τ w / ρ. 1. Boundary-layer properties at five selected downstream stations. In the table, Ma τ = u τ /c w is the friction Mach number, B TR q = q TR w /(ρ w u τ h w ) and B V q = q V w /(ρ w u τ h w ) are the roto-translational and vibrational dimensionless heat fluxes. Lastly, Δx + , Δy + w , Δy + δ and Δz + denote the grid sizes in inner variables in the x-direction, y-direction at the wall and at the boundary-layer edge and in the z-direction, respectively.

Results

Global flow properties

The streamwise evolution of selected quantities at the wall is first discussed. Figures 2(a) and 2(b) report the distributions of the skin friction coefficient C f and the heat flux coefficients C q and C V q , defined as

C f = 2τ w ρ e u 2 e , C TR q = q TR w ρ e u 3 e , C V q = q V w ρ e u 3 e , (3.1a-c)
the net total heat flux being given by the sum of the two contributions. A ramp up starting at (xx 0 )/δ in ≈ 1100 leads to a sharp overshoot in the C f profile, with a peak of the wall shear stress almost five times larger than its corresponding laminar value. The flow achieves a turbulent regime starting from (xx 0 )/δ in ≈ 1800, with a smoothly decreasing C f up to the end of the computational domain. The roto-translational contribution of the wall heat flux C TR q is shown to be largely predominant with respect to the vibrational counterpart C V q (which is multiplied by ten to match the range of figure 2b), mainly because of the vibrational thermal conductivity λ V , whose values are approximately one order of magnitude smaller than λ TR . While C TR q closely follows the skin friction coefficient profile, C V q exhibits a different evolution. It rapidly decreases after the forcing strip and stays approximately constant up to the breakdown to turbulence; afterwards, the value in the turbulent region almost doubles the one in the laminar regime and keeps a constant value up to the end of the domain. The different trend of the heat fluxes is related to the evolution of vibrational temperature gradients across the boundary layer, as discussed later in § 3.3. In order to isolate the contributions of the mean and fluctuating field to the skin friction, the decomposition of Renard & Deck (2016) (later extended to compressible flows by [START_REF] Li | Decomposition of the mean skin-friction drag in compressible turbulent channel flows[END_REF] has been computed; of note, the general hypotheses under which it is derived allows a straightforward application even in presence of thermal and chemical non-equilibrium effects. The skin friction coefficient can then be rewritten under the

C q V ×10 C f (×10 -3 ) (×10 -4 ) (a) ( b) C q TR Figure 2.
Evolution of (a) the skin friction coefficient C f (black line) and (b) the heat flux coefficients C TR q and C V q (black and red lines, respectively). The filled symbols in (a) denote the streamwise location of the five stations selected in table 1.

following form:

C f = 2 ρ ∞ u 3 ∞ δ 0 τ xy ∂ ũ ∂y dy C f ,1 + 2 ρ ∞ u 3 ∞ δ 0 -ρu v ∂ ũ ∂y dy C f ,2 + 2 ρ ∞ u 3 ∞ δ 0 (ũ -u ∞ ) ρ ũ ∂ ũ ∂x + ṽ ∂ ũ ∂y - ∂ ∂x τ xx -ρ u u -p dy C f ,3 , (3.2)
where C f ,1 , C f ,2 and C f ,3 represent the mean-field molecular dissipation, the turbulent dissipation and the effects related to boundary-layer spatial growth, respectively. The sum of the three terms and their separate contributions are displayed in figure 3. An excellent agreement with respect to the direct C f computation is observed in the pseudo-laminar and fully turbulent regions, with only minor deviations in the transition region. The C f ,3 contribution is negligible everywhere but at the breakdown-to-turbulence location, where it takes negative values due to the abrupt thickening of the boundary layer. The decrease of C f ,3 is counterbalanced by a large increase of the Reynolds stress-related term, whereas the growth of the mean-field contribution, C f ,1 , is slightly delayed with respect to the other two. In the turbulent region, C f ,1 and C f ,2 are almost superposed, differently from what observed by [START_REF] Passiatore | Finite-rate chemistry effects in turbulent hypersonic boundary layers: a direct numerical simulation study[END_REF] where the former was shown to be predominant. Such a discrepancy can be ascribed to the much larger friction Reynolds numbers reached in the current configuration, leading to an increased turbulent contribution [START_REF] Fan | Decomposition of the mean friction drag in zero-pressure-gradient turbulent boundary layers[END_REF].

Mean flow analysis

Different scalings for the profiles of the averaged streamwise velocity have been tested. First, the transformations of [START_REF] Van Driest | The Problem of Aerodynamic Heating[END_REF] and [START_REF] Trettel | Mean velocity scaling for compressible wall turbulence with heat transfer[END_REF] high speeds [START_REF] Passiatore | Finite-rate chemistry effects in turbulent hypersonic boundary layers: a direct numerical simulation study[END_REF]. On the other hand, it becomes inaccurate when highly cooled configurations are considered, whether they be boundary layers [START_REF] Zhang | Direct numerical simulation database for supersonic and hypersonic turbulent boundary layers[END_REF][START_REF] Hirschfelder | Simulation and modeling of cold-wall hypersonic turbulent boundary layers on flat plate[END_REF], pipe flows [START_REF] Ghosh | Compressible turbulent channel and pipe flow: similarities and differences[END_REF] or channel flows [START_REF] Modesti | Reynolds and mach number effects in compressible turbulent channel flow[END_REF][START_REF] Sciacovelli | Direct numerical simulations of supersonic turbulent channel flows of dense gases[END_REF].

u VD = 1 u τ ū 0 ρ ρ w du, u TL = ū 0 ρ ρ w 1 + 1 2 1 ρ d ρ dy y - 1 μ d μ dy y du (3.
.0 C f 1.2 (×10 -3 ) (×10 -3 ) 1000 2000 3000 -0.2 0 0.2 0.4 0.6 (x -x 0 )/δ in (x -x 0 )/δ in C f,1 C f,2 C f,3 (b) (a)
u GFM = δ 0 1 μ + ∂u + ∂y * 1 + 1 μ + ∂u + ∂y * -μ + ∂u + ∂y + dy * (3.4)
where μ + denotes normalization of the mean viscosity with respect to its wall value. Such a scaling has been shown to successfully collapse channels, pipe flows and boundary-layer configurations, even at large Mach numbers. A comparison between u TL and u GFM as a function of y is shown in figure 4(c) for a velocity profile extracted at Re θ = 6200. In the inner layer, viscous stresses are predominant and roughly correspond to the total shear stress resulting in a very good collapse for both scalings. On the contrary, turbulent shear stresses become dominant in the logarithmic region, which effect is most correctly taken into account by the total-stress-based scaling. This results in a better collapse of u GFM onto the universal logarithmic profile with respect to u TL , for which the slope of the logarithmic region is largely overestimated. Yet, the u GFM transformation still predicts a higher slope and intercept in the logarithmic region compared with the classical incompressible values, in accordance with the recent study of [START_REF] Lee | New velocity transformation based on compressibility-corrected total stress[END_REF]. Finally, we verified that the transformation based on the constant-stress-layer assumption (i.e. τ/τ w ≈ 1 across the boundary layer, leading to (3.4)) and the one based on the actual total shear stress (extracted from DNS data) do not exhibit any discernible differences, confirming the Figure 5 shows the Reynolds stress profiles at the turbulent station Re θ = 6200 normalized with respect to the semi-local friction velocity u τ ; results are compared with low-enthalpy configurations extracted from [START_REF] Zhang | Direct numerical simulation database for supersonic and hypersonic turbulent boundary layers[END_REF] (cases M14Tw018 and M8Tw048) and [START_REF] Xu | Compressibility effect in hypersonic boundary layer with isothermal wall condition[END_REF] (case M8T1). Despite the important differences in the values of the free-stream Mach numbers, friction Reynolds numbers, absolute wall temperatures and wall-cooling rates, the turbulent intensity profiles are comparable and do not exhibit any marked influence attributable to TCNE effects. Consistently with previous observations (Duan et al. 2011;[START_REF] Lagha | A numerical study of compressible turbulent boundary layers[END_REF][START_REF] Zhang | Direct numerical simulation database for supersonic and hypersonic turbulent boundary layers[END_REF], the larger streamwise component and the smaller cross-flow one may indicate the presence of strong compressibility effects, as discussed later in § 3.6. Another scaling often investigated in wall-bounded compressible turbulence is the one that relates the mean velocity to the mean temperature profile for zero-pressure-gradient boundary layers. The modified Crocco relation derived by [START_REF] Walz | Boundary Layers of Flow and Temperature[END_REF] writes

T T e = T w T e + T aw -T w T e ũ u e + T e -T aw T e ũ u e 2 , ( 3.5)
with T aw /T e = 1 + r((γ -1)/2)M 2 e and r the recovery factor set equal to 0.9; the relation is shown in figure 6(a). In previous studies with high Mach numbers, wall-cooled, high-enthalpy boundary layers (Duan & Martín 2011b), it was found that such a relation deviates from the exact T/T e profile extracted from DNS data; a significant discrepancy is indeed shown in the range 0.2 < ũ/u e < 0.7. This should not be surprising since relation (3.5) was derived under calorically perfect gas hypotheses. In order to remove the explicit dependence on thermal and chemical models, Duan & Martín (2011b) where h aw = h e + 1 2 ru 2 e . The function f (ũ/u e ) has to be independent of free-stream conditions, wall temperature and surface catalysis (if any). For calorically perfect gases, D. Passiatore, L. Sciacovelli, P. Cinnella and G. Pascazio

Thermal non-equilibrium effects

The presence of thermal non-equilibrium is analysed by comparing the vibrational relaxation time with characteristic time scales of the dynamic field. For that purpose, the following vibrational Damköhler numbers are considered:

Da VS v = μ w /τ w t m , Da LE v = δ /u τ t m , Da RT v =
x/u e t m , (3.8a-c)

with Da

(•) v = O(1)
denoting the regime in which non-equilibrium effects are relevant. The numerator of each Damköhler number is fixed for a given streamwise station and equal for each of the three molecules; the different trends of these parameters, for each of the three diatomic species, are therefore entirely dictated by the species molecular relaxation times t m , which depend upon the local state of the gas. The evolution of these dimensionless numbers is displayed in figure 7. The first ratio in (3.8a-c) relates the characteristic time of fluid motion in the viscous sublayer (VS) to the vibrational relaxation time of molecule m. Figure 7(a,d,g) shows the evolution of this quantity for the three molecules present in the mixture. It is shown that Da VS v 1 everywhere for N 2 and O 2 whereas it is of the order of unity only for NO, whose concentration is, however, almost negligible (see § 3.5). The decoupling of the two time scales indicates that there exists a layer of thickness l v in the inner (colder) region of the boundary layer in which the molecules are vibrationally frozen. The higher values of Da VS v 1 observed at the laminar position are due to the predominance of viscous stresses. Of major interest is Da LE v , representing the ratio between the large-eddy (LE) turnover time and the relaxation time, constructed by means of the displacement thickness δ and u τ . The orders of magnitude of the two time scales are comparable for molecular nitrogen (figure 7b), for which the peak of Da LE v ranges from 0.4 to 0.7 at y ≈ 10 as the flow evolves in the turbulent region. On the other hand, figures 7(e) and 7(h) show that Da LE v 1 for O 2 and NO, implying that the two molecules reach a vibrational equilibrium with respect to the characteristic time of the large-scale turbulent motions much faster. Similar considerations hold for Da RT v , which compares the vibrational relaxation time with the local flow residence time (RT). As the flow evolves along the plate, there is enough time for O 2 and NO to achieve thermal equilibrium, whereas N 2 remains significantly out of equilibrium at all stations. This result is a consequence of the high vibrational temperature (and slower relaxation rate) characteristic of molecular nitrogen. Of note, both Da LE v and Da RT v exhibit much higher values at the three turbulent stations due to the contribution of turbulent mixing. Wall-normal profiles of the normalized mean translational and vibrational temperatures are reported in figures 8(a) and 8(b), respectively. The translational temperature reaches a peak value of approximately [2.3 -2.5]T w , at a wall-normal location y ≈ 10, corresponding to the peak of turbulent intensity. The vibrational temperature T V is 20 % to 30 % smaller than T across the boundary layer, proving the existence of a vibrational-relaxation layer for the thermodynamic conditions under investigation, and also exhibits a peak located between the buffer layer and the lower logarithmic region. The major gap is visible at Re θ = 2960, which corresponds to the C f peak. Additionally, the wall gradient ∂ T V /∂y is little sensitive to the peak's magnitude and location, explaining the nearly constant trend of C V q previously observed in figure 2(b). The normalized temperature difference ( T -T V )/T w is reported in figure 9(a,b) as a function of y and y/δ, respectively. Such a quantity increases in the transitional region, reaches a maximum of the order of ≈2000 K where C f peaks, and then tends to decrease moving towards the turbulent region. The wall-normal maximum value moves from the outer part of the boundary layer (y/δ ≈ 0.6) to the buffer layer, where it stabilizes at y ≈ 8. Thermodynamic non-equilibrium conditions are kept until the end of the computational domain, as already confirmed by inspection of the vibrational Damköhler number values. Positive T -T V values are observed from the VS (we recall that T = T V = 1800 K at the wall) up to the first portion of the logarithmic region. However, as the flow evolves along the flat plate, a region characterized by T -T V < 0 appears, meaning that the flow becomes vibrationally over-excited. The height at which the sign changes smoothly moves from the edge of the boundary layer up to y/δ ≈ 0.45, due to the boundary-layer thickening. This is clearly visible in figure 10, where we report a snapshot of the temperature difference along a longitudinal cut plane (top) and the isocontours of the second invariant of the velocity gradient tensor coloured by (T -T V ). [START_REF] Fiévet | Numerical investigation of the coupling of vibrational nonequilibrium and turbulent mixing using state-specific description[END_REF], who conducted an analysis on the coupling between vibrational non-equilibrium and turbulent mixing, observed that the majority of non-equilibrium states occur in under-excited conditions. The rare regions in which there exists an over-excited non-equilibrium state are dictated by the link between the roto-translational and mechanical energies. The static temperature is coupled with the other variables through the equation of state; hence, a decrease of turbulent intensity (moving towards the edge of the boundary layer) produces a decrease of the translational temperature. At the same time, the vibrational energy remains nearly constant because it relaxes slowly, leading to an over-excitation in the outer layer, until the free-stream equilibrium conditions are reached by both temperatures. Normalized root mean square (r.m.s.) of the Favre fluctuations of translational and vibrational temperatures across the boundary layer are displayed in figure 11. In figure 11(a), a peak emerges on each side of the Tmax location (y ≈ 10), as classically observed for wall-cooled TBLs [START_REF] Duan | Direct numerical simulation of hypersonic turbulent boundary layers. Part 2. Effect of wall temperature[END_REF][START_REF] Zhang | Direct numerical simulation database for supersonic and hypersonic turbulent boundary layers[END_REF]Di Renzo & Urzay 2021). The predominance of the outer r.m.s. peak in the turbulent region is mainly caused by the decrease of T toward the outer boundary, the absolute temperature fluctuations being comparable in the inner and outer regions. In contrast, the vibrational temperature fluctuations exhibit no local minimum (figure 11b), and increase monotonically up to the boundary-layer edge, where they drop abruptly. This is in contrast to profiles typically observed in the laminar region (see the profile at Re θ = 1600 reported in the figure for comparison), and can be ascribed to thermal non-equilibrium: the vibrational relaxation times being slower than the translational ones and comparable two temperatures conditioned on sweep events (characterized by u > 0 and v < 0) and ejection events (u < 0 and v > 0) at Re θ = 6200 and for several wall-normal locations. Within the VS, sweeps entrain hot fluid from the overlying regions. Consequently, they are characterized by higher mean temperatures with respect to the total average in the inner region. The opposite is true for ejection events. The largest differences between sweeps and ejections conditionally averaged temperatures are registered for the roto-translational modes due to the sharper peak observed for T. An opposite trend is observed at normal locations beyond the temperature peak (y ≈ 10). In this case, sweeps bring colder air from the outer layers towards the hotter, inner region; here, the translational modes rapidly reach equilibrium, whereas a longer relaxation time is needed for the vibrational ones. As a result, T V is always smaller than T. On the other hand, ejection motions bring hotter, thermally out-of-equilibrium mixture towards the outer region of the boundary layer, where vibrational relaxation times rapidly increase because of the colder environment; therefore, in this region T V > T. The largest temperature difference between sweeps and ejections events occurs in the first portion of the logarithmic region (up to y ≈ 100), with values as large as 800 K for vibrational modes and 1000 K for the roto-translational ones. Further evidence of this phenomenon is provided by the conditional probability density functions (p.d.f.s) of the two temperatures, shown in figure 14. Three wall-normal stations T ˜ej.

μ ¯w/τ ¯w t N 2 t N 2 t O 2 t O 2 t O 2 t NO t NO t NO δ /u τ δ /u τ δ /u τ x/
(h) (g) ( i) (d ) ( f )
T V T V SW. T V ej.
Figure 13. Mean values of T and T V conditioned on sweep and ejection events, at Re θ = 6200 at various wall-normal heights. Blue lines: roto-translational modes; red lines: vibrational modes; dashed lines, diamonds: average over sweep events; dashed-dotted lines, triangles: average over ejection events; solid lines, circles: global average. are considered: the location where the maximum temperature difference is registered (panels a and b, y ≈ 10), the logarithmic region (panels c and d, y ≈ 100) and the boundary layer edge (panels e and f, y ≈ 3000). Close to the boundary layer edge, a heavy right tail is observed for the p.d.f.s of both events, with T V values slightly larger than T ones. Moving inside the boundary layer, at y ≈ 100, both p.d.f.s tend to become symmetric. The most remarkable differences are observed in the buffer layer with left-skewed distributions for T and symmetric ones for T V , their shape being related to the presence of a local minimum for T 2 and a small plateau for T 2 V as shown in figure 11.

Thermal modelling and Reynolds analogies

For thermally out-of-equilibrium flows, a new unclosed term appears by Favre averaging the vibrational energy equation (2.5), namely, the turbulent transport of vibrational energy ρu j e V . Such a term is the analogue of the turbulent transport term appearing in the total energy equation. The latter is usually modelled by introducing a turbulent Prandtl number Pr t , defined as

Pr t = ρu v ∂ T/∂y ρv T ∂ ũ/∂y . (3.9)
According to the classical strong Reynolds analogy (SRA), Pr t is expected to take values close to one throughout the flow [START_REF] Morkovin | Effects of compressibility on turbulent flows[END_REF]. Figure 15(a) (displayed in external units for better clarity) shows that the estimation fails at the peak of T, located at y/δ ≈ 0.03, since the coexistence of null temperature gradient and turbulent heat flux causes an overshoot. Such a behaviour was previously observed by Duan & Martín (2011b) and Di Renzo & Urzay (2021) (for a discussion about the validity of this assumption for adiabatic hypersonic boundary layers, see [START_REF] Passiatore | Finite-rate chemistry effects in turbulent hypersonic boundary layers: a direct numerical simulation study[END_REF]. By analogy, the turbulent flux of vibrational energy can be modelled by introducing a vibrational turbulent Prandtl number

Pr V t = ρu v ∂ T V /∂y ρv T V ∂ ũ/∂y . (3.10)
Here, Pr V t (shown in figure 15b) takes values close to unity in a relatively large region corresponding to 0.2 y/δ 0.7, where turbulent transport dominates, and it diverges at the boundary-layer edge (where ∂ ũ/∂y tends to vanish) and close to the wall, similarly to Pr t . We can get further insights by reporting, in figure 16, correlations of u and v fluctuations with T and T V , respectively. As observed in both low-and high-enthalpy boundary layers [START_REF] Zhang | Direct numerical simulation database for supersonic and hypersonic turbulent boundary layers[END_REF]Di Renzo & Urzay 2021), the logarithmic and outer regions are characterized by a negative correlation between u and T (a) and positive correlation for v and T (b). This comes from (i) sweeps (u > 0, v < 0) entraining colder (T < 0) fluid towards the buffer layer, and (ii) ejections (u < 0, v > 0) entraining hotter (T > 0) fluid towards the edge, resulting in R u T < 0 and R v T > 0. Below the buffer layer, both correlations change their sign for the same reason, consistently with the trend previously discussed in figure 13. Similar considerations hold for R u T V and R v T V . However, while streamwise velocity fluctuations close to the wall are perfectly correlated with the translational temperature ones (R u T ≈ 1), the correlation coefficient with the vibrational temperature (c) is weaker (R u T V is between 0.7 and 0.85). This effect is attributed once again to the larger relaxation time of vibrational modes, T V needing more time than T for adapting to the surrounding temperature. Together with the evolution of Pr t and the temperature-velocity fluctuation correlations, the SRA hypothesis also predicts that

T / Trms (γ -1)M 2 (u /ũ rms ) ≈ 1. (3.11)
As already observed in the literature, the relation in (3.11) exhibits large deviations from the SRA estimation (not shown) since it fails to take into account the diabatic behaviour of the wall. Several modified forms have been proposed, one of the most successful being that of Huang, Coleman & Bradshaw (1995)

T / Trms (γ -1)M 2 (u /ũ rms ) ≈ 1 Pr t (1 -(∂ T t /∂ T)) , ( 3.12)
usually referred to as HSRA, with T t the stagnation temperature. By removing the assumption of calorically perfect gas (Duan & Martín 2011b), one can write a generalized HSRA (GHSRA)

T rms ≈ - 1 Pr t ∂ T ∂ ũ u rms . (3.13)
The wall-normal evolution of the ratio between the right-and the left-hand sides of (3.13) is plotted in figure 17(a), where a satisfactory agreement is observed for most of the boundary-layer thickness (at least in regions where both velocity and temperature fluctuations are not nullified). By using the same procedure followed for Pr t , the GHSRA can also be extended to the vibrational temperature such that

T V,rms ≈ - 1 Pr V t ∂ T V ∂ ũ u rms . (3.14)
The results shown in figure 17(b) prove that the GHSRA maintains its validity under thermal non-equilibrium conditions and holds also for vibrational energy.

In addition to turbulent transport, further attention should be paid to the modelling of the nonlinear turbulence-thermal non-equilibrium interactions induced by the source terms of the vibrational energy equation, which are associated with translational-vibrational energy exchanges (Q TV ) and chemical dissociation/recombination processes ( m ωm e Vm ). While thermal non-equilibrium can be reasonably quantified by the magnitude of T -T V , the analysis of the source terms allows one to isolate the physical processes originating such non-equilibrium conditions in the first place. In this specific configuration, chemical activity is extremely weak and therefore m ωm e Vm is almost negligible; as a consequence, changes in T -T V are tightly related to Q TV values. Positive values indicate vibrational under-excitation since T > T V (and the opposite for negative values). The wall-normal profiles of Q TV , shown in figure 18(a) as a function of y and y/δ, closely follow those shown for ( T -T V ) in figures 9(a) and 9(b), albeit that the negative values visible in the zoomed view of the outer region are much less pronounced (two orders of magnitude lower) than those registered in the inner zone. Such a large difference is due to the concurrent effects of lower absolute temperature values and (to a lesser extent) smaller ( T -T V ) differences, causing the denominator to increase and the numerator to decrease in (2.15), respectively. Although the physical significance of ( T -T V ) and Q TV is similar in the present case, the interest in analysing the latter quantity stems from a modelling standpoint; that is, only an appropriate modelling of the source term can ensure a correct description of the non-equilibrium state. In this respect, following the technique used by Duan & Martín (2011a) for assessing turbulence/chemistry interactions and with the aim of quantifying the coupling between turbulence and thermal relaxation, we built an indicator defined as

Q I TV = Q TV (T, T V , ρ, p, Y n ) -Q TV ( T, T V , ρ, p, Y n ) Q TV ,max , (3.15)
where Q TV ,max denotes the maximum wall-normal value of Q TV at the selected station.

Deviations from null values represent a measure of vibrational energy variations associated with turbulent fluctuations. Figure 18(b) shows that computing Q TV from averaged quantities would underestimate its value across the whole boundary layer, with discrepancies larger than 15 % of Q TV ,max in the inner layer. These data further confirm the tight coupling between thermal non-equilibrium and turbulence discussed in the current section and underline that attention must be paid to the development of suitable turbulence closures, especially in the context of RANS modelling.

Chemical activity

In this section we investigate chemical activity triggered in the high-enthalpy boundary layer by friction heating. The profiles of Favre-averaged mass fractions of species present in the mixture are reported in figure 19. At the flow conditions under investigation, molecular oxygen is the species which dissociates most, leading to the formation of atomic oxygen. The maximum mass fraction of atomic oxygen in the mixture is approximately 1 %. The maximum temperatures in the boundary layer are too low to allow nitrogen dissociation in any significant amount, even at the temperature peak, and consequently the quantity of atomic nitrogen is negligible compared with the other species. Nevertheless, the present conditions are such that a small amount of nitric oxide (approximately 0.7 %) is produced by the forward shuffle reaction R4 and by the backward shuffle reaction R5. The rather low chemical activity is partly due to the strong wall cooling, limiting the maximum temperature to values just above the threshold for oxygen dissociation (T > 2000 K). Thermal non-equilibrium further limits chemical reactions. Highly vibrationally excited molecules are meant to dissociate more rapidly, since they require less energy exchange during the collisions. On the contrary, in under-excitation conditions, dissociation is less likely to occur. This mechanism is numerically taken into account in the computation of the reaction rate coefficients (for the calculation of which a geometric average of the two temperatures is considered, T 0.7 T 0.3 V ); they turn out to be lower when strong non-equilibrium effects are present. On the other hand, the endothermic chemical reactions tend to drain roto-translational energy, although this effect is partly compensated by the relaxation of vibrational modes, which supply energy and partly counteract temperature reduction. The chemical Damköhler number based on the LE turnover time (δ /u τ )/(ρ/ ωn ) is at most of the order of 10 -2 for all species, denoting little interaction between chemical reactions and turbulent motions and an almost frozen-chemistry behaviour. A detailed discussion can be found in [START_REF] Passiatore | Finite-rate chemistry effects in turbulent hypersonic boundary layers: a direct numerical simulation study[END_REF] for an adiabatic turbulent chemically out-of-equilibrium boundary layer; despite the different conditions, similar qualitative observations can be made in the present work.

Compressibility effects

The present analysis is completed by discussing the role of compressibility effects and their interaction with thermal non-equilibrium conditions. In figure 20 we report the evolutions of the turbulent and fluctuating Mach numbers, defined as M t = u rms /c and M rms = (u/c) rms , respectively. Despite the high external Mach number, M t is shown to be smaller than one everywhere, with maximum values of approximately 0.8-0.9 registered at the location where the turbulent kinetic energy production peaks. Such values are, however, large enough to enable the formation of eddy shocklets, as also proven by the location of the sonic line (identified by the three vertical lines in figure 20a) within the VS. Accordingly, the fluctuating Mach number takes values larger than one starting from y 10. The outer M rms peaks at y ∈ [1000,2000] (whose values are compatible with those obtained by [START_REF] Zhang | Direct numerical simulation database for supersonic and hypersonic turbulent boundary layers[END_REF] for high-speed boundary layers of calorically perfect gases) are linked to the large fluctuations of thermodynamic properties resulting from the presence of turbulent/non-turbulent interfaces along the boundary-layer edge. These structures are clearly visible in figure 21, displaying the normalized density gradient magnitude on a streamwise slice of the computational domain. Such large fluctuations of the thermodynamic variables are not accompanied by an intense activity of the dynamic field, as testified by the values of the semi-locally scaled fluctuating velocity divergence ϑ μ/τ w (with ϑ = ∂u k /∂x k ) shown in figure 22 in the last portion of the domain. The stronger dilatation events are registered close to the wall and their magnitude monotonically decreases with the wall distance, as shown in figure 23(a), where results from [START_REF] Xu | Compressibility effect in hypersonic boundary layer with isothermal wall condition[END_REF] are reported for comparison. Figure 23(b) displays two p.d.f.s of ϑ μ/τ w computed at wall-normal positions corresponding to y ≈ 10 and 1000. Both the profiles are slightly skewed toward negative values, highlighting a mild predominance of the occurrence of strong compressions over strong expansions, although such events are much stronger in the inner region. However, joint p.d.f.s of (ϑ ( μ/τ w ), |∇ρ|/(ρ e /δ 99,end )) (not shown) did not exhibit any preferential correlations; this amounts to saying that, even if shocklets might exist inside the boundary layer, their number and/or intensity is not large enough to engender any remarkable effect on thermodynamic and dynamic quantities.

The interplay between compressibility and thermal non-equilibrium effects has previously been investigated in mixing layers [START_REF] Fiévet | Numerical investigation of the coupling of vibrational nonequilibrium and turbulent mixing using state-specific description[END_REF]) and homogeneous isotropic turbulence [START_REF] Zheng | Vibrational relaxation in compressible isotropic turbulence with thermal nonequilibrium[END_REF] configurations, in which it was observed that vibrational non-equilibrium is somehow enhanced by strong dilatation motions (whether they be expansions or compressions). Impulsive changes of thermodynamic quantities, as in the case of strong density gradients, can increase the difference between the two temperatures because of their different relaxation times. Here, high-|∇ρ| regions are mainly found near the boundary-layer edge where the average temperature values are 18. All the profiles are strongly right skewed because of the predominance (in magnitude and number) of under-excited states, as explained in § 3.3. At both locations, strong compressions are associated with slightly higher Q TV /Q TV,rms values than strong expansions; the wider differences between the distributions in panel (b) are caused by the much smaller source term fluctuations observed in the outer region. Globally, the differences are shown to be marginal, validating our previous statement.

Conclusions

The present work illustrates the results of the numerical simulation of a hypersonic, high-enthalpy TBL. The selected edge conditions are such that the flow undergoes TCNE conditions. The interaction between high-temperature effects and turbulent quantities is then investigated, with specific insights on the vibrational-relaxation dynamics. The forcing strip located right after the inlet similarity solution induces an abrupt transition to turbulence, characterized by a steep increment of the wall shear stress in the transitional region and an asymptotic decay in the fully turbulent portion. The skin friction decomposition of [START_REF] Poinsot | A theoretical decomposition of mean skin friction generation into physical phenomena across the boundary layer[END_REF] shows that the mean field and turbulent terms contribute equally to the C f value, whereas the term related to the boundary-layer growth is non-negligible only in the region of breakdown to turbulence. With regard to the wall heat flux, it is observed that its roto-translational contribution follows the same trend of C f , whereas the vibrational counterpart is nearly constant across the turbulent region, with values being one order of magnitude smaller. The van Driest and Trettel & Larsson transformations perform poorly in collapsing the streamwise velocity profiles onto the universal logarithmic law. On the contrary, the Trettel & Larsson transformation and the new total-stress-based transformation of [START_REF] Griffin | Velocity transformation for compressible wall-bounded turbulent flows with and without heat transfer[END_REF] provide a perfect match with the linear law, the latter improving the collapse in the logarithmic region as well. Furthermore, the relation between the mean velocity and mean temperature profiles derived by [START_REF] Walz | Boundary Layers of Flow and Temperature[END_REF] resulted to be inappropriate for calorically imperfect gases; an improvement was found with the analogous enthalpy-based equation proposed by Duan & Martín (2011b), despite the presence of thermochemical out-of-equilibrium effects.

With the aim of quantifying the thermal non-equilibrium, three vibrational Damköhler numbers are defined based on different characteristic flow time scales. The results illustrate that, aside from the near-wall region where the flow is vibrationally frozen, strong thermal non-equilibrium conditions persist up to the edge of the boundary layer. Molecular nitrogen is responsible for thermal non-equilibrium, since its characteristic temperature is high enough to cause slow relaxation rates. The vibrational excitation is largely sustained by turbulent mixing, as observed by the tight coupling with sweep and ejection events inside the boundary layer. Sweeps convey colder air from the outer regions towards the flat-plate wall, whereas ejections bring hotter air far from the inner layer. The translational modes rapidly adapt to such large temperature fluctuations, contrary to the vibrational modes which re-adjust slowly. At the macroscopic level, the direct consequence is that vibrational temperature lags behind the translational one up to the first half of the logarithmic layer; from there onwards, the flow transitions towards a vibrationally over-excited state which is kept until the edge of the boundary layer. Turbulence/thermal relaxation interactions, quantified by means of a parameter measuring the vibrational energy variations due to turbulent fluctuations, are shown to be substantial across the buffer layer, confirming that turbulent transport plays a major role. Following usual RANS modelling strategies, a vibrational turbulent Prandtl number Pr V t can be defined to close the turbulent vibrational energy flux and it assumes values close to unity for most of the boundary layer. Similarly to the classical Pr t , Pr V t deviates from unity in the inner region and at the boundary-layer edge, where ∂ ũ/∂y vanishes. The remaining SRA relations exhibit an acceptable agreement with the estimation and the GHSRA was found to work reasonably well even when extended to vibrational temperature.

Second-order statistics show classical self-similar trends for the Reynolds stresses, whose values are similar to those observed at low-enthalpy conditions for comparable values of the free-stream Mach numbers. The isothermal wall condition leads to a double-peak-shaped profile for the temperature fluctuations, the outer maximum being predominant. On the contrary, vibrational temperature fluctuations exhibit monotonously increasing profiles due to the marked non-equilibrium conditions existing in the buffer layer.

Gas dissociation at the wall is almost negligible since the imposed wall temperatures (T = T V = 1800 K) are not high enough to trigger significant chemical activity. Little interaction is observed between chemical reactions and turbulent motions, the Damköhler numbers being at most of the order of 10 -2 for all species. Statistical analyses on compressibility effects, measured by means of the semi-locally scaled velocity divergence fluctuations, have shown that strong dilatation motions are mainly found close to the wall, with strong compressions slightly more probable than strong expansions. With the aim of relating compressibility effects and thermal non-equilibrium, p.d.f.s of the translational-vibrational energy exchange source term Q TV conditioned on the local divergence values are computed. It is observed that strong compressions are associated with slightly higher values of Q TV , as observed in previous studies about thermally out-of-equilibrium mixing layers and isotropic turbulence configurations. However, this correlation is shown to be much weaker for boundary-layer configurations, since the strongest density and temperature gradients are registered in the outer portion of the boundary layer where the flow exhibits a vibrationally frozen behaviour due to the cold temperature.

The present DNS represents a first contribution towards the knowledge of wall-bounded turbulence in presence of thermochemical effects. Many aspects require future research efforts. For instance, the behaviour of the wall in terms of participation to the chemical Table 3. Constants for diffusion coefficient curve fits from [START_REF] Gupta | A review of reaction rates and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30 000[END_REF]. Of note, coefficient A 1 is identically zero.

processes is of the utmost importance: the hypothesis of non-catalytic surface should not be taken for granted, since the materials adopted for the thermal protection systems become catalytic above certain degrees of dissociation and temperature values. At the same time, ablation causes defects on the surface that could induce different transition-to-turbulence scenarios, to be considered in addition to the effect of wall roughness. Moreover, the different combination of operating conditions, such as Mach number, stagnation enthalpy and altitude, may trigger one phenomenon more intensely than another. Turbulent regimes are more likely to occur at lower altitudes, whereas thermal-relaxation phenomena would be predominant higher in the atmosphere. At the same time, higher densities and kinetic energy contents promote molecular collisions, which may result in more intense chemical activity. Spanning a range of control parameters as wide as possible would enable further insights on high-enthalpy hypersonic flows.

MFigure 1 .

 1 Figure 1. Sketch of the configuration under investigation.

  3a,b) are shown for the three turbulent stations in figures 4(a) and 4(b), respectively. The former scaling has been shown to work reasonably well for adiabatic boundary layers, even at

Figure 3 .

 3 Figure 3. Evolution of (a) the skin friction coefficient C f (black line) and Renard-Deck decomposition (red symbols) and (b) separate contribution of each term of (3.2).

  Figure4(a) confirms such a trend, both the linear and logarithmic regions being offset with respect to the analytical laws. Of note, in figure 4 the logarithmic region is described by (1/κ) log y + + C and (1/κ) log y * + C, with κ = 0.41 and C = 5.2. The semi-local scaling of Trettel & Larsson, as expected, improves the near-wall prediction since it explicitly accounts for the stress-balance condition within the entire inner layer. A large scatter is, however, observed in the logarithmic region which shows a Re-dependence similar to the van Driest scaling. Although such a scaling works reasonably well for internal flows, it is not as good for external configurations, most likely because of an interaction with the wake region which is shown to be over-stretched in figure 4(b). Recently,[START_REF] Griffin | Velocity transformation for compressible wall-bounded turbulent flows with and without heat transfer[END_REF] proposed a new total-stress-based transformation (called hereafter Griffin-Fu-Moin scaling, u GFM ) with a constant-stress-layer assumption, which reads

Figure 4 .

 4 Figure 4. Wall-normal profiles of (a) the van Driest-transformed streamwise velocity, (b) Trettel & Larsson's transformation and (c) comparison between Trettel & Larsson transformation and total-stress-based scaling of Griffin et al. (2021) at Re θ = 6200.

Figure 5 .Figure 6 .

 56 Figure 5. Reynolds stresses in the (a) streamwise, (b) wall-normal and (c) spanwise directions, and (d) Reynolds shear stress. The M14Tw018 and M8Tw048 profiles are extracted from Zhang et al. (2018), whereas M8T1 is from Xu et al. (2021).

Figure 7 .Figure 8 .Figure 9 .

 789 Figure 7. Vibrational Damköhler numbers for N 2 (a-c), O 2 (d-f ) and NO (g-i). Damköhler number based on inner units (a)-(d)-(g), large-eddy turnover times (b)-(e)-(h) and flow residence time (c)-( f )-(i).

Figure 10 .Figure 11 .Figure 12 .

 101112 Figure 10. Instantaneous isocontours of T -T V in a xy-plane and second invariant of the velocity gradient tensor, Q, coloured by T -T V .

Figure 14 .

 14 Figure14. Conditional probability density function (p.d.f.) of translational and vibrational temperatures at Re θ = 6200 for y approximately equal to 10, 100 and 3000, respectively; sweeps events (a,c,e) and ejections events (b,d, f ).

Figure 16 .

 16 Figure 16. Correlation coefficients between Favre fluctuations of vibrational temperature and streamwise velocity (a) and wall-normal velocity (b) and of translational temperature and streamwise velocity (c) and wall-normal velocity (d).

Figure 17 .

 17 Figure 17. Classical GHSRA expressed in (3.13) (a) and GHSRA based on T V , reported in (3.14) (b); -•-•-•• SRA estimation.

Figure 18 .

 18 Figure 18. Wall-normal profiles of (a) Q TV [J/(m 3 s)] and (b) Q ITV , defined in (3.15). In the inset of (a), a zoom of the profile is shown in outer scale; the dashed line denotes the null value.

Figure 19 .

 19 Figure 19. Mean mass fraction of molecular nitrogen (a), molecular oxygen (b), nitric oxide (c), atomic oxygen (d) and atomic nitrogen (e).

Figure 20 .

 20 Figure 20. Profiles of turbulent Mach number (a) and r.m.s. Mach number (b).

Figure 21 .Figure 22 .Figure 23 .

 212223 Figure 21. Instantaneous snapshot of the normalized density gradient magnitude along a slice of the computational domain.

Figure 24 .

 24 Figure 24. Conditional p.d.f. of normalized translational-vibrational energy exchange at Re θ = 6200 for y ≈ 10 (a) and y ≈ 1000 (b). The p.d.f.s are conditioned on strong events, defined by |ϑ/ϑ rms | ≥ 2.

  ). The p.d.f.s are conditioned on strong events, defined by |ϑ/ϑ rms | ≥ 2.relatively low; as a consequence, compressibility/thermal non-equilibrium interactions are expected to be less marked for boundary-layer configurations. For this purpose, we show in figure24the p.d.f.s of Q TV /Q TV,rms at the last turbulent station for the same wall-normal positions shown in figure23(b). The p.d.f.s are conditioned on strong dilatation motions, defined as |ϑ/ϑ rms | ≥ 2. The p.d.f.s in figures 24(a) and 24(b) peak at positive (negative) values of Q TV /Q TV,rms , coherently with the local mean under-excitation (over-excitation) state depicted in figure

Table 2 .

 2 Coefficients to compute pure species viscosity from[START_REF] Blottner | Chemically reacting viscous flow program for multi-component gas mixtures[END_REF].

		A n		B n	C n
	N 2	0.0268142	0.3177838	-11.3155513
	O 2	0.0449290	-0.0826158	-9.2019475
	NO	0.0436378	-0.0335511	-9.5767430
	O	0.0203144	0.4294404	-11.6021403
	N	0 . 0 1 1 5 5 7 2	0 .6031679	-12.4327495
			A 2		A 3	A 4
	N 2 -N 2		0.0112	1.6182	-11.3091
	N 2 -O 2		0.0465	0.9271	-8.1137
	N 2 -NO		0.0291	1.2676	-9.6878
	N 2 -O		0 .0140	1.5824	-10.8891
	N 2 -N		0 .0195	1.4880	-10.3654
	O 2 -O 2		0.0410	1.0023	-8.3597
	O 2 -NO		0.0438	0.9647	-8.2380
	O 2 -O		0 .0226	1.3700	-9.6631
	O 2 -N		0 .0169	1.4848	-10.2810
	NO -NO	0.0364	1.1176	-8.9695
	NO -O		0 .0179	1.4848	-10.3155
	NO -N		0 .0185	1.4882	-10.3301
	N -N		0 .0033	1.5572	-11.1616
	N -O		-0.0048	1.9195	-11.9261
	O -O		0 .0034	1.5572	-11.1729

Table 4 .

 4 Parameters needed for the five-species chemical model of[START_REF] Park | Nonequilibrium Hypersonic Aerothermodynamics[END_REF]. Pre-exponential C,

																			• for
	C a T act A 1 A 2 A 3 A 4 A 5	× 10 21 -1.6 113 200 1.4766 1.6291 1.2153 -11.4570 -0.009444	× 10 21 -1.6 113 200 1.4766 1.6291 1.2153 -11.4570 -0.009444	× 10 21 -1.6 113 200 1.4766 1.6291 1.2153 -11.4570 -0.009444	× 10 22 -1.6 113 200 1.4766 1.6291 1.2153 -11.4570 -0.009444	× 10 22 -1.6 113 200 1.4766 1.6291 1.2153 -11.4570 -0.009444	× 10 21 -1.5 59 500 0.50989 2.4773 1.7132 -6.5441 0.029591	× 10 21 -1.5 59 500 0.50989 2.4773 1.7132 -6.5441 0.029591	× 10 21 -1.5 59 500 0.50989 2.4773 1.7132 -6.5441 0.029591	× 10 22 -1.5 59 500 0.50989 2.4773 1.7132 -6.5441 0.029591	× 10 22 -1.5 59 500 0.50989 2.4773 1.7132 -6.5441 0.029591	× 10 15 0 75 500 0.50765 0.73575 0.48042 -7.4979 -0.016247	× 10 15 0 75 500 0.50765 0.73575 0.48042 -7.4979 -0.016247	× 10 15 0 75 500 0.50765 0.73575 0.48042 -7.4979 -0.016247	× 10 17 0 75 500 0.50765 0.73575 0.48042 -7.4979 -0.016247	× 10 17 0 75 500 0.50765 0.73575 0.48042 -7.4979 -0.016247	× 10 17 -1.0 38 400 0.96921 0.89329 0.73531 -3.9596 0.006818	× 10 12 0 19 400 -0.95365 -0.04585 -0.002428 -1.7415 -1.2331	act and coefficients A activation temperature T exponent a,	computing the equilibrium constant K eq .
		⇐⇒ 2N + N 2 7.0	⇐⇒ 2N + O 2 7.0	⇐⇒ 2N + NO 7.0	⇐⇒ 2N + O 3.0	⇐⇒ 2N + N 3.0	⇐⇒ 2O + N 2 2.0	⇐⇒ 2O + O 2 2.0	⇐⇒ 2O + NO 2.0	⇐⇒ 2O + O 1.0	⇐⇒ 2O + N 1.0	⇐⇒ N + O + N 2 5.0	⇐⇒ N + O + O 2 5.0	⇐⇒ N + O + NO 5.0	⇐⇒ N + O + O 1.0	⇐⇒ N + O + N 1.0	⇐⇒ NO + N 6.4	⇐⇒ N + O 2 8.4
	Reaction	N 2 + N 2	N 2 + O 2	N 2 + NO	N 2 + O	N 2 + N	O 2 + N 2	O 2 + O 2	O 2 + NO	O 2 + O	O 2 + N	NO + N 2	NO + O 2	NO + NO	NO + O	NO + N	N 2 + O	NO + O

Table 5 .

 5 Coefficients for computing the relaxation times t m from[START_REF] Park | Review of chemical-kinetic problems of future NASA missions. I-Earth entries[END_REF].

Acknowledgements. This work was granted access to the HPC resources of IDRIS and TGCC under the allocation A0092B10947 made by GENCI (Grand Equipement National de Calcul Intensif). We also acknowledge CINECA for awarding access to the Galileo supercomputer under the allocation HP10CRTQH1. Funding. D. Passiatore and G. Pascazio were partially supported by the Italian Ministry of Education, University and Research under the Program Department of Excellence Legge 232/2016 (Grant No. CUP -D94I18000260001).

Declaration of interests. The authors report no conflict of interest.

Author ORCIDs.

D. Passiatore https://orcid.org/0000-0002-6255-3079; L. Sciacovelli https://orcid.org/0000-0002-2463-4193; P. Cinnella https://orcid.org/0000-0001-9979-0783; G. Pascazio https://orcid.org/0000-0002-6684-4113.

Appendix A. Thermochemical model and transport properties

Pure species' viscosity and thermal conductivities are computed using curve fits by [START_REF] Blottner | Chemically reacting viscous flow program for multi-component gas mixtures[END_REF] and Eucken's relations [START_REF] Hirschfelder | Simulation and modeling of cold-wall hypersonic turbulent boundary layers on flat plate[END_REF])

where A n , B n and C n are fitted parameters reported in table 2. The corresponding mixture properties are evaluated by means of the mixing rules of Wilke ( 1950)

with f being either μ, λ TR or λ V , and

As for mass diffusion processes, the equivalent coefficient D n of (2.12) is computed as in [START_REF] Hirschfelder | Simulation and modeling of cold-wall hypersonic turbulent boundary layers on flat plate[END_REF])

, withD mn = 1 p exp (A 4,mn )T A 1,mn (ln T) 2 +A 2,mn ln T+A 3,mn , (A5)

where D mn is the binary diffusion coefficient of species m into species n, and A 1,mn , . . . , A 4,mn are curve-fitted coefficients computed as in [START_REF] Gupta | A review of reaction rates and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30 000[END_REF]; the values are reported in table 3.

Concerning the chemical model of [START_REF] Park | Nonequilibrium Hypersonic Aerothermodynamics[END_REF], the forward and backward reaction rates k f ,r and k b,r shown in (2.14) are modelled by means of Arrhenius' law

Here, T act,r is the activation temperature of reaction r, whose values are reported in table 4, together with the values of the pre-exponential parameter C r and of the exponent a r . In (A6a,b), K eq,r is the equilibrium constant of the rth reaction, for which Park employed a curve fit of the form

with Z = 10 000/T avg , and the constants A •r given in table 4. Lastly, the coefficients a mn and b mn , used to compute the relaxation time in (2.16), are reported in table 5.