Enora Belz
email: enora.belz@univ-rennes1.fr

Estimating Inequality Measures from Quantile Data

Keywords: inequalities, income, distribution, aggregated data, Lorenz curve, Gini, Pietra, Theil

This article focuses on the problem of dealing with aggregate data. It proposes an innovative method for modelling Lorenz curves and estimating inequality indices on small populations, when (only) quantiles are available. When dealing with small population areas and due to privacy restrictions, individual or income share data are often not available and only quantiles are reported. The method is based on conditional expectation in order to find the different income shares and thus model a Lorenz curve with the functional forms already proposed in the literature. From this Lorenz curve, inequality indices (Gini, Pietra, Theil indices) can be derived. A simulation study is performed to evaluate this method and compare it with the other methods used. An example based on real Parisian data is presented to illustrate the method.

Introduction

A suitable and widely used approach to depicting economic inequalities is to provide indices to measure the degree of inequality. These measures enable comparisons of living standards between different countries, regions or among time. Descriptive measures could be useful in understanding economic relationship [START_REF] Kaplow | Why measure inequality?[END_REF]. In this regard, several indices have been described in the literature, including the Gini coefficient, the Pietra index or the Theil indexes. Similarly, the Lorenz curve is a relevant indicator of income distribution and most indices are related to it. The modelling of a Lorenz curve then provides a useful tool. Two main strategies have been developed for modelling a Lorenz curve, either by approximating the empirical Lorenz curve or either by modelling an income distribution and deriving the Lorenz curve from it.

The first is focused on the parametric approximation of the size distribution of income. Based on this estimation, a Lorenz curve and the several indices can be derived (see [START_REF] Champernowne | Economic inequality and income distribution[END_REF] for a survey). Several relevant functional forms have been used to describe the income distribution (see [START_REF] Chotikapanich | Modeling Income Distributions and Lorenz Curves[END_REF] for an exhaustive list). The most popular forms are the so-called generalized distribution of the second kind (GB2) proposed by [START_REF] Mcdonald | Some Generalized Functions for the Size Distribution of Income[END_REF] and its special and limited cases. Special and limited cases include Pareto distributions, Log-normal distributions, Gamma-type size distributions (Gamma, Generalized Gamma, Weibull) and Beta-type size distribution (GB2, Singh-Maddala, Dagum, Fisk, Beta distribution of the first and the second kind) (see Kleiber and Kotz (2003) for a survey). Alternatively, an approach is to apply on shares the generalized non-parametric Pareto interpolation technique developed by [START_REF] Blanchet | Generalized pareto curves: theory and applications[END_REF].

The second literature examines the parametric approximation of a Lorenz curve. In the same way, several functional forms have been developed to fit Lorenz curves. The pioneers are [START_REF] Kakwani | On the estimation of Lorenz curves from grouped observations[END_REF] who proposed a functional form consistent with the distribution of Australian data. Many others authors suggested functional forms, in particular, [START_REF] Rasche | Functional forms for estimating the Lorenz curve: comment[END_REF], [START_REF] Pakes | On Income Distributions and their Lorenz Curves[END_REF], [START_REF] Gupta | Functional form for estimating the Lorenz curve[END_REF], [START_REF] Arnold | A class of hyperbolic Lorenz curves[END_REF], [START_REF] Villaseñor | Elliptical Lorenz curves[END_REF], [START_REF] Basmann | A general functional form for approximating the Lorenz curve[END_REF], [START_REF] Ortega | A new functional form for estimating Lorenz curves[END_REF], [START_REF] Chotikapanich | A comparison of alternative functional forms for the Lorenz curve[END_REF], [START_REF] Sarabia | A hierarchy of Lorenz curves based on the generalized Tukey's lambda distribution[END_REF] or [START_REF] Rohde | An alternative functional form for estimating the Lorenz curve[END_REF]. Several methods allow the estimation of the functional form parameters. The earlier models are based on linear or non-linear least squares. [START_REF] Blanchet | Generalized pareto curves: theory and applications[END_REF] and [START_REF] Sarabia | An ordered family of Lorenz curves[END_REF] proposed an alternative approach by considering the median or least median square. More recently, [START_REF] Chotikapanich | Estimating Lorenz curves using a Dirichlet distribution[END_REF] developed a maximum likelihood estimator based on a Dirichlet distribution to capture the cumulative nature of the shares. The selection of the best functional form is subsequently based on an adjustment criterion. [START_REF] Chotikapanich | Averaging lorenz curves[END_REF] suggested an alternative way by averaging functional forms using a Bayesian model averaging approach.

However, problems arise due to the lack of income data at the individual level. Data are usually aggregated and provide less information than individual data. Institutions sometimes report inequality indices with aggregate data, although most only the Gini index. In such cases, an interesting approach is to reconstruct the proposed index and provide alternative indices. Analysis is straightforward when individual data are accessible. The income distribution and the Lorenz curve can be estimated directly from their empirical forms. However, personal income data are not widely available. They are usually aggregated, whether census or survey data. Census data appear to be more relevant to provide an overview over time or across sub-regions. Nevertheless, the institutions provide several different forms of available data. At the national or regional scale, data are mostly reported in the form of income shares. Institutions may provide data as class mean income and deciles or quintiles group shares like on World Income Inequality Database (WIID), World Inequality Database (WID) or even World Bank. The available information thus depicts points on the Lorenz curve. For example, the poorest 10% hold 3.6% of total income in France in 2016 (WIID). When dealing with smaller areas, the available data may be only quantiles. The income quartiles, quintiles or deciles can be provided but not income shares or class mean income. The available information does not depicts points on the Lorenz curve. For example, the poorest 10% receive less than e10,860 in France in 2016 (INSEE). Unfortunately, their share in the total income is not available. The information given is therefore nearby but not precisely the same. In other words, when focusing on a high geographical level (region or country), income shares can be found, although when dealing with smaller areas such as cities, the data provided are frequently reported as income quantiles.

The purpose of this paper is to determine a methodology for estimating a Lorenz curve and associated indices with only quantile data. The methods need to be adjusted to obtain a Lorenz curve from these quantiles. Quantiles can easily be transformed into tabulated data. The classic method with tabulated data is to assume the midpoint of each class as the class mean. This assumption is rather strong. For this reason, we develop an innovative method based on conditional expectations to compute class means. We will illustrate this method by measuring inequalities within the city of Paris. French municipalities can be subdivided into narrower areas called iris (Ilots Regroupés pour l'Information Statistique). The iris generally have between 1800 and 5000 inhabitants and are built in relation to large sections of the urban network. The narrowness of the areas limits the data available due to the privacy of personal data and INSEE (Institut National de la Statistique et des Études Économiques) provides only quartiles and deciles of income. This methodology is also applicable to other quantile data such as ZIP code income data in the United States.

In Section 2, we will present the proposed methodology for modelling a Lorenz curve from quantile data. Section 3 will illustrate the method with an application on simulated income data and Parisian iris data. A R vignette with codes is available from https://github.com/EnoraBelz/Inequality.

Methods

Definition of a Lorenz curve

The Lorenz curve, introduced by [START_REF] Lorenz | Methods of measuring the concentration of wealth[END_REF], is the curve defined by the points (p, L(p)) where p is the cumulative proportion of the income-receiving units sorted from the poorest to the richest and L(p) the cumulative proportion of income received by these units. [START_REF] Gastwirth | A general definition of the Lorenz curve[END_REF] proposed a general definition of the Lorenz curve such that, if X is the income of a member of the population and assumed to be a random variable with cumulative distribution function F (x), quantile function F -1 (x) and mean µ = xdF (x), then the Lorenz curve is the mapping

p → L(p) = 1 µ p 0 F -1 (t)dt (1)
A Lorenz curve will have some expected properties: (1) if p = 0 then L(p) = 0, (2) if p = 1 then L(p) = 1, (3) L(p) ≤ p and (4) L(p) is continuous and differentiable and the slope of the curve increases monotonically.

Several authors developed functional forms to fit a Lorenz curve (see Table 1).

p α e -β(1-p) β > 0, α ≥ 1 Rasche et al. (1980) (1 -(1 -p) α ) β β ≥ 1, 0 ≤ α ≤ 1 Arnold (1986) p(1 + (α -1)p) 1 + (α -1) + β(1 -p) α, β > 0, α -β < 1 Ortega et al. (1991) p α (1 -(1 -p) β ) α ≥ 0, 0 < β < 1 Chotikapanich (1993) e kp -1 e k -1 k > 0 Sarabia (1997) π 1 p + π 2 p α1 + (1 -π 1 -π 2 )(1 -(1 -p) α2 ) 0 ≤ π 1 , π 2 ≤ 1, α 1 ≥ 1, 0 < α 2 < 1 Rohde (2009) p( β -1 β -p ) β > 1

From quantile data to tabulated data

These parametric forms require grouped or individual data which allow the derivation of an empirical Lorenz curve. Quantiles are not sufficient in their current form. Information on the cumulative proportion of the population is given, but the cumulative proportion of income is unknown. Indeed, a quantile indicates the income level such that p% is below this threshold. The information concerning the total income received by p% is not available and therefore the share in the total income of p% is undefined. Nevertheless, quantiles can be transformed into tabulated data. Quantiles become the boundaries of income ranges and the size of the subset becomes the share of individuals. In the case of deciles, the boundaries are the deciles and the proportion of individuals in each range is 10 percent. The use of both quartiles, deciles or other q-quantiles is also feasible. The proportion of individuals depends on being between two same subdivisions or between two different subdivisions. Suppose X a random variable of income and different associated quantiles Q k where 0 < k < K is the rank of the quantile and s k the population share below this quantile so that P(X

< Q k ) = s k . The ranges are therefore [0, Q 1 [ with proportion s 1 for the first range, [Q k , Q k+1 [
with proportion s k+1 -s k from the second to the penultimate range and [Q K , +∞[ with proportion 1 -s K for the last range (see Table 2 for an example).

Table 2. From quantile data to tabulated data with quartiles and deciles From quantile data...

D 1 D 2 Q 1 D 3 D 4 Q 2 /D 5 D 6 D 7 Q 3 D 8 D 9 10% 20% 25% 30% 40% 50% 60% 70% 75% 80% 90%
...to tabulated data

[0, D 1 [ [D 1 , D 2 [ [D 2 , Q 1 [ [Q 1 , D 3 [ [D 3 , D 4 [ [D 4 , Q 2 [ [Q 2 , D 6 [ [D 6 , D 7 [ [D 7 , Q 3 [ [Q 3 , D 8 [ [D 8 , D 9 [ [D 9 , ∞[ 10% 10% 5% 5% 10% 10% 10% 10% 5% 5% 10% 10%

Income shares

The first step is to determine the total income received by each income bin and thereby the shares of income received. This requires determining the mean income within each income bin. In many studies, the mean income of each class is assumed to be the midpoint of the interval and a Pareto-tail for the open-ended interval (Midpoint method). In other words, it means that, somehow, we assume that the income is uniformly distributed in the ranges and consequently the center of mass (i.e. the mean) equals the center of the interval.

This assumption is rather strong. Indeed, this assumes that each individual in a group has the same income level and therefore neglects intra-group variations. In addition, some income levels are more predominant in a population (e.g. minimum income, agreement income).

It can be noticed that quantiles enable to fit an income distribution function. Several different functional forms have been developed to model the distribution of income. [START_REF] Mcdonald | Some Generalized Functions for the Size Distribution of Income[END_REF] suggested the Generalized Beta distribution of the second kind (GB2) which includes most of the previous functional forms as limited or special cases. For this reason, a GB2 distribution is assumed to describe the size of income and the four parameters are optimizied on quantile data. Suppose X a random variable of income with known density f (x), the conditional expectation of being in a range [a, b[ can be determined as

E(X|x ∈ [a, b[) = b a xf (x)dx b a f (x)dx (2)
where f (x) is the density of X. Therefore, knowing the size distribution of income allows to determine the means of the income per bins. The conditional means between each range can be computed by plugging-in the (four) estimated parameters of the GB2 distribution (Conditional Expectation method). Using the class mean income and the class population share , the income share of each bin in the total income can then be determined and consequently the cumulative income shares.

Functional form optimization

Quantile data become similar to grouped data and classical functional forms of Lorenz curves can be applied. In this way, seven different forms are used : [START_REF] Kakwani | On the estimation of Lorenz curves from grouped observations[END_REF] , [START_REF] Rasche | Functional forms for estimating the Lorenz curve: comment[END_REF], [START_REF] Arnold | A class of hyperbolic Lorenz curves[END_REF], [START_REF] Ortega | A new functional form for estimating Lorenz curves[END_REF], [START_REF] Chotikapanich | A comparison of alternative functional forms for the Lorenz curve[END_REF], [START_REF] Sarabia | A hierarchy of Lorenz curves based on the generalized Tukey's lambda distribution[END_REF] and [START_REF] Rohde | An alternative functional form for estimating the Lorenz curve[END_REF]. The parameters are optimized by non-linear least squares (NLS) estimator as

θ = arg min θ K k=1 (F (p k , θ) -s k ) 2 (3)
where F (., θ) is the functional form to be optimized according to the parameter(s) θ, p k is the population share of the k th income range and s k is the income share of the k th income range. Thereafter, the different functional forms are compared according to a goodness-of-fit measure, especially the value of the Chi-squared statistic

χ 2 = K k=1 (s k -F (p k , θ)) 2 F (p k , θ) (4)
where F (., θ) is the functional form with the optimal parameter(s) θ, p k is the population share of the k th income range and s k is the income share of the k th income range.

The algorithm is therefore as follows:

• Convert quantiles into tabulated data

• Estimate a GB2 family distribution on quantiles using ML estimation techniques

• Compute conditional expectations of each bin using the fitted distribution

• Calculate income shares and cumulative income shares

• Estimate a functional form on shares using a NLS estimator

• Compare the functional forms with goodness-of-fit measures.

Inequality measures

The Lorenz curve enables to determine the inequality indices. Gini coefficient G [START_REF] Gini | Sulla Misura della Concentrazione e della Variabilità dei Caratteri[END_REF] can be evaluated as

G = E(|Y -X|) 2E(X) (5)
where X and Y are i.i.d with common distribution F (x). Alternatively, the Gini coefficient can be determined in terms of the Lorenz curve, corresponding to two times the area between L(p) and the egalitarian line as

G = 2 1 0 (p -L(p))dp (6)
As a result, Gini coefficients can be derived from the functional forms as a function of parameters (Table 3).

Table 3. Gini coefficient of the functional forms

Gini Index [START_REF] Kakwani | On the estimation of Lorenz curves from grouped observations[END_REF] 1

-2e -β 1 + α 1 F 1 (1 + α; 2 + α; β) Rasche et al. (1980) 1 -2 α B( 1 α ; β + 1) Arnold (1986)      β β -α + 1 + 2αβ (β -α + 1) 2 [1 + β + 1 β -α + 1 log( α β + 1 )] if β -α + 1 = 0 β 3(1 + β) if β -α + 1 = 0 Ortega et al. (1991) α -1 α + 1 + 2B(α + 1; β + 1) Chotikapanich (1993) (k -2)e k + (k + 2) k(e k -1) Sarabia (1997) π 2 (1 -2 1 + α 1 ) -(1 -π 1 -π 2 )(1 -2 1 + α 2 ) Rohde (2009) 2β[(β -1) log( β -1 β ) + 1] - 1
Note: B is the Beta function. 1F1 is the confluent hyper-geometric function.

The Pietra index P [START_REF] Pietra | Nuovi contributi alla metodologia degli indici di variabilita e di concentrazione[END_REF] is also widely used to measure inequality, it can be defined as

P = E(|X -E(X)|) 2E(X) (7)
The Pietra index can also be defined in terms of the Lorenz curve, it corresponds to the maximum deviation between L(p) and p as

P = max p∈[0,1] {p -L(p)} (8)
In addition, it corresponds to two times the area of the largest triangle that can be inscribed between L(p) and the equalitarian line [START_REF] Arnold | Pareto and generalized pareto distributions[END_REF].

In his book, Theil (1967) also suggests inequality indices, Theil's indexes T L (low) and T H (high), based on the generalized entropy (see [START_REF] Cowell | Measuring inequality[END_REF] for a discussion on generalized entropy measures).

T L = GE 0 = -E(log( X µ
)) ( 9)

T H = GE 1 = E( X µ log( X µ
)) ( 10) [START_REF] Rohde | Lorenz curves and generalised entropy inequality measures[END_REF] relates the concept of Generalized Entropy (GE) to the Lorenz curve and provides mathematical expressions to define the GE measures in terms of the Lorenz curve. 11)

T L = - 1 0 log(L (p))dp (
T H = 1 0 L (p) log(L (p))dp (12)
where L (p) is the first derivative of the Lorenz Curve.

Applications

In a first part, we will simulate samples to evaluate the performance of the methodology and compare it to other methods, then, we will use it to get local inequality indices in Paris.

Simulated data

Simulations are performed both to evaluate the performance of the method and to compare with the midpoint method. The different methods are applied on income distribution simulations. Several known distributions are selected to perform the simulations such as GB2, Log-normal and Singh-Maddala. The parameters of the distributions are determined to ensure realistic income levels. For each simulation, the income of 2000 individuals is simulated according to a defined distribution. Two areas considered, one with a low Gini index and one with a high index. Three distributions2 are used: for the High Gini Index case, LN (10.6, 1.01), GB2(40000, 1.7, 0.98, 1.02) and SM (30000, 1.9, 0.7), respectively with Gini index 0.52, 0.58 and 0.69, and for the Low Gini Index case, LN (10.5, 0.7), GB2(35000, 2.5, 0.95, 1.02) and SM (50000, 2.2, 1.8), respectively with Gini index 0.37, 0.39 and 0.35. For each of them, (i) we generate a sample of size 2000 and calculate Gini index from those individual observations directly -called Ind.

(ii) we can get the shares and partial information (quantiles), and then compute Gini index -called Shares we can use only quantiles -and estimate the shares (iii) the first step of estimating a GB2 distribution allows to compute a Gini index -called First Step (iv) with the conditional expectation method -called Cond. Expectation (v) with the midpoint method -called Midpoint.

Figure 1 provides a comparison of the different methods in terms of estimated Gini index. Method (i) is the one with detailed data, with the other four are based on aggregated quantities. The distribution of the Gini estimate based on shares data (ii) is very close to that observed on the individual data (i). Then, when data are in the form of shares, having aggregate or individual data gives a fairly close result. However, if we make as if the shares are not available (as in the dataset we have) and only quantiles are available, the lack of information leads to estimating a Gini index less close to the true value.

Admittedly, the first step of the Conditional Expectation method (iii) also allows to compute a Gini index. This consists of fitting a GB2 distribution on the quantiles to determine conditional expectations, however, indices and the Lorenz curve can be computed using the GB2 distribution. Nevertheless, the underlying income distribution is not necessarily GB2 and, in addition, the GB2 distribution is a fourparameter distribution estimated on only a few points, and can therefore be misspecified, especially on the distribution tail. The second step (iv) should adjust the method by improving the results. The results suggest that the Conditional Expectation method (iv) computes a Gini index close than the GB2 one (iii). When the underlying distribution is GB2, both methods provide the same value close to the true value. When the The red line corresponds to the true value of the Gini index of the simulated distribution (obtained by 1, 000, 000 simulations). The parameters of the simulated distributions are: High LN (10.6, 1.01), GB2(40000, 1.7, 0.98, 1.02) and SM (30000, 1.9, 0.7) and Low LN (10.5, 0.7), GB2(35000, 2.5, 0.95, 1.02) and SM (50000, 2.2, 1.8).

distribution is not GB2, the results indicate an overestimation of Gini coefficient (or underestimation in the case of Singh-Maddala High Gini Index), however, the Gini index of the second step is more accurate. As a result, the Conditional Expectation method adjusts the estimate of Gini index.

Conditional expectation (iv) and Midpoint (v) methods can also be compared. The Midpoint method underestimates the true values in all cases, while the conditional expectation method only underestimates the true values for Singh-Maddala High Gini Index and overestimates for Singh-Maddala Low Gini Index and Log-normal. If the distribution is GB2 or Singh-Maddala, the Conditional Expectation method significantly outperforms the Midpoint method. When the distribution is Log-normal, the result is less clear. The Conditional Expectation method overestimates the true value while the midpoint method underestimates it. The performance of both methods is then related to the underlying distribution.

Parisian income data

This paper focuses on the measure of inequality within Paris. The Institut National de la Statistique et des Études Économiques (INSEE) provides income data at a very narrow scale named iris (see Figure 2 for a visualisation of median income in the iris). The deciles and quartiles of income3 in these areas are available. Paris includes 966 iris which are either living (861), activity (88) or miscellaneous (17). For some iris with insufficient population (especially activity or miscellaneous iris), income information is not available. Thus, 865 iris are reported with income data. Using these data, a Lorenz curve can then be modelled, the Gini index can be determined and compared with that provided by INSEE, and the Pietra and Theil indices not included can be calculated. The first step is to determine the parameters of the GB2 distribution of each area by MLE (see dotted line in Figure 3). The parameters are used to identify the density and the distribution of income in the area. Conditional means of income bins can then be calculated (see blue lines in Figure 3). Those means are used to compute the cumulative shares of income and the cumulative shares of population are the shares between two quantiles. Table 4 displays an example of the data obtained for an iris ("Bel Air 5" in the 12th arrondissement). In this area, the first decile is at e10,570, implying that 10% of the population have an income below e10,570. In this range, the conditional mean is found to be at e7,403. The cumulative share of income can be derived from the total income (weighted sum of the conditional means). In this area, the poorest 10% of the population receive 2% of the total income. By plotting the cumulative share of income in relation to the cumulative share of the population, the points of the empirical Lorenz curve are observed (Figure 4). To approximate the empirical Lorenz curve, all alternative functional forms can be applied on the cumulative shares of income and population. The parameters are optimized by NLS estimator. For each functional form, the most optimal parameters and Lorenz curve are obtained. The different forms could also be compared with each other using goodness-of-fit measures (see column 2 of the Table 5). Figure 5 represents the best functional form according to the chi-squared of each Parisian iris. 58% of the iris tend to be in favour of the [START_REF] Sarabia | A hierarchy of Lorenz curves based on the generalized Tukey's lambda distribution[END_REF] form, 38% for [START_REF] Ortega | A new functional form for estimating Lorenz curves[END_REF] form and 3% for [START_REF] Rasche | Functional forms for estimating the Lorenz curve: comment[END_REF] form. A rather notable pattern appears, the iris of the center and 16th arrondissements are generally in favour of an [START_REF] Ortega | A new functional form for estimating Lorenz curves[END_REF] form, while the iris of the outlying arrondissements are in favour of [START_REF] Sarabia | A hierarchy of Lorenz curves based on the generalized Tukey's lambda distribution[END_REF] form. This spatial pattern is similar to the one of inequality indices, and hence the functional form choice appears to be related to the income inequalities of the area. The Lorenz curve is used to compute inequality indices (see Table 5 for Bel Air 5 iris). The Gini index can be defined as a function of each functional form parameters. The iris database provides the Gini index calculated from individual data. Then, a comparison between the one found with the functional form and the individual one should be considered. Figure 6 depicts the Gini derived from functional forms as a function of the one given by INSEE and Figure 7 the difference between them for each functional form. The method tends to slightly overestimate the Gini index. The difference remains quite small, in most cases ranging from 0 to 0.02. In addition, inequality indices not provided by INSEE, such as the Pietra or Theil indices, can be calculated according to the best functional form (see Figures 8 for Gini index, 10 for Pietra index, 9 for Theil L index and 11 for Theil H index).

Conclusion

The objective of this paper is to propose an innovative method to model Lorenz curves and estimate inequality indices on small populations, when only quantiles are available. The method is based on conditional expectation in order to find the different income shares and thus model a Lorenz curve with the functional forms already proposed in the literature. Real and simulated data are used to evaluate the proposed method and compare it to other methods used. We note from simulated data that it is more difficult to estimate a Gini index when shares are not available and only quantiles are available. However, the proposed Conditional Expectation method outperforms the traditional Midpoint method. Similarly, the method applied to the Parisian iris data provides a Gini index very similar to the true value. Finally, the proposed methodology enables to model a Lorenz curve and hence to estimate inequality indices with quantile data. Therefore, this method is useful for measuring inequalities when data are limited. This approach can be applied on income data in quantile form. However, it can also be used for data in class form with an underlying distribution to find conditional expectations.

Figure 1 .

 1 Figure 1. Comparison of the different methods in terms of estimated Gini

Figure 2 .

 2 Figure 2. Median income of the Parisian iris

Figure 4 .

 4 Figure 4. Empirical Lorenz curve of Bel Air 5 iris

Figure 5 .

 5 Figure 5. Best choice of functional forms obtained for each Parisian iris

Figure 6 .Figure 11 .

 611 Figure 6. Comparison between estimated and observed Gini index by functional form

Table 1 .

 1 Lorenz curve functional forms

Functional form L(p)

Parameter constraints

[START_REF] Kakwani | On the estimation of Lorenz curves from grouped observations[END_REF]

Table 4 .

 4 Tabulated data obtained using the conditional expectation method for Bel Air 5 iris

		Quantile	Income	Proportion	Mean	Income	Population
		(e)	Range (e)		Income (e)	share	share
	D 1	10, 570	Below 10,570	0.10	7, 403.64	0.02	0.10
	D 2	18, 696	10,570-18,696	0.10	14, 947.50	0.06	0.20
	Q 1	21, 558	18,696-21,558	0.05	20, 135.63	0.08	0.25
	D 3	24, 344	21,558-24,344	0.05	22, 948.68	0.12	0.30
	D 4	28, 298	24,344-28,298	0.10	26, 296.94	0.18	0.40
	Q 2	32, 626	28,298-32,626	0.10	30, 412.82	0.26	0.50
	D 6	37, 782	32,626-37,782	0.10	35, 113.38	0.36	0.60
	D 7	43, 444	37,782-43,444	0.10	40, 488.39	0.46	0.70
	Q 3	46, 088	43,444-46,088	0.05	44, 737.59	0.52	0.75
	D 8	50, 926	46,088-50,926	0.05	48, 410.77	0.58	0.80
	D 9	61, 920	50,926-61,920	0.10	55, 936.68	0.73	0.90
			61,920 and over	0.10	102, 703.30	1	1

Table 5 .

 5 Inequality measures for the different functional forms of the Bel Air 5 iris

	Functional form	χ 2	Rank	Gini	Pietra	T L	T H
	Kakwani and Podder (1973)	0.01028	4	0.345	0.261	0.205	0.185
	Rasche et al. (1980)	0.00154	3	0.351	0.248	0.221	0.214
	Arnold (1986)	0.01493	6	0.342	0.259	0.189	0.184
	Chotikapanich (1993)	0.01102	5	0.344	0.261	0.199	0.184
	Sarabia (1997)	0.00036	1	0.356	0.244	0.24	0.287
	Ortega et al. (1991)	0.00112	2	0.352	0.247	0.225	0.221
	Rohde (2009)	0.02149	7	0.340	0.259	0.183	0.183

Lognormal distribution LN (µ, σ), GB2 distribution GB2(µ, σ, ν, τ ) and Singh-Maddala distribution SM (µ, σ, ν).

Quantiles of declared household income per consumption unit for the year 2014. The used scale (OECD scale) has the following weighting: 1 UC for the first adult in the household, 0.5 UC for other persons aged 14 or over and 0.3 UC for children under 14 years.

Acknowledgements

I would like to thank the Brittany region and the ACTINFO chair for their support. I also thank Olivier L'Haridon, Benoît Tarroux, Arthur Charpentier, Emmanuel Flachaire and the participants of the CNRS Thematic School on Public Policy Evaluation (ETEPP), the annual conference of the Société Canadienne de Science Économique (SCSE) and the R conference in Quebec 2019 for their constructive suggestions and discussions.


T Maquart

T Elguedj
email: thomas.elguedj@insa-lyon.fr

A Gravouil

M Rochette

A note on topological properties of volumes constructed from surfaces

Converting surfaces into a volume has a long interest in several communities, e.g. the computational mechanics community. This process involves having specific surfaces which can be converted into a solid, i.e., a volume. This paper presents in a clear and brief manner the topological properties conserved during surface to volume transformation. We state the limits of this approach if a specific volume structure is required. Volume structures can be a coarse volume organization or meshes. For that purpose, surface manifolds are mathematically turned into a volume manifold. Topological tools are presented to understand which properties are transmitted to the volume and which ones are unset. Developments are submitted both for continuous and discrete manifolds using CW-complexes.

Introduction

Problematic of surface to volume transformation finds its interest in various fields of application. For many years, different approaches have been used to solve this conversion. Despite this particular interest, it currently not exists generic methods to generate a volume with a prescribed simplicial complex connectivity from arbitrary surfaces. Such notes are given to better understand topological aspects relative to this issue. This investigation will serves application fields who apply physics on solid domains [START_REF] Maquart | Trivariate models generation from unstructured surface manifolds for isogeometric analysis -application to reduced order modeling with geometric parameters[END_REF].

Current challenges in computational mechanics and more globally in the physics field are the generation of volume domains which are structured. For instance, pure hexahedral meshes have always been greatly appreciated. They stand for a more homogeneous solution field with less local effects and a more stable computation during large deformations. Structured mesh generation goes with topology. Isogeometric analysis [HCB05, [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF][START_REF] Occelli | Explicit Dynamics Isogeometric Analysis. LR B-Splines Implementation in the Radioss solver[END_REF] is a new method that represents the geometry better than standard meshes. Despite this benefit, isogeometric meshes also need structured domains. Attempts to compute fields with a specific symmetry embedded in volumes are on current research [START_REF] Chemin | Representing three-dimensional cross fields using 4th order tensors[END_REF] to produce hexahedral meshes [RS15, VCD + 16, SRUL16, SVB17].

Furthermore, today's artificial intelligence applications need relevant ordered data. In other words, the discretization of the geometry should be the same in order to learn important information from identical structures [START_REF] Maquart | Trivariate models generation from unstructured surface manifolds for isogeometric analysis -application to reduced order modeling with geometric parameters[END_REF][START_REF] Lu | Adaptive sparse grid based hopgd: Toward a nonintrusive strategy for constructing space-time welding computational vademecum[END_REF].

We present mathematical formulations on topological properties of 2, 3-dimensional manifolds embedded in R 3 . In particular with a goal to convert many 2-dimensional manifolds into a 3dimensional manifold. Using the well known generalized Euler characteristic and references in the literature, our developments will be securely anchored. Close links between surfaces and its associated volume are then determined.

Presented concepts are given without prerequisites, references in the topology field are recommended [HLS00, Hat01, Jub09, SOA07]. In the following we consider only finite and connected volumes constructed from one or more surfaces.

Surfaces to volume conversion

We strive to find topological properties of 3-dimensional manifolds V that respect some constraints inherited from topological and geometrical properties of a surface M , and in the following from many surfaces M . We seek for constraints provided from the surface to be mapped into the volume manifold. Into a volume manifold V , the number of volumes k 3 is defined by equation (1a). Naturally, the considered surface M becomes the boundary of V as formulated in equation (1b).

k 3 = -χ(V ) + k 0 -k 1 + k 2 , (1a)
s.t. {k 0 M -k 1 M + k 2 M = χ(∂V ) = χ(M c ) & |∂V | = 1}. (1b)
χ(V ) is the Euler characteristic of the volume V and ∂V its boundary. k i M are the CWcomplex or simplicial complex entities of the closed surface M c . Likewise, k i are related to V . M c is given with no boundaries, i.e., without 1-dimensional manifolds and is provided from a surface M possibly with boundaries. Since equation (1a) has overabundant unknown terms, and in general not enough constraints in equation (1b), closed surface M c can refer to many associated volumes V with different structures. The constraint |∂V | = 1 if the volume manifold is built from only one compact, connected and oriented 2-dimensional manifold. Hence, for one given closed surface M c , there exists many different volume structures composed by more or less k 3 entities that satisfying above constraints. Depending on the surface features, choosing the best volume structure is done by analyzing geometry. It leads us to the following statement in proposition 1.

Proposition 1 Non-isomorphic associated volumes. Let M c be a closed compact oriented 2-dimensional manifold M . If we want to turn M c into a volume manifold V , there exists a set U V with an infinite size of non-isomorphic volume structures that satisfies constraints provided in equation (1b):

M c ⇒ V i ≡ iso V j , ∀i, j < |U V | = ∞.
(2)

Where ≡ iso denotes the isomorphic equivalence in terms of simplicial complex entities, i.e., related to the number of entities of each dimension. Transformations with a discontinuous function are permitted.

Non-isomorphism between different volumes V in equation (2) differ just by number of volumes k 3 , number of surfaces k 2 and number of edges k 1 for a same χ(V ). In other words, two volumes are isomorphic if they have the same number of simplices of each dimension n, n ∈ {0, 1, 2, 3}, if a mapping with a discontinuous function is allowed. Homeomorphic volumes are isomorphic volumes holding exactly the same simplicial complex connectivity. The set U V has an infinite size as explained in proof 1.

Proof 1 Infinite size of the set U V . The set of non-isomorphic volumes U V of an associated closed surface M c has an infinite size because of the properties of the 3-simplex subdivision. It is always possible to subdivide a n-simplex in the way that its (n -1)-simplex entities remains topologically and geometrically the same; and the created manifold is also embedded in the same n-dimensional topological space, so:

|U V | = ∞, ∀V. (3)
Taking into account previous developments and constraints, we now seek to find which volume Euler characteristic χ(V ) is interesting. In other words, we expect a 3-dimensional manifold embedded in R 3 with specific properties inherited from a surface M or several surfaces M . Furthermore, a volume which lies in R 3 with one ore more (n-1)-dimensional compact, connected and orientable entities is desired. However, we can not accept any 1-dimensional entity, i.e., boundaries of 2-dimensional manifolds. Finally, equation (4) and equation (5) in proposition 2 state our interesting volume Euler characteristics.

Proposition 2 Equation 2χ(V ) = χ(M c ). In R 3 , a 3-dimensional manifold V constructed from a compact, connected and oriented 2-dimensional manifold M c is restricted to the following volume Euler characteristic χ(V ):

χ(V ) = χ(M c ) 2 = 1 -g. ( 4
)
Using g as the genus of the surface M or the closed surface M c , i.e., g(M ) or g(M c ). Indeed, adding 1-dimensional manifold entities to a 2-dimensional one does not change its genus g. If several compact, connected and oriented 2-dimensional manifolds are used to create the volume, the disjoint union of Euler characteristics implies:

χ(V ) = χ(M 1 c M 2 c ... M ns c ) 2 = ns i=1 χ(M i c ) 2 = ns i=1 1 -g(M i c ). ( 5
)
With n s the number of given surface manifolds which serve to generate the volume V . An example is shown in figure 1.

This result can be demonstrated using CW-complex properties and Betti numbers in a topological space R 3 [START_REF] Damiand | Contributions aux cartes combinatoires et cartes généralisées: Simplification, modèles, invariants topologiques et applications[END_REF], see equation (9) in proof 2. Notice that meshes are simplicial complexes due to their combinatorial form. Remember that CW-complexes are themselves both simplicial complexes and continuous entities on which homotopy theory is applicable.

Proof 2 Equation 2χ(V ) = χ(M c ). To turn any 3-dimensional volume manifold V embedded in R 3 into a topological object homeomorphic to the volume ball B 3 , we have to add CWcomplex objects to V . We first create g surfaces k 2 , i.e., for each 1-dimensional tunnel, see e.g. figure 2 (A). Secondly, we create two edges k 1 and one surface k 2 for each 2-dimensional cavity, see e.g. figure 2 (B):

χ(V ) = k 0 -(k 1 + 2b 2 ) + (k 2 + b 1 + b 2 ) -k 3 , (6)
= k 0 S -k 1 S + k 2 S -2b 2 + b 1 + b 2 -k 3 , (7)
= χ(S) + b 1 -b 2 -1, (8)
= χ(S) -χ(V ) ⇒ 2χ(V ) = χ(S). (9)
The Betti number b 3 = 0, because the topological space R 3 is considered. S is the set of closed surfaces,

S = M i c . k 0 = k 0 S , k 1 = k 1 S and k 2 = k 2 S . S remains unchanged. k 3 = 1 implies only one connected volume.
Developments leading to equation (9) establish the topological relation between the associated volume V of many closed surfaces M i c . We can now express in equation (10) the Betti numbers of V using the previous results.

    b 0 = 1, b 1 = b 0 + b 2 -χ(Mc) 2 or b 1 = b 0 + b 2 - ns i=1 χ(M i c ) 2 , b 2 = |M |. (10)
With |M | the number of connected 2-dimensional manifolds used to define V minus one. If only one surface manifold is used to generate the volume, |M | = 0. We can now re-write equation (1) to include interesting volume Euler characteristics and needed topological space constraints inherited from several surfaces. This idea is formulated in equation (11).

k 3 = ns i=1 1 -g(M i c ) + k 0 -k 1 + k 2 , (11a) s
.t. { ns i=1 (k i 0 M -k i 1 M + k i 2 M ) = χ(∂V ) = ns i=1 χ(M i c )} and (11b) {|∂V | = n s & V ∈ R 3 & ∂V ≡ M i c }. (11c)
Where ≡ denotes the classic geometric equivalence: if the n s surfaces M i c are meshes, boundaries of V are. ∂V and surfaces M i c are topologically and geometrically the same. ∂V replicates perfectly the closed surfaces M i c . Notice that, in equation (11b) and equation (11c), constraints are not independent. Since ∂V is provided from the surfaces M i c , |∂V | and χ(∂V ) are already established.

Conclusion

We have presented topological aspects in order to determine properties of volumes constructed from surfaces. We have shown that topology concepts are fundamental features in the transformation of boundary surfaces to a volume. In spite of mathematical tools given above, the problem of converting surfaces into a structured volume is still hard. Indeed, there is many possibilities of volume structures satisfying input information during this kind of conversion. It is easy to see that the interior volume is not defined. A lack of constraints and unset properties are clearly observed. Additional constraints should be added to tackle this uncertainty on the solid structure.

Within the computational mechanics community, researchers are trying to turn scanned meshes into volume meshes. For instance, several techniques exist to convert boundary surfaces into solids, in order to apply physics on them. Presented mathematics state the limits of this kind of approaches.

Figure 1 :

 1 Figure 1: A 3-dimensional manifold V . Three compact, connected and oriented 2-dimensional manifolds M 1 c , M 2 c and M 3 c are used to determine V . According to Betti numbers b i , χ(V ) = 1 -3 + 2 = 0. Using equation (5) the same result is found.

Figure 2 :

 2 Figure 2: 3-dimensional manifold V to volume ball B 3 conversion using CW-complexes. (A) g surfaces k 2 are created for a genus-g surface (red surface), i.e., for each tunnel. (B) Two edges k 1 and one surface k 2 are constructed (red curves and surface) for a cavity. Black and blue objects are not created beacause they are lying on a lower dimensional object than a volume.



Acknowledgments

T. Maquart was partially supported by a CIFRE fellowship of the french Association Nationale de la Recherche et de la Technologie and ANSYS research. T. Maquart, T. Elguedj and A. Gravouil were also partially supported by a research contract with ANSYS. These supports are gratefully acknowledged.


domain: spi phys

text:
Lorenzo Lanzillotta

Olivier Léon
email: olivier.leon@onera.fr

David Donjat
email: david.donjat@onera.fr

Guy Le Besnerais

3D density reconstruction of a screeching supersonic jet by synchronized multi-camera Background Oriented Schlieren

This work describes the application of 3D Background Oriented Schlieren (3DBOS) to the experimental analysis of the mean and instantaneous 3D structures of supersonic under-expanded screeching jets. These jets feature various dynamical states as the jet nozzle pressure ratio (NPR) is varied. Hence, six flow conditions were considered using acoustic measurements and 3D density reconstructions. Relying on the latter and on Proper Orthogonal Decomposition (POD) to isolate the main flow features, clear visualizations of the instantaneous 3D state and the reconstructed average dynamics of modes A1 and C, that are respectively axisymmetric and helical, were obtained.

Introduction

Supersonic axisymmetric jets are commonly found in aerospace and industrial applications, e.g. aircraft engines, space launchers or pressurized tank leakage. Such jets are often non-ideally expanded, such that the flow static pressure at the nozzle exit, p 0 , is greater than the ambient one, p a . The operating condition of these under-expanded jets is defined by the Nozzle Pressure Ratio NPR ≡ p t0 /p a where p t0 is the stagnation pressure at the nozzle exit. Equivalently, the ideally expanded Mach number can be used and is defined as

M j ≡ 2 γ -1 NPR (γ-1) /γ -1 1 /2 ( 1
)
where γ is the heat capacity ratio taken equal to 1.4. The physics of such jets is complex: in order for the flow static pressure to match the ambient one, quasi-periodic shock cell structures are generated, leading to expansion and compression waves trapped into the jet plume as illustrated by a Schlieren visualization in Fig. 1. As the NPR is increased, different regions in the jet may be identified such as Mach disks, barrel shocks or reflected shocks. Such under-expanded jets have been extensively studied in the literature with a particular effort made toward understanding the flow field average characteristics [START_REF] Franquet | Free underexpanded jets in a quiescent medium: A review[END_REF] and dynamics. [START_REF] Daniel Edgington-Mitchell | Coherent structure and sound production in the helical mode of a screeching axisymmetric jet[END_REF] A common feature of underexpanded jets is the presence of an aeroacoustic feedback mechanism yielding acoustic screech tones and strong flow oscillations. As initially described by Powell,[START_REF] Powell | On the mechanism of choked jet noise[END_REF] these tones are the result of the following feedback loop: the coherent structures associated with the jet shear-layer convective instabilities propagate downstream and interact with the shock cells, generating upstream-propagating waves (either freestream acoustic waves or neutral upstream modes inside the jet, see Edgington-Mitchell 8 for a recent review) that, in return, excite the shear layer at the nozzle-lip, a highly receptive region. For an axisymmetric jet, four different oscillation modes were identified by Powell 22 and are commonly referred to as modes A, B, C and D. Using Schlieren visualizations and azimuthal microphone measurements, the general structure and dynamics of such jet plume global instabilities have been discussed in the literature, these four modes being generally associated with axisymmetric toroidal, rotating/precessing flapping, helical, and flapping motions respectively.

Despite decades of research effort and progress, a precise understanding of the mechanisms behind all these oscillation modes is still lacking. [START_REF] Edgington-Mitchell | Aeroacoustic resonance and self-excitation in screeching and impinging supersonic jets -a review[END_REF] For example, open questions remain regarding the exact nature of the screech noise production process (that is the interaction of the downstream propagating instabilities and the shock cells), the possibility of coexistence of two screech modes, or the nature of the upstream waves. Regarding the latter, recent studies conducted by Gojon et al [START_REF] Gojon | Oscillation modes in screeching jets[END_REF] and Edgington-Mitchell et al [START_REF] Daniel Edgington-Mitchell | Upstream-travelling acoustic jet modes as a closure mechanism for screech[END_REF] suggest for example that modes A and C can be satisfactorily modeled and explained relying on a simple vortex-sheet model where the screech feedback loop is sustained by upstream neutral acoustic wave modes, whereas the most commonly acknowledged mechanism relies on upstream near-field acoustic waves.

These questions are difficult to address because of the complexity of the different hydrodynamic and acoustic processes brought into play: hydrodynamic instabilities, shock waves disturbances, acoustic leakage and propagation, receptivity at the nozzle... each of these topics requiring deep investigations on their own. Furthermore, experimentally studying under-expanded jets to provide answers to these questions is usually challenging and studies have mainly relied on Schlieren visualizations (providing qualitative images of the flow in the form of integrated density gradients along light rays, even though quantitative information may also be retrieved [START_REF] André | Shock-tracking procedure for studying screechinduced oscillations[END_REF][START_REF] Mercier | Experimental characterisation of the screech feedback loop in underexpanded round jets[END_REF] ) and on acoustic measurements [START_REF] Michael | Acoustic study of b helical mode for choked axisymmetric nozzle[END_REF] (giving access to only one part of the screech feedback loop). While these techniques have already led to a better understanding of the screech process over the past fifty years, it is expected that more physical insight could be obtained relying on quantitative measurements of the whole hydrodynamic field. Steps in this directions have been made using more advanced measurement techniques such as Particle Image Velocimetry (PIV) [START_REF] André | Investigation of the mixing layer of underexpanded supersonic jets by particle image velocimetry[END_REF][START_REF] Daniel Edgington-Mitchell | Coherent structure and sound production in the helical mode of a screeching axisymmetric jet[END_REF] and Rayleigh scattering. [START_REF] Panda | Measurement of shock and shock-vortex interaction in underexpanded jets using rayleigh scattering[END_REF] The present work aims at going further in this direction and investigates the use of 3D Background Oriented Schlieren (3DBOS) to analyze the hydrodynamic features of the various screech modes in an underexpanded jet. Indeed, this non-intrusive measurement technique provides instantaneous and average 3D density fields, a quantity that is usually hardly accessible while being of prime importance in the study of supersonic flows.

The manuscript is organized as follows. First, descriptions of the jet facility, the 3DBOS setup, the acoustic measurements and the post-processing methodology are provided in Sec. 2. Then, near-field acoustic measurements are discussed in Sec. 3 in order to compare the properties of our screeching jet with the ones found in the literature. The first 3DBOS results are given in Sec. 4 where mean density fields are discussed. It is then detailed in Sec. 5 how phase-averaged density fields representative of the main plume instability were extracted from the 3DBOS measurements. The link between these average structures and screech noise is also discussed in this section. Finally, conclusions drawn from this preliminary work focusing on the use of 3DBOS to analyze screeching jets are given in Sec. 6.

Methodology

Jet facility

The jet studied was issued from a contoured convergent nozzle of exit diameter D = 22 mm and designed to ensure a straight sonic condition at the exit. A schematic of this nozzle is given in Fig. 1. The flow was regulated in mass-flow rate and temperature, such that the total temperature was set to T t0 = 293 K. This stagnation temperature and the Nozzle Pressure Ratio (NPR) were monitored throughout the experiment using respectively a thermocouple and a total pressure probe located 0.7 m upstream of the convergent nozzle (equivalent to ≈ 32D).

As illustrated in Fig. 1, one particular feature of the nozzle considered in this work lies in its large lip thickness t = 6.9 cm = 3.1D. This geometry has a significant effect on the screech properties of the jet, this point being further developed in Sec. 3.

3D Background Oriented Schlieren setup

This study builds on the previous work performed by Nicolas et al, [START_REF] Nicolas | 3d reconstruction of a compressible flow by synchronized multi-camera BOS[END_REF] with slight improvements brought to the 3DBOS setup in terms of spatial resolution. We briefly recall that the principle of BOS is to evaluate the apparent displacement of a pattern observed by a camera and induced by aero-optical effects associated with the presence of a gradient of index in-between (hot plume, compressibility effects, etc.). Pattern displacement fields and thus light ray deviation maps obtained from a given number of point of views can then be combined to estimate the 3D density field causing these observed displacements by solving an inverse problem.

BOS images were acquired using 8 cameras (JAI BM-500GE) placed on an azimuthal arc of 1 m radius centered on the jet axis and homogeneously distributed over an azimuthal angle of 135°, as schematically shown in Fig. 2. These 5MPx cameras have a 2/3 inch sensor with a pixel size of 3.45 µm and were equipped with 50 mm focal length lenses (Schneider). With this optical setup, each camera observed an area of about 16 cm × 14 cm in a plane containing the jet axis. Four background plates were placed on the opposite side, at a distance of 0.4 m from the jet axis. The BOS backgrounds were designed using a semi-random dot pattern printed with opaque ink on a glass beads retroreflective layer. These dots had a diameter of 0.35 mm, yielding a size of 3.5 px in the images. Unlike usual printed backgrounds that exhibit Lambertian reflectance, these retroreflective backgrounds reflect the light back to its source with minimum scattering: the quantity of light reaching the camera is larger, enabling smaller apertures and thus finer spatial resolutions.

Illumination of the backgrounds was achieved using a 532 nm double-pulse laser (Quantel EverGreen 200). Only one laser pulse was used during the cameras exposure time, leading to an effective exposure time of about 10 ns, the pulse duration. This time scale ensured that the turbulent structures of the flow imaged by the cameras did not move more than about 0.1 px during each acquisition. This was regarded as sufficient in this work to consider the acquisitions as instantaneous images of the flow. Because of the high directionality of the retroreflective backgrounds reflection, each camera needed its own light source. Consequently, the laser beam was split into eight beams using a separation table made of seven 50:50 beamsplitter plates. Eight liquid guides equipped with diverging lenses were then used to direct this light toward the backgrounds. With such a setup, the amount of light collected by each camera was sufficient to operate at an optimal f-number f # = 16. [START_REF] Nicolas | 3d reconstruction of a compressible flow by synchronized multi-camera BOS[END_REF] Under these conditions, the circle of confusion of the optical system in the region of the flow had an approximate diameter of 0.9 mm. This value is similar to the size of the interrogation window used in the post-processing algorithm (15 px) back-projected to the flow region. It is then expected that the global spatial resolution of the present BOS measurements is about 1 mm, corresponding to about 0.05D.

The 8 cameras and the laser were synchronized using a TTL pulse generator, ensuring a synchronous observation of the flow by the cameras. During the acquisition, the camera images were transferred to the RAM of a computer and stored on a SSD disk. In order to obtain well converged first and second statistical moments of the BOS displacement fields, 6000 images were acquired for the 6 aerodynamic conditions tested, corresponding to NPR ∈ {2.1, 2.3, 2.7, 4.0, 5.0, 6.6} as explained in Sec. 3. A maximum acquisition rate of 5 Hz was achieved, yielding statistically independent realizations but no relevant information about the temporal dynamics of the flow, the typical fundamental screech frequencies measured in the present experiments ranging between 3.3 and 11.2 kHz.

3D BOS processing

3DBOS measurements first require an accurate geometrical calibration of the multicamera system. This calibration process consists in identifying all cameras internal parameters (focal length, principal points and distortion parameters) and all cameras poses (position and orientation) in the same frame of reference. To do so, we relied on the multicamera calibration methodology specifically developed for 3DBOS by Le Sant et al [START_REF] Le Sant | Multi-camera calibration for 3dbos[END_REF] that consists in recording several images of a 2D calibration body randomly moved inside the common field of view of all the cameras. The projection errors of all visible dots of the calibration body are then minimized through an optimization on the calibration parameters. Note that with the cameras arranged in almost half a ring as in the present work, the calibration body cannot be observed by all the cameras at once. As a consequence, the calibration process needs to handle a sequence of incomplete but overlapping views of the body to yield a global set of calibration parameters.

BOS displacement fields were then evaluated by correlating the acquired images with a reference background image that was recorded prior to the experiment. This reference image is the result of an average of 50 reference images in order to decrease the influence of the measurement noise in the correlation process of individual images. The correlation process was performed using the ONERA software FOLKI-PIV, [START_REF] Champagnat | Fast and accurate PIV computation using highly parallel iterative correlation maximization[END_REF] implemented on GPU (Graphic Processing Unit) and relying on an optical flow approach. Interrogation windows with a size of 15 px × 15 px were used. The deviation maps were then evaluated from the estimated displacement fields relying on the approximation that light ray deviations occur at the mid-point of the investigated volume. [START_REF] Atcheson | Time-resolved 3d capture of non-stationary gas flows[END_REF] Using these camera calibration results and (instantaneous or average) deviation maps, 3D reconstructions of the density fields were performed. The software developed at ONERA for this purpose (MIRAGE) relies on a direct tomographic method, discussed in details by Nicolas et al, [START_REF] Nicolas | A direct approach for instantaneous 3d density field reconstruction from background-oriented schlieren (BOS) measurements[END_REF] aiming at inverting the equation

= G n ext ray ∇ρ ds = T (D(ρ)) = A(ρ) , (2)
where is the deviation angle of a ray of light, ρ is the density, s the curvilinear coordinate along the ray and G the Gladstone-Dale constant. The observation operator A combines the tomographic operator T and the spatial derivation operator D of the density field ρ. Inverting this problem is not straightforward since A is ill-conditioned, requiring the use of a regularization method. We consider here a simple Tikhonov regularization [START_REF] Tikhonov | Solutions of ill-posed problems[END_REF] and we thus seek a minimization of a cost function J including a least-square data-fidelity term and a regularization function based on the density gradient such that

J (ρ) = || -A(ρ)|| 2 + λ||∇ρ|| 2 . ( 3
)
The balance between the two terms in the r.h.s. is controlled by the regularization parameter λ that was identified using a L-curve method. [START_REF] Per | Analysis of discrete ill-posed problems by means of the l-curve[END_REF] Minimization of J is achieved by a conjugate gradient method. We note that other regularization methods that better account for the presence of discontinuities could have been considered, but that, overall, the present approach was observed to yield shock waves that are not excessively smoothed. Nonetheless, future work may consider the use of other regularization methods to improve the sharpness of the solution. A 3D mask in the form of a truncated cone diverging downstream was designed to reduce the number of voxels considered in the reconstruction process and to limit the diffusion of the information out of the region of interest. The use of such a mask additionally helps to reduce the computational cost. [START_REF] Atcheson | Time-resolved 3d capture of non-stationary gas flows[END_REF][START_REF] Nicolas | A direct approach for instantaneous 3d density field reconstruction from background-oriented schlieren (BOS) measurements[END_REF] Free boundary conditions were imposed on both inlet and outlet boundaries to account for the jet flow, while a (constant) external density, measured during the experiment, was enforced on the remaining lateral boundary of the 3D mask. Finally, we note that the code is implemented on a GPU architecture to reduce the computation time and to allow the reconstruction of large fields. Here the reconstruction volume was about 13 cm × 6.5 cm × 6.5 cm and made of about 4.5 × 10 6 square voxels having a size of 0.5 mm.

Acoustic measurement setup

The temporal dynamics of the jet flow was indirectly scrutinized through acoustic measurements. This was performed using 1/4 inch microphones (GRAS 46DB) having a dynamic range upper limit of 166 dB and a frequency range of 70 kHz. Two microphone setups were considered. In the first configuration, four microphones were mounted along an arc of radius r = 20D, centered on the jet axis at x = 1D and with polar angles measured from downstream θ = 30°, 50°, 70°and 90°. This setup allowed studying the jet noise directivity and was employed to analyze the screech properties of the jet presented in the following Sec. 3. In the second setup that is schematically represented in Fig. 2, eight microphones were evenly mounted (every 45°) on a ring centered on the jet axis in order to analyze the azimuthal distribution of the jet acoustic field. This ring was positioned at x = 21D and the microphones were located at r = 13D. These acoustic measurements were not performed in anechoic conditions, but care was taken to attenuate the acoustic reflections on the floor and the walls of the facility as well as on the BOS structure itself using acoustic foam panels. These muffled conditions were observed to be sufficient for the present analysis. Indeed, very limited reflection artifacts were observed on the acoustic spectra. Furthermore, it should be noted that this work is primarily focusing on the analysis of the jet screech effects that occur at very specific frequencies and with large sound pressure levels that easily emerge from the background noise. Finally, it was verified that the presence of the 3DBOS bench did not modify the frequencies and the levels of the screech noise measured, ensuring the relative non-intrusivity of our setup and measurement technique.

Screech properties of the jet studied

Preliminary acoustic measurements were performed to characterize the jet screech frequencies obtained with this setup. As mentioned in Sec. 2.1, the jet studied is issuing from a nozzle designed with a very thick lip, a geometrical feature that needs to be considered when comparing the present results with the literature. Indeed, it was shown in previous works that the lip thickness is an important parameter in the overall screech feedback process. More precisely, Ponton and Seiner 21 and Shen and Tam 24 observed a delayed screech cessation (in terms of NPR or M j values) as the lip thickness was increased, together with a significant increase in screech noise levels (also reported by Norum [START_REF] Norum | Screech suppression in supersonic jets[END_REF] and Raman [START_REF] Raman | Cessation of screech in underexpanded jets[END_REF] ). Furthermore, nozzle lip thickness was reported to induce slight shifts in screech frequencies by Aoki et al. [START_REF] Aoki | An experimental study of the nozzle lip thickness effect on supersonic jet screech tones[END_REF] Finally, this geometrical parameter was observed to play a determinant role in the transition between screech modes: [START_REF] Lima De Assunção | Experimental study of underexpanded round jets: nozzle lip thickness effects and screech closure mechanisms investigation[END_REF][START_REF] Ponton | The effects of nozzle exit lip thickness on plume resonance[END_REF] first, in thick-lip configurations, an earlier transition in terms of NPR between modes B and C was reported compared to thinlip cases; second, the range of NPR over which intermittent switching between these two modes occurs was observed to be reduced; third, mode C was observed to extend over a larger range of NPR values, meaning that the transition to mode D should occur at larger NPR values compared to thin-lip cases.

Fig. 3 was obtained using the first acoustic measurement setup presented in Sec. 2.4 using only the most downstream microphone (at θ = 30°from the jet axis). This figure provides contours of Sound Pressure Level (SPL) as a function of the NPR (or equivalently M j using Eq. ( 1)) and the Strouhal number St j ≡ f D j /U j where f is the frequency, D j the fully expanded diameter and U j the fully expanded jet velocity. Screech tones are observed in this figure as thin spectral bands of large SPL decreasing as the NPR increases. The main modes reported in the literature are labeled next to the fundamental spectral bands while the first harmonics are left aside. Together with these acoustic measurements are plotted some classical results reported in the literature. The solid black line corresponds to the semi-empirical law proposed by Tam et al [START_REF] Tam | Proposed relationship between broadband shock associated noise and screech tones[END_REF] that reads

St j = 0.67(M 2 j -1) -1 /2 1 + 0.7M j 1 + γ -1 2 M 2 j -1 /2 T a T t0 -1 /2 -1 , (4)
where T a /T t0 is the ratio of the ambient to the total temperature of the jet. As one can observe, this smooth function provides overall satisfactory estimates of the screech frequencies for modes A, C and D. Note however, that, by design, it cannot account for the modes staging phenomenon and for the frequencies of mode B. In addition to this curve, experimental results obtained by Powell [START_REF] Powell | On the mechanism of choked jet noise[END_REF] and Ponton and Seiner [START_REF] Ponton | The effects of nozzle exit lip thickness on plume resonance[END_REF] for thinner lip configurations are superposed. A satisfactory agreement in obtained for NPR values lower than 4. For NPR > 5, in agreement with the literature mentioned previously, the mode C from our thick-lip screeching jet extends over a larger range, up to NPR = 6.1. In addition, no screech cessation was observed up to the maximum NPR = 7 reached.

Based on these results, six values of NPR were selected to further investigate the different screeching modes using 3DBOS: 2.1 (mode A1), 2.3 (mode A2), 2.7 (mode B), 4.0 and 5.0 (mode C) and 6.6 (mode D). These values were chosen to lie in the middle of each mode spectral band in order to ensure stable conditions and avoid at best modes intermittency and switching.

Mean density fields results

This section details the properties of the 3D mean density fields obtained by 3DBOS for the six NPR conditions identified in Sec. 3 using acoustic measurements. The results are gathered in Fig. 4, displaying 3D density iso-surfaces and density contours in a longitudinal plane. Examining the 3D iso-surface plots, axisymmetric mean density fields are well recovered for all the conditions, which is to be expected. Shock cells appear to be furthermore well captured and one can clearly visualize the increase of the spacing between the cells as the NPR is increased. These shock cells are more readily observed in the density contours extracted from the y = 0 plane, which also illustrate the presence of Mach disks for the last three NPR conditions.

In order to compare the mean properties of the present jet with the ones reported in the literature, quantitative data regarding the topology of these jets were extracted from these 3DBOS reconstructions. First, the mean shock spacing L s was evaluated for each NPR. These results are compared in Fig. 5a with the Prandtl-Pack formula that reads

L s D = 1.306(M 2 j -1) 1 /2 . (5)
A satisfactory agreement is observed, displaying an increase of the shock spacing for increasing values of the NPR. the black diamonds are experimental data from Addy; 1 the blue plus markers were obtained from 3DBOS by Nicolas et al; [START_REF] Nicolas | 3d reconstruction of a compressible flow by synchronized multi-camera BOS[END_REF] the green crosses correspond to Schlieren visualizations; the red inverted triangle is from RANS simulation; the magenta dots correspond to the results of the present experiment. Note that the last four data sets were obtained with the same jet facility.

Second, the average diameter of the first Mach disk D M D appearing for a NPR ≈ 4 was evaluated based on the gradient amplitude of the reconstructed density fields for the three highest NPRs studied. These results are compared in Fig. 5b with experimental results reported by Nicolas et al [START_REF] Nicolas | 3d reconstruction of a compressible flow by synchronized multi-camera BOS[END_REF] and by Addy 1 for a contoured nozzle similar to the one here studied. Furthermore, estimates of D M D were also extracted from Schlieren visualizations performed under similar conditions and from a RANS simulation performed by the authors. The overall agreement with all the results is satisfactory, especially with Addy's, suggesting that the present facility and jet main characteristics do not feature significant geometrical or aerodynamic imperfections. We finally note that the main source of uncertainty in the estimation of D M D from our 3D density results is expected to lie in the spatial resolution of the technique. This uncertainty remains nonetheless sufficiently low to have confidence in the accuracy of the spatial features of the flow evidenced by the 3DBOS technique. Particularly, one can observe an improvement on the estimation of D M D compared to the previous 3DBOS setup used by Nicolas et al [START_REF] Nicolas | 3d reconstruction of a compressible flow by synchronized multi-camera BOS[END_REF] at NPR = 5.

3D density fields of screech modes

The low acquisition rate of the 3DBOS system prevents us from performing a spatio-temporal analysis of the jet dynamics associated with screech. Nonetheless, the absolute instability at the root of screech is expected to significantly alter the jet density field in the form of large coherent oscillations in the jet plume: its coherent imprints should be observed on the instantaneous BOS displacement fields obtained from each camera, as shown in Sec. 5.2. Following the work of Nicolas et al, [START_REF] Nicolas | 3d reconstruction of a compressible flow by synchronized multi-camera BOS[END_REF] we thus propose to rely on a Proper Orthogonal Decomposition (POD) of the displacement fields (and not the 3DBOS instantaneous density fields for a computational cost reason) in an attempt to isolate the main coherent features associated with the plume instability. A similar approach was also considered by Edgington-Mitchell et al [START_REF] Daniel Edgington-Mitchell | Coherent structure and sound production in the helical mode of a screeching axisymmetric jet[END_REF] to study the dynamics of the screech mode C relying on PIV data.

In the present work, the displacement field d obtained with only one camera among the eight available was considered since it was observed to provide sufficiently reliable data for this first investigation. A consequence of this choice is that complex screech dynamics such as flapping were not satisfactorily captured (that is for modes B and D), and only axisymmetric and purely helical features (modes A and C) were observed to be efficiently isolated. Hence, only these two last modes A and C (that is the two cases where NPR equals 2.1 and 4.0) are further analyzed in this section. Future work will however consider the use of all the cameras to thoroughly isolate the dominant structures at all the screech modes. Note that the choice of the camera considered for the following analysis does not modify the results: each camera having a similar point of view on the jet by rotational symmetry around the jet axis, similar decompositions will be obtained.

Proper Orthogonal Decomposition results

POD is performed using the Snapshot POD method proposed by Sirovich 25 since the number of images acquired per camera (N = 6000) is much lower than the number of displacement vectors evaluated per image (typically around 20 000). We briefly recall that this approach seeks a decomposition of the displacement fluctuation field d (x, t j ) in the form of orthonormal spatial modes φ i that optimally capture the variance of the data. We have

d (x, t j ) = N i=1 a i (t j )φ i (x) (6)
where a i (t j ) are random variables of the time t j . Note that since the present experiment is not timeresolved, these a i coefficients are functions of the discrete time variable t j and do not represent the temporal dynamics of the modes. These spatial modes φ i and expansion coefficients a i are obtained following this classical procedure: first, the (temporal) correlation matrix R associated with the displacement vectors d is evaluated using the L 2 inner product; second, an eigenvalue decomposition of R provides "temporal" modes a i and eigenvalues λ i that represent the energy content of each mode (usually sorted in descending order); finally, spatial POD modes are retrieved by the projection

φ i (x) = 1 N λ i N j=1 a i (t j )d (x, t j ) . (7)
Applying this approach to the displacement fields obtained for the two NPRs 2.1 and 4.0 yielded the distribution of eigenvalues λ i given in Fig. 6. One can clearly observe the emergence of the first two eigenvalues in both cases. For the first case (NPR = 2.1), they contribute to about 30% of the total variance, while for the second case (NPR = 4.0) to about 25%. This observation suggests that we may rely on such two modes to analyze the main coherent organization of the flow in both cases. It appears however for now very uncertain to claim here that they are exactly representative of the spatio-temporal dynamics of screech modes. This point will be further discussed in Sec. 5.3.

The spatial POD modes φ i = (φ

(x) i , φ (z)
i ) T associated with these first two eigenvalues for the two NPR cases are displayed in Fig. 7, showing both the horizontal (along x) and vertical (along z) displacement components. Examining these spatial POD modes, it can be observed that for the two aerodynamic conditions, the first two modes selected display a similar spatial structure, with comparable spatial axial wavelength (shifted streamwise by about a quarter-wavelength between the two modes). For the case NPR = 2.1, we observe a wave-packet-like perturbation that is symmetric along the axis x for the axial component and anti-symmetric for the vertical component. Such a structure is reminiscent of what would be obtained with a Kelvin-Helmholtz-type instability. [START_REF] Daniel Edgington-Mitchell | Upstream-travelling acoustic jet modes as a closure mechanism for screech[END_REF] For the second case NPR = 4.0 we observe that the two spatial modes still display a wave-packet-like structure outside of the jet (that might be the imprint of screech noise radiation) but with also a significant contribution of fluctuations from the shock cells inside the jet. The symmetry along the axis x is different compared to the case NPR = 2.1: the axial component of displacement is antisymmetric along x while the vertical component is symmetric. Such a structure is reminiscent of what would be obtained with an helical instability with unit azimuthal wavenumber. [START_REF] Daniel Edgington-Mitchell | Coherent structure and sound production in the helical mode of a screeching axisymmetric jet[END_REF] These observations are thus consistent with the general physical behavior of screech modes reported in the literature for the modes A and C here studied.

The previous observation that the first two dominant POD modes of each NPR case have similar structures suggests that they could be coupled and representative of a similar convective and oscillatory feature like the global instability associated with screech. One way to investigate such a coupling is to plot the state of each BOS acquisition in the subspace generated by the two dominant modes (φ 1 , φ 2 ), yielding the scatter plot (a 1 (t j ), a 2 (t j )). If the measurements were temporally resolved, this would lead to a Lissajous figure (or a phase-portrait in terms of dynamic systems) and a purely oscillating dynamics would generate a circle in such a phase-space. The results for the two NPR values considered are displayed in Fig. 8. One can clearly observe that both scatter plots statistically draw a circle with some radial dispersion, supporting a phase-coupling between the two modes that are thus likely representative (statistically) of a convective oscillatory dynamics linked to the plume instability. Further support on this point will be provided in Sec. 5.3.

The plots given in Fig. 8 finally offer a way to sort all the BOS images obtained with all the cameras: dividing the two plots in twelve angular sectors, we can define for each BOS acquisition a given "phase" relative to the first two dominant POD modes identified with one camera. This procedure is illustrated in this figure by filling the markers with 12 different colors. Evaluating the associated phase-average of the displacement fields for each camera based on this sorting procedure yields 12 mean displacement fields for each camera and finally 12 phase-averaged 3D density fields after reconstruction. These 3D results are displayed in Fig. 9. As one can observe, relying on this methodology leads to phase-averaged 3D density fields that are axisymmetric for NPR = 2.1 and helical for NPR = 4.0. These results are in line with the accepted topology of the plume instabilities associated with screech modes A1 and C. This is however the first time to the knowledge of the authors that such average 3D structures are actually reconstructed, assuming that the present methodology provides a representative image of the instabilities associated with screech.

We conclude this section by noting that these phase-averaged reconstructed density fields allow us to qualitatively observe the phase-averaged dynamics of the shock cells for the two NPRs considered. For NPR = 2.1 (mode A1), the first shock cells appear to be significantly disturbed by the perturbation, with forward and backward oscillations. For NPR = 4.0 (mode C), it is observed that the second, third and fourth shock cells remain at the same axial locations but appear to sustain a precession around the x axis, moving in a circular motion within a plane that tilts slightly and periodically with respect to the (yz) plane.

Instantaneous 3D density fields

To highlight the fact that 3D structures identified by the previous methodology relying on POD does not provide completely artificial results, Fig. 10 provides instantaneous 3D density reconstructions obtained for the case NPR = 4.0. These four instantaneous fields were arbitrarily chosen from the data shown in the scatter plot of Fig. 8b in the angular sectors corresponding to ψ = 0°, 90°, 180°and 270°. While these reconstructions do not provide as smooth density fields as the phase-averaged results, a similar helical organization of the flow can clearly be observed, confirming the physical relevance of the 3D structures previously identified.

On the link between POD results and screech modes

Finally, a central question that one may ask is whether the 3D coherent structures previously identified using (spatial-only) POD are relevant in the description of the (spatio-temporal) screech modes or not. To investigate the possible connection between POD modes and screech noise, 3DBOS measurements were performed together with acoustic measurements using the second acoustic setup presented in Sec. 2.4, that is with an azimuthal array of microphones shown in Fig. 2. For each BOS acquisition, an oscilloscope acquired the laser Q-switch signal, defining the exact time at which BOS images were recorded, and the acoustic signal of six microphones over at least ten acoustic periods associated with the screech frequency. These acoustic measurements were then post-processed using a narrow digital 100 Hz band-pass filter centered on the screech frequency (with zero phase-shift) and a Hilbert transform to provide a phase estimate of each BOS acquisition with respect to screech. Note that in this process, one actually needs to account for the phase-shift induced by the location of the microphones with respect to the source of screech noise, a distance that is not precise and likely stochastically fluctuates due to the motion of the shock waves.

To illustrate the results obtained following this approach, let us consider the case NPR = 4.0. Similarly to what was performed with POD, BOS acquisitions were divided into 12 "screech phase" classes. An easy way to visualize the results in comparison with the ones obtained by POD is to consider the scatter plot of POD coefficients in Fig. 8b and to set the colors of the markers as a function of the "screech phase" here evaluated. If the two approaches were equivalent, the colored version of Fig. 8b would be exactly retrieved. The plot obtained is shown in Fig. 11: a dispersion of the markers of similar color can be observed, suggesting that some differences exist. Nonetheless, one can notice that this dispersion is not uniform, but that it is centered around each phase class considered, over an angle of about 90°. This Fig. 11 appears as a dispersed version of Fig. 8b maintaining the overall phase class order previously identified. We thus suggest that, while the details of the two approaches are different, they capture on average the same physics. In other words, we believe that these results support the idea that the POD analysis performed previously is actually relevant in isolating the average modes associated with screech noise. The dispersion observed in the present results is likely the consequence of either some jittering of the screech modes (we observed for example slight variations of its frequency over a complete acquisition sequence) or some spatial fluctuation of the source of screech with respect to the microphone locations; an other possible origin may be the uncertainty in the phase estimation process using the Hilbert transform.

Conclusion

In this work, the non-intrusive, seedless 3DBOS technique was applied to measure the 3D density fields of under-expanded screeching jets. Compared to the previous work of Nicolas et al, [START_REF] Nicolas | 3d reconstruction of a compressible flow by synchronized multi-camera BOS[END_REF] a significant improvement on the spatial resolution of the technique was achieved using retro-reflective backgrounds, yielding smaller lens apertures and reduced astigmatism effects. As a consequence, satisfactory mean and instantaneous 3D density reconstructions were obtained, enabling a future deeper analysis of the structure of the 3D shock cells and their impact on the downstream instabilities associated with screech. Furthermore, it was shown that the fluctuating density field associated with the main (average) instability wave driving screech for modes A1 and C could be satisfactorily isolated from the BOS measurements relying on the first two POD modes of the BOS displacement fields. The relevance of these isolated coherent structures in the description of the screech process was evidenced using acoustic measurements. A clear 3D visualization of the two modes A1 and C was obtained, displaying axisymmetric and helical structures. The other modes proved more difficult to analyze following this methodology and future work will address this point.

Figure 1 :

 1 Figure 1: Schematic view of the contoured nozzle geometry studied issuing an under-expanded supersonic jet visualized using Schlieren photography; the field of view represented corresponds to the approximate one obtained with the 3DBOS setup.

Figure 2 :

 2 Figure 2: Schematic view of the 3DBOS setup mounted around the supersonic jet nozzle axis; the 8 cameras are placed on an arc covering an azimuthal angle of 135°and the 8 microphones are distributed on a circle every 45°.

Figure 3 :

 3 Figure 3: Contours of sound pressure level (SPL) as function of jet operating conditions (NPR and Mach number M j ). The solid black line is the semi-empirical relation proposed by Tam, 26 the blue and the red dots are representative of measurements obtained by Ponton and Seiner 21 and Powell 22 respectively.

  NPR = 6.6

Figure 4 :

 4 Figure 4: Mean 3DBOS results obtained for the six NPR selected; left column: 3D iso-surfaces of mean density; right column: contours of mean density in the plane y = 0.

Figure 5 :

 5 Figure 5: (a) Mean shock spacing L s /D as a function of jet operating conditions: the solid line is the Prandtl-Pack formula (Eq. (5)) and red plus markers correspond to the present mean 3DBOS results; (b) First Mach disk diameter D M D /D as a function of the jet NPR: the blue solid line is the relation proposed by Addy; 1the black diamonds are experimental data from Addy; 1 the blue plus markers were obtained from 3DBOS by Nicolas et al;[START_REF] Nicolas | 3d reconstruction of a compressible flow by synchronized multi-camera BOS[END_REF] the green crosses correspond to Schlieren visualizations; the red inverted triangle is from RANS simulation; the magenta dots correspond to the results of the present experiment. Note that the last four data sets were obtained with the same jet facility.

2 (

 2 b) NPR = 4.0 (mode C)

Figure 6 :

 6 Figure 6: Normalized POD eigenvalues λ i /σ 2 obtained on the BOS displacement fields of one camera for two NPR values. These eigenvalues are made dimensionless using the total displacement variance σ 2 . The cumulative sum of these eigenvalues (right vertical axis) indicates the fraction of energy captured with an increasing number of modes.

Figure 7 :

 7 Figure 7: Axial (φ (x) i ) and vertical (φ (z) i ) components of the first two POD spatial modes (i ∈ [1..2]) obtained by considering the BOS displacement fields of one camera for 2 values of NPR.

a 1 a 1 Figure 8 :Figure 9 :

 1189 Figure 8: Scatter plots (a 1 (t j ), a 2 (t j )) giving the state of all BOS acquisitions (of one camera, j ∈ [1..N ]) in the subspace generated by the first two POD spatial modes φ i (x) i∈[1..2] for the two NPR conditions corresponding to screech modes A1 and C. Colors highlight the division of the scatter plot to define 12 phase classes.

a 1 Figure 11 :

 111 Figure 11: Scatter plot similar to Fig. 8b, but with marker colored by the phase value determined by acoustic measurements.

Acknowledgments

This work has benefited from fruitful discussions with F. Nicolas, F. Champagnat and J. Troyes. F. Micheli is warmly acknowledged for his help during the 3DBOS measurement campaign. LL was supported by the fellowship from the Région Occitanie, France.