INTRODUCTION

Gaseous combustion inside a combustion chamber has been studied for a long time due to its importance in industrial applications. However, spray combustion necessitates specific attention because injecting a liquid fuel inside a combustion chamber involves many different phenomena, like liquid atomization, evaporation, and the formation of a film on the chamber walls, which are not present in the purely gaseous case. Among the different questions, a point of interest is the determination of the global burning of a spray and the possible collective effects between the droplets. In both theoretical and numerical approaches as presented in the next paragraph, it is usually addressed by considering that a representative length scale of the droplet interdistance is equal to the mean droplet density at one-third power. While such an estimate is always valuable to build models, for instance, a more precise value of this scale-and the associated nearestneighbor interdroplet distance distribution and spray spatial distribution-could be beneficial for modeling spray evaporation and combustion.

Starting from the works of [START_REF] Godsave | Studies of the Combustion of Drops in a Fuel Spray-The Burning of Single Drops of Fuel[END_REF] or [START_REF] Spalding | The Combustion of Liquid Fuels, 4th Symposium on Combustion[END_REF] on the combustion of an isolated droplet of fuel, efforts have been made in the past on the theoretical, numerical, or experimental grounds to answer this question [START_REF] Annamalai | Interactive Processes in Gasification and Combustion. Part I: Liquid Drop Arrays and Clouds[END_REF]. In the 1980s, various research groups developed theoretical models describing spray combustion [START_REF] Chiu | Internal Group Combustion of Liquid Droplets, 19th Symposium on Combustion[END_REF][START_REF] Correa | The Group Combustion of a Spherical Cloud of Monodisperse Fuel Droplets[END_REF][START_REF] Kerstein | Percolation in Combusting Sprays I: Transition from Cluster Combustion to Percolate Combustion in[END_REF]. For instance, Chiu and his coworkers proposed the concept of droplet group combustion, considering a droplet cloud with monodisperse droplets arrayed at uniform intervals. The distance is deduced from the average density number of the cloud n. They proposed a droplet group combustion diagram with four spray combustion regimes: single droplet combustion, internal group combustion, external group combustion, and external sheath combustion [see Fig. 1 of [START_REF] Chiu | Internal Group Combustion of Liquid Droplets, 19th Symposium on Combustion[END_REF]]. The regimes were classified through a group combustion number G [START_REF] Jiang | Advanced Modelling of Spray Combustion Processes in Airbreathing Propulsion Combustors[END_REF]:

G = 1.5 Le 1 + 0.276 Sc 1/3 Re 1/2 N 2/3 d D i ( 1
)
where Le is the Lewis number, Sc is the Schmidt number, Re is the Reynolds number based on the droplet slip velocity, N is the total number of droplets contained in the cluster, d is the droplet diameter, and D i is the mean distance between droplet centers. Since then, various group combustion diagrams have been proposed (Borghi and Champion, 2020;[START_REF] Candel | [END_REF][START_REF] Reveillon | Analysis of Weakly Turbulent Dilute-Spray Flames and Spray Combustion Regimes[END_REF]. Borghi proposed a description of spray combustion in a premixed flow based on the Kerstein and Law model. The value n 1/3 r F , a ratio between the flame radius r F around a single droplet and the uniform distance given by the droplet density number n, depends on the ratio d/e L , where e L represents the flame thickness. For values n 1/3 r F < 0.41, droplets burn in groups. For values n 1/3 r F > 0.73, pockets of gas are surrounded by flames.

The intermediate case is a combination of the first two, in addition to the presence of a diffusion flame. Reveillon and Vervisch based their combustion diagram on some numerical simulations of a turbulent dilute-spray jet flame. Direct numerical simulation was used for the gaseous phase and a Lagrangian description was used for the spray. Particular attention was paid to the influence of the equivalence ratio and the dilute-spray density. Concerning the theoretical/analytical models, refer to the review paper by [START_REF] Annamalai | Interactive Processes in Gasification and Combustion. Part I: Liquid Drop Arrays and Clouds[END_REF], which describes several interaction models, from two droplets to arrays of droplets, and involves many different methods.

One of the most interesting approaches is that proposed by [START_REF] Elperin | Analysis of Evaporation and Combustion of Random Clusters of Droplets by a Modified Method of Expansion into Irreducible Multipoles[END_REF]. They considered a random cluster of polydisperse droplets in evaporation or combustion, using an irreducible multipole method. However, analytical models such as the previous one usually rest on some simplifying assumptions, like a uniform distance between droplets, gas at rest, or the constancy of the Lewis number Le = 1, which are hardly ever encountered in real situations. On the numerical side, spray combustion in most industrial solvers is generally based on the isolated droplet model under evaporation using the point-source approximation. Although some works have been done to describe the behavior of an isolated droplet in combustion more precisely [START_REF] Beck | Investigation of the Effect of Incomplete Droplet Prevaporization on NOx Emissions in LDI Combustion Systems[END_REF][START_REF] Beck | Identification of Droplet Burning Modes in Lean, Partially Prevaporized Swirl-Stabilized Spray Flames[END_REF] or to provide results regarding the influence of droplet density on vaporization and combustion in the case of droplet arrays (Imaoka and Sirignano, 2005a,b), we discuss only spray models. Chiu and some of his coworkers [START_REF] Jiang | Advanced Modelling of Spray Combustion Processes in Airbreathing Propulsion Combustors[END_REF][START_REF] Jiang | Multistate Behavior of a Droplet in Dilute Sprays[END_REF][START_REF] Jiang | Advanced Modelling of Spray Combustion Processes in Airbreathing Propulsion Combustors[END_REF] proposed several different approaches. In the most complex case [START_REF] Jiang | Multistate Behavior of a Droplet in Dilute Sprays[END_REF], the modeling refers to a droplet multistate behavior accounting for its ignition, the presence of a wake flame or an envelope flame, and the flame extinction. The droplet state is chosen by comparing several specific Damkholer numbers and the local Damkholer number. The ignition and extinction numbers are based on the asymptotic theory set forth by [START_REF] Law | Asymptotic Theory for Ignition and Extinction of Droplet Burning[END_REF] and [START_REF] Mawid | Chemical Kinetics Effects on the Ignition of a Fuel Droplet[END_REF]. As pointed out by [START_REF] Law | An Ignition Criterion for Droplets in Sprays[END_REF], improvement of the droplet ignition criterion would necessitate defining a characteristic length representative of the droplet distance. On their part, Wang and his colleagues [START_REF] Wang | A Two-Phase Turbulent Combustion Model and Its Validation for Spray Flames[END_REF] proposed to extend the evaporation model by considering four criteria: comparison of the evaporation and ignition delay time scales, respect of the flammability limit by the local fuel-air ratio, a minimal local temperature value, and a minimal droplet diameter value. The values of the last two criteria are respectively fixed at 1200 K and 1 micrometer. Assessment of the model is based on a large eddy simulation of an annular configuration with flow air for a methanol air spray flame. Finally, [START_REF] Paulhiac | Modélisation de la Combustion d'un Spray dans un Brûleur Aéronautique[END_REF] and [START_REF] Paulhiac | Analysis of the Spray Flame Structure in a Lab-Scale Burner Using Large Eddy Simulations and Discrete Particle Simulation[END_REF] also proposed a multistate model, called Mustard, for spray combustion. Compared to the previous model, Paulhiac introduces an ignition time delay and considers only an envelope flame for a burning droplet. In addition, the first step of the model determines the possible state of the droplet; that is, the vaporization state or envelope flame state. In order to do so, the interdroplet distance is estimated from the local liquid volume fraction or, equivalently, the mean local droplet density n. Such a model has been applied to large real-scale simulations.

On the experimental side, McDonnell et al. (1992) examined the structures of methanol sprays obtained by air-assisted or pressure-swirl atomizers under nonreacting and reacting conditions. Two-phase flow characterization is based on phase Doppler anemometry (PDA) and the infrared extinction/scattering method. In the case of swirling air-assisted spray, the analysis focuses on the differences between gas velocities, droplet velocities and diameters, and vapor concentration for both kinds of conditions. On their part, [START_REF] Akamatsu | Measurement of the Local Group Combustion Number of Droplet Clusters in a Premixed Spray Stream[END_REF] and [START_REF] Hwang | Evaluation of Combustion Mechanism of Droplet Cluster by Simultaneous Time-Series Measurement in Premixed Spray Flame[END_REF] worked on a premixed kerosene spray flowing through an annular burner. They performed simultaneous measurements of OH-and CH-light emissions, droplet Mie scattering emission, and droplet size and velocity using PDA. From their experimental data, Akamatsu and coworkers concluded that droplet clusters are present, but the local aspect of some measurements makes this issue questionable with regard to its extent or its presence. The authors also note that for many clusters, the group combustion mode does not correspond to Chiu's classification. Indeed, some assumptions are not supported by the experimental spray, such as the sphericity of the cluster or the flow uniformity.

More recently, [START_REF] Beck | Investigation of the Effect of Incomplete Droplet Prevaporization on NOx Emissions in LDI Combustion Systems[END_REF][START_REF] Beck | Identification of Droplet Burning Modes in Lean, Partially Prevaporized Swirl-Stabilized Spray Flames[END_REF] studied the effect of incomplete droplet prevaporization on NOx emissions in lean direct injection systems. In this experiment, two liquid fuel injectors are used to control the degree of prevaporization. The first one is far upstream of the combustion chamber and generates a fully prevaporized premixed incoming mixture. The second injector is located at the combustor inlet inside a swirl system, enabling the aerodynamic stabilization of the flame. Spray characterization has been performed using PDA and simultaneous CH-planar fluorescence imaging (PLIF)/Mie scattering. Isolated droplets burn in the very dilute part of the flow under two possible regimes: an envelope flame surrounding the droplet or a wake flame behind the droplet. According to their numerical analysis, these regimes have a significant influence on NOx emissions.

Two-phase flow combustion has also been studied by [START_REF] Lovett | On the Influence of Fuel Distribution on the Flame Structure of Bluff-Body Stabilized Flames[END_REF], within the framework of bluff body stabilized flames for a non-premixed, jet-in-crossflow fuel injection. Results on the flame structure are discussed in relation to classical combustion reaction zone regimes. Experimental analysis is mainly based on chemiluminescence and PLIF of OH radicals, to understand the development of the reactive zones with respect to the fuel injection or the upstream flow conditions. In all of these works, the characterization of the droplet spray in terms of spatial distribution or interdroplet distance is not addressed.

Finally, the work of Sahu, Hardalupas, and Taylor [START_REF] Sahu | Droplet-Turbulence Interaction in a Confined Polydispersed Spray: Effect of Droplet Size and Flow Length Scales on Spatial Droplet-Gas Velocity Correlations[END_REF][START_REF] Sahu | Droplet-Turbulence Interaction in a Confined Polydispersed Spray: Effect of Turbulence on Droplet Dispersion[END_REF][START_REF] Sahu | Interaction of Droplet Dispersion and Evaporation in a Polydispersed Spray[END_REF] provides one of the most comprehensive descriptions of the interaction between a spray and a turbulent flow. The experimental apparatus consists of a vertical test rig. Water or acetone is injected through an air-assisted nozzle inside a quiescent air flow characterized by large r m.s. velocities at ambient pressure and temperature. Nonevaporative water spray is used as a point of comparison with the evaporating acetone spray. Simultaneous coupling between interferometric laser imaging for droplet sizing (ILIDS) and PLIF measurements can describe the behavior of the spray. Thorough analysis is based on some correlations like the one between the fluctuations of droplet number densities and vapor mass fractions. Sahu and his colleagues also investigated the occurence of clusters of droplets and described their presence with respect to Chiu's theory thanks to a group evaporation number G. Such an approach should be extended to more realistic configurations with a swirling flow and under reacting conditions. This short review of the experimental studies shows that, although various aspects of spray combustion have been-and are still being-addressed, constituting an exhaustive experimental database could be of primary importance for increasing our knowledge of such flows and for assessing numerical approaches. Thus a new experimental device, the Prométhée-LACOM test rig, has been developed to study reacting and nonreacting two-phase flow [START_REF] Rouzaud | Experimental Analysis of Droplet Spatial Distribution in a Spray Burner[END_REF][START_REF] Vicentini | Régimes de Combustion dans les Foyers Aéronautiques: Mise en Place Expérimentale de l' Écoulement Diphasique, 14ème Congrès Francophone de Techniques Laser[END_REF][START_REF] Vicentini | Experimental Investigation of Spray Combustion Regimes in Aero-Engines Combustors, 6th European Conf. for Aeronautics and Space Sciences[END_REF]. The configuration corresponds to a turbulent air flow circulating around a bluff body in a square-section channel. Liquid fuel is injected at the rear face of the bluff body through a flat-fan atomizer. The obstacle also acts as a flameholder. The setup has been designed to facilitate experimental observations and to reproduce, at least partially, the operating conditions of an aero-engine combustion chamber. The former constraint has been achieved by retaining a two-dimensional geometrical configuration. Characterization of the nonreacting flow has been presented elsewhere [START_REF] Vicentini | Régimes de Combustion dans les Foyers Aéronautiques: Mise en Place Expérimentale de l' Écoulement Diphasique, 14ème Congrès Francophone de Techniques Laser[END_REF] and the scope of the paper is thus limited to the reacting case. More particularly, emphasis is laid on the droplet spatial distribution; that is, the interdroplet spacing and spatial distribution law. To the best of our knowledge, such a characterization has not been done for burning sprays, but realized in other fields like metallography [START_REF] Bansal | Average Nearest-Neighbor Distances between Uniformly Distributed Finite Particles[END_REF], particle transport (Neumann and Umhauer, 1991) and meteorology [START_REF] Kostinski | On the Spatial Distribution of Cloud Particles[END_REF][START_REF] Kostinski | Scale-Dependent Droplet Clustering in Turbulent Clouds[END_REF]. A usual analysis considers the clustering of the particles with several available treatments like the Cluster Index [START_REF] Monchaux | Analyzing Preferential Concentration and Clustering of Inertial Particles in Turbulence[END_REF]. For instance, Kostinski et al. demonstrated that the spatial distribution of particles inside a cloud is described rather by a Poisson mixture process than by a Poisson process (uniformly random distribution; [START_REF] Chandrasekhar | Stochastic Problems in Physics and Astronomy[END_REF]. In our case, the spray analysis is based on a straightforward analysis of droplet Mie scattering data.

After a description of the experimental apparatus and operating conditions in Section 4, the methodology employed to process the droplet Mie scattering data is presented in Section 3. Such a procedure was used to estimate, as locally as possible, values of the spray density number and of the average nearest-neighbor interdroplet distance in the flow based on images. Since the analysis of the experimental data has also been extended to the spatial distribution law, this section also contains a short reminder of two theoretical spray distributions, the uniform and uniformly random distributions, which will be used for comparison with our data, and a description of an in-house Monte-Carlo solver developed to further support the analysis. This reproduces a spray with uniformly randomly distributed droplets. Results are presented in Section 4. At first, the method was used with the experimental droplet Mie scattering data. Six different locations in the flow were chosen in order to describe the spray evolution inside the channel. It constitutes the main objective of this paper and brings answers to the question addressed at the end of the first paragraph, the description of a nearest-neighbor interdroplet distance and spray spatial distributions. The method was also applied to Direct Numerical Simulation-Direct Particle Simulation (DNS/DPS) with an isotropic homogeneous turbulence field. Six different cases were considered with different Stokes particle numbers to assess the interest of the present procedure and tentatively the influence of the clustering effect on the analysis.

EXPERIMENTAL SETUP AND MEASUREMENTS

Test Setup

A new air-breathing propulsion test setup, Prométhée-LACOM, has been recently developed at ONERA (see Fig. 1) to study spray combustion. The air feeding system is made up of a large pressure vessel, gas pressure regulators, and an electric heater (1 MW). The air mass flow rate is measured with a sonic nozzle located downstream of the electric heater (± 3.3% accuracy). Subsequently, the preheated air stream passes through a succession of flow conditioners and turbulence grid. The test chamber consists of a square internal section (120 mm × 120 mm). The confined turbulent air stream flows around a fuel injection system prior to the combustion chamber. The injection system consists of a bluff body in which a liquid fuel nozzle is mounted. The bluff body spans over the entire width of the test section. The fuel feeding line is equipped with a Coriolis flowmeter (± 0.2% accuracy). The spray nozzle is centered and screwed into the rear face of the obstacle. A spray of droplets is generated downstream directly into the combustion chamber. The fuel injector corresponds to a flat fan nozzle that produces an elliptical-shaped polydisperse spray. The spray opening angles of the major and minor axis are respectively 150°a nd 25°. The wide angle is set up to spray toward the upper and lower walls of the combustor. Note that in this study, the injector provides a nonsymmetrical spray with a larger number of droplets in the lower part of the flow. Such a feature let us consider different droplet density numbers in the upper and lower parts of the flow. The water-cooled combustor is equipped with three UV-transparent windows that can perform optical measurements downstream of the injection system. Ignition of the air-fuel mixture is triggered by a hydrogen-oxygen torch. The burnt gases are ejected into the exhaust pipe. This setup partially reproduces the main features of the flow inside an actual aero-engine combustor. The air flow is confined, turbulent, and exhibits instabilities and large-scale structures. Besides, the flame is anchored aerodynamically by creating a main recirculation zone. The fuel-air equivalence ratio replicates the operating lean-burn conditions during a flight phase at cruising speed. A previous sizing study has been carried out on a prototype in 1:1 scale to choose the bluff body geometry [START_REF] Vicentini | Régimes de Combustion dans les Foyers Aéronautiques: Mise en Place Expérimentale de l' Écoulement Diphasique, 14ème Congrès Francophone de Techniques Laser[END_REF] under nonreacting conditions. For a 6 m/s inert bulk airflow velocity at 298 K, a vortex-shedding phenomenon and a large recirculation zone are observed downstream of the bluff body. The criterion of choice for the bluff body geometry and dimensions was based on the intensity of the pressure signal measured on the bluff body upper and lower faces with respect to the noise. As a result, a trapezoidal bluff body with a 42% ratio blockage, defined as the ratio of the projected area of the obstacle upon the unobstructed cross section of the channel, has been retained (Fig. 1). Note that in reacting conditions, the von Kármán vortex street is no longer observed while the recirculation zone is still present. A large number of droplets is expected to be located around the vertical mid-plane of the combustor. Since some of the experimental diagnostics used are planar (imaging devices), the configuration promotes the measurement of many relevant flow data (droplets positions, sizes, velocities, and flame position) compared to an axisymmetric one.

The n-decane is selected as a surrogate fuel because it has chemical properties close to kerosene ones (density, surface tension, auto-ignition temperature). Furthermore, evaporation or chemical kinetics modeling are much simpler for a monocomponent fuel than for a mixture, should this configuration be calculated. Under nonreactive conditions, the droplet spray is polydisperse with droplet sizes from 10 to 150 µm [START_REF] Vicentini | Régimes de Combustion dans les Foyers Aéronautiques: Mise en Place Expérimentale de l' Écoulement Diphasique, 14ème Congrès Francophone de Techniques Laser[END_REF].

Operating Conditions

Although nonreacting and reacting two-phase flows have been investigated, we only consider reacting flow in this paper. The nominal operating point is set at an air mass flow rate equal to 64 g/s. At the combustor inlet, the air flow is at standard atmospheric pressure and 450 K. The corrected mass flow rate corresponds to 1.36 kg/s/bar/K 1/2 . Since the combustion section is 120 mm × 120 mm, the bulk airflow velocity is close to 6 m/s. Consequently, the Reynolds number based on the hydraulic diameter is about 22,000. The liquid fuel temperature is measured just upstream from the nozzle and is equal to 330 K. The n-decane fuel (95% purity) is injected at 1 g/s to provide 44 kW thermal power. Therefore the fuel-air equivalence ratio is around 0.24. The water cooling of the combustion chamber protects the inner walls (T wall < 650 K) and we can perform tests lasting several tens of minutes.

Implementation of Measurements

The main objective of the experimental program is to build a spray combustion database for a better understanding of droplet combustion and for the assessment of numerical approaches. In order to do so, gaseous and liquid phases have been studied in nonreacting and reacting conditions. Concerning the reacting conditions, we employ OH-PLIF imaging, droplet Mie imaging, and phase Doppler interferometry (PDI) techniques. Chemiluminescence and OH-PLIF imaging are employed to study the hot regions of the flow. The first technique uses spontaneous emission of the OH flame marker while in the second technique, excitation of the marker is produced using of a laser. Droplets illuminated by a laser sheet generate light through the Mie scaterring process. The light is captured by a camera in order to evaluate the positions of droplets. Simultaneous OH-PLIF and Mie scaterring measurements provide instantaneous information about the droplet spatial distribution with respect to the reacting zones. Finally, application of a twocomponent PDI measures the instantaneous droplet velocities (longitudinal, transverse) and their sizes. Except for the OH-chimiluminescence, all the measurements have been performed in the vertical mid-plane of the combustor.

Figure 2 presents the optical system for both OH-PLIF and Mie scattering techniques. Due to the scope of the paper, we only give details on the Mie scattering setup. To light droplets, a millimetric thin laser sheet is produced with a high-frequency Nd:YLF laser Quantronix Darwin (λ = 527 nm) using a cylindrical lens. The laser sheet crosses the combustion chamber from the bottom window. The frequency of the system is fixed to 1 kHz with a laser pulse duration around 200 ns. A high-speed camera Phantom V341 is positioned aside to capture the Mie scatterred light from the droplets through a lateral window of the combustor. Image resolution is 1792 pixels × 1600 pixels. In the sequel, analysis of Mie scaterring is performed in six specific regions of the flow, considering only small parts of the overall Mie image (Table 1). Such a choice allows a relatively local characterization of the spray. Figure 3 The droplets are in white and their sizes have been enlarged to make them visible on the image. The significance of the Roman numbers will be explained later.

SPRAY ANALYSIS

Analysis of the Interdroplet Distance

Since the experimental work is based on Mie scattering images, we have at our disposal a certain number N Mie of them. Several steps are needed before performing a statistical analysis on the interdroplet distance D i . As a preliminary remark, we have to mention that the droplets are located in a volume of which the dimensions correspond to the size of the Mie scattering image times the thickness of the laser sheet ǫ (close to 1 millimeter in this case). It means that the interdroplet distance is based on the orthogonal projection of the droplets from the measurement volume to the image surface (projection process, Fig. 4). At first, each image R has been processed to suppress the soot luminance and the particles that are out of the depth of field. To meet the first objective, a smoothing method, based on a sliding-window algorithm (16 pixels × 16 pixels), is applied. This algorithm acts as a linear lowpass filter on the luminance intensity. Concerning the out-of-depth of field particles, we define a threshold value for the light intensity (value of 85 out of 4095 levels) and retain only the droplets with scattered-light intensity above this level.

FIG. 4: Description of the data processing for the error analysis

The second step consists of determining the mean density number of particles n over all the images. By denoting N R the number of particles in the current image R, the associated density number per unit surface n R is straightforwardly given by

n R = N R Window area (2)
leading to the definition of the mean density number

n = 1 N Mie N Mie R = 1 n R (3)
Determination of the interdroplet distance constitutes the next step of the data processing. For each droplet S in an image R, its center of gravity is deduced from the first-order moment integral of its pixel coordinates according to

       X R,S = 1 A S A X dXdY Y R,S = 1 A S A Y dXdY (4)
where A S is the droplet area. The interdroplet distance D R,S i is the minimal value of the Euclidean distance between the Sth droplet and all other droplets S * belonging to the same image,

D R,S i = Min ∀S * = S (X R,S -X R,S * ) 2 + (Y R,S -Y R,S * ) 2 (5)
Once all the images have been processed, the mean value and the standard deviation of the interdroplet distance are evaluated as

E [D i ] = 1 N Mie N R N Mie R = 1 N R S = 1 D R,S i (6)
and

SD [D i ] = 1 N Mie N R N Mie R = 1 N R S = 1 D R,S i -E [D i ] 2 (7)
An additional quantity is the equivalent regular distance D unif i written as

D unif i = 1 N Mie N Mie R = 1 n R -1/2 (8)
Such a distance is often used to underline the neighborhood effects between droplets in evaporation or in combustion ratio, through the distance parameter C, between this distance and the droplet diameter or the flame diameter around a droplet. For instance, according to [START_REF] Crowe | Multiphase Flows with Droplets and Particles[END_REF] and [START_REF] Virepinte | Droplet Spacing on Drag Measurement and Burning Rate for Isothermal and Reacting Conditions[END_REF], there is no interaction between two droplets in terms of evaporation for a distance parameter from 10 to 15. Simulations of Imaoka and Sirignano (2005a,b) on droplet arrays confirm that the evaporation/combustion interaction is quite reduced over a certain value of the distance parameter. In most aeronautical combustion chambers, the value of this parameter is between 5 and 40 [START_REF] Candel | [END_REF][START_REF] Correa | The Group Combustion of a Spherical Cloud of Monodisperse Fuel Droplets[END_REF]. At this point, the experimental procedure to get the quantities of interest from the Mie scattering images has been established. However, two additional remarks have to be made to fully describe the experimental results that will be presented in the next section. First, only restricted parts of the images, the so-called analysis windows in the sequel, are taken into account during the restriction process (Fig. 4). These windows are square-shaped with an edge length l. The reason for such a choice is that, since the spray is not at all uniformly distributed inside the test section, averaging over all the droplets of all the images will produce quite inaccurate values of interest. On the contrary, by carefully choosing the sizes of the analysis windows, we may expect to get local values of the quantities of interest with a sufficient number of droplets to perform a relevant statistical analysis on these values. In this analysis, six analysis windows have been located in different parts of the Mie scattering images in order to describe the spatial evolution of the spray. Their locations have been chosen empirically using the Mie images to describe different positions in the two-phase flow and different mean density numbers. Their sizes vary from 35 × 35 to 200 × 200 pixels depending on their locations inside the flow.

Second, the values provided by Eq. ( 5) may be biased by at least two sources of errors (Fig. 4). The first one, called the projection error, is due to the orthogonal projection of the droplets from the probe volume to the image surface previously mentioned. Most of the time, it leads to an underestimation of the distance between two droplets when they are not in a plane perpendicular to the projection direction. Occasionally, the distance is conserved when the two droplets are inside a plane perpendicular to the projection direction. The second type of error comes from the restriction process, described below and illustrated in Fig. 4. When considering a droplet inside an analysis window, it may happen that its nearest-neighbor droplet is not in the analysis window itself but just outside of it. In that case, the droplet distance is systematically overestimated and its associated error is called the restriction error or border zone effect as mentioned by Neumann and Umhauer (1991). A last comment arises about a third source of error associated with the presence of preferential segregation in the spray. This question will be partially addressed in Section 4.2.

Finally, although this procedure has been developed primarily for Mie scattering images, its extension to a three-dimensional field of droplets, like the ones obtained by a two-phase flow simulation, is quite straightforward. We have to replace the window area and the density number per unit surface in Eqs. ( 2) and (3) by a window volume and mean density number per unit volume

n R 3D = N R Window volume (9) n3D = 1 N Mie N Mie R = 1 n R 3D (10)
An additional feature is to account for the third coordinate of the barycenter droplet Z R,S and the distance between a pair of droplets:

D R,S i = Min ∀S * = S (X R,S -X R,S * ) 2 + (Y R,S -Y R,S * ) 2 + (Z R,S -Z R,S * ) 2 (11)
In the next section, such a procedure will be applied to a DNS/DPS simulation to extract data on the droplet spatial distribution, starting from the second step of the procedure, thus avoiding the projection step.

Analysis of the Spray Spatial Distribution

To get further insight in the spray organization, we can also try to characterize the spatial distribution of the droplets in each analysis window by comparison with the uniform (droplets arrayed at uniform intervals) and the uniformly random distributions. In the two-dimensional case, given a mean surfacic density number n, the interdroplet distance between nearest neighbor writes, in the case of the uniform distribution, as

D unif i = n-1/2 (12)
while in the three-dimensional case, we have

D unif i = n-1/3 3D ( 13
)
where n3D represents the usual mean density number. Note that the standard deviation is null for this type of distribution in all cases. The Hertz-Chandrasekhar distribution refers to a distribution where the droplets or the particles are uniformly randomly dispersed. Using the Hertz-Chandrasekhar theory (see Appendix), we can demonstrate that, in the two-dimensional case, the mean distance and its standard deviation are expressed as

E D H-C i = 0.5 n-1/2 , SD D H-C i ≈ 0.2613 n-1/2 (14)
while the three-dimensional counterparts [START_REF] Chandrasekhar | Stochastic Problems in Physics and Astronomy[END_REF] read as

E D H-C i ≈ 0.5539 n-1/3 3D , SD D H-C i ≈ 0.2013 n-1/3 3D (15)
Thus, for both laws, the interdroplet distance varies linearly with respect to either the inverse of the square-root of the surfacic density number n or the inverse of the cubic-root of the volumic density number n3D , but with different slope coefficients. In the sequel, we will compare the behavior of the experimental data with the two theoretical laws. Two main questions drive this comparison. The first one addresses the relevancy of the uniform distribution to describe a spray. Such a choice is often done to estimate the interdroplet distance and, possibly, used for modeling. The second one is to verify if a linear relationship between the interdroplet distance and the density number n exists in our experiment such as

E [D i ] = α n-1/2 (16)
and

SD [D i ] = β n-1/2 (17)
Finally, to further support our analysis, we have developed a Monte-Carlo solver to compute the interdroplet distance for a set of point particles uniformly randomly dispersed. Although this configuration is not exactly similar to a real spray (the droplets are not pointwise), it is expected that such a solver will contribute to explain some of the tendancies of our results. The locations of the particles are randomly chosen into a three-dimensional volume of which the dimensions are provided by the user. The mean density number nl is also specified but must be chosen according to the case under consideration. For instance, referring to the analysis of the Mie scattering images, we can associate a surfacic mean density number n to an analysis window. To build a three-dimensional droplet distribution corresponding to n, the mean number of particle per unit volume n3D writes as n3D ǫ = n (18)

where ǫ represents the thickness of the three-dimensional slice. The previous relationship simply states that the number of particles is exactly the same in the three-dimensional volume and in the two-dimensional window. If the droplets are uniformly randomly dispersed, we can expect to approximate the three-dimensional field of droplets prior to the projection process, and to reproduce the projection and restriction steps of the experimental procedure.

RESULTS AND DISCUSSION

Spray Description

Figure 5 presents the system of coordinates and the locations of the analysis windows. The (X, Y ) coordinates correspond respectively to the longitudinal and transversal axis, and are in the vertical mid-plane of the combustor. The origin of the frame of reference is on the lower wall of the channel. It is vertically aligned with the rear face of the bluff body. The center of the injection system is at coordinates (0, 60) in millimeters. Downstream of the obstacle, we can see two OH-PLIF images taken in the upper part and in the lower part of the flow but not at the same instant. They have been assembled for a better understanding of the flow. Figure 5 also indicates the position of the six windows of analysis that have been considered to study the droplet spatial distribution using Mie scattering images. The coordinates of the center of the six windows are given in Table 1. Two points are located in the upper part of the channel above Y = 60 mm and the others in the lower part. Their locations have been chosen in order to account for the spray density, the distance to the injection system, and the flame presence. The table also specifies the sizes of the analysis windows in micrometers. These sizes FIG. 5: Positions of the analysis window with respect to the obstacle have been considered in order to provide statistical data on the spray distribution as locally and as accurately as possible. The local constraint means that the size of the window should not be too large to define values representative of the local spray distribution. On the accuracy part, a sufficient number of samples (droplet positions) is needed to perform reliable statistics. The next paragraph provides more details on this aspect.

The experimental work presented hereafter is split in two parts. The first part of the work deals with the experimental analysis of the center-to-center interdroplet distance to the nearest neighbor since this value is commonly considered in evaporating or reacting sprays to characterize droplet interactions. In the second part, the spray distribution law is investigated and compared to two existing laws, Hertz-Chandrasekhar and one meaning a cubic arrangement of droplets.

Droplet Spacing

Experimental Analysis of the Interdroplet Distance

For each analysis window, we have at our disposal N Mie = 2000 images. Table 2 summarizes the main results obtained at the six points under consideration. Starting from the second column, the data yield the number of samples, the mean density number n, the mean interdroplet distance E [D i ] and its associated standard deviation SD [D i ], the so-called uniform distance D unif i , and the distance parameter C. The last value corresponds to the ratio of the interdroplet distance to the local mean arithmetic droplet diameter d 10 . This latter quantity is based on PDI measurements done inside each analysis window. The regular distance D unif i is calculated as if the droplets were uniformly arrayed inside the window.

A first comment bears on the evolution of the mean density number and the mean interdroplet distance. According to the data, the mean density number varies from a larger value close to the flat fan injection (point I) to smaller values far away (points III and VI). Such an evolution is consistent with the presence of the flame involving high-temperature and a more rapid evaporation of the smaller droplets compared to the nonreactive case [START_REF] Vicentini | Mise en Évidence Expérimentale et Modélisation des Régimes de Combustion Diphasique Présents dans les Foyers Aéronautiques[END_REF]. Indeed, there is a factor 20 in the density number at point III between the nonreacting and the reacting cases as observed by Vicentini. Note also that the quite dense spray at point I precludes the flame penetration close to the injection system as revealed by the Mie/OH-PLIF superimposed images [START_REF] Vicentini | Mise en Évidence Expérimentale et Modélisation des Régimes de Combustion Diphasique Présents dans les Foyers Aéronautiques[END_REF]. The mean interdroplet distance values are consistent with the mean density number behavior. In particular, it increases with the distance to the injection system. Referring to the distance parameter C in the last row, we note that the spray evolves between TABLE 2: Outcomes of the experimental distributions of droplets

Window Nb. samples n (1/cm 2 ) E [D i ], SD [D i ]
(µm, µm) droplet interactions and isolated droplet behaviors according to the admitted description [START_REF] Crowe | Multiphase Flows with Droplets and Particles[END_REF]. Second, we compare the interdroplet distance to the regular one. In all cases, the experimental values are smaller than the regular values derived from the experimental density number with a factor varying between 1.3 and 2. Thus, in this experiment, the regular distance is not the most representative distance for droplet interactions.

D unif i (µm) C = E [D i ]/

Size of the Analysis Windows

Since we intend to determine values of the droplet density number as local as possible, the size of the windows should be small with respect to the overall Mie image. On the other side, performing a statistical analysis on the interdroplet distance D i requires considering a sufficient number of samples (droplets). Due to the spray inhomogeneity, the trade-off between these two constraints leads to choosing different sizes according to the window locations.

In the case of window I, Table 3 gathers the mean density number, the mean nearest-neighbor distance, and the associated standard deviation for three different sizes of the analysis window from simple to double. To choose the size of analysis window I, the main guidelines were to get a sufficient number of samples (arbitrarily fixed above 25,000) and to check that an increase or decrease of its size would provide limited differences as seen in Table 1. For the other analysis windows, the size of the analysis window has been approximately defined by the ratio of their density number to the density number of window I and we have checked that decreasing or increasing its size provided consistent results.

Droplet Spatial Distribution Law

To get further insight in the spray organization, we now try to characterize the spatial distribution of the droplets in each analysis window by comparison with droplets arrayed at uniform intervals and uniformly random distributions. Since the experimental values of the density number depend on Mie scattering images taken at different instants, we have decided to retain only the images of which the density number is close to the mean value n in order to limit the too-large density number fluctuations and to be as close as possible to the perfect conditions of the two theoretical laws. In this work, the selection criterion is fixed to ± 10%, meaning that we retain images for which n ∈ [0.9n, 1.1n]. When compared to the results accounting for all the images, such a choice leads to a relative error on the mean interdroplet distance around 10% with a reduced standard deviation.

The mean distance and standard deviation of both theoretical and experimental distributions are plotted respectively in Fig. 6(left) and Fig. 6(right) for the six windows of analysis. The extremal values of the mean density number correspond to the windows for which the spray is respectively less or more dense. As mentioned in Section 2.1, the nonsymmetrical injection

Error Analysis

To estimate the accuracy of the experimental results on the droplet spacing or the spatial distribution, the evaluation of two different sources of error have been performed. The first one is associated with the identification of the droplet barycenters on each Mie scaterring image. The resolution of the images is 35.2 pixels per millimeter (one pixel corresponds to 28 micrometers).

Thus the error committed on the barycenter coordinates of a spherical droplet is expected to be around 20 µm. Since the coordinates of two droplets are used to calculate their interdroplet distance, the error on this distance is, in the worst case, around 40 µm. Referring to Table 2, the relative error with respect to the mean distance E [D i ] varies between 4.5% and 22%.

Consider now the errors associated with the projection and restriction steps (Fig. 4). To assess their importance on the interdroplet distance, the Monte-Carlo numerical approach has been employed. Numerical droplets are randomly dispersed according to the Hertz-Chandrasekhar law inside a slab volume of which the dimensions are equivalent to the dimensions of the experimental Mie scattering images. The statistical analysis of the Monte-Carlo data follows the same procedure than the experimental one. The projection step consists of projecting the droplets on the (X, Y ) plane, discarding the third coordinate. The interdroplet distance is then calculated for both three-dimensional and two-dimensional distributions, and the comparisons of these two values provide an estimation of the projection error.

In the next step-the restriction one-the large two-dimensional numerical image is split into a set of smaller images of which the sizes are comparable with the sizes of the analysis windows. Comparison of the interdroplet distance for the large image and the set of the smaller images provides an estimation of the restriction error. As previously mentioned, this difference is due to the fact that, for some droplets located in the border zone of a small image, its nearestneighbor droplet is out of the image. Such a border zone effect (Neumann and Umhauer, 1991) is all the more important as the size of the image is smaller.

Table 4 summarizes the projection and restriction errors for the six analysis windows. The total error corresponds to their product. As expected, the projection step decreases the interdroplet distance (minus sign). The projection error is all the more important as the density number inside the measurement volume is large. Possible improvements would be to reduce the laser sheet thickness. On the other side, the restriction effect overestimates the interdroplet distance because of the border zone effect. Minimizing the error would be possible by adding an additional treatment for the droplets of which the distance to the image boundaries is smaller than its calculated interdroplet distance. In that case, we should either discard these droplets from the data or calculate their real interdroplet distance considering outer droplets.

Finally, we can also question the point particle approximation influence in the estimation of the nearest-neighbor distance. By point particle approximation, we mean that each particle is described by the coordinates of its barycenter and that its size is not taken into account in the distance analysis. [START_REF] Bansal | Average Nearest-Neighbor Distances between Uniformly Distributed Finite Particles[END_REF] have considered such an influence in three different cases: the nearest-neighbor spacing in a three-dimensional problem, the nearest-neighbor spacing of finite parallel cylinders in a two-dimensional problem, and the nearest-neighbor spacing in a plane section through a three-dimensional array of finite spheres. This last problem is very similar to our case. For volume fractions α l lower than 0.01, [START_REF] Bansal | Average Nearest-Neighbor Distances between Uniformly Distributed Finite Particles[END_REF] demonstrate that the point particle approximation provides very accurate results. Referring to the analysis windows, the largest value of the volume fraction is associated with window I and is equal to 0.005. Thus the point particle approximation yields consistent results.

Analysis of the Clustering Effect on Nearest-Neighbor Interdroplet Distance from DNS/DPS Simulations

For two-phase flows experiencing preferential segregation, the spray distribution will be modified compared to the uniformly random one. Therefore it is interesting to describe as far as possible the effect of some clustering effect on the values of the coefficients α and β. Such a description could also serve to complete the analysis of the present experimental data. In order to do so, numerical DNS/DPS have been used. [START_REF] Fede | Monte-Carlo Simulation of Colliding Particles or Coalescing Droplets Transported by a Turbulent Flow in the Framework of a Joint Fluid-Particle PDF Approach[END_REF] realized DNS/DPS calculations of solid particles. The particles are spherical and the spray is monodisperse with a particle diameter equal to 600 micrometers. The particles are submitted to an isotropic homogeneous turbulence inside a cubic box of 128 millimeters long containing 450,000 particles corresponding to a mean density number of n = 215 per cubic centimeter. Six calculations have been done using six different particle densities and the Stokes number:

St = τ F f p τ t f @p F
defined as the ratio of the fluid-particle relaxation time scale τ f p to the turbulent Lagrangian time scale τ t f@p , varies from 0.22 up to 4.58 (0.22, 0.43, 0.85, 1.70, 3.25, 4.58). Preferential concentration behavior is clearly exhibited for Stokes numbers lower than 1.70 as shown in [START_REF] Fede | Monte-Carlo Simulation of Colliding Particles or Coalescing Droplets Transported by a Turbulent Flow in the Framework of a Joint Fluid-Particle PDF Approach[END_REF] (Fig. 13).

Previous methodology used for Mie scattering images processing has been applied to these fields, but in a three-dimensional configuration (see end of Section 3.1). More precisely, the computational domain has been cut in smaller cubic boxes with two different sampling lengths respectively equal to 30 and 45 millimeters. The average number o f p all el er y domains are respectively close to 6000 and 20,000. The numbers N , N equal to (6000, 64) and (20,000, 8) for the smallest length and the largest one. In a second step, droplet distances have been calculated and statistics performed to estimate the dimensionless mean and standard deviation α and β defined in Eqs. ( 16) and ( 17) but for a three-dimensional case. Their respective theoretical values according to the Hertz-Chandrasekhar theory are equal to 0.5539 and 0.2013. In that case, the projection error is null since we consider the real distance between droplets while the restriction error is still present.

The results are gathered in Fig. 7 presenting the α and β coefficients versus the Stokes number. At first, we note that, referring to the two plots, the α and β coefficients for the mean minimal distance and its standard deviation tend respectively to their theoretical values 0.5539 and 0.2013 with an increasing Stokes number. This tendancy is observed for both sizes of the sampling box. A second point of concern is the evolution of the α and β coefficients when the Stokes number FIG. 7: Clustering effect decreases, leading to a decrease of the interdroplet distance and an increase for the standard deviation. A possible explanation of this behavior is the enhancement of the segregation phenomenon, which becomes more and more important as the Stokes number decreases. This last effect is well known in the literature [START_REF] Simonin | Connection between Two Statistical Approaches for the Modelling of Particle Velocity and Concentration Distributions in Turbulent Flow: The Mesoscopic Eulerian Formalism and the Two-Point Probability Density Function Method[END_REF][START_REF] Squires | Preferential Concentration of Particles by Turbulence[END_REF] although characterized differently. In a recent work, Boutsikakis (2020) performed several additional DNS simulations with a larger range of the Stokes number and observed the same tendencies on the coefficients evolutions. He also noticed that, below a certain value of the particle Stokes number, the α and β coefficients respectively increase and decrease, and seem to tend toward the values of an uniformly random distribution as the Stokes number tends to zero. However that may be, using the α and β coefficients seems appropriate to account for the clustering effect when considering a three-dimensional field of droplets. A last comment is about the sampling lengths, which provide similar results. Although no other sampling length has been tested, it is our belief that there is a range of sampling lengths that provide valuable results, while using too small or too large lengths may produce poor or inconsistent results [START_REF] Monchaux | Analyzing Preferential Concentration and Clustering of Inertial Particles in Turbulence[END_REF].

Finally, if we refer to Fig. 6, we observe that the experimental data is above the Monte-Carlo simulation with the present procedure. Consider now two three-dimensional droplet fields with the same density number but associated respectively with a preferential segregation effect and the uniformly random distribution. Based on the analysis of the previous paragraph, the α coefficient of the first case with preferential segregation is lower than the α coefficient. Applying our experimental procedure to these two fields and assuming that the restriction errors is similar in both cases, the projection error will decrease their respective mean minimal distances. Thus we can expect that the mean minimal value as well as the α coefficient in the first case will still be lower than their counterparts in the uniformly random case. This is in contradiction with the plot of Fig. 6 and tend to prove that the clustering effect does not seem to play an important role in our experiment.

CONCLUSIONS

A new experimental setup (Prométhée-LACOM) has been designed to study spray combustion and provide an experimental database to acquire a better knowledge of such flows and for further assessment of the numerical approaches. The setup consists of a trapezoidal bluff body located inside a square-channel duct. The upstream air flow is turbulent and liquid n-decane is injected at the rear face of the bluff body with a flat-fan atomizer. Operating conditions reproduce partially the conditions met inside an aero-engine combustor.

In this work, a method is proposed to estimate the interdroplet between nearest-neighbor particles from Mie scattering images and to characterize the spray spatial distribution, with a linear relationship between the square root of the mean density number and either the mean interdroplet distance or its standard deviation. The whole process has been applied to DNS/DPS simulations on one side, and to experimental data from Prométhée-LACOM setup on the other side.

Measurements have been done for both nonreacting and reacting cases to describe the gaseous and the liquid phases. In the paper, emphasis is laid on the spray distribution in the reacting case. At first, using droplet Mie scattering helped to estimate as local as possible values of the droplet density number in six regions of the flow and to derive values of the interdroplet distance between nearest-neighbor particles. This distance is believed to be an important parameter since it is often used to characterize droplet interactions by comparison with the droplet diameter or the droplet flame diameter. Further data analysis using Monte-Carlo simulation also shows that droplets are randomly dispersed according to the uniformly random distribution. Comparison with the regular distribution also demontrates that this distribution, although currently used in deriving some models, is not appropriate. Such a result could be useful to specify a characteristic length scale related to the interdroplet distance and more representative of the spray distribution. Ongoing work is done to use the available nonreacting data, further assess these initial results, and extend the experimental database. The six DNS/DPS simulations consider monodisperse solid particles in an isotropic homogeneous turbulence field with six different particle Stokes numbers St. Depending on its value, the particles may or may not experiment preferential concentration effects. Spray analysis can exhibit a dependency between the spray distribution and the Stokes number due to the preferential concentration phenomenon. It also proves that the preferential segregation effect is not important in the present experimental data, thus confirming that the uniformly random distribution is appropriate to describe the flow under consideration.

APPENDIX A.

Consider a two-dimensional spray with a mean droplet number density n and assume that its spatial distribution is randomly uniform. It obeys the Poisson law:

P (λ, k) = λ k k! exp(-λ)
where the parameter k corresponds to the number of droplets enclosed in a surface and λ is the product of the mean density and the surface. Considering a droplet, the probability H(D i ) that its nearest neighbor is located at a distance D i is equivalent to the product of the probability of having no droplet enclosed in the surface πD 2 i times the probability of having a droplet between D i and D i + dD i . Thus we have

H(D i ) = 2πnD i exp -πnD 2 i
By definition, the mean value of the nearest neighbor is

E(D i ) = ∞ o D i 2πnD i exp -πnD 2 i dD i
Using the intermediate variable u = πnD 2 i and after some algebraic operations, the mean value expression is

E(D i ) = 1 √ π n-1/2 Γ 3 2
By definition, the variance of this variable is

Var(D i ) = ∞ o [D i -E(D i )] 2 2πnD i exp -πnD 2 i dD i = 1 πn Γ(2) -Γ 3 2 2 Γ(1)
or, equivalently, using the standard deviation,

σ(D i ) = 1 √ πn Γ(2) -Γ 3 2 2 Γ(1) 1/2
These formulas use different values of the Gamma function: Γ(1) = 1, Γ(3/2) = √ π/2, Γ(2) = 1. By replacing these values inside the relationships, we finally obtain Eq. ( 14).

For the three-dimensional case, refer to [START_REF] Chandrasekhar | Stochastic Problems in Physics and Astronomy[END_REF].

FIG. 1 :

 1 FIG. 1: View of the test setup

  FIG. 2: Setting and digital timing diagram for the synchronization of OH-PLIF and Mie scattering imaging

FIG. 3 :

 3 FIG. 3: Simultaneous PLIF-OH and Mie images

TABLE 1 :

 1 Centers and sizes of the analysis windows

Window Coordinates (mm) Window size (µm 2 )

	I II III IV V VI	(17, 52) (31, 63) (49, 76) (19, 39) (31, 41) (49, 28)	993 × 993 3319 × 3319 5673 × 5673 2837 × 2837 1986 × 1986 5673 × 5673

TABLE 3 :

 3 Influence of the size of an analysis window: window I

	Window size (µm 2 )	Nb. samples	n (cm 2 )	E [D i ] (µm)	SD [D i ] (µm)
	709 × 709 993 × 993 1390 × 1390	20,270 38,000 74,200	2016 1929 1920	183 178 171	64 61 58

TABLE 4 :

 4 Estimation of the interdroplet distance errors

Window Projection error Restriction error Total error estimation

	I	-45%	+10%	-39%
	II	-22%	+8%	-16%
	III	-4%	+21%	+16%
	IV	-20%	+11%	-12%
	V	-31%	+10%	-23%
	VI	-12%	+7%	-6%

ACKNOWLEDGMENTS

The work has been realized thanks to the partial financial support of the Région Midi-Pyrénées. The authors would also like to thank F. Bigot for his invaluable involvement in the manufacturing of the Prométhé-LACOM setup.


Jun Jin

Shan Jiang

Hongping Hu

Lamin Zhan

Xiaohong Wang

Vincent Laude

Affiliations

Acousto-optic cavity

published or not. The documents may come   L'archive ouverte pluridisciplinaire

2

A prerequisite for the simultaneous localization of light and sound is the presence of both photonic and phononic band gaps, ideally complete, within which the propagation of light and sound is forbidden whatever the polarization and the wave vector. Compared with their 3D counterparts, 7 1D and 2D PxCs have received the most attention since they are easier to model as well as to manufacture. We focus on 2D infinite PxCs in the following.

As observed in previous related works, [START_REF] Joannopoulos | Photonic crystals: molding the flow of light[END_REF]16 for the photonic band gaps in the 2D case, connectivity of high-ε region (ε denotes dielectric permittivity) is conducive to transverse electric (TE) band gaps, whereas isolated patches of high-ε regions lead to transverse magnetic (TM) band gaps. For instance, TM band gaps prefer lattices of air holes in a dielectric matrix and TE band gaps favor lattices of dielectric columns in air. Therefore, the complete photonic band gap is the result of a compromise between two conflicting objectives. On the other hand, for holey PnCs, heavy masses and slim connections are generally necessary for the unit cell to open a wide band gap for both in-plane and out-ofplane modes. [17][18][19][20] As a whole, the optimal strategy of PxCs design is to find out a structure that possesses all these geometrical attributes, or at least achieves a good balance between them.

One should also be aware of some obvious distinctions between PtCs and PnCs. For example, the triangular lattice of circular air columns in a dielectric matrix is not appropriate for the generation of phononic band gaps, whereas it is the most widely used in PtC devices. As another example, the width of photonic band gaps is conditioned by the contrast in the refractive index between the air holes and the matrix. The choice of material is hence vital for PtCs. The case of PnCs is different, however, since the presence of air can be neglected and elastic contrast is given by the traction-free boundaries of the holes. Such boundary conditions are independent of the choice of the material.

During the past decades, many interesting PtCs and PnCs have been reported for efficiently modulating electromagnetic, acoustic and elastic waves. The simultaneous existence of photonic and phononic band gaps have also been disclosed in a variety of PxCs with different materials, 10,[21][22][23][24][25] periodicity, 14,23 dimensions, 7,9,15,[26][27][28] and topologies 8,29- 32 . However, even though abundant types of PxCs have been considered in the literature, it is still a challenge to broaden simultaneously photonic and phononic band gaps according to practical application demands.

Because of technological constraints of fabrication at the micro-or nanoscale, the unit cell is quite often limited to a topology with simple shapes. A too complicated topology means that it could be difficult to fabricate even though the band gaps are very wide. Moreover, other parameters like the tolerable amount of feature rounding and fillet sizes are more important. Hence the smallest feature width is extremely challenging even for cutting edge nano-fabrication techniques.

For two-dimensional (2D) PxCs, several works have reported the simultaneous existence of photonic and phononic band gaps. Maldovan and Thomas 21 reported that phoxonic band gaps can be obtained in a 2D square or hexagonal lattice crystals made of air holes in a silicon matrix.

Sadat-Saleh et al. 22 and Bria et al. 23 reported the influence of the choice of lattice and filling fraction on band gaps in 2D systems, considering LiNbO3 (lithium niobate) for the former and sapphire/silicon for the later. Recently, Dong et al. 30 investigated the opening of phoxonic band gaps by topology optimization. On the one hand, simple conventional topologies, like circular holes in a square lattice, [21][22][23] can hardly open large simultaneous band gaps. On the other hand, topologies obtained from topology optimization are always complex and quite difficult to manufacture, even though they theoretically achieve large simultaneous band gaps.

In a recent work, we proposed a convex-concave combination in a 2D square lattice PnC with a cross-like hole in the center of the unit cell surrounded by four square holes at the corners, 20 whose most notable advantage is the generation of ultra-wide band gaps in such a simple topology. In this work, we further investigate the optical characteristics of the convex-concave structure. A comprehensive search of complete dual band gaps is conducted. In addition to the broad phononic band gaps, large photonic band gaps are also identified. The optimal gap-tomidgap ratios for the complete photonic and phononic band gaps are 11.5% and 90.7%, respectively. In summary, the following advantages of the proposed PxC are demonstrated: (i) simultaneous broad complete photonic and phononic band gaps; (ii) broad band gaps that can be obtained in a wide range of geometrical parameters which provide good tolerance for the technological realization and flexibility in designing.

These traits are beneficial to improve feasibility, applicability and reliability, making the proposed PxC as a potential candidate for practical applications.

Furthermore, with large complete band gaps, the structure is suitable for AO devices due to the following two reasons: (i) More pure defect modes can be generated within the large band gaps; (ii) Photonic defect modes are unlikely to move out of the band gap even under strong AO coupling. Numerous approaches have been proposed to study AO interaction in phoxonic crystals with defects or cavities. Transfer matrix method, 33 and layer-multiple-scattering method, 34 were employed to investigate 1D PxCs with defects respectively. Transmittance and reflection of optical and acoustic defect modes were calculated to search the shift of optical defect frequency under AO coupling. 35 1D PxCs with defects were also studied by a Born series approach with a Green's function to derive the expression of the electric field. 35 It is worth noting that the first-order and second-order effects can be determined by firstorder and second-order Born series, respectively. For 2D and 3D PxCs with cavities, the linear AO coupling rates of the photo-elastic (PE) effect and the moving interface (MI) effect were obtained from first order perturbation theory. 6,[36][37][38] The equation to obtain the quadratic AO coupling was given 39 and the case of degenerate modes was further included. 40 In optomechanics, the shift of the optical frequency is given by the product of the AO coupling rate and the zero-point fluctuation amplitude of the mechanical oscillator xzpf that depends on both the effective mass and the acoustic frequency. 41 Obviously, the effective mass depends on the thickness of the structure. Hence the results of the 2D model have to be compensated for by assuming a given thickness. Furthermore, in order to ensure the safety of the structure, a maximum strain in the cavity was limited to 1% or 0.1%. This can be achieved by setting the maximum modal displacement. Certainly, for PE and MI effects, the maximum mode displacement needs to be set the same in the calculation of optical frequency shift.

In this work, owing to the large complete band gaps, the AO cavity coupling in the 2D PxC with combined convex and concave holes is further investigated. The AO couplings for the PE and the MI effects are studied by the first-and second-order perturbation methods and the finite element method (FEM). An appropriate normalization on displacement is performed to estimate numerically the PE and the MI effects. Whether the acousto-optic coupling belongs to linear coupling or quadratic regime is further distinguished by FEM analysis. Moreover, optomechanical coupling rates relative to the zero-point motion are studied by 2D and 3D PxC models.

II. PXC MODEL AND THEORETICAL BASIS

As described in Fig. 1 (b), a 2D square lattice PxC with a convexconcave topology is proposed by considering a cross-like hole in the center of the unit cell and 4 square convex holes at the corners.

Alternatively, the structure can also be regarded as a central square surrounded by 4 L-shaped holes at the corners, as depicted in Fig. 1 (c).

In either case, the unit cell is formed by four lumps connected by four Lshaped connectors. Three independent parameters (b, c, d or their derivatives) are required to determine the structure. The cross-like hole in the center is determined by b and c, while the side length of the square hole at each corner is d/2. The width of the connections is thus h = (a-c-d)/2. All these parameters are normalized against the lattice constant a in the following calculations. The chosen material is single crystal silicon, with refractive index n = 3.6, 30 where  is the gradient operator, ω denotes the angular frequency, and For PnCs the forming mechanism of large band gaps with combined convex and concave holes was investigated previously. 20 Regarding

III. OPTIMIZED STRUCTURE FOR LARGE COMPLETE PHOXONIC BAND GAPS

PtCs, TM and TE polarizations need to be discussed separately. For the TM polarization, the magnetic energy of mode M1 is mainly confined within the vacuum holes (low-ε) regions. On the contrary, mode M2 is mainly concentrated in the solid dielectric (high-ε) regions. Therefore, for the adjacent modes, the dramatic change of the magnetic energy in the high-ε regions results in a large TM band gap. For the TE polarization, the electric energy of both modes E3 and E4 distributes mostly in solid dielectric (high-ε) regions. However, mode E3 concentrates on lumps whereas E4 concentrates on L-shape connections.

Thus, the difference in the distribution of electric energy in the high-ε 0 6 regions can explain the formation of a large TE band gap. As a result, the upheaval of magnetic or electric fields in the high-ε regions between consecutive bands forms a band gap. Moreover, a larger dielectric difference between regions opens a wider band gap. This result is in agreement with the previous idea that connectivity of high-ε regions is conducive to TM gaps whereas isolated patches of high-ε regions lead to TE gaps.

A thorough geometrical optimization was carried out to reveal the relationship between band gaps and geometrical parameters. Generally, it is preferable that the topology involves several geometrical parameters, offering more possibilities to tune band gaps. For the proposed convexconcave PxC, three independent parameters need to be considered. We take b, c and w (w = a-2h) as those independent parameters. Parameter w is introduced as the void ratio in the x or y direction and is not shown in

IV. PHOXONIC CRYSTAL CAVITY DEFECT MODES

Once large complete phononic and photonic band gaps have been achieved simultaneously, a phoxonic defect cavity can be introduced.

The lattice constant a is fixed as 875 nm, for optical operation around λ = 1550 nm in silicon. A 7×7 super-cell provides a good compromise between accuracy of the results and reasonable computation time. 6 As shown in Fig. 6 center of the structure, and the photonic modal shapes are displayed using displacements and electromagnetic fields in Fig. 6(e), (f) and (g), respectively. The displacements are highly localized in the cavity region.

The electromagnetic fields also show a high confinement inside the cavity. We notice that among the eleven phononic modes, three pairs are degenerate (b and c, g and h, m and n). The same observation is also made for photonic TE modes α and β. The occurrence of degenerate modes can be explained by the cubic symmetry of silicon combined with the spatial symmetry of the structure.

The symmetry of the phononic modes is indeed one of the key factors governing the acousto-optic coupling strength. Following Ref. 5 , modes are classified as odd (o) or even (e) with respect to the two symmetry axes P and P' shown in Fig. 6(d). Both modes d and f have ee symmetry. In details, the vibration of mode d is a stretch-contract motion where the cavity stretches along symmetry axis P while it contracts along symmetry axis P', and the extension-contraction state alternates once during every acoustic period. Mode f is a breathing mode which displays a breathing motion without too much distortion of the cavity shape.

Besides, the photonic TM modes α' and β' have ee symmetry, whereas modes γ and γ' have oo symmetry. The degenerate photonic TE modes α and β have neither ee nor oo symmetry since they are a mixture; they can be obtained from one another by 90° clockwise and anticlockwise rotations, respectively. Besides, compared with the other symmetries of phononic modes, ee symmetric phononic modes have a stronger effect on the AO coupling. Hence, we focus on the AO coupling with modes d and f in the following.

V. ACOUSTO-OPTIC COUPLING STRENGTH

The frequencies of the photonic TE and TM defect modes sustain a modulation around their initial values with a strength given by the magnitude of the acoustic strain inside the cavity, when that cavity is submitted to a periodic acoustic deformation. The frequency modulation results from the sinusoidal variation of the permittivity that is induced by the acoustic wave perturbation. Two mechanisms are responsible for it:

(1) The moving interface (MI) effect 44 or permittivity variation caused by the moving of the interfaces and (2) The photoelastic (PE) effect [START_REF] Yariv | Optical waves in crystals[END_REF] or bulk permittivity variation induced by strain field. Both effects are taken into consideration in the following calculations.

For Pockels effects, the variation of the relative permittivity Δεij is determined by the strain distribution considered as frozen in time.

         - +  - +  - +  =  = -  =  =  =  = (2)
where pij denotes the photo-elastic constants, Si the acoustic strains written in Voigt notation. For silicon, p11 = -0.1, p12 = 0.01 and p44 = -0.051. 47 Thus, the symmetry of the strain field of the acoustic modes are inherited by the permittivity variation.

A. Coupling with phononic breathing mode f

Generally, the acousto-optic coupling strength is measured by the normalized optical frequency shift induced by a realistic acoustic wave perturbation. The maximum of the deformation umax induced by acoustic wave perturbation is fixed as 1% of the lattice constant a, i.e. umax = 0.01a, which is a reasonable choice that does not exceed strength limit of the material. Fig. 7 displays the optical frequency modulations induced by phononic breathing mode f during one acoustic period (0<Ωt<2π)

with the acoustic angular frequency Ω and time t, which is obtained by FEM. The contributions of PE effect and MI effect are both illustrated.

Whether for PE effect or MI effect, the modulations all display a sine function behavior representing linear coupling, in other words, it corresponds to one-phonon absorption and emission process by a photon.

Specially, the degenerate photonic TE modes α and β display a complete synchronized modulation. The contribution of PE and MI can be inphase and out-of-phase. In most case, the MI effect is dominant except for the phononic TM mode γ'. The coupling for modes α, β, α' and β' is obviously larger than for the other two modes γ and γ'. The reason is that for the former, the electromagnetic field mainly distributes in the silicon area; whereas for the latter, it distributes in both silicon and hole areas.

Besides, the deformation induced by acoustic wave only affects silicon.

The fact that the electromagnetic field and the displacement field highly overlap certainly enhances the coupling strength.

The acousto-optic coupling can also be evaluated from perturbation theory. The maximum deformation is chosen as the perturbation parameter. Based on first-order perturbation theory, the equation for nondegenerate phononic modes corresponding to one-phonon processes is derived as 44 ( ) ( )

PE MI 1 0 2 EE EE      +   = - (3)
   = -   + + + +    ( 4
)
where n is the refractive index. Re{} denotes taking real part. Superscript * represents conjugation. For the TE polarization, ( )

2 4 PE 0 12 d z xx yy E E n E p S S S   = - +  (5)
The MI contribution of the PxC is given by ( )(

)

2 1 2 MI || d E E l    - ⊥  =    -   U n E D (6)
where U is the modal displacement field. A maximum of the deformation ( ) max max u = U is defined as above. n is the outward normal unit vector of the surface. E and D denote electric field and electric displacement field, respectively. The subscripts and ⊥ represent the components parallel and perpendicular to the interface between two dielectric materials, respectively. In our case,

   +    = -   - +     (7)
where

ij ij i i j j EE E E E E    = (8)
This explicit expression can be derived based on Eqs.( 4), ( 5) and ( 6). E1 and E2 denote the electric fields of the degenerate modes, and Δω (1)+ and Δω (1)-represent their frequency shifts.

The normalized frequency shifts (Δωa/2πc0) obtained by perturbation theory are listed in Table II and agree well with the amplitude of the modulations in Fig. 7. A symmetry analysis further clarifies the result of perturbation theory. Because the deformation of phononic breathing mode f in x and y directions is in-phase, the superposition of the integral in both directions is enhanced whatever the symmetry of the photonic mode. Linear coupling is absolutely dominant.

For photonic degenerate modes α and β, Γ11 = Γ22 and Γij = 0 (i ≠ j) because of rotational symmetry, which results in a synchronized modulation according to Eq. ( 7). It is noteworthy that the degenerate perturbation theory can be reduced to non-degenerate perturbation theory when Γij = 0 (i ≠ j).

B. Coupling with phononic extension-contraction mode d

The optical frequency modulations of the six photonic defect modes coupled with phononic extension-contraction mode d are shown in Fig. 8. The modulation of the degenerate mode pair α and β still shows a sine function behavior, however, it is not synchronized, just as listed in Table II. This fact can be explained according to Eq. ( 7): we have Γ11 = -Γ22 and Γij ≠ 0 (i ≠ j), therefore, the modulations have opposite sign.

Nevertheless, the other four photonic modes, γ, α', β' and γ', show squared sine function behavior. The reasons are analyzed as follows: The phononic extension-contraction mode d has out-of-phase deformation in x and y directions. These four photonic modes all have either ee or oo symmetry so that the term |E| 2 always has an ee symmetry. The integrals in the x and y directions cancel each other according to Eqs.( 4), ( 5) and (6). As a result, between each one of these four photonic modes and phononic mode d, linear coupling vanishes, and quadratic coupling is dominant over higher-order nonlinear effect. It is worth noting that quadratic coupling though weak has been proposed as a means for realizing quantum nondemolition (QND) measurement of the phonon number because of the square dependence with the motion displacement. [48][49][50] In order to distinguish whether coupling is linear or quadratic, we further investigate the dependence of acousto-optic cavity coupling on the acoustic displacement by FEM. The normalized frequency shift of PE and MI as a function of the maximum of the deformation umax is depicted in Fig. 9. Normalized frequencies of modes α and β both shift proportionally to the maximum of the deformation. Thus, these acoustooptic couplings are linear. Furthermore, for umax/u0 = 1, the linear coupling hypothesis agrees with the conclusion drawn by the perturbation theory, which is shown by optical frequency modulation amplitude in Fig. 8(a) and is also listed in Table II. In contrast, for the PE effect of modes γ and β', and the MI effect of α', β', normalized frequencies shift in proportion to the square of the maximum of the deformation. Therefore, these acousto-optic couplings are quadratic.

However, a very fine mesh needs to be applied in FEM analysis and the MI effect of modes γ and γ' is still unstable even under a very small maximum of the deformation. Thus, we infer that a higher order acousto-optic coupling exists in the MI effect of photonic modes γ and γ' induced by phononic extension-contraction mode d.

Based on second-order perturbation theory, the equation for nondegenerate phononic modes corresponding to two-phonon processes is derived as

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 3
       =  +    +  - -  j ii ii ii d j k i i j i k ji j j i i EE EE EE E E E E (9)
Whereas for twofold degenerated modes,

+    +  - -  j ii ii ii d j k i i j i k ji j j i i EE EE EE E E E E (10)
Furthermore, Table III lists the normalized frequency shifts of quadratic coupling obtained by FEM and 2nd-order perturbation theory.

It is worth noting that 2nd-order perturbation theory has a large error compared with FEM in this case.

C. Coupling rates relative to zero-point motion

In a quantum mechanical approach, the acousto-optic coupling strength is quantified by the frequency shift imparted by the zero-point motion of the mechanical resonator. 1,2,[51][52][53][54] For linear coupling, g (1) = g'xzpf. For quadratic coupling, ( )

2 2 zpf g g x  =
where g and g  denote the coupling coefficients and , where the thickness dt is fixed as 220 nm. The coupling rates for six photonic defect modes coupled with phononic modes d and f obtained by FEM are shown in Table . IV.

Actually, these results can be obtained through a linear or quadratic transformation from the results of the linear or quadratic effect in Subsections A and B of Section V, respectively. Owing to its optical field being highly confined in silicon and overlapping efficiently with the displacement field, mode β' has both maximum linear and quadratic coupling rates. The largest linear coupling rate g (1) /2π is about 760 kHz for TM mode β' coupling with phononic breathing mode f. The largest quadratic coupling rate g (2) /2π is about 0.01 Hz for TM mode β' coupling with stretch-contract phononic mode d. The corresponding coupling coefficients g' and g'' are 181 GHz/nm and 643 MHz/nm 2 , respectively. phononic modes and photonic modes are obtained. Considering both the photoelastic and the moving interface mechanisms, the acousto-optic coupling was evaluated and the effect of symmetry on acousto-optic coupling was analyzed. The coupling strength was also investigated by FEM and first-order perturbation theory. The results of both methods are in good agreement for linear coupling, while second-order perturbation theory needs a correction to estimate quadratic coupling. Finally, the optomechanical coupling rate relative to zero-point motion was investigated by FEM. A photonic TM mode is found to have both the largest linear and quadratic coupling rates, in accordance with the optical field being strongly confined in the central silicon defect.

  the mass density ρ = 2320 kg/m 3 and the elastic constants c11 = 165.7 GPa, c12 = 64.1 GPa and c44 = 79.6 GPa. Transverse and longitudinal speeds of sound are ct=5840 m/s and cl=8433 m/s, respectively. The refractive index n of air in the holes is set to 1.For electromagnetic waves propagating within the considered 2D infinite PxCs in the xy-plane, Maxwell's equations can be decoupled into two scalar wave equations for transverse electric (TE) and transverse magnetic (TM) polarizations[START_REF] Born | Principles of optics: electromagnetic theory of propagation, interference and diffraction of light, Seventh[END_REF][START_REF] Yariv | Photonics: optical electronics in modern communications[END_REF]

FIG. 1 .

 1 FIG. 1. Cross-section of the proposed 2D square-lattice PxC crystal: (a) an extended unit cell shows the arrangement of concave and convex holes; (b) and (c) show two possible equivalent primitive unit cells, as outlined by the blue and the red-dotted lines in (a). In (b), the primitive cell is centered on the concave cross-like hole surrounded by four square holes placed at the corners. In (c), the primitive cell is centered on the convex square hole surrounded by four L-shaped holes placed at the corners. In (b) and (c), the red and blue areas both represent lumps and L-shaped connections.

  c0 represents the speed of light in a vacuum. The distribution of the dielectric constant is periodic in the xy plane and uniform in the z direction. It satisfies ε(x) = ε(x + a), where x= (x, y) is the position vector and a is a lattice vector. The electric and magnetic fields in the z-direction Ez and Hz satisfy the Floquet-Bloch conditions Ez(x+a) = exp(ik•a)Ek(x) and Hz(x+a) = exp(ik•a)Hk(x). Ek and Hk are cell periodic fields and k is the Bloch wave vector restricted to the first Brillouin zone. 'i' is the imaginary unit. Band structures are calculated by the finite element method (FEM). Floquet periodic boundary conditions are imposed on pairs of opposite external boundaries of the unit cell. All solutions are obtained by sweeping the wave vector k along the edges of the first irreducible Brillouin zone. In the band structures, normalized frequencies (ɷa/2πc0 or ɷa/2πct) are functions of the reduced wavenumber (ka/2π). The width of a band gap can be measured using the gap-to-midgap ratio BG% = 200% (ftop-fbot)/(ftop+fbot), a non-dimensional parameter that avoids frequency dependence, where ftop and fbot are the bounding frequencies of the band gap. Generally, a higher BG% value means a better performance.

Figure. 2

 2 Figure. 2(a) and (b) depict phononic and photonic band structures for the proposed PxC, with geometrical parameters b/a = 0.58, c/a = 0.32 and d/a = 0.62. Large complete phononic (90.7%) and photonic (11.5%) band gaps are visible. The complete photonic band gap for normalized electromagnetic frequencies ɷa/2πc0 ranges from 0.549 to0.616, and the complete Phononic band gap for normalized elastic frequencies ɷa/2πct is between 0.358 and 0.952. It is remarkable that although photonic and phononic band gaps in PxCs occur at similar wavelengths of the order of twice the crystal periodicity a, the frequencies of the photonic and phononic gaps are very different since light and sound propagation

  velocities in solids differ by several orders of the magnitude. The dimensionless representation of the dispersion curve makes it easy to choose the lattice constant based on wavelength and frequency. For instance, if a PxC displays a photonic band gap with midgap frequency at about 200 THz (telecom wavelength λ = 1550 nm), the lattice constant should be a = 875 nm. The midgap frequency of the phononic band gap is then 4.39 GHz. A comparison of dual complete band gaps with published 2D PxCs is summarized in Table I to illustrate the advantages of the proposed topology. For each case, Table I lists the filling fraction of solid material, as well as the thickness of the connections h which is the smallest feature limiting the fabrication. For clarify, except for the data directly quoted from Refs. 21,28 , values are obtained considering silicon to facilitate comparison. Obviously, the proposed model has significant advantages over the first four kinds of PxCs regarding complete band gaps.

FIG. 2 .

 2 FIG. 2. (a) PnC and (b) PtC band structures of the proposed PxC and the corresponding eigenmodes at edges of the band gaps, for parameters b/a = 0.58, c/a = 0.32 and d/a = 0.62. In (a), the color of bands represents the kinetic energetic ratio ex of the x-polarization of the elastic wave. The blue curves represent bands of out-of-plane mode, and the rest curves are bands of in-plane mode. In (b), the blue and red lines denote the TM and TE polarizations, respectively.

Fig. 1 .FIG. 3 .FIG. 4 .FIG. 5 .

 1345 Fig. 1. The width of L-shaped connections h is critical for the generation of band gaps and is the smallest feature size limiting fabrication. However, in order to display the geometric parameter range corresponding to the maximum bandwidth on the front view of Fig. 3 (a) and (b), the ordinate represents w instead of h.

Figure 5

 5 Figure 5 illustrates the formation of large complete phoxonic band gaps. Fig. 5(a) and (b) shows how the locations of phononic and photonic band gaps vary with d/a, when b/a = 0.58 and c/a = 0.32 are fixed. In-plane (or TM) and out-of-plane (or TE) band gaps are represented by blue and magenta regions, respectively. Complete band gaps appear when they overlap. Comparatively, Fig. 5(b) shows a much more barren terrain.

Fig. 6 .FIG. 6 .

 66 photonic and phononic defect modes in the band gap range are shown in Fig.6.The phononic and photonic dispersion relations are depicted only along the ΓΧ direction of the irreducible Brillouin zone, since the bands associated with defect flat modes are flat and independent of the wave vector. One can notice that some new flat bands are born inside both phononic and photonic band gaps after introduction of the defect cavity into the super-cell. The flat in-plane defect bands include eleven

FIG. 7 .

 7 FIG. 7. Optical frequency modulation of (a) TE modes α, β and γ, (b) TM modes α´, β´ and γ´, with reduced frequency during one period of PnC breathing mode f. The moving interfaces (MI, red-dashed lines), the photoelastic (PE, green-dotted lines), and the full acousto-optic coupling (PE+MI, blue-solid lines) are represented separately.

  degenerate modes α and β, the frequency shift of a twofold degeneracy is obtained by degenerate perturbation theory.40

FIG. 8 .

 8 FIG. 8. Optical frequency modulation of (a) TE modes α, β and γ, (b) TM modes α´, β´ and γ´, with normalized frequency during one period of PnC extensioncontraction mode d. The moving interfaces (MI, red-dashed lines), the

  represents the effective motional mass, Ω is the angular frequency of the acoustic mode. Since the PnC defect modes are all in-plane modes in our 2D model, the effective motional mass can be rewritten as

FIG. 9 .U

 9 FIG. 9. Optical normalized frequencies shift as a function of the maximum of the deformation induced by acoustic extension-contraction mode d. (a) PE, (b) MI. The results are calculated by FEM, where u0 = 0.01a is fixed. For phononic breathing mode f, the effective motional mass and zeropoint motion are 129 fg and 4.21 fm, respectively. As for phononic extension-contraction mode d, the corresponding values are 178 fg and 4.04 fm, respectively. Comparing with a full 3D model with

TABLE I .

 I Comparison of maximum relative band gap widths for the proposed model and some published structures. h denotes the smallest size and Vf is the filling fraction of the solid material.

	Topology optimization 30	n.a.	94.2%		14.7%	>0.5
	This work	0.03	90.7%		11.5%	0.347
	PxCs	h/a	PnC	Complete BG%	PtC	Vf
	Square lattice of circular air holes 21	1/40	44.5%	No complete band gaps 0.292
	Triangular lattice 11	0.01	No complete band gaps	3.93%	0.360
	Honeycomb lattice 22	<0.01	69.9%	No complete band gaps	n.a.
	Ref. 23	1/25	38.1%		8.32%	0.336
	Dielectric circles connected by veins 29	0.03	83.8%		11.3%	0.307

TABLE II .

 II Normalized frequency shift (×10 -3 ) of linear coupling obtained by first-order perturbation theory.

	PtC modes	0 c PE  2 a  	0 c MI  2 a  		0 c PE&MI 2  	a
		Induced by breathing mode f
	TE α	-0.28861	3.3206	3.0320
	TE β	-0.28865	3.3206	3.0320
	TE γ	-0.30014	1.1909	0.89073
	TM α'	1.2041	3.3670	4.5711
	TM β'	0.85570	3.6800	4.5357
	TM γ'	2.4520	-0.30605	2.1460
		Induced by extension-contraction mode d
	TE α	-0.057357	3.2933	3.3506
	TE β	0.057357	-3.2933	-3.3506

TABLE III .

 III Comparison of the normalized frequency shifts of quadratic coupling obtained by either FEM or 2nd-order perturbation theory.

	Nromalized frequency shift				PE					MI			PE&MI
	Δωa/2πc0(×10 -3 )		2 	0		2 	0	2 	0
	PnC-PtC	FEM	2nd perturbation	FEM	2nd perturbation	FEM	2nd perturbation
	d-γ	-0.01686					-0.00741	0.081031			0.00878	0.09568	0.00765
	d-α'	0.005591					-0.065	-0.2159			-0.555	-0.08995	-0.22058
	d-β'	0.1052					0.083	0.4098			0.692	0.142992	0.27526
	d-γ'	-0.05493					-0.00171	-0.00489			-0.00135	-0.05602	-0.00055
						c	a				c	a	c	a

Table V

 V lists maximum linear and quadratic coupling coefficients of some published structures and this work. It can be noted that the structure proposed in this work achieves the largest linear coupling coefficient and the second largest quadratic coupling coefficient. For optical modes coupled with acoustic modes of a few GHz, the coupling rate seems considerable. Nevertheless, the lower acoustic frequency, the larger the zero-point motion amplitude, and the larger the coupling rate, especially for quadratic coupling. A convex-concave 2D square-lattice PxC was designed on the basis of a large phononic bang gap. Based on FEM and geometrical optimization, a comprehensive investigation of band gaps has been performed. The results indicate that the convex-concave PxC can open simultaneously large complete photonic and phononic band gaps, while at the same time preserving a simple topology. Moreover, large complete band gaps can be achieved for a wide range of parameters. Maximum complete band gaps of up to 90.7% for PnC and 11.5% for PtC are achieved after optimization. As a combination of convex and concave holes, the topology constituted by L-shaped connections and lumps is critical to induce dual band gaps. The present study starts a new path for the design of PxCs. Moreover, with large complete band gaps, the structure is a suitable choice to design a PxC cavity. Several defect

	VI. CONCLUSION

TABLE IV .

 IV Linear and quadratic coupling rates g(1) /2π and g(2) /2π.

					PnC-PtC	Linear coupling rate g (1) /2π (kHz)
					f-α	-47	526	478
					f-β	-49	519	470
					f-γ	-52	191	139
					f-α'	174	555	729
					f-β'	150	610	760
					f-γ'	424	-18	404
					d-α	9	-499	-490
					d-β	-14	544	530
						Quadratic coupling rate g (2) /2π (Hz)
					d-γ	-0.0013	0.0064	0.0070
					d-α'	0.0004	-0.0158	-0.0066
					d-β'	0.0077	0.0300	0.0105
	Mode pair	PE	MI	PE&MI	d-γ'	-0.0040	-0.0004	-0.0041

TABLE V .

 V Comparison of maximum linear and quadratic coupling coefficients for the proposed model and some published structures.

		Coupling coefficients
	Phoxonic crystal cavities	Linear coupling g'	Quadratic coupling g''	Methods
		(GHz/nm)	(MHz/nm 2 )
	tilted membrane within an optical cavity 55	0.0028	4.46	FEM and experiment
	Fiber cavity 56	3	20000	experiment
	Snowflake cavity 32	17	n.a.	experiment
	waveguide with air slot 57	36	n.a.	FEM
	paddle cavity 39	n.a.	400	FEM
	This work	181	643	FEM

ACKNOWLEDGMENTS

National Natural Science Foundation of China (NSFC) (11872186); Fundamental Research Funds for the Central Universities (HUST: 2016JCTD114).

DATA AVAILABILITY.

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

DISCLOSURES.

The authors state no conflict of interest.