Stéphane Lagarde

Anthony Meilland

Bruno Lopez

Romain G Petrov

Fatmé Allouche

Alexis Matter

Jean-Charles Augereau

Phillippe Berio

Felix Bettonvil

Roy Van Boekel

Philippe Berio

Paul Bristow

Pierre Cruzalèbes

William C Danchi

Andreas Glindemann

Massinissa Hadjara

Matthias Heininger

Michiel Hogerheijde

Walter Jaffe

James Leftley

Michael Lehmitz

Florentin Millour

Claudia Paladini

Eric Pantin

Miguel Riquelme

Thomas Rivinius

Sylvie Robbe-Dubois

Markus Schoeller

Jozsef Varga

Gerd Weigelt

Julien Woillez

Gérard Zins

MATISSE, the VLTI mid-infrared imaging spectro-interferometer at the time of the first published astrophysical results

Keywords: Interferometry, mid-infrared, astrophysical results, VLTI

INTRODUCTION

The Multi AperTure mid-Infrared SpectroScopic Experiment, MATISSE, is the mid-infrared spectrograph and imager of the VLTI. This second generation interferometric instrument was designed with prime science objectives focusing on the Young Stellar Objects (YSOs) and on their associated stellar and planet formation processes, on the Active Galactic Nuclei (AGNs) and their active central engine, on the Stellar Physics, the surface structure and mass loss of stars at different evolutionary stages and in frequent binary interactions.

In particular, MATISSE is perfectly designed to observe the small spatial scales in YSO environments at wavelengths at which the dust material that represents the elementary bricks initiating the growth of the planet core embryos, emits efficiently. MATISSE will also be of importance to probe the inner dust structure in AGNs by testing the validity of the AGN unified model. That model is supposedly edge-on and is hiding the central source in type-II AGN and it is face-on and dominated by the inner dust sublimation radius in type-I AGN. Clumpy structures with important polar extensions are expected to trace the dust driven wind in the brightest AGNs.

MATISSE gives access to the L (2.8-4.2) , M (4.5-5.0 m) and N bands (8.0-13.0 m) band, various spectral resolutions between R 30 and R 3500, which allow to study very specific gas and dust material spectral signatures (see Table 1). MATISSE combines four beams from the VLT's unit or auxiliary telescopes and can thus provide mid-infrared images (closure-phase aperture synthesis imaging) down to an angular resolution of about 3 mas in L band. Images are now routinely produced. Even in the case where image reconstruction are not performed through aperture synthesis and a good uv coverage approach, the access to closure phases is of considerable value since they bring new unexplored constraints in this specific wavelength domain.

Through the MATISSE project, ESO and our consortium are contributing to the new generation of midinfrared instrumentation becoming available to the astronomical community. MATISSE is an essential piece completing a worldwide ensemble of astronomical instruments. From spectral considerations and regarding the near-infrared domain (NIR), MATISSE extends the operating wavelength range of existing instruments such as PIONIER [START_REF] Zins | PIONIER: A Four-telescope Instrument for the VLTI[END_REF] and GRAVITY, 2 which are both also sensitive to IR radiation but at a shorter wavelength. Moreover MATISSE complements the (sub)millimetre domain in which high angular resolution observations are covered by the radio interferometer ALMA. [START_REF] Carpenter | The ALMA Development Program: Roadmap to 2030[END_REF] MATISSE gives access to long baseline lengths compared to the midinfrared interferometric instrument of the Large Binocular Telescope. [START_REF] Sallum | ELT Imaging of MWC 297 from the 23 m LBTI: Complex Disk Structure and a Companion Candidate[END_REF] Note that MATISSE is thus a bridge in wavelengths as well as an exploratory path for high angular resolution in its spectral bands. That is of essential importance with the advent of the JWST 5 and the installation of METIS, the ELT mid-infrared instrument, [START_REF] Brandl | The Mid-IR ELT Imager and Spectrograph (METIS) and its Science Goals in the Context of AKARI[END_REF] both of which will operate in the same wavelength domain as MATISSE. In terms of sensitivity, these instruments will have advantages over MATISSE. JWST will benefit from the absence of the atmosphere and low background observing conditions due to its location in space, while METIS will benefit from the large collecting area of the ELT. However, from an imaging point of view, MATISSE has a resolving power 20× and 4-5× greater than that of the JWST and METIS respectively. The MATISSE instrument (Fig. 1) consists of a warm optical system and two cold optical benches (COBs) housed in two separate cryostats, one for the L-M band, and one for the N band. The WOP operates at an ambient temperature. It is fed by four beams from the VLTI, which are propagated via dichroic beam splitters to the two COBs (Fig. 2).

OVERVIEW OF THE MAIN INSTRUMENT CHARACTERISTICS

The MATISSE instrument contains several specific functions of importance : -the Beam Commuting Device(s): MATISSE features a specific module called Beam Commuting Device (BCD), to calibrate closure and differential phases. Two BCDs are installed at the entrance of MATISSE, one for each pair of beams. By commuting beams 1 and 2, and/or 3 and 4, we invert all the phase contributions before the instrument BCD (i.e., including the contributions of the astrophysical source, the atmosphere, and the VLTI beam trains) while keeping the instrumental contributions after the BCD device fixed. By commuting the beams, 2 differential phases and the 4 closure phases of the astrophysical object can be obtained independently of the phase produced by the instrument but not by the VLTI beam train or the atmosphere. On the closure phase, the BCDs allow a reduction of the phase residuals from several degrees to less than one degree.

-the OPD modulator(s): The thermal background contribution is huge in the mid-infrared and dominates the astrophysical signal by several orders of magnitude in N-band. In principle, it only affects the low-frequency peak of the Fourier Transform of the interferogram. However, the finite detector window and spatial variations of the detector gain can induce a crosstalk between the different peaks. Therefore, the high-frequency peaks (i.e. the fringe peaks) may be contaminated by the thermal background emission. To reduce such a contamination , a temporal modulation is performed. This modulation consists in applying an artificial Optical Path Difference (OPD) to each fringe signal with a periodic sequence during the atmospheric coherence time. It is a generalisation of the natural 0-π scheme existing with a co-axial combination. This function is provided by a set of piezo actuators. The modulation is made in ten steps of λ /10 during the coherence time in which the background remains almost constant. The resulting 10 signals are then numerically re-phased and summed. The lower frequency peak, which is not affected by the OPD modulation, and its contamination to the adjacent fringe peaks is then eliminated and only the high frequency peaks remain.

-the beam shaper and combiner: The MATISSE beam combination is multi-axial. All the telescope beams are mixed together to produce a common interference pattern. A non-redundant beam configuration arrangement with separations between beams, respectively equal to 3D, 9D, and 6D (where D is the beam diameter), allows the fringe peaks to be separated in the Fourier space. It produces dispersed fringes on two different detectors simultaneously, the HAWAII-2RG from Teledyne Technologies for the L&M bands, and the AQUARIUS from Raytheon Company for the N band.

In the L band, the detector integration time is driven by the seeing coherence time or by the fringe tracking time (if an external fringe tracker is used). In the N band, it is driven by the high thermal background level, which imposes the individual detector integration time to be short enough to avoid saturation. Indeed the maximum exposure time in N is less than 30 ms while in L, several seconds of integration time are allowed.

The instrument has been available in open time since April 2019. Its integration at Paranal took place from November 2017 to February 2018. A two-year commissioning period followed and allowed for the instrument performance to be assessed in all modes with ATs and UTs, and within an upgraded VLTI infrastructure [START_REF] Woillez | VLTI: First Light for the Second Generation[END_REF] involving the implementation of the ATs' adaptive optics, NAOMI, [START_REF] Woillez | NAOMI: the adaptive optics system of the Auxiliary Telescopes of the VLTI[END_REF] and of the fringe tracker mode, GRA4MAT (Woillez et al. in preparation).

MAIN PERFORMANCES OF MATISSE WITH RECENT PERFORMANCE ASSESSMENTS RESULTING FROM 2022 COMMISSIONING SESSION

MATISSE is a highly sensitive spectro-interferometer 9 which benefits from the large telescope apertures of the VLT. The ATs and the UTs are large-sized telescopes when compared to other worldwide interferometric array. This is important for mid-infrared observations which requires strong source signal at the detector level to beat the photons and background fluctuation noise generated in this spectral domain by both the atmosphere and the optical train emission.

The different contributions to the errors, which affect the measurements, are:

• The fundamental noise errors that affect the precision of all measurements.

• The broad-band seeing errors that affect the absolute visibility and closure phase independently from the source brightness.

• The broad-band photometric errors likely caused by imperfect thermal background subtraction, which produce a source flux dependent error in the absolute visibilities.

• The group delay estimation error that introduces a bias on the coherent flux, in the case of coherent integration, for faint targets.

• The chromatic OPD that has a strong impact on the differential phase Tables 2 and 3 combine all these error sources and shows our updated estimates without fringe tracking recently reassessed considering the recent and final results of the commissioning. This limiting magnitudes correspond to the ones offered by ESO to the community of users. They are set by the following precision criteria, to be achieved per spectral channel after one minute of observation:

• Visibility: σ V = 0.1 • Closure phase: σ ψ = 5 • • Differential phase: σ φ = 4 •
Table 2: MATISSE standalone (i.e. without Fringe Tracker) sensitivity limits in Jy in L and M, as defined for the following precision criteria: σ V = 0.1, σ ψ = 5 • and σ φ = 4 • , to be achieved per spectral channel after one minute of observation. The symbol "-" means that the combination resolution-versus-band does not exist. Magnitude 0 in L (at 3.5µm) and M (at 4.75µm) corresponds to 290Jy and 165Jy respectively.

Telescopes Resolution V ψ φ L M L M L M ATs LOW

A SELECTION OF RECENTLY PUBLISHED RESULTS

MATISSE was designed as a successor to MIDI, a prolific instrument which produced numerous results on protoplanetary disks, AGNs, and evolved stars. Beyond its huge success, MIDI was nevertheless limited to two telescopes. Going from a two-beam to a four-beam combiner increases the number of simultaneous visibility measurements by a factor 6 and gives access to closure phase measurements. Moreover, MATISSE also opens two new spectral windows at VLTI, the L and M bands, filling the gap of the first generation VLTI instruments. Its very extended spectral domain allows to probe dust, from the sublimation temperature down to about 200K, and free-free emission from hot gas. Moreover, by using its various spectral resolutions, solid-state features, molecular bands and atomic lines can be spectrally resolved, allowing to spatially disentangle gas and dust emissions, constrain their physics, infer the dust mineralogy and even probe the gas kinematics.

In Fig. 3 we present three of the first reconstructed images from MATISSE: the circumstellar disk around the unclassified B[e] star FS CMa, 10 the dusty torus around the AGN NGC 1068, [START_REF] Rosas | Thermal imaging of dust hiding the black hole in NGC 1068[END_REF] and the spectro-imaging in the Brα emission line of the wind collision zone around the massive binary star η Car. [START_REF] Weigelt | VLTI-MATISSE chromatic aperture-synthesis imaging of η Carinae's stellar wind across the Brα line. Periastron passage observations in February 2020[END_REF] These images not only demonstrate the capability of MATISSE to image environment at a few milli-arcseconds resolution, but also its ability to produce spectrally-resolved images either in a large spectral domain (from 3 to 13 µm) in the case of FS CMa and NGC 1068 or at a resolution high enough to study kinematics in spectral lines in the case of η Car (R≳1000). [START_REF] Leftley | New Evidence for the Dusty Wind Model: Polar Dust and a Hot Core in the Type-1 Seyfert ESO 323-G77[END_REF] . The 2-telescope MIDI instrument was unable to make images and could only give angular sizes as a function of the position angle. It was limited to the N band where it is difficult to discriminate the effects of the temperature and composition of the emitting dust from the foreground absorbing medium. MATISSE solves these two problems by making images in the N band where dust emission dominates. The inner dust structure in AGN is the tracer of the exchanges between the central engine, its surrounding broad line region (BLR), and the host galaxy. Observing it down to the pc scale can reveal the mechanisms of co-evolution of the super massive black holes and galaxies. Moreover, understanding the geometry of the nuclei well can allow independent distance estimates in combination with infrared (IR) reverberation mapping [START_REF] Hönig | A dust-parallax distance of 19 megaparsecs to the supermassive black hole in NGC 4151[END_REF] . Finally, the recent breakthrough made by GRAVITY on the spectro-astrometry of BLRs 16 makes this MATISSE information of an even higher value.

In the case of YSOs and protoplanetary disks, reaching the mas scale in the mid-infrared allows to probe au or even sub-au scale for the closest star-forming regions (∼ 150 pc). Notably, accessing that spatial scale is required to resolve the water ice line in disks, which is thought to represent the outer boundary of the telluric planet forming region. In addition, mapping disks at such a spatial scale could allow to detect au-scale gaps or radial-inhomogeneties that could be linked to planet formation. A survey is being conducted by the MATISSE consortium and more than 60 young stars have been observed as of July 2022. one example of star observed in the frame of this survey is the Herbig Ae star HD 163296 17 . MATISSE obser-vations revealed a significant brightness asymmetry in the L-band emission at a very small scale, i.e., 0.33±0.01 au. The authors proposed that this asymmetric structure, located in or near the inner rim of the dusty disk, orbits the star. To find the physical origin of the asymmetry, they tested the hypothesis of a vortex created by a Rossby wave instability. They found that a unique large-scale vortex is compatible with the MATISSE data. Such vortices are typically proposed as a mechanism to trap dust particles and facilitate planetesimal growth in proto-planetary discs. Another important matter regarding YSOs is to trace carbon in their protoplanetary disks. While detected in the outer disk regions, little is known about the solid carbon species in the inner disk regions and key questions are whether there is a significant reservoir of C-based solid species in the inner 10 au region of disks, in what form are those C-based species and how can solid carbon be incorporated into planetesimals and ultimately planets. Thanks to its large wavelength coverage (3-13µm) and using its medium spectral resolution mode (R>500) MATISSE can study various carbonaceous molecules and nanograins (PAHs, nano-diamonds, CO) down to the sub-au scale. In [START_REF] Kokoulina | First MATISSE L-band observations of HD 179218. Is the inner 10 au region rich in carbon dust particles?[END_REF] , MATISSE observations allowed to detect the presence of two separated carbonaceous dust populations in the disk around the Herbig star HD 179218: very hot nano-grains that dominate the near-to mid-IR emission in the inner optically thin 10 au region and colder and larger grains that dominate the emission outward. The survival of nano-grains in the highly irradiated inner 10 au region brings new insights in the different photo-induced processes acting on such small grains, and their regeneration/replenishment: accretion could maintain a significant inflow of material to maintain a stable nano-grain population there as well as local photo-fragmentation of larger carbonaceous grains.

In parallel to the two main drivers that are YSOs and AGNs, MATISSE has started to produce many results in the field of evolved stars. MATISSE can image dusty and gaseous circumstellar environments around hot stars such as the ones presented in 3 but it can also probe the surface of giant and supergiants stars as shown in [START_REF] Chiavassa | The extended atmosphere and circumstellar environment of the cool evolved star VX Sagittarii as seen by MATISSE[END_REF] in the case of the enigmatic evolved star VX Sgr. In that case MATISSE observations allowed to reconstruct multi-band images of the stellar surface and compare them to outputs from radiative transfer models of both RSG and AGB star, showing that this object belong most likely to the later class of objects. Another example is the case of the AGB star R Scl. By modeling its extended atmosphere, 20 managed to disentangle the emission from different chemical species (acetylene, Hydrogen cyanide, and amorphous carbon and silicon carbide dust). They showed that the molecules are located closer to the star than the dust, as shown in Fig 4, which presents the spectro-radial profiles of R Scl derived from the MATISSE data. Such pioneering work will help to improve our understanding of the chemical process aroung AGBs star and dust formation.

FORESEEN VLTI UPGRADES FOR FRONTLINE SCIENCE

The MATISSE instrument is now routinely operating at Paranal for science observations and will soon become officially property of ESO. Indeed, the PAC (Provisional Acceptance in Chile) phase, which represents the final delivery step of the Consortium, is being concluded. Beyond the PAC, the already engaged commissioning activity on GRA4MAT will be pursued during 2022 and future upgrades are foreseen. They concern the MATISSE-Wide mode, which is the use of GRAVITY as a wide field-of-view off-axis fringe tracker for MATISSE, and the proposition of a new generation fringe tracking of the HFT-type (Petrov et al. 2022, this conference). The MATISSE-wide upgrade consists in using the so-called beam B for MATISSE and IRIS while the set of beam A are used to feed GRA4MAT (see Figure 5). The proposition of a new generation fringe tracking of the HFT-type consists in developing a dedicated fringe tracker for MATISSE and the VLTI. This was foreseen years ago by ESO at the beginning of the MATISSE Phase A. Such a new development would be to take advantage of the state-of-Art knowledge in this field by using in particular several spectral bands in once to fringe sense with the highest performance the signal (Petrov et al. 2022, this conference).

These upgrades are science driven by a major case which is the access to a much larger number of AGNs and Quasars 21 through a multispectral approach at the VLTI. The mid-infrared domain covered by MATISSE offers crucial constraints to determine temperature maps and morphologies, and will contribute in particular by this mean to better constraints on the distance estimates through reverberation mapping methods. It is the Through the MATISSE-wide upgrade, it is not only MATISSE perfomance which will be improved for faint object observations but it is globally the VLTI multi-spectral capability which will be increased to an homogeneous level allowing the community to continue to use the ensemble of instruments PIONIER, GRAVITY and MATISSE for the multi-spectral analysis of the individual sources.

  MATISSE operates in three atmospheric transmission windows, the L, M, and N bands of the mid-infrared. It combines four beams coming from either the 4 UTs or the 4 ATs of the VLTI. The offered spectral bands and resolutions are presented in Tab. 1.

Figure 1 :

 1 Figure 1: Top: Layout of MATISSE inside the VLTI focal laboratory. In front of MATISSE, the feeding optics table reflects towards MATISSE either four UT or AT beams, after their transportation through the VLTI optical train. Bottom: Location of MATISSE inside the VLTI laboratory with the Warm OPtics table called the WOP (and its cover) in the foreground.

Figure 2 :

 2 Figure 2: Sketch of the instrument with its main functions.

Figure 3 :

 3 Figure 3: A few of the first Images reconstructed from MATISSE observations. Top left: the circumstellar Disk around the unclassified B[e] star FS CMa. 10 Top Right: NGC1068 dusty torus. 11 Bottom : spectro-imaging in the Brα emission line of the wind collision zone around the massive binary star η Car.12

Figure 4 :

 4 Figure 4: R Scl spectro-radial maps in L-and-N-bands 20 obtained by plotting the best intensity Hankel profile normalized at one for each observed wavelength. The modelling reveals the location of different species detected in the upper atmosphere and various shells around the star: hydrogen cyanide (HCN), acetylene (C 2 H 2 ), and dust

Figure 5 :

 5 Figure 5: Sketch of the MATISSE feeding optics setup allowing the off-axis use of GRA4MAT.

Table 1 :

 1 Spectral resolutions, calculated with a pinhole of 1.5λ/D in L and M band, and 2λ/D in N band.

	Spectral Mode Spectral bandwidth (µm) Resolution
	LOW-N	8.0-13.0	31.5
	HIGH-N	8.0-13.0	218
	LOW-LM	2.8-4.2; 4.5-5.0	31.5
	MED-LM	2.8-4.2; 4.5-5.0	499
	HIGH-L	2.8-4.2	979
	VHIGH-L	3.95-4.2	3370
	VHIGH-M	4.5-5.0	3370

Table 3 :

 3 MATISSE standalone (i.e. without Fringe Tracker) sensitivity limits in Jy in N1 (8.5µm), N2(10.5µm), N3(11.5µm) as defined for the following precision criteria: σ V = 0.1, σ ψ = 5 • and σ φ = 4 • , to be achieved per spectral channel after one minute of observation. Magnitude 0 in N (at 8.5µm) corresponds to 50Jy.

				1.1		2.1	0.4		1.9	0.3	1.1
			MED	3.8 13.2 3.3 11.5 2.4	8.4
			HIGH	20.1	-	14.7	-	10.7	-
		UTs	LOW	0.3		0.4 0.07 0.15 0.06 0.15
			MED	1.1		1.1	0.8		0.9	0.6	0.7
			HIGH	2.4		-	1.7		-	1.2	-
	Telescopes Resolution		V				ψ		φ
			N1	N2	N3		N1	N2		N3	N1	N2	N3
	ATs	LOW	16.8 31.4 44.0 9.4 23.6 36.7 2.9	7.5 11.7
		HIGH	30.3 51.3 83.7 29.9			25.3 45.9 58.9
	UTs	LOW	0.9	1.4	1.9		0.3	0.8		1.3	0.2 0.55 0.6
		HIGH	1.6	3.4	4.1		1.5			1.1	3.5	3.7

* * MATISSE was designed, funded, and built in close collaboration with ESO by a consortium composed of institutes in France (J.-L.


Jolan Reynaud
email: jolan.reynaud@onera.fr

Pierre-Elie Weiss
email: pierre-elie.weiss@onera.fr

Sébastien Deck
email: sebastien.deck@onera.fr

Sbastien Deck

On the use of bi-species ZDES for multi-nozzle space launcher configurations

Keywords: multi-nozzle afterbodies, supersonic jets, turbulence modelling, bi-species flow, RANS/LES, ZDES

Space launchers experience high fluctuating pressure levels during the ascent (for both expendable and reusable launch vehicles) and return phases (for reusable launch vehicles only). To simulate the fluctuating pressure field occurring on such configurations accurately, a numerical workflow combining ZDES Mode 2 (2020) and a hybrid scheme ensuring robustness in shock wave regions and low dissipation levels in vortical regions is used in the framework of bi-species inert flows. The assessment of the performance of this numerical strategy is based on the simulation of a four-nozzle launcher model previously studied experimentally by Musial and Ward [1]. A comparison of pressure coefficients shows that ZDES gives an improvement in the capability of predicting the base pressure over standard RANS models. Spectral analysis of the fluctuating pressure at the wall shows that the flow is dominated by the antisymmetric mode m = 1 contributing to up to 80% of the energy at a dimensionless frequency StD = 0.2.

Introduction

CFD is playing an expanding role in the aerospace industry and enhancing its capabilities to simulate complex unsteady turbulent flows is considered as a major objective in the field [START_REF] Sagaut | Large Eddy Simulation for Aerodynamics: Status and perspectives[END_REF][START_REF] Deck | High-fidelity simulations of unsteady civil aircraft aerodynamics: stakes and perspectives. application of zonal detached eddy simulation[END_REF]. As an example, to support the development of multi-nozzles for reusable launcher designs, the accurate prediction of afterbody flows occurring On the use of bi-species ZDES for multi-nozzle space launcher configurations behind the base and around the nozzles during flight is of primary importance to assess the aerodynamic performance and the integrity of the vehicles. Indeed, the reactive multi-gas interactions occurring between the hot propulsive jets and the base flow can induce adverse phenomena such as base drag, due to the low pressure recirculation zone forming behind the base. The pressure loads in this area can damage mechanical structures. Moreover, base heating attributable to the recirculation of hot gases towards the base may arise. As a result of this complexity, the aerodynamic design of new generation launcher afterbodies remains a challenge for both experimental and numerical tools.

As far as numerical approaches are concerned, two main challenges are classically acknowledged in the literature. First, the need to resolve the large-scale turbulent eddies driving the dynamics of base flows in order to get accurate predictions of the wake flow topology and wall pressure levels has been well documented for generic space launcher configurations [START_REF] Reynaud | Numerical workflow for scale-resolving computations of space launcher afterbody flows with and without jets[END_REF] [START_REF] Reynaud | Numerical workflow for scale-resolving computations of space launcher afterbody flows with and without jets[END_REF] and Statnikov et al. 2016 [5]). RANS approaches have indeed shown limited capacities to reproduce such massively separated flows whereas RANS/LES approaches, combined with low dissipative numerical schemes, have displayed a good cost/accuracy compromise [START_REF] Pain | Zonal Detached Eddy Simulation of theFlow Around a Simplified Launcher Afterbody[END_REF] [START_REF] Pain | Zonal Detached Eddy Simulation of theFlow Around a Simplified Launcher Afterbody[END_REF], Weiss and Deck [START_REF] Weiss | On the coupling of a zonal body-fitted/immersed boundary method with ZDES: Application to the interactions on a realistic space launcher afterbody flow[END_REF]). Then, authors such as Gusman et al. 2011 [START_REF] Gusman | Best practices for cfd simulations of launch vehicle ascent with plumes-overflow perspective[END_REF] have advocated the need to take into account the thermodynamic properties of the hot reactive propulsive jets to increase the representativeness of space launcher simulations. As the computational cost of both RANS/LES approaches and multi-species reactive flow models is important, authors often have to establish a compromise between an accurate treatment of turbulence and the representation of chemistry effects.

Regarding wind tunnel experiments, only few studies have been dedicated to multi-nozzles space launcher configurations with hot propulsive jets (Mehta et al. 2013 [9]). Among them, the case of a four-JP4/LOX-nozzle space launcher configuration studied experimentally in Musial and Ward 1961 [START_REF] Musial | Base flow characteristics for several four-clustered rocket configurations at Mach numbers from 2[END_REF] is considered as a reference case and has been studied numerically by several authors. Among these studies, Mehta et al. 2013 [START_REF] Mehta | Base heating sensitivity study for a 4-cluster rocket motor configuration in supersonic freestream[END_REF] and Pu and Jiang 2019 [START_REF] Pu | Assessing turbulence models on the simulation of launch vehicle base heating[END_REF] employed advanced thermodynamic models but only used RANS models for turbulence modelling. Though these studies have provided useful analyses of the flow field, noticeable errors in the evaluation of wall pressure levels have been reported as well.

To contribute to the study of launcher base flow prediction methods, this article presents a numerical study of a four-nozzle space launcher configuration using a scale resolving framework based on the use of a simple two-species physical model and hybrid numerical methods in order to be able to do scale-resolving simulations at an affordable cost. The hybrid model used is the ZDES mode 2 (2020) model by Deck and Renard 2020 [START_REF] Deck | Towards an enhanced protection of attached boundary layers in hybrid RANS/LES methods[END_REF]). The results will be compared with a two-species RANS computation and with RANS computations including chemistry effects from the literature in order to investigate the benefits of a scale resolving approach for multi-nozzle afterbody computations.

In the following, the test case as well as the salient features of the numerical framework used for the bi-species ZDES and RANS computations are first briefly described before going into the instantaneous and Reynolds averaged properties of the flow. In particular, the focus is put on the resulting wall pressure distribution on the base and on the nozzles of the launcher. Finally, the salient features of the fluctuating pressure field are investigated with second-order statistics and spectral analyses.

Test case and mesh setup

The present study is based on a wind tunnel experiment provided by Musial and Ward [START_REF] Musial | Base flow characteristics for several four-clustered rocket configurations at Mach numbers from 2[END_REF], the geometry of the model for this study is shown in Fig. 1. It consists of a four-nozzle subscale rocket model defined by the following parameters r b = 152.4 mm, L n /D e = 1.53, and D s /D e = 1.67. The four rocket engines operate with a propulsive mixture combining JP-4 as a fuel and LOX as an oxidizer. The selected flow conditions for pressure correspond to a flight altitude of 14.9 km and are summarized in Table 1 together with the total pressure P c and total temperature T c inside the combustion chamber. These test conditions correspond to underexpanded jets with an exit Jet Pressure Ratio (JPR=3.4) (Mehta et al. 2013 [9]). The wind tunnel model was instrumented with pressure transducers and temperature probes on its base.

As shown in Fig. 2b a structured mesh has been designed using a secondary O-H topology for each of the nozzles and a main O-H topology built around the launcher body. The radial point distribution shown in Fig. 2a is set such that ∆y + ≈ 1 is ensured for every attached boundary layer with exponential expansion ratios lower than 1.05 and to cluster points in the mixing layer forming behind the main body base. The computational domain extends over 5 times the length L c of the main body (1L c upstream and 3L c downstream) and has a global radius that is equal to 5 times the radius of the base. The domain ends with a sponge region and non-reflective boundary conditions to avoid any spurious numerical waves. With 240 points in the azimuthal direction (∆θ = 1.5°) the mesh contains 55.6×10 3 Computational set-up

Two-species model

To study the 4-nozzle launcher configuration, an implicit finite volume formulation of the two-species Navier-Stokes equations, presented in Reynaud et al. 2022 [START_REF] Reynaud | A comprehensive framework for high fidelity computations of two-species compressible turbulent flows[END_REF] is used. This approach enables to define the propulsive gas mixture and the surrounding air as two perfect gases. In the following, all relevant quantities for the jet gas are indicated by the subscript j. To provide adequate properties for the propulsive mixture, JP-4/LOX equilibrium flow data from Huff et al. 1956 [13] are used to define constant values for the Sutherland's law constants [µ 0j ,S j ,T 0j ], the specific heat ratio γ j , the volume constant heat ratio c vj , the Prandtl number P r j of the jet gas and a turbulent Prandtl number P r t (see Tab.2). The turbulent Schmidt number is chosen as 0.5 for the mixing between air and the propulsive mixture. [START_REF] Weiss | On the coupling of a zonal body-fitted/immersed boundary method with ZDES: Application to the interactions on a realistic space launcher afterbody flow[END_REF]). In these last references, the numerical results are thoroughly compared with the available experimental data of both Reynolds averaged quantities as well as pressure and velocity fluctuations including spectral and second-order analysis.

The time integration is carried out by means of the second-order-accurate backward scheme of Gear. Newton-like inner iterations (typically 4) are used to reach second-order time accuracy. A decrease of the inner-residuals of at least one order is obtained. Besides, the time step needs to be chosen in order to correctly describe the important physical phenomena (such as mixing layer and wake dynamics) which are being simulated. The Courant-Friedrichs-Levy (CFL) condition relates the maximal allowable time step to a local velocity and local mesh spacing. The acoustic CF L a number based on maximum acoustic velocity U + a (where U and a denote respectively the local velocity magnitude and speed of sound) and the grid step taken as the minimum grid extension spacing in each direction, namely:

CF L a = (U + a)∆t CF D min(∆x, ∆y, ∆z) (1)
In the present study, the time-step for the ZDES computations is set to ∆t CF D = 2.10 -7 s yielding CF L a ≤ 6 in the base flow mixing layer. The highest values are observed in the very early stages of the mixing layer due to the fine grid cells inherited from the attached boundary layer in the frame of a structured mesh. One can note that the corresponding CFL value based only on the maximum hydrodynamic velocity U is lower than one (i.e. CF L = U ∆t CF D min(∆x,∆y,∆z) ≤ 1) in the whole massively separated flow region of interest.

The choice of the spatial scheme deserves particular attention when scale resolving simulations in a supersonic framework are considered since the scheme has to challenge conflicting demands between accuracy in LES regions (i.e. low numerical dissipation) and robustness for shocks. Hence, Reynaud et al. [START_REF] Reynaud | Numerical workflow for scale-resolving computations of space launcher afterbody flows with and without jets[END_REF] developed a hybrid numerical framework based on Ducros's sensor [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF] and designed to switch from a low-dissipation formulation in presence of vortical structures to a robust formulation around high gradients. In practice, the scheme relies on a salient combination of the AUSM+ (Liou 1996 [27]) for shocks and an adapted low-dissipative AUSM+(P-LD) (Mary & Sagaut 2002 [17]) in the resolved turbulence region. The numerical dissipation of the spatial scheme is dynamically optimized based on local properties, especially for turbulent flows containing shocks. Further details concerning the numerical method and implementation of turbulence models can be found in Pechier, Guillen & Cayzac (2001) [START_REF] Pechier | Magnus effect over finned projectiles[END_REF], [START_REF] Deck | Development and application of spalart-allmaras one equation turbulence model to three-dimensional supersonic complex configurations[END_REF] [START_REF] Deck | Development and application of spalart-allmaras one equation turbulence model to three-dimensional supersonic complex configurations[END_REF] while the numerical framework dedicated to the resolution of bi-species flows is thoroughly detailed in the paper by Reynaud et al. [START_REF] Reynaud | A comprehensive framework for high fidelity computations of two-species compressible turbulent flows[END_REF].

Turbulence modelling

The ZDES approach developed by ONERA [START_REF] Deck | Zonal-Detached-Eddy Simulation of the flow around a high-lift configuration[END_REF][START_REF] Deck | Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF] belongs to the family of multiresolution approaches which covers the full range of modeling from RANS to Wall-Modelled LES (WMLES) [START_REF] Deck | A rapid and low noise switch from rans to wmles on curvilinear grids with compressible flow solvers[END_REF][START_REF] Deck | Zonal Detached Eddy Simulation (ZDES) of a spatially developing flat plate turbulent boundary layer over the reynolds number range 3 150 ≤ Re θ ≤ 14 000[END_REF][START_REF] Vaquero | Outer layer turbulence dynamics in a highreynolds-number boundary layer up to re θ ≈ 24,000 recovering from mild separation[END_REF]] and aims at treating in a simple model all classes of flow problems indicated in Fig. 3. The zonal nature of ZDES permits the use of various model formulations within the same computation (see Deck & Laraufie [START_REF] Deck | Numerical investigation of the flow dynamics past a three-element aerofoil[END_REF] for an example). The current study is solely based on ZDES mode 2 (2020) [START_REF] Deck | Towards an enhanced protection of attached boundary layers in hybrid RANS/LES methods[END_REF], a fully automatic model where the model sets dynamically by itself the RANS and LES branches (i.e. the RANS and LES zones are not set in advance by the user). In the frame of the underlying Spalart-Allmaras model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF], dmode 2 ZDES replaces the distance to the wall d w involved in the SA model according to:

dmode 2 ZDES = f ∆, d w , U i,j , ν, ν, ∂ ν ∂n , ∂ω ∂n (2)
with ∆(∆x, ∆y, ∆z, ω) denoting the subgrid length scale entering dmode 2 ZDES that depends on the grid cell size, U i,j the velocity gradients, d w the distance to the wall, ω the magnitude of vorticity, ν and ν the kinematic and pseudo eddy viscosity as well as ∂/∂n the derivative in the wall-normal direction. The definition of dmode 2 ZDES is not a minor adjustment in the DES framework since the modified lengthscales depend not only on the grid but also on the velocity and eddy viscosity fields and their wall-normal derivatives. ZDES mode 2(2020) is to the authors' best knowledge the only published formulation providing a safe RANS shielding of the whole turbulent boundary layer profile for any grid cell size (including infinite mesh refinement) and pressure gradient even when the separation moves in time (see [START_REF] Deck | Towards an enhanced protection of attached boundary layers in hybrid RANS/LES methods[END_REF] for an example of transonic buffet over a supercritical airfoil).

Figure 4 displays an instantaneous snapshot of the eddy viscosity in a plane of symmetry of the current configuration. The attached turbulent boundary layer treated in RANS mode is clearly identified. The eddy viscosity in the base region (i.e. LES region) is not destructed as fast as in ZDES mode 1 calculation (see [START_REF] Pain | Large scale dynamics of a high Reynolds number axisymmetric separating/reattaching flow[END_REF] for example), but the impact in the resolved turbulence seen in Figure 6a is marginal. Let us be reminded that within ZDES mode 2 (2020) the RANS and LES domains are not set in advance by the user (unlike as in mode 1) so that mode 2 appears as a good compromise between efficient shielding of attached boundary layers and destruction of eddy viscosity in the free shear layer. (c) Instantaneous pressure distribution and Q-criterion (black iso-lines) in the wake of the launcher in a plane between the nozzles Fig. 6: Instantaneous flow topology pressure gradient on the outer side of the jet boundaries. Looking at the underexpanded propulsive jets, we observe that in accordance with results from Mehta et al. 2013 [START_REF] Mehta | Base heating sensitivity study for a 4-cluster rocket motor configuration in supersonic freestream[END_REF] and Pu and Jiang 2019 [START_REF] Pu | Assessing turbulence models on the simulation of launch vehicle base heating[END_REF], the jet boundaries collide near the axis of the On the use of bi-species ZDES for multi-nozzle space launcher configurations (a) Distribution of longitudinal velocity and streamtraces in a plane containing two nozzles.

(b) Distribution of longitudinal velocity and streamtraces in a plane between the nozzles Fig. 7: Average dimensionless velocity distribution in the wake of the launcher launcher and create a high pressure zone on the inner side of the jet boundaries. These high pressure gradient zones together with the normal shocks visible inside the jets justify the use of hybrid numerical methods able to capture shocks while maintaining a low numerical dissipation in the rest of the flow to resolve turbulent structures. One can indeed note in Fig. 6b that in addition to the azimuthal vortices issuing from the base, the mixing process between the jets yield the formation of a wide variety of turbulent structures in the inter-nozzle area. To complete the description of the organisation of the flow in this area, a cutting plane located between the nozzles is presented in Fig. 6c. On both sides of the main high-pressure zone corresponding to the confluence of the four jets, two secondary high-presure zones can be observed. Analysing the turbulent structures forming in the afterbody area, one can note that interactions occur between the structures from the main annular mixing layer and the structures from the interaction of the jets.

Mean flow topology

The Reynolds averaged field is computed on the fly during the unsteady ZDES calculation. As the vortex shedding period related to the base diameter was estimated to be T s = 2 * r b 0.2 * U∞ ≈ 1.6 ms, a period of 32 ms= 20 × T s was adopted to clear the flow from transient effect before using a 80 ms ≈ 50 × T s period to compute statistics. As an example, figures 7a and 7b display the mean velocity field as well as pseudo streamlines for both RANS and ZDES calculations. First, looking at a longitudinal plane containing two nozzles (Fig. 7a), one can observe the formation of a recirculation zone on the outer side of the nozzles. This recirculation zone interacts with the propulsive jets which entrain the external flow through a mixing layer. This entrainment effect appears stronger in the RANS computation as the recirculation bubble reattaches closer to the nozzle exit in the RANS computation ( Xr Ln = 1.10) than in the ZDES computation ( Xr Ln = 1.35). This first flow pattern corresponds to the one usually observed for axisymmetric base flow with a single nozzle (Reynaud et al. 2021 [START_REF] Reynaud | Numerical workflow for scale-resolving computations of space launcher afterbody flows with and without jets[END_REF], Statnikov et al. 2016 [START_REF] Statnikov | Experimental and numerical investigations of the turbulent wake flow of a generic space launcher at M ∞ = 3 and M ∞ = 6[END_REF]). Due to the expansion of the jets, one can notice that the area delimited by the launcher base, the nozzle wall and the jet boundary displays similarities with a supersonic cavity geometry with a L c = X r length (corresponding the location of the reattachement) and a L c /De c ≈ 2 depth ratio with De c its depth. Such an analogy, used experimentally by Wong et al. [START_REF] Wong | Experimental and theoretical investigation of base-flow buffeting on ariane5 launch vehicles[END_REF], is described in Fig. 8 presenting parameters L c and De c . In the present case, the cavity analogy concerns the geometrical area formed by the base and a nozzle. Fig. 8: Analogy between the afterbody area and a supersonic cavity, adapted from Wong et al. [START_REF] Wong | Experimental and theoretical investigation of base-flow buffeting on ariane5 launch vehicles[END_REF].

In the present case with four propulsive nozzles, additional flow interactions occur. Indeed, due to the impingement between the supersonic propulsive jets, the fluid located between the nozzles (see Fig. 7b) is subject to a strong adverse pressure gradient that leads to a reverse flow heading towards the base. We note that the RANS computation provides a higher reverse velocity than what is observed in the results from the ZDES simulation similarly to what is observed for single axisymmetric base flows (see Simon et al. [START_REF] Simon | Zonal-Detached-Eddy Simulation of Projectiles in the subsonic and transonic regimes[END_REF]). This reverse flow impacts the base and induces the formation of radial wall jets. In a plane containing nozzles (see Fig. 7a), these wall jets will be deflected as they go around a nozzle and part of the reverse flow will regain a positive longitudinal velocity. In a plane located between the nozzles (see Fig. 7b), these wall jets are pushed towards the base outer radius and feed the external recirculation zone. In such a plane, the prediction of the recirculation zone topology again differs between the two computations as the S-A model predicts a center position of the recirculation closer to the base than ZDES. Despite some local differences, the two bi-species computations are able to reproduce the main flow features of a 4-nozzle launcher wake flow (Mehta et al. 2013 [9]).

Gas mixing

To investigate the interaction between the propulsive jets and the base flow further, the average distribution of the propulsive mixture mass fraction (Y f uel ) in the afterbody area is displayed in Fig. 9. One can note that the RANS computation predicts a much more intense mixing between the two gases than the ZDES computation. This is shown by the greater predicted growth rate of the jet external mixing layer and by Y f uel levels in the inter-nozzle area nearly twice as high as those found with ZDES. Such discrepancies in the amount of hot propulsive gases advected back to the base area can be attributed to the strong influence of the turbulent Schmidt number on RANS mixing prediction (Reynaud et al. 2021 [START_REF] Reynaud | Towards Binary Gas Mixture ZDES for Space Launcher Base Flow Prediction[END_REF]) and is expected to lead to differences in base temperature predictions. This influence is strongly minimized thanks to the use of ZDES where turbulence in the base region is mostly resolved.

Base pressure distribution

In order to guide the design of future space launchers, numerical methods have to predict the base drag and the mechanical loads resulting from the previously described flow topology accurately.

As an example, Mehta et al. 2013 [START_REF] Mehta | Base heating sensitivity study for a 4-cluster rocket motor configuration in supersonic freestream[END_REF] performed computations with different levels of physical modelling for the propulsive jets; a frozen flow computation (Me1) where the propulsive mixture is considered as a perfect gas with a fixed γ j =1.15 (value corresponding to chamber conditions); a variable γ j computation (Me2) where γ j is a function of temperature and a reactive computation (Me3) where the propulsive gas is modelled as a 10-species mixture and chemical reactions are taken into account. Then, in Pu and Jiang 2019 [START_REF] Pu | Assessing turbulence models on the simulation of launch vehicle base heating[END_REF] several RANS models were tested to perform reactive computations.

Looking at the base pressure coefficient distributions displayed in Fig. 11, results from Metha et al. 2013 [START_REF] Mehta | Base heating sensitivity study for a 4-cluster rocket motor configuration in supersonic freestream[END_REF] showed that, with the Menter BSL (Baseline) RANS model, the use of the simplest physical model (refered to as Me1) surprisingly provided the best predictions for base pressure levels as both the variable γ j computation (refered to as Me2) and the reactive computation (refered to as Me3) lead to a significant underestimation of base pressure levels in comparison with experimental data. The same magnitude of error is found in Pu and Jiang 2019 [START_REF] Pu | Assessing turbulence models on the simulation of launch vehicle base heating[END_REF] as shown here with the predictions obtained with the RNG (Renormalization Group) turbulence model (Pu RNG), which provided the best results in their study. One can furthermore observe that all of the mentioned RANS computations are unable to reproduce the almost flat pressure profile found experimentally as they induce an overestimation of the pressure difference occurring between the center of the base and the outer radius. These discrepancies motivate the use of more advanced turbulence modelling such as the ZDES approach used here. Focusing on base pressure coefficient distributions obtained with the present computations, one can note that the use of an equilibrium flow hypothesis in the nozzles (with γ j =1.224) provides predictions coherent with numerical results from the literature, in between the frozen flow computation and the computations including chemistry effects. The use of ZDES improves predictions in comparison with RANS as it recovers a flat pressure profile and higher pressure levels near the outer radius. These disparities can be linked to the previously observed differences in the topology of the recirculation zone (i.e. characteristic sizes and maximum backflow velocity).

Nozzle pressure distribution

As reminded in the introduction, the accurate prediction of the engine bay environment and especially on the base is of major importance. As an example, the mean wall pressure (C p = P -P∞ For each simulation, similar pressure fields are predicted on the four nozzles. Notable differences are however observed when comparing the two computations for both the pressure levels and the spatial variation of C p .

To analyse these differences, the evolution of the pressure coefficient (Cp) has been studied along generating lines located at different angular locations (outer position (OUT), inner position (IN) and side positions (SIDE)) on a nozzle external wall as displayed in Fig. 13. The obtained pressure profiles are displayed in Fig. 14 and reveal the deviations between the RANS and the ZDES computations. Indeed, the use of ZDES provides relatively flat Cp profiles (close to Cp = -0.065) on every position around the nozzle whereas the RANS calculation predicts noticeable pressure axial variations for each location around the nozzle and a substantial pressure gap between the inner position and the outer position. Each turbulence modelling approach would thus lead to different estimations of the fluctuating pressure field used to study the mechanical behaviour of the nozzles. Based on previous numerical results obtained on launcher configurations with a single nozzle (Pain et al. 2014 [6], Reynaud et al. 2021 [START_REF] Reynaud | Numerical workflow for scale-resolving computations of space launcher afterbody flows with and without jets[END_REF]) one can note that the ZDES predictions are used to be more representative than the RANS ones.

Pressure fluctuations in the engine bay

To get a better understanding of the unsteady properties of the base pressure field, the unsteady pressure data obtained from the ZDES computation are investigated in the engine bay. First, as a consequence of the unsteady interactions occurring between the different flow features described in the previous section, the instantaneous distribution of wall pressure coefficient (Cp) shown in Fig. 15a, displays no symmetry and can be used to estimate the risk of nozzle vibrations due to side-loads. The rms value of pressure fluctuations namely the Cp rms (Cp rms = prms q∞ with

q ∞ = 1 2 γp ∞ M 2 ∞
the dynamic pressure) distribution on the launcher walls is displayed in Fig. 15b. The pressure fluctuations levels observed on the base lie in the range [0.001, 0.0025], such values are lower than the ones reported in [START_REF] Kawai | Time-series and time-averaged characteristics of subsonic to supersonic base flows[END_REF] (Cp rms ≈ 0.004) for a blunt base flow with a similar Mach number. The presence of the nozzles Looking at the propulsive nozzles, the maximum fluctuation levels (Cp rms > 0.01) are located on the edge of every nozzles near the launcher axis. To analyse the pressure fluctuations further, the next sections describe their evolution in selected areas: a ring located on the base r/r b = 0.85 shown in Fig. 16 as well as the lines along the nozzle walls presented in Fig13.

The azimuthal distribution of Cp rms on the ring r/r b = 0.85 (Figure 17) displays a quasi-periodic variation of the pressure fluctuation intensity with peaks located between the nozzles (45°, 135°, 225°, 315°) and an important decrease around locations aligned with the center of the nozzles (0°, 90°, 180°, 270°). The amplitude of Cp rms is not exactly the same and can partly be attributed to the limited Fig. 14: Axial evolution of Cp at different angular locations along a nozzle duration of the calculation. Nevertheless one can note that discrepancies between the extrema (e.g. between peak 1 and peak 2) correspond to small levels typically P rms ≈ 0.023% of the dynamic pressure.

These observations confirm the stabilizing effect of the nozzle wall and jet as reported by Deck and Thorigny [START_REF] Deck | Unsteadiness of an axisymmetric separating-reattaching flow: Numerical investigation[END_REF].

The evolution of the fluctuating pressure coefficient along the nozzles is then presented in Fig. 18. The overall increase of fluctuating pressure levels in the direction of the flow observed with the present computation is in line with previous RANS/LES results and measurements described in Statnikov et al. 2016 [START_REF] Statnikov | Experimental and numerical investigations of the turbulent wake flow of a generic space launcher at M ∞ = 3 and M ∞ = 6[END_REF] on a launcher with one nozzle placed in an external flow at M ∞ = 3. Furthermore, the computed fluctuation levels are in the same range as the ones reported by Statnikov et al. 2016 [START_REF] Statnikov | Experimental and numerical investigations of the turbulent wake flow of a generic space launcher at M ∞ = 3 and M ∞ = 6[END_REF] with Cp rms ∈ [0.0015, 0.01]. Important differences are however observed depending on the angular position. Indeed, the inner side of the nozzle is subject to pressure fluctuations twice as high as the outer side.

To complete this description, one can note that the predicted fluctuation levels for this supersonic configuration are an order of magnitude lower than the ones computed in Reynaud et al. 2021 [START_REF] Reynaud | Numerical workflow for scale-resolving computations of space launcher afterbody flows with and without jets[END_REF] for a transonic configuration. Indeed, as described in [START_REF] Reynaud | Numerical workflow for scale-resolving computations of space launcher afterbody flows with and without jets[END_REF], the transonic regime is more critical in terms of unsteady loads for a space launcher. At the altitudes corresponding to a transonic regime, a lower Jet Pressure Ratio (JPR) could furthermore correspond to a flight configuration without jet collision [START_REF] Faheem | Experimental study on the mean flow characteristics of a supersonic multiple jet configuration[END_REF] and thus with a flow topology notably different from the one studied here. Further computations using the present numerical tools to treat a transonic multinozzle configuration with hot propulsive jets could thus be the next step to guide the design of space launcher for this critical flight regime. The realisation of experimental studies on such a configuration would of course be necessary to assess the accuracy of the numerical studies. Such detailed databases are unfortunately currently missing in the open literature. The additional information given by the analysis of the spectral content of the pressure fluctuations is useful for the mechanical design of the afterbody components and can help identify the fluid phenomena driving the flow.

Spectral analysis

One-point analysis

The one-sided Power Spectral Density (PSD) function of pressure fluctuations, named G(f ) and expressed in Pa 2 /Hz describes how the rms value of pressure fluctuations p rms previously described is distributed in frequency according to:

p 2 rms = ∞ 0 G(f )d(f ) (3)
These spectra are computed using the classical Welch periodogram [START_REF] Welch | The use of fast fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms[END_REF] on a signal of total duration T = 62.5T S , where T S is the vortex shedding period (i.e. 100 ms in physical time). The PSD fluctuations are expressed as a function of the Strouhal number based on base diameter St D = f ×D b U∞ . For uni-nozzle configuration, one must remember that an important contribution to pressure fluctuations has been identified as an anti-symmetric (m = 1) mode occurring at a Strouhal number St D ≈ 0.2 [START_REF] Deck | Unsteadiness of an axisymmetric separating-reattaching flow: Numerical investigation[END_REF]. In the following, we investigate if this characteristic feature for canonical flows is still showed in this complex engine bay afterbody flow.

To begin with, the spectral content of wall pressure signals is obtained by data sampling on the r/r b =0.85 ring of the launcher base (see Figure 16 for a reminder of the base geometry). We first focus on the PSD for a probe located at θ = 225°, i.e., in a plane between two nozzles (Figure 19b). This spectra displays several features namely a peak around St D ≈ 0.2 (classically attributed to a vortex shedding phenomenon), a succession of peaks located in the range St D ∈ [0.5, 2], that will be discussed in the following, and a broadband contribution in the range St D ∈ [4, 10] due to small-scale turbulence. The probe located at θ = 180°(Figure 19a) shows similar properties but the broadband contribution (St D ∈ [4; 10]) seems strongly attenuated. According to the results presented in [START_REF] Weiss | Simulation numérique et analyse physique d'un écoulement d'arrière-corps axisymétrique et application au contrôle des charges latérales[END_REF], the broadband contribution in the range St D ∈ [4, 10] can be associated with fluctuations induced by Kelvin-Helmholtz instabilities in the initial zone of the annular mixing layer.

To get further insight into the physical origin of the PSD peaks observed in the frequency range St D ∈ [0.5, 2], one can follow the analogy between the afterbody flow topology and a supersonic cavity proposed by Wong et al. 2007 [START_REF] Wong | Experimental and theoretical investigation of base-flow buffeting on ariane5 launch vehicles[END_REF] (see Figure 8). For such a cavity, the unsteady behaviour of the flow is driven by Rossiter modes (Rm) observed experimentally [START_REF] Rossiter | Wind tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds[END_REF][START_REF] Zhuang | Supersonic cavity flows and their control[END_REF] and numerically [START_REF] Larchevêque | Large-eddy simulation of a compressible flow in a three-dimensional open cavity at high Reynolds number[END_REF][START_REF] Larchevêque | Large-eddy simulation of a compressible flow past a deep cavity[END_REF] in the literature. These modes are generated by a feedback loop linked to the convection of hydrodynamic instabilities in the mixing layer. These instabilities are generated by the separation at the upstream wall and travelled at the speed k c * U ∞ towards the downstream wall (in the present case: the propulsive jet boundary) where their impact generates pressure waves which propagate in the upstream direction (with a dephasing factor r c ) [START_REF] Larchevêque | Large-eddy simulation of a compressible flow in a three-dimensional open cavity at high Reynolds number[END_REF] in the subsonic part of the base flow. To estimate the Strouhal number St Db|Rm of these modes in the present case, an adaptation of the Rossiter [START_REF] Rossiter | Wind tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds[END_REF] formula is proposed to take into account the temperature of the flow in the recirculation zone:

St Db|Rm = D b L c m -r c U ∞ /a c + 1/k c (4)
The chosen parameters are k c = 0.57 to estimate the instabilities convection velocity and a dephasing factor r c =0.25 (both recommended for cavities with dimension between L c /D c =1 and 2 [START_REF] Rossiter | Wind tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds[END_REF]), a mean sound velocity a c ≈ 760 m.s -1 in the recirculation zone and L c = X r = 1, 35 × L n = 0.15525 m according to ZDES results (see figure 7a). The formula then provides the following values for the first three cavity modes: St Db|R1 = 0.49, St Db|R2 = 1.14, St Db|R3 = 1.8. These frequencies obtained for a supersonic cavity with M ∞ = 2.75 are coherent with the peaks observed in the range St Db ∈ [0.5, 2]. However, these frequencies do not corresponds exactly to the peaks noticed on Figure 19b. However, one can note that the under-predicted base pressure is associated to the prediction of the temperature and thus the sound speed in the recirculation zone. As a consequence, this observation can have an impact on the use of such an expression for the characteristic modes, in analogy to the Rossiter formula. Similar analyses on less complex configurations could be realised in a future work to support this hypothesis.

The vortex shedding frequency St D ≈ 0.2 therefore still plays a significant role in the dynamics of the flow as discussed in the next section.

Two-point analysis

To get a better knowledge of the spatial organisation of the fluctuating pressure field, one can consider the azimuthal coherence between two pressure transducers [START_REF] Reynaud | Numerical workflow for scale-resolving computations of space launcher afterbody flows with and without jets[END_REF][START_REF] Weiss | On the dynamics of axisymmetric turbulent separating/reattaching flows[END_REF][START_REF] Deck | Unsteadiness of an axisymmetric separating-reattaching flow: Numerical investigation[END_REF]. Indeed, as the base pressure is 2π-periodic, one can decompose the pressure field into azimuthal Fourier modes. As an example, the coherence C(f, ∆φ) between two probes located on the r/r b =0.85 ring and separated by an angle ∆φ can be decomposed in azimuthal modes as follows:

C r (f, ∆φ) = ∞ m=0 C r,m (f )cos(m∆φ) (5)
Figure 20 shows the C r,m spectrum for the first three azimuthal modes for a ring of sensors located at r/r b = 0.85. The axisymmetric mode (m = 0) is confined at low frequencies. The contribution of the ovalization mode (m = 2) is not significant over the whole frequency range. Of great interest is the spectrum of C r,1 which displays a peak near St Db ≈ 0.2 and shows that up to 80% of the pressure fluctuations at this frequency are due to the antisymmetric mode m = 1. Let us be reminded that the m = 1 mode corresponds to an antiphase relationship between probes facing each other yielding to the so-called buffet loads [START_REF] Weiss | On the coupling of a zonal body-fitted/immersed boundary method with ZDES: Application to the interactions on a realistic space launcher afterbody flow[END_REF][START_REF] Deck | Unsteadiness of an axisymmetric separating-reattaching flow: Numerical investigation[END_REF]. This important feature is also observed in more canonical flows [START_REF] Weiss | On the dynamics of axisymmetric turbulent separating/reattaching flows[END_REF][START_REF] Pain | Large scale dynamics of a high Reynolds number axisymmetric separating/reattaching flow[END_REF][START_REF] Marié | From pressure fluctuations to dynamic loads on axisymmetric step flows with minimal number of kulites[END_REF][START_REF] Weiss | Control of the antisymmetric mode (m = 1) for high Reynolds axisymmetric turbulent separating/reattaching flows[END_REF] showing the robustness of the dynamics of the antisymmetric mode in launcher base flows.

Conclusion

The flow dynamics of a four-nozzle afterbody has been investigated numerically with a scale-resolving approach. Assuming a two-species flow, an advanced numerical simulation has been performed taking advantage of both the automatic mode of ZDES namely Mode 2 (2020) which allows the RANS/LES interface to be set dynamically and a recently developed hybrid scheme that ensures robustness when shock waves occur and low dissipation levels in vortical regions. The numerical results are compared with standard RANS calculations as well as with other results published in the open literature. Concerning the mean pressure base, it is shown experimentally that an important characteristic is reproduced namely a constant level at the base of the afterbody. This important aspect is not recovered with a RANS approach while ZDES mode 2 (2020) simulates this salient feature. The analysis of the pressure fluctuations (Cp rms ) revealed high level of fluctuations that can reach up to 1% of the dynamic pressure in the central region of the base (i.e. between the four nozzles). Furthermore, one-point spectral analyses reveal the occurrence of a large scale dynamics related to a Strouhal number St Db = 0.2 together with other peaks at higher frequencies. An analogy is discussed with the occurrence of these peaks with those observed for supersonic cavity flows as suggested by Wong et al. [START_REF] Wong | Experimental and theoretical investigation of base-flow buffeting on ariane5 launch vehicles[END_REF].

The dynamics of this complex engine bay flow shares similar properties with its axisymmetric counterpart. Indeed, a two-point spectral analysis reveals that up to 80% of the energy of the pressure fluctuations are due to the antisymmetric mode m = 1 at a dimensionless frequency St Db = 0.2. On the use of bi-species ZDES for multi-nozzle space launcher configurations This shows the robustness of mode m = 1 which is observed for both simple and launcher base flows. Accurate validation of space launcher engine systems in extreme aerothermodynamics conditions will also require experimental databases which are still missing in the open literature.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Schematic view of the 4-nozzle launcher geometry (from Pu and Jiang 2019 [10])

Fig. 3 :

 3 Fig.3: Classification of typical flow problems, associated ZDES modes and examples of applications. mode 1: separation fixed by the geometry, mode 2: separation induced by a pressure gradient on a curved surface, mode 3: separation strongly influenced by the dynamics of the incoming boundary layer. RANS (mode 0) is the default mode (adapted from Deck[START_REF] Deck | Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF]).

Fig. 4 :

 4 Fig. 4: Instantaneous eddy viscosity field.

  topology 4.1.1 Instantaneous flow topology

Fig. 5 :

 5 Fig. 5: Instantaneous Mach number distribution

  (a) Mach number distribution and iso-surface of Q-criterion in the wake of the afterbody (b) Instantaneous pressure distribution and Q-criterion (black iso-lines) in the wake of the launcher in a plane containing two nozzles.

Fig. 9 :

 9 Fig. 9: Distribution of the fuel mass fraction in a plane containing the nozzles.

Fig. 10 :

 10 Fig. 10: Pressure coefficient iso-contours at the afterbody base (Top: RANS, Bottom: ZDES).

1 2 γairP∞M 2 ∞

 2 ) distributions predicted on the base and on the

Fig. 11 :

 11 Fig. 11: Base pressure distribution along an inter-nozzle radial line

Fig. 12 :

 12 Fig. 12: Distribution of mean wall pressure coefficient predicted on the external walls of the launcher.

Fig. 13 :

 13 Fig. 13: Angular positions of the generating lines along the nozzles.

  (a) Instantaneous wall pressure coefficient. (b) Fluctuating wall pressure coefficient.

Fig. 15 :

 15 Fig. 15: Iso-contours of wall pressure fluctuations in the base area.

Fig. 16 :

 16 Fig. 16: Description of the base geometry for the study of pressure fluctuations.

Fig. 17 :

 17 Fig. 17: Angular distribution of Cp rms on the r/r b = 0.85 ring. Vertical bars represent 90, 180, and 270 degrees.

Fig. 18 :

 18 Fig. 18: Cp rms distribution along the nozzle generating lines.

  (a) θ = 225°(b) θ = 180°F ig. 19: Power spectral density G(f )[Pa 2 .Hz -1 ] of pressure fluctuations on the base of the launcher. Green lines represent the calculated Rossiter modes.

Fig. 20 :

 20 Fig. 20: Contribution of the first three azimuthal modes to pressure fluctuations on the r/r b = 0.85 ring (initial point at 0°, aligned with a nozzle).

Table 1 :

 1 6 cells overall. Low altitude case initial flow conditions (Musial & Ward 1961[START_REF] Musial | Base flow characteristics for several four-clustered rocket configurations at Mach numbers from 2[END_REF]). • c refers to the nozzle chamber conditions.

	Parameters	Values
	Pc	41.37 bar
	Tc	3469.8 K
	P∞	12209.5 Pa
	T∞	297.7 K
	M∞	2.75

Table 2 :

 2 Sc t on RANS results. However, Reynaud et al. 2021[START_REF] Reynaud | Towards Binary Gas Mixture ZDES for Space Launcher Base Flow Prediction[END_REF] have shown that conversely to RANS methods, a weak dependence to Sc t values is observed for scale resolving simulations. Thanks to the reduced influence of unresolved turbulent fluctuations, these authors have shown that their numerical predictions with the scale resolving ZDES computations appear almost independent of Sc t . This weak dependence of the numerical assessments with regards to Sc t (similarly with P r t ) is precisely one of the interest to resorting to scale resolving simulations since the large-scale mixing processes are mostly resolved rather than modelled. Properties of the perfect gas simulating a JP4/LOX mixture

	µ 0j (N.s/m 2 )	1.85.10 -5
	T 0j (K)	370
	S j (K)	168
	γ j	1.224
	cv j (J.kg -1 .K -1 )	1997.5
	P r j	0.78
	P r t	0.90
	3.2 Numerical set-up
	3.2.1 Numerical method
	The FLU3M code developed by ONERA solves the compressible Navier-Stokes
	equations on multiblock structured grids. The accuracy of the solver for DNS, LES
	and hybrid RANS/LES purposes has been assessed in various applications including
	transitional flows (Mary & Sagaut 2002 [17]), wall-bounded turbulent flows (Pamies
	et al. 2009 [18]; Gand et al. 2010 [19]; Deck et al. 2014 [20], Deck et al. 2018 [21]),
	massively separated flows (Simon et al. 2007 [22]; Larcheveque et al. 2004 [23];
	Weiss et al. 2009 [24]; Weiss & Deck 2013 [25], Weiss & Deck 2018

Many authors (e.g. Baurle & Edwards 2010 [14], Clifton & Cutler 2007 [15]) have shown the dramatic importance of the value of the turbulent Schmidt number

(a) RANS (b) ZDES

Acknowledgments. The authors wish to thank the Centre National d' Études Spatiales (CNES) for financial support. The Ph.D. work of J. Reynaud is funded by CNES and ONERA. In addition, all the people involved in the evolution of the FLU3M code are warmly thanked.

Declarations

• The thesis of Jolan Reynaud is funded by CNES and ONERA.